Here's a cross section of the situation, with two oriented planes. If you're rotating the red/green plane in the direction of its green side, does that mean clockwise or counterclockwise?
So I don't think you can use the orientation of the planes to fix an interpretation of the angle. However, you write that you have a preferred orientation of the line of intersection, and that sounds more promising.
Namely, you can decide to turn the first plane counterclockwise when looking in the preferred direction along the intersection line.
Somewhat tedious formula-churning implementation of this: Fix some point on the intersection line and consider the plane that is normal to it at that point. Equip that plane with an orthonormal coordinate system that creates a right-handed system together with the preferred direction of the intersection. The orthonormal coordinates create a parametric representation of the normal plane; compose that with the equations for the two original planes; that gives 2D equations within the normal plane for its intersection with the original plane. Now find the angle in the normal plane using standard 2D techniques (e.g as the difference of atan2 values).
Computationally slicker: You can find $|\cos\theta|$ by the dot product or normals in 3D without needing to project into a specific cross section; the only remaining information you need is whether to use $\theta_0$, $\pi-\theta_0$, $\theta_0-\pi$ or $-\theta_0$, where $\theta_0$ is the arccosine of the dot product. But you can actually distinguish between those four cases simply by the combination of (a) the sign of the dot product and (b) whether the cross product of the normals is parallel or antiparallel to the chosen direction of the intersection line. I'd probably make some embarassing mistake if I tried to work our exactly what the connection should be, but there are only 4! possible ways to do it anyway, so just figure it out with trial and error...
The general equation of a plane in 3-D is given by $$\mathbf{(p-p_0).n}=0$$
where $\mathbf{p}$ is any general point on the plane, and $\mathbf{p_0}$ is any known point on the plane. $\mathbf{n}$ is a vector normal to the plane.
The equation of a line is given by $$\mathbf{p = l_0+}t\mathbf{l} $$
where $\mathbf{l_0}$ is any point on the line.
If the line lies in the plane, it must satisfy two conditions-
It must be perpendicular to the normal to the plane i.e. $\mathbf{l.n}=0$
$\mathbf{l_0}$ must lie in the plane i.e. satisfy the plane's equation. So, $\mathbf{(l_0-p_0).n} = 0$
You can calculate $\mathbf{p_0,l_0}$ very easily with the information you have.
$\mathbf{p_0}$ can be calculated by choosing any two of $x,y,z$ and finding the third to satisfy the equation of the plane. $\mathbf{l_0}$ is precisely $\vec O$ that you already know. You can use them to cross-check whether the fit is good or not. But they do not depend on $\mathbf{l}$. So, they are secondary, and may be used as a sanity check later on.
Suppose you are given the equation of the line as $ax+by+cx+d = 0$. You can recast it as $[(x,y,z)-\mathbf{p_0}].(a,b,c)$. So, $(a,b,c)$ is your normal vector.
Suppose you have $m$ planes with normals $\mathbf{n_1,n_2,\ldots,n_m}$. Then, your overall constraints are $$\mathbf{l.n_1} = 0 \\ \mathbf{l.n_2} = 0 \\ \ldots \\\mathbf{l.n_m}=0$$.
It is a system of linear equations with 3 variables and $m$ equations. You need to find out $\mathbf{l}$, which can be found up to a constant factor by the standard least squares method, provided $m>3$, which should not be a worry for you. If you code carefully enough, you can implement all of it in matrix terms.
Hope it helps.
Best Answer
The normal of a plane, and the direction of a line, are each unique up to non-zero scalar multiplication.
The angle between two lines, between two planes, and between a line & a plane are each (defined to be) between $0$ and $\frac{\pi}2$ (inclusive); on the other hand, the angle between two vectors is (defined to be) between $0$ and $\pi$ (inclusive).