Geometry

Collaboration diagram for Geometry:

Classes
struct	FitPoint
struct	fit_point2D
class	Bspline_model

Geometrical operations
void	Uproject_to_plane (const Matrix1D< double > &point, const Matrix1D< double > &direction, double distance, Matrix1D< double > &result)
void	Uproject_to_plane (const Matrix1D< double > &r, double rot, double tilt, double psi, Matrix1D< double > &result)
void	Uproject_to_plane (const Matrix1D< double > &r, const Matrix2D< double > &euler, Matrix1D< double > &result)
double	spherical_distance (const Matrix1D< double > &r1, const Matrix1D< double > &r2)
double	point_line_distance_3D (const Matrix1D< double > &p, const Matrix1D< double > &a, const Matrix1D< double > &v)
double	point_plane_distance_3D (const Matrix1D< double > &p, const Matrix1D< double > &a, const Matrix1D< double > &v)
void	least_squares_plane_fit (FitPoint *IN_points, int Npoints, double &plane_A, double &plane_B, double &plane_C)
void	least_squares_plane_fit_All_Points (const MultidimArray< double > &Image, double &plane_A, double &plane_B, double &plane_C)
void	least_squares_line_fit (const std::vector< fit_point2D > &IN_points, double &line_A, double &line_B)
void	Bspline_model_fitting (const std::vector< FitPoint > &IN_points, int SplineDegree, int l0, int lF, int m0, int mF, double h_x, double h_y, double x0, double y0, Bspline_model &result)
void	rectangle_enclosing (const Matrix1D< double > &v0, const Matrix1D< double > &vF, const Matrix2D< double > &V, Matrix1D< double > &corner1, Matrix1D< double > &corner2)
void	box_enclosing (const Matrix1D< double > &v0, const Matrix1D< double > &vF, const Matrix2D< double > &V, Matrix1D< double > &corner1, Matrix1D< double > &corner2)
bool	point_inside_polygon (const std::vector< Matrix1D< double > > &polygon, const Matrix1D< double > &point)
void	def_affinity (double u1x, double u1y, double u2x, double u2y, double u3x, double u3y, double t1x, double t1y, double t2x, double t2y, double t3x, double t3y, Matrix2D< double > &A, Matrix1D< double > &T, Matrix2D< double > &invW)
double	triangle_area (double x1, double y1, double x2, double y2, double x3, double y3)
int	line_plane_intersection (const Matrix1D< double > normal_plane, const Matrix1D< double > vector_line, Matrix1D< double > &intersection_point, const Matrix1D< double > point_line, double point_plane_at_x_y_zero=0.)

Euler operations
template<typename T >
void	Euler_angles2matrix (T a, T b, T g, Matrix2D< T > &A, bool homogeneous=false)
void	Euler_anglesZXZ2matrix (double a, double b, double g, Matrix2D< double > &A, bool homogeneous=false)
double	Euler_distanceBetweenMatrices (const Matrix2D< double > &E1, const Matrix2D< double > &E2)
template<typename T >
T	Euler_distanceBetweenAngleSets (T rot1, T tilt1, T psi1, T rot2, T tilt2, T psi2, bool only_projdir)
template<typename T >
T	Euler_distanceBetweenAngleSets_fast (const Matrix2D< T > &E1, T rot2, T tilt2, T psi2, bool only_projdir, Matrix2D< T > &E2)
void	Euler_Angles_after_compresion (const double rot, double tilt, double psi, double &new_rot, double &new_tilt, double &new_psi, Matrix2D< double > &D)
void	Euler_direction (double alpha, double beta, double gamma, Matrix1D< double > &v)
void	Euler_direction2angles (Matrix1D< double > &v, double &alpha, double &beta, double &gamma)
void	Euler_matrix2angles (const Matrix2D< double > &A, double &alpha, double &beta, double &gamma, bool homogeneous=false)
void	Euler_up_down (double rot, double tilt, double psi, double &newrot, double &newtilt, double &newpsi)
void	Euler_another_set (double rot, double tilt, double psi, double &newrot, double &newtilt, double &newpsi)
void	Euler_mirrorY (double rot, double tilt, double psi, double &newrot, double &newtilt, double &newpsi)
void	Euler_mirrorX (double rot, double tilt, double psi, double &newrot, double &newtilt, double &newpsi)
void	Euler_mirrorXY (double rot, double tilt, double psi, double &newrot, double &newtilt, double &newpsi)
void	Euler_apply_transf (const Matrix2D< double > &L, const Matrix2D< double > &R, double rot, double tilt, double psi, double &newrot, double &newtilt, double &newpsi)
void	Euler_rotate (const MultidimArray< double > &V, double rot, double tilt, double psi, MultidimArray< double > &result)
void	Euler_rotate (const MultidimArrayGeneric &V, double rot, double tilt, double psi, MultidimArray< double > &result)
void	computeCircleAroundE (const Matrix2D< double > &E, double angCircle, double angStep, std::vector< double > &outputEulerAngles)
#define	EULER_CLIPPING(rot, tilt, psi)
#define	EULER_CLIPPING_RAD(rot, tilt, psi)

Intersections
double	intersection_unit_sphere (const Matrix1D< double > &u, const Matrix1D< double > &r)
double	intersection_unit_cylinder (const Matrix1D< double > &u, const Matrix1D< double > &r)
double	intersection_unit_cube (const Matrix1D< double > &u, const Matrix1D< double > &r)

Detailed Description

Macro Definition Documentation

EULER_CLIPPING

#define EULER_CLIPPING	(		rot,
			tilt,
			psi
	)

Value:

rot = realWRAP(rot, 0, 360); \
    tilt = realWRAP(tilt, 0, 360); \
    psi = realWRAP(psi, 0, 360);

Getting the Euler angles to a range (0-360). No direction equivalence is applied, ie, there is no correction of the direction of projection making use that a view from the top is the same as a view from the bottom but reversed ... Just a wrapping of the angles is done until the angles fall in the specified ranges. The angles given must be variables and they are modified with the new values.

Definition at line 471 of file geometry.h.

EULER_CLIPPING_RAD

#define EULER_CLIPPING_RAD	(		rot,
			tilt,
			psi
	)

Value:

rot = realWRAP(rot, 0, 2.0*PI); \
    tilt = realWRAP(tilt, 0, 2.0*PI); \
    psi = realWRAP(psi, 0, 2.0*PI);

Getting the Euler angles to a range (0-2*PI).

The same as before but the angles are expressed in radians.

Definition at line 480 of file geometry.h.

Function Documentation

box_enclosing()

void box_enclosing	(	const Matrix1D< double > &	v0,
		const Matrix1D< double > &	vF,
		const Matrix2D< double > &	V,
		Matrix1D< double > &	corner1,
		Matrix1D< double > &	corner2
	)

Box which encloses a deformed box

Given a box characterized by the top-left corner (most negative) and the right-bottom (most positive) corner, and given a matrix after which the box is deformed. Which is the minimum box which encloses the preceding one? This function is useful for stablishing for loops which will cover for sure the deformed box. All vectors are supposed to be 3x1 and the deformation matrix is 3x3. The corner (x0,y0,z0) goes to V*(x0,y0,z0)' and (xF,yF,zF) to V*(xf,yF,zF)'. After that you can make a loop from corner1 to corner2.

The v0 and vF vectors can be reused as outputs.

Definition at line 366 of file geometry.cpp.

{
    SPEED_UP_temps012;
    Matrix1D<double> v(3);
    corner1.resize(3);
    corner2.resize(3);

    // Store values for reusing input as output vectors
    double XX_v0 = XX(v0);
    double YY_v0 = YY(v0);
    double ZZ_v0 = ZZ(v0);
    double XX_vF = XX(vF);
    double YY_vF = YY(vF);
    double ZZ_vF = ZZ(vF);

    VECTOR_R3(v, XX_v0, YY_v0, ZZ_v0);
    M3x3_BY_V3x1(v, V, v);
    XX(corner1) = XX(v);
    XX(corner2) = XX(v);
    YY(corner1) = YY(v);
    YY(corner2) = YY(v);
    ZZ(corner1) = ZZ(v);
    ZZ(corner2) = ZZ(v);

#define DEFORM_AND_CHOOSE_CORNERS3D \
    M3x3_BY_V3x1(v,V,v); \
    XX(corner1)=XMIPP_MIN(XX(corner1),XX(v)); \
    XX(corner2)=XMIPP_MAX(XX(corner2),XX(v)); \
    YY(corner1)=XMIPP_MIN(YY(corner1),YY(v)); \
    YY(corner2)=XMIPP_MAX(YY(corner2),YY(v)); \
    ZZ(corner1)=XMIPP_MIN(ZZ(corner1),ZZ(v)); \
    ZZ(corner2)=XMIPP_MAX(ZZ(corner2),ZZ(v));

    VECTOR_R3(v, XX_vF, YY_v0, ZZ_v0);
    DEFORM_AND_CHOOSE_CORNERS3D;
    VECTOR_R3(v, XX_v0, YY_vF, ZZ_v0);
    DEFORM_AND_CHOOSE_CORNERS3D;
    VECTOR_R3(v, XX_vF, YY_vF, ZZ_v0);
    DEFORM_AND_CHOOSE_CORNERS3D;
    VECTOR_R3(v, XX_v0, YY_v0, ZZ_vF);
    DEFORM_AND_CHOOSE_CORNERS3D;
    VECTOR_R3(v, XX_vF, YY_v0, ZZ_vF);
    DEFORM_AND_CHOOSE_CORNERS3D;
    VECTOR_R3(v, XX_v0, YY_vF, ZZ_vF);
    DEFORM_AND_CHOOSE_CORNERS3D;
    VECTOR_R3(v, XX_vF, YY_vF, ZZ_vF);
    DEFORM_AND_CHOOSE_CORNERS3D;
}

Bspline_model_fitting()

void Bspline_model_fitting	(	const std::vector< FitPoint > &	IN_points,
		int	SplineDegree,
		int	l0,
		int	lF,
		int	m0,
		int	mF,
		double	h_x,
		double	h_y,
		double	x0,
		double	y0,
		Bspline_model &	result
	)

Least-squares fit of a B-spline 2D model

For fitting a set of values that are distributed between (x0,y0) and (xF,yF) with a cubic Bspline centered on each corner, the right call is

Bspline_model model;
Bspline_model_fitting(list_of_points, 3, -1, 2, -1, 2, xF-x0, yF-y0, x0,
    y0, model);

Once the model is returned you can evaluate it at any point simply by

model.evaluate(x,y);

Definition at line 243 of file geometry.cpp.

{
    // Initialize model
    result.l0 = l0;
    result.lF = lF;
    result.m0 = m0;
    result.mF = mF;
    result.x0 = x0;
    result.y0 = y0;
    result.SplineDegree = SplineDegree;
    result.h_x = h_x;
    result.h_y = h_y;
    result.c_ml.initZeros(mF - m0 + 1, lF - l0 + 1);
    STARTINGY(result.c_ml) = m0;
    STARTINGX(result.c_ml) = l0;

    // Modify the list of points to include the weight
    int Npoints = IN_points.size();
    std::vector<FitPoint> AUX_points = IN_points;
    for (int i = 0; i < Npoints; ++i)
    {
        double sqrt_w = sqrt(AUX_points[i].w);
        AUX_points[i].x *= sqrt_w;
        AUX_points[i].y *= sqrt_w;
        AUX_points[i].z *= sqrt_w;
    }

    // Now solve the normal linear regression problem
    // Ax=B
    // A=system matrix
    // x=B-spline coefficients
    // B=vector of measured values
    int Ncoeff = YSIZE(result.c_ml) * XSIZE(result.c_ml);
    Matrix2D<double> A(Npoints, Ncoeff);
    Matrix1D<double> B(Npoints);
    for (int i = 0; i < Npoints; ++i)
    {
        B(i) = AUX_points[i].z;
        double xarg = (AUX_points[i].x - x0) / h_x;
        double yarg = (AUX_points[i].y - y0) / h_y;
        for (int m = m0; m <= mF; ++m)
            for (int l = l0; l <= lF; ++l)
            {
                double coeff=0.;
                switch (SplineDegree)
                {
                case 2:
                    coeff = Bspline02(xarg - l) * Bspline02(yarg - m);
                    break;
                case 3:
                    coeff = Bspline03(xarg - l) * Bspline03(yarg - m);
                    break;
                case 4:
                    coeff = Bspline04(xarg - l) * Bspline04(yarg - m);
                    break;
                case 5:
                    coeff = Bspline05(xarg - l) * Bspline05(yarg - m);
                    break;
                case 6:
                    coeff = Bspline06(xarg - l) * Bspline06(yarg - m);
                    break;
                case 7:
                    coeff = Bspline07(xarg - l) * Bspline07(yarg - m);
                    break;
                case 8:
                    coeff = Bspline08(xarg - l) * Bspline08(yarg - m);
                    break;
                case 9:
                    coeff = Bspline09(xarg - l) * Bspline09(yarg - m);
                    break;
                }
                A(i, (m - m0)*XSIZE(result.c_ml) + l - l0) = coeff;
            }
    }

    Matrix1D<double> x = (A.transpose() * A).inv() * (A.transpose() * B);
    for (int m = m0; m <= mF; ++m)
        for (int l = l0; l <= lF; ++l)
            result.c_ml(m, l) = x((m - m0) * XSIZE(result.c_ml) + l - l0);
}

computeCircleAroundE()

void computeCircleAroundE	(	const Matrix2D< double > &	E,
		double	angCircle,
		double	angStep,
		std::vector< double > &	outputEulerAngles
	)

Compute circle around Euler matrix

Given an input Euler matrix, this function returns a set of Euler angles such that they sample a circle around the original projection direction (a sample every angStep). The projection directions in the circle are separated by angCircle.

The output is in outputEulerAngles whose structure is (newrot1,newtilt1,newpsi1,newrot2,newtilt2,newpsi2,...)

Definition at line 1078 of file geometry.cpp.

{
    outputEulerAngles.clear();

    // Get the projection direction and a perpendicular direction
    Matrix1D<double> projectionDirection, perpendicular;
    E.getRow(1,perpendicular);
    E.getRow(2,projectionDirection);
    Matrix2D<double> newEt;
    newEt = E.transpose();
    Matrix2D<double> rotStep, sampling;
    rotation3DMatrix(angCircle,perpendicular,rotStep,false);
    rotation3DMatrix(angStep,projectionDirection,sampling,false);

    // Now rotate
    newEt = rotStep*newEt;
    for (double i = 0; i < 360; i += angStep)
    {
        newEt=sampling*newEt;

        // Normalize
        for (int c=0; c<3; c++)
        {
            Matrix1D<double> aux;
            newEt.getCol(c,aux);
            aux/=aux.module();
            newEt.setCol(c,aux);
        }
        Matrix2D<double> newE=newEt.transpose();

        double newrot, newtilt, newpsi;
        Euler_matrix2angles(newE,newrot,newtilt,newpsi);
        outputEulerAngles.push_back(newrot);
        outputEulerAngles.push_back(newtilt);
        outputEulerAngles.push_back(newpsi);
    }
}

def_affinity()

void def_affinity	(	double	u1x,
		double	u1y,
		double	u2x,
		double	u2y,
		double	u3x,
		double	u3y,
		double	t1x,
		double	t1y,
		double	t2x,
		double	t2y,
		double	t3x,
		double	t3y,
		Matrix2D< double > &	A,
		Matrix1D< double > &	T,
		Matrix2D< double > &	invW
	)

Affine transformation

Definition at line 441 of file geometry.cpp.

{
    double den = (u1x*u2y - u1y*u2x - u1x*u3y + u1y*u3x + u2x*u3y - u2y*u3x);

    //std::cout << "den" <<k << std::endl;

    invW.initZeros(6,6);
    MAT_ELEM(invW,0,0) = MAT_ELEM(invW,2,3) = (u2y - u3y)/den;
    MAT_ELEM(invW,0,1) = MAT_ELEM(invW,2,4) = -(u1y - u3y)/den;
    MAT_ELEM(invW,0,2) = MAT_ELEM(invW,2,5) = (u1y - u2y)/den;

    MAT_ELEM(invW,1,0) = MAT_ELEM(invW,3,3) = -(u2x - u3x)/den;
    MAT_ELEM(invW,1,1) = MAT_ELEM(invW,3,4) = (u1x - u3x)/den;
    MAT_ELEM(invW,1,2) = MAT_ELEM(invW,3,5) = -(u1x - u2x)/den;

    MAT_ELEM(invW,4,0) = MAT_ELEM(invW,5,3) = (u2x*u3y - u2y*u3x)/den;
    MAT_ELEM(invW,4,1) = MAT_ELEM(invW,5,4) = -(u1x*u3y - u1y*u3x)/den;
    MAT_ELEM(invW,4,2) = MAT_ELEM(invW,5,5) = (u1x*u2y - u1y*u2x)/den;



    Matrix1D<double> t_vec;
    t_vec.initZeros(6);
    VEC_ELEM(t_vec,0) = t1x;
    VEC_ELEM(t_vec,1) = t2x;
    VEC_ELEM(t_vec,2) = t3x;
    VEC_ELEM(t_vec,3) = t1y;
    VEC_ELEM(t_vec,4) = t2y;
    VEC_ELEM(t_vec,5) = t3y;

    double dettt = invW.det();
    //std::cout << "determinant" << dettt << std::endl;

    //Matrix1D<double> sol = invW*t_vec;
    //std::cout << "sol = " << sol << std::endl;

    if (fabs(dettt) <  DBL_EPSILON )
    {
        //std::cout << "I'm in if" << std::endl;
        A.initZeros(2,2);
        T.initZeros(2);
        VEC_ELEM(T,0) = DBL_MAX ;
        VEC_ELEM(T,1) = DBL_MAX ;
    }
    else
    {
        //std::cout << "I'm in else" << std::endl;
        Matrix1D<double> sol = invW*t_vec;
        //std::cout << "sol = " << sol << std::endl;

        A.initZeros(2,2);
        T.initZeros(2);
        //std::cout << A << std::endl;
        MAT_ELEM(A,0,0) = VEC_ELEM(sol,0);
        MAT_ELEM(A,0,1) = VEC_ELEM(sol,1);
        MAT_ELEM(A,1,0) = VEC_ELEM(sol,2);
        MAT_ELEM(A,1,1) = VEC_ELEM(sol,3);
        VEC_ELEM(T,0) = VEC_ELEM(sol,4);
        VEC_ELEM(T,1) = VEC_ELEM(sol,5);

        //std::cout << "A ==" << sol << std::endl;
    }
}

Euler_angles2matrix()

template<typename T >

void Euler_angles2matrix	(	T	a,
		T	b,
		T	g,
		Matrix2D< T > &	A,
		bool	homogeneous = false
	)

Euler angles –> Euler matrix.

This function returns the transformation matrix associated to the 3 given Euler angles (in degrees).

As an implementation note you might like to know that this function calls always to Matrix2D::resize

See http://xmipp.cnb.csic.es/twiki/bin/view/Xmipp/EulerAngles for a description of the Euler angles.

Definition at line 624 of file geometry.cpp.

{
    static_assert(std::is_floating_point<T>::value, "Only float and double are allowed as template parameters");

    if (homogeneous)
    {
        A.initZeros(4,4);
        MAT_ELEM(A,3,3)=1;
    }
    else
        if (MAT_XSIZE(A) != 3 || MAT_YSIZE(A) != 3)
            A.resizeNoCopy(3, 3);

    T ca = std::cos(DEG2RAD(alpha));
    T sa = std::sin(DEG2RAD(alpha));
    T cb = std::cos(DEG2RAD(beta));
    T sb = std::sin(DEG2RAD(beta));
    T cg = std::cos(DEG2RAD(gamma));
    T sg = std::sin(DEG2RAD(gamma));

    T cc = cb * ca;
    T cs = cb * sa;
    T sc = sb * ca;
    T ss = sb * sa;

    MAT_ELEM(A, 0, 0) =  cg * cc - sg * sa;
    MAT_ELEM(A, 0, 1) =  cg * cs + sg * ca;
    MAT_ELEM(A, 0, 2) = -cg * sb;
    MAT_ELEM(A, 1, 0) = -sg * cc - cg * sa;
    MAT_ELEM(A, 1, 1) = -sg * cs + cg * ca;
    MAT_ELEM(A, 1, 2) = sg * sb;
    MAT_ELEM(A, 2, 0) =  sc;
    MAT_ELEM(A, 2, 1) =  ss;
    MAT_ELEM(A, 2, 2) = cb;
}

Euler_Angles_after_compresion()

void Euler_Angles_after_compresion	(	const double	rot,
		double	tilt,
		double	psi,
		double &	new_rot,
		double &	new_tilt,
		double &	new_psi,
		Matrix2D< double > &	D
	)

Angles after compresion

Let be two volumes f and g related by g(x,y,z) = f(D(x,y,z)) (where D is a lineal transformation) then the projection direction parallel to the vector w in f is going to be related with the projection direction parallel to the vector w_prime in g. Given the w Euler angles this routine provide the w_prime angles

Definition at line 955 of file geometry.cpp.

{
    Matrix1D<double> w(3);
    Matrix1D<double> new_w(3);
    Matrix2D<double> D_1(3, 3);

    //if D has not inverse we are not in business
    D_1 = D.inv();

    Euler_direction(rot, tilt, psi, w);
    if (fabs(w(2)) > 0.999847695)/*cos one degree */
    {
        Euler_direction(rot, 10., psi, w);
        new_w = (Matrix1D<double>)(D_1 * w) / ((D_1 * w).module());
        Euler_direction2angles(new_w, new_rot, new_tilt, new_psi);

        Euler_direction(rot, tilt, psi, w);
        new_w = (Matrix1D<double>)((D_1 * w) / ((D_1 * w).module()));
        new_tilt = SGN(new_tilt) * fabs(ACOSD(new_w(2)));
        new_psi = psi;

        // so, for small tilt the value of the rot is not realiable
        // doubleo overcome this problem I first calculate the rot for
        // any arbitrary large tilt angle and the right rotation
        // and then I calculate the new tilt.
        // Please notice that the new_rotation is not a funcion of
        // the old tilt angle so I can use any arbitrary tilt angle
    }
    else
    {
        new_w = (Matrix1D<double>)(D_1 * w) / ((D_1 * w).module());
        Euler_direction2angles(new_w, new_rot, new_tilt, new_psi);
        new_psi = psi;
    }
}

Euler_anglesZXZ2matrix()

void Euler_anglesZXZ2matrix	(	double	a,
		double	b,
		double	g,
		Matrix2D< double > &	A,
		bool	homogeneous = false
	)

Euler angles –> Euler matrix.

This function returns the transformation matrix associated to the 3 given Euler angles (in degrees).

Definition at line 661 of file geometry.cpp.

{
    Matrix2D<double> RZ1, RX2, RZ3;
    rotation3DMatrix(a,'Z',RZ1,homogeneous);
    rotation3DMatrix(b,'X',RX2,homogeneous);
    rotation3DMatrix(g,'Z',RZ3,homogeneous);
    A=RZ3*RX2*RZ1;
}

Euler_another_set()

void Euler_another_set	(	double	rot,
		double	tilt,
		double	psi,
		double &	newrot,
		double &	newtilt,
		double &	newpsi
	)

The same view but differently expressed

As you know a projection view from a point can be expressed with different sets of Euler angles. This function gives you another expression of the Euler angles for this point of view. Exactly the operation performed is:

newrot = rot + 180;
newtilt = -tilt;
newpsi = -180 + psi;

Euler_another_set(rot, tilt, psi, newrot, newtilt, newpsi);

Definition at line 1002 of file geometry.cpp.

{
    newrot  = rot + 180;
    newtilt = -tilt;
    newpsi  = -180 + psi;
}

Euler_apply_transf()

void Euler_apply_transf	(	const Matrix2D< double > &	L,
		const Matrix2D< double > &	R,
		double	rot,
		double	tilt,
		double	psi,
		double &	newrot,
		double &	newtilt,
		double &	newpsi
	)

Apply a geometrical transformation

The idea behind this function is the following. 3 Euler angles define a point of view for a projection, but also a coordinate system. You might apply a geometrical transformation to this system, and then compute back what the Euler angles for the new system are. This could be used to "mirror" points of view, rotate them and all the stuff. The transformation matrix must be 3x3 but it must transform R3 vectors into R3 vectors (that is a normal 3D transformation matrix when vector coordinates are not homogeneous) and it will be applied in the sense:

New Euler matrix = L * Old Euler matrix * R

where you know that the Euler matrix rows represent the different system axes. See Euler_angles2matrix for more information about the Euler coordinate system.

Matrix2D< double > R60 = rotation3DMatrix(60, 'Z');
R60.resize(3, 3); // Get rid of homogeneous part
Matrix2D< double > I(3, 3);
I.initIdentity();
Euler_apply_transf(I, R60, rot, tilt, psi, newrot, newtilt, newpsi);

Definition at line 1038 of file geometry.cpp.

{

    Matrix2D<double> euler(3, 3), temp;
    Euler_angles2matrix(rot, tilt, psi, euler);
    temp = L * euler * R;
    Euler_matrix2angles(temp, newrot, newtilt, newpsi);
}

Euler_direction()

void Euler_direction	(	double	alpha,
		double	beta,
		double	gamma,
		Matrix1D< double > &	v
	)

Euler direction

This function returns a vector parallel to the projection direction. Resizes v if needed

Definition at line 721 of file geometry.cpp.

{
    double ca, sa, cb, sb;
    double sc, ss;

    v.resize(3);
    alpha = DEG2RAD(alpha);
    beta  = DEG2RAD(beta);

    ca = cos(alpha);
    cb = cos(beta);
    sa = sin(alpha);
    sb = sin(beta);
    sc = sb * ca;
    ss = sb * sa;

    v(0) = sc;
    v(1) = ss;
    v(2) = cb;
}

Euler_direction2angles()

void Euler_direction2angles	(	Matrix1D< double > &	v,
		double &	alpha,
		double &	beta,
		double &	gamma
	)

Euler direction2angles

This function returns the 3 Euler angles associated to the direction given by the vector v. The 3rd Euler angle is set always to 0

Definition at line 746 of file geometry.cpp.

{
    double abs_ca, sb, cb;
    double aux_alpha;
    double aux_beta;
    double error, newerror;
    Matrix1D<double> v_aux;
    Matrix1D<double> v;

    //if not normalized do it so
    v.resize(3);
    v = v0;
    v.selfNormalize();

    v_aux.resize(3);
    cb = v(2);

    if (fabs((cb)) > 0.999847695)/*one degree */
    {
        std::cerr << "\nWARNING: Routine Euler_direction2angles is not reliable\n"
        "for small tilt angles. Up to 0.001 deg it should be OK\n"
        "for most applications but you never know";
    }

    if (fabs((cb - 1.)) < FLT_EPSILON)
    {
        alpha = 0.;
        beta = 0.;
    }
    else
    {/*1*/

        aux_beta = acos(cb); /* beta between 0 and PI */


        sb = sin(aux_beta);

        abs_ca = fabs(v(0)) / sb;
        if (fabs((abs_ca - 1.)) < FLT_EPSILON)
            aux_alpha = 0.;
        else
            aux_alpha = acos(abs_ca);

        v_aux(0) = sin(aux_beta) * cos(aux_alpha);
        v_aux(1) = sin(aux_beta) * sin(aux_alpha);
        v_aux(2) = cos(aux_beta);

        error = fabs(dotProduct(v, v_aux) - 1.);
        alpha = aux_alpha;
        beta = aux_beta;

        v_aux(0) = sin(aux_beta) * cos(-1. * aux_alpha);
        v_aux(1) = sin(aux_beta) * sin(-1. * aux_alpha);
        v_aux(2) = cos(aux_beta);
        newerror = fabs(dotProduct(v, v_aux) - 1.);
        if (error > newerror)
        {
            alpha = -1. * aux_alpha;
            beta  = aux_beta;
            error = newerror;
        }

        v_aux(0) = sin(-aux_beta) * cos(-1. * aux_alpha);
        v_aux(1) = sin(-aux_beta) * sin(-1. * aux_alpha);
        v_aux(2) = cos(-aux_beta);
        newerror = fabs(dotProduct(v, v_aux) - 1.);
        if (error > newerror)
        {
            alpha = -1. * aux_alpha;
            beta  = -1. * aux_beta;
            error = newerror;
        }

        v_aux(0) = sin(-aux_beta) * cos(aux_alpha);
        v_aux(1) = sin(-aux_beta) * sin(aux_alpha);
        v_aux(2) = cos(-aux_beta);
        newerror = fabs(dotProduct(v, v_aux) - 1.);

        if (error > newerror)
        {
            alpha = aux_alpha;
            beta  = -1. * aux_beta;
        }
    }/*else 1 end*/
    gamma = 0.;
    beta  = RAD2DEG(beta);
    alpha = RAD2DEG(alpha);
}/*Eulerdirection2angles end*/

Euler_distanceBetweenAngleSets()

template<typename T >

T Euler_distanceBetweenAngleSets	(	T	rot1,
		T	tilt1,
		T	psi1,
		T	rot2,
		T	tilt2,
		T	psi2,
		bool	only_projdir
	)

Average distance between two angle sets. If the only_projdir is set, then only the projection direction is considered.

Definition at line 681 of file geometry.cpp.

{
    static_assert(std::is_floating_point<T>::value, "Only float and double are allowed as template parameters");

    Matrix2D<T> E1, E2;
    Euler_angles2matrix(rot1, tilt1, psi1, E1, false);
    return Euler_distanceBetweenAngleSets_fast(E1,rot2,tilt2,psi2,only_projdir,E2);
}

Euler_distanceBetweenAngleSets_fast()

template<typename T >

T Euler_distanceBetweenAngleSets_fast	(	const Matrix2D< T > &	E1,
		T	rot2,
		T	tilt2,
		T	psi2,
		bool	only_projdir,
		Matrix2D< T > &	E2
	)

Average distance between two angle sets. E1 must contain the Euler matrix corresponding to set1, E2 is used as an auxiliary variable for storing the second Euler matrix.

Definition at line 693 of file geometry.cpp.

{
    static_assert(std::is_floating_point<T>::value, "Only float and double are allowed as template parameters");

    Euler_angles2matrix(rot2, tilt2, psi2, E2, false);
    T aux=MAT_ELEM(E1,2,0)*MAT_ELEM(E2,2,0)+
               MAT_ELEM(E1,2,1)*MAT_ELEM(E2,2,1)+
               MAT_ELEM(E1,2,2)*MAT_ELEM(E2,2,2);
    T axes_dist=std::acos(CLIP(aux, -1, 1));
    if (!only_projdir)
    {
        for (int i = 0; i < 2; i++)
        {
            T aux=MAT_ELEM(E1,i,0)*MAT_ELEM(E2,i,0)+
                       MAT_ELEM(E1,i,1)*MAT_ELEM(E2,i,1)+
                       MAT_ELEM(E1,i,2)*MAT_ELEM(E2,i,2);
            T dist=acos(CLIP(aux, -1, 1));
            axes_dist += dist;
        }
        axes_dist /= 3.0;
    }
    return RAD2DEG(axes_dist);
}

Euler_distanceBetweenMatrices()

double Euler_distanceBetweenMatrices	(	const Matrix2D< double > &	E1,
		const Matrix2D< double > &	E2
	)

Distance between two Euler matrices.

The distance is defined as 1/3*(X1.X2 + Y1.Y2 + Z1.Z2)

Definition at line 671 of file geometry.cpp.

{
    double retval=0;
    FOR_ALL_ELEMENTS_IN_MATRIX2D(E1)
    retval+=MAT_ELEM(E1,i,j)*MAT_ELEM(E2,i,j);
    return retval/3.0;
}

Euler_matrix2angles()

void Euler_matrix2angles	(	const Matrix2D< double > &	A,
		double &	alpha,
		double &	beta,
		double &	gamma,
		bool	homogeneous = false
	)

"Euler" matrix –> angles

This function compute a set of Euler angles which result in an "Euler" matrix as the one given. See Euler_angles2matrix to know more about how this matrix is computed and what each row means. The result angles are in degrees. Alpha, beta and gamma are respectively the first, second and third rotation angles. If the input matrix is not 3x3 then an exception is thrown, the function doesn't check that the Euler matrix is truly representing a coordinate system.

Euler_matrix2angles(Euler, alpha, beta, gamma);

Definition at line 839 of file geometry.cpp.

{
    double abs_sb, sign_sb;

    if (homogeneous)
        if (MAT_XSIZE(A) != 4 || MAT_YSIZE(A) != 4)
            REPORT_ERROR(ERR_MATRIX_SIZE, "Euler_matrix2angles: The Euler matrix is not 4x4");
    else
        if (MAT_XSIZE(A) != 3 || MAT_YSIZE(A) != 3)
            REPORT_ERROR(ERR_MATRIX_SIZE, "Euler_matrix2angles: The Euler matrix is not 3x3");

    abs_sb = sqrt(A(0, 2) * A(0, 2) + A(1, 2) * A(1, 2));
    if (abs_sb > 16*FLT_EPSILON)
    {
        gamma = atan2(A(1, 2), -A(0, 2));
        alpha = atan2(A(2, 1), A(2, 0));
        if (ABS(sin(gamma)) < FLT_EPSILON)
            sign_sb = SGN(-A(0, 2) / cos(gamma));
        // if (sin(alpha)<FLT_EPSILON) sign_sb=SGN(-A(0,2)/cos(gamma));
        // else sign_sb=(sin(alpha)>0) ? SGN(A(2,1)):-SGN(A(2,1));
        else
            sign_sb = (sin(gamma) > 0) ? SGN(A(1, 2)) : -SGN(A(1, 2));
        beta  = atan2(sign_sb * abs_sb, A(2, 2));
    }
    else
    {
        if (SGN(A(2, 2)) > 0)
        {
            // Let's consider the matrix as a rotation around Z
            alpha = 0;
            beta  = 0;
            gamma = atan2(-A(1, 0), A(0, 0));
        }
        else
        {
            alpha = 0;
            beta  = PI;
            gamma = atan2(A(1, 0), -A(0, 0));
        }
    }

    gamma = RAD2DEG(gamma);
    beta  = RAD2DEG(beta);
    alpha = RAD2DEG(alpha);

#ifdef double

    Matrix2D<double> Ap;
    Euler_angles2matrix(alpha, beta, gamma, Ap);
    if (A != Ap)
    {
        std::cout << "---\n";
        std::cout << "Euler_matrix2angles: I have computed angles "
        " which doesn't match with the original matrix\n";
        std::cout << "Original matrix\n" << A;
        std::cout << "Computed angles alpha=" << alpha << " beta=" << beta
        << " gamma=" << gamma << std::endl;
        std::cout << "New matrix\n" << Ap;
        std::cout << "---\n";
    }
#endif

#ifdef DEBUG
    std::cout << "abs_sb " << abs_sb << std::endl;
    std::cout << "A(1,2) " << A(1, 2) << " A(0,2) " << A(0, 2) << " gamma "
    << gamma << std::endl;
    std::cout << "A(2,1) " << A(2, 1) << " A(2,0) " << A(2, 0) << " alpha "
    << alpha << std::endl;
    std::cout << "sign sb " << sign_sb << " A(2,2) " << A(2, 2)
    << " beta " << beta << std::endl;
#endif
}

Euler_mirrorX()

void Euler_mirrorX	(	double	rot,
		double	tilt,
		double	psi,
		double &	newrot,
		double &	newtilt,
		double &	newpsi
	)

Mirror over X axis

Given a set of Euler angles this function returns a new set which define a mirrored (over X axis) version of the former projection.

-----> X               Y
|                       ^
|                       |
|               ======> |
v                       |
Y                        -----> X

The operation performed is

newrot = rot;
newtilt = tilt + 180;
newpsi = 180 - psi;

Euler_mirrorX(rot, tilt, psi, newrot, newtilt, newpsi);

Definition at line 1020 of file geometry.cpp.

{
    newrot  = rot;
    newtilt = tilt + 180;
    newpsi  = 180 - psi;
}

Euler_mirrorXY()

void Euler_mirrorXY	(	double	rot,
		double	tilt,
		double	psi,
		double &	newrot,
		double &	newtilt,
		double &	newpsi
	)

Mirror over X and Y axes

Given a set of Euler angles this function returns a new set which define a mirrored (over X and Y axes at the same time) version of the former projection.

-----> X                       Y
|                               ^
|                               |
|               ======>         |
v                               |
Y                        X<-----

The operation performed is

newrot = rot;
newtilt = tilt;
newpsi = 180 + psi;

Euler_mirrorX(rot, tilt, psi, newrot, newtilt, newpsi);

Definition at line 1029 of file geometry.cpp.

{
    newrot  = rot;
    newtilt = tilt;
    newpsi  = 180 + psi;
}

Euler_mirrorY()

void Euler_mirrorY	(	double	rot,
		double	tilt,
		double	psi,
		double &	newrot,
		double &	newtilt,
		double &	newpsi
	)

Mirror over Y axis

Given a set of Euler angles this function returns a new set which define a mirrored (over Y axis) version of the former projection.

-----> X               X<------
|                              |
|                              |
|               ======>        |
v                              v
Y                             Y

The operation performed is

newrot = rot;
newtilt = tilt + 180;
newpsi = -psi;

Euler_mirrorY(rot, tilt, psi, newrot, newtilt, newpsi);

Definition at line 1011 of file geometry.cpp.

{
    newrot  = rot;
    newtilt = tilt + 180;
    newpsi  = -psi;
}

Euler_rotate() [1/2]

void Euler_rotate	(	const MultidimArray< double > &	V,
		double	rot,
		double	tilt,
		double	psi,
		MultidimArray< double > &	result
	)

Rotate a volume after 3 Euler angles

Input and output volumes cannot be the same one.

Definition at line 1061 of file geometry.cpp.

{
    Matrix2D<double> R;
    Euler_angles2matrix(rot, tilt, psi, R, true);
    applyGeometry(1, result, V, R, xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP);
}

Euler_rotate() [2/2]

void Euler_rotate	(	const MultidimArrayGeneric &	V,
		double	rot,
		double	tilt,
		double	psi,
		MultidimArray< double > &	result
	)

Rotate a volume after 3 Euler angles

Input and output volumes cannot be the same one.

Definition at line 1068 of file geometry.cpp.

{
  Matrix2D<double> R;
  Euler_angles2matrix(rot, tilt, psi, R, true);
#define APPLYGEO(type) applyGeometry(1, result, *((MultidimArray<type> *)V.im), R, xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP);
  SWITCHDATATYPE(V.datatype, APPLYGEO)
#undef APPLYGEO
}

Euler_up_down()

void Euler_up_down	(	double	rot,
		double	tilt,
		double	psi,
		double &	newrot,
		double &	newtilt,
		double &	newpsi
	)

Up-Down projection equivalence

As you know a projection view from a point has got its homologous from its diametrized point in the projection sphere. This function takes a projection defined by its 3 Euler angles and computes an equivalent set of Euler angles from which the view is exactly the same but in the other part of the sphere (if the projection is taken from the bottom then the new projection from the top, and viceversa). The defined projections are exactly the same except for a flip over X axis, ie, an up-down inversion. Exactly the correction performed is:

newrot = rot;
newtilt = tilt + 180;
newpsi = -(180 + psi);

Euler_up_down(rot, tilt, psi, newrot, newtilt, newpsi);

Definition at line 993 of file geometry.cpp.

{
    newrot  = rot;
    newtilt = tilt + 180;
    newpsi  = -(180 + psi);
}

intersection_unit_cube()

double intersection_unit_cube	(	const Matrix1D< double > &	u,
		const Matrix1D< double > &	r
	)

Intersection of a ray with a unit cube

The cube is centered at (0,0,0) and has got unit size length in all directions, i.e., the cube goes from (-0.5, -0.5, -0.5) to (0.5, 0.5, 0.5). The ray is defined by its direction (u) and a passing point (r). See Cube to know how you can intersect any ray with any cube.

Definition at line 1190 of file geometry.cpp.

{
    double t1=0., t2=0., t;
    int found_t = 0;

#define ASSIGN_IF_GOOD_ONE \
    if (fabs(XX(r)+t*XX(u))-XMIPP_EQUAL_ACCURACY<=0.5 && \
        fabs(YY(r)+t*YY(u))-XMIPP_EQUAL_ACCURACY<=0.5 && \
        fabs(ZZ(r)+t*ZZ(u))-XMIPP_EQUAL_ACCURACY<=0.5) {\
        if      (found_t==0) {found_t++; t1=t;} \
        else if (found_t==1) {found_t++; t2=t;} \
    }

    // Intersect with x=0.5 and x=-0.5
    if (XX(u) != 0)
    {
        t = (0.5 - XX(r)) / XX(u);
        ASSIGN_IF_GOOD_ONE;
        t = (-0.5 - XX(r)) / XX(u);
        ASSIGN_IF_GOOD_ONE;
    }

    // Intersect with y=0.5 and y=-0.5
    if (YY(u) != 0 && found_t != 2)
    {
        t = (0.5 - YY(r)) / YY(u);
        ASSIGN_IF_GOOD_ONE;
        t = (-0.5 - YY(r)) / YY(u);
        ASSIGN_IF_GOOD_ONE;
    }

    // Intersect with z=0.5 and z=-0.5
    if (ZZ(u) != 0 && found_t != 2)
    {
        t = (0.5 - ZZ(r)) / ZZ(u);
        ASSIGN_IF_GOOD_ONE;
        t = (-0.5 - ZZ(r)) / ZZ(u);
        ASSIGN_IF_GOOD_ONE;
    }

    if (found_t == 2)
        return fabs(t1 -t2);
    else
        return 0;
}

intersection_unit_cylinder()

double intersection_unit_cylinder	(	const Matrix1D< double > &	u,
		const Matrix1D< double > &	r
	)

Intersection of a ray with a unit cylinder

The cylinder is centered at (0,0,0), has got unit radius on the plane XY, and in Z goes from -h/2 to +h/2. The ray is defined by its direction (u) and a passing point (r). If the ray belongs to the lateral circular wall of the cylinder the length returned is h (h is computed as 1/ZZ(u), this is so because it is supposed that this intersection is computed after a coordinate transformation process from any cylinder to a unit one).

See Cylinder to know how you can intersect any ray with any cylinder.

Definition at line 1151 of file geometry.cpp.

{
    // Intersect with an infinite cylinder of radius=ry
    double A = XX(u) * XX(u) + YY(u) * YY(u);
    double B = XX(r) * XX(u) + YY(r) * YY(u);
    double C = XX(r) * XX(r) + YY(r) * YY(r) - 1;

    double B2_AC = B * B - A * C;
    if (A == 0)
    {
        if (C > 0)
            return 0;       // Parallel ray outside the cylinder
        else
            return 1 / ZZ(u); // return height
    }
    else if (B2_AC < 0)
        return 0;
    B2_AC = sqrt(B2_AC);

    // Points at intersection
    double t1 = (-B - B2_AC) / A;
    double t2 = (-B + B2_AC) / A;
    double z1 = ZZ(r) + t1 * ZZ(u);
    double z2 = ZZ(r) + t2 * ZZ(u);

    // Check position of the intersecting points with respect to
    // the finite cylinder, if any is outside correct it to the
    // right place in the top or bottom of the cylinder
    if (ABS(z1) >= 0.5)
        t1 = (SGN(z1) * 0.5 - ZZ(r)) / ZZ(u);
    if (ABS(z2) >= 0.5)
        t2 = (SGN(z2) * 0.5 - ZZ(r)) / ZZ(u);

    return ABS(t1 -t2);
}

intersection_unit_sphere()

double intersection_unit_sphere	(	const Matrix1D< double > &	u,
		const Matrix1D< double > &	r
	)

Intersection of a ray with a unit sphere

The sphere is centered at (0,0,0) and has got unit radius. The ray is defined by its direction (u) and a passing point (r). The function returns the length of the intersection. If the ray is tangent to the sphere the length is 0. See Ellipsoid to know how you can intersect any ray with any ellipsoid/sphere.

Definition at line 1122 of file geometry.cpp.

{

    // Some useful constants
    double A = XX(u) * XX(u) + YY(u) * YY(u) + ZZ(u) * ZZ(u);
    double B = XX(r) * XX(u) + YY(r) * YY(u) + ZZ(r) * ZZ(u);
    double C = XX(r) * XX(r) + YY(r) * YY(r) + ZZ(r) * ZZ(r) - 1.0;
    double B2_AC = B * B - A * C;

    // A degenerate case?
    if (A == 0)
    {
        if (B == 0)
            return -1; // The ellipsoid doesn't intersect
        return 0;            // The ellipsoid is tangent at t=-C/2B
    }
    if (B2_AC < 0)
        return -1;

    // A normal intersection
    B2_AC = sqrt(B2_AC);
    double t1 = (-B - B2_AC) / A;  // The two parameters within the line for
    double t2 = (-B + B2_AC) / A;  // the solution
    return ABS(t2 -t1);
}

least_squares_line_fit()

void least_squares_line_fit	(	const std::vector< fit_point2D > &	IN_points,
		double &	line_A,
		double &	line_B
	)

Least-squares-fit a line to an arbitrary number of (x,y) points

Plane described as Ax + B = y

Points are defined using the struct

struct fit_point2D
{
    double x;
    double y;
    double w;
};

where w is a weighting factor. Set it to 1 if you do not want to use it

Definition at line 212 of file geometry.cpp.

{

    double  sumx = 0.;
    double  sumy = 0.;
    double  sumxy = 0.;
    double  sumxx = 0.;
    double  sumw = 0.;
    double  W2;
    const fit_point2D * point;

    int n = IN_points.size();

    for (int i = 0; i < n; i++)
    {
        point = &(IN_points[i]);
        W2 = point->w * point->w;
        sumx  += point->x * point->w ;
        sumy  += point->y * point->w ;
        sumxx += point->x * point->x * W2 ;
        sumxy += point->x * point->y * W2 ;
        sumw  += point->w ;
    }
    line_a = (sumx*sumy - sumw*sumxy) / (sumx*sumx - sumw*sumxx) ;
    line_b = (sumy - line_a*sumx) / sumw ;
}

least_squares_plane_fit()

void least_squares_plane_fit	(	FitPoint *	IN_points,
		int	Npoints,
		double &	plane_A,
		double &	plane_B,
		double &	plane_C
	)

Least-squares-fit a plane to an arbitrary number of (x,y,z) points

Plane described as Ax + By + C = z

Points are defined using the struct

struct fit_point
{
    double x;
    double y;
    double z;
    double w;
};

where w is a weighting factor. Set it to 1 if you do not want to use it

Definition at line 101 of file geometry.cpp.

{
    double  D = 0;
    double  E = 0;
    double  F = 0;
    double  G = 0;
    double  H = 0;
    double  I = 0;
    double  J = 0;
    double  K = 0;
    double  L = 0;
    double  W2 = 0;
    double  denom = 0;
    const FitPoint * point;

    for (int i = 0; i < Npoints; i++)
    {
        point = &(IN_points[i]);//Can I copy just the address?
        W2 = point->w * point->w;
        D += point->x * point->x * W2 ;
        E += point->x * point->y * W2 ;
        F += point->x * W2 ;
        G += point->y * point->y * W2 ;
        H += point->y * W2 ;
        I += 1 * W2 ;
        J += point->x * point->z * W2 ;
        K += point->y * point->z * W2 ;
        L += point->z * W2 ;
    }

    denom = F * F * G - 2 * E * F * H + D * H * H + E * E * I - D * G * I;

    // X axis slope
    plane_a = (H * H * J - G * I * J + E * I * K + F * G * L - H * (F * K + E * L)) / denom;
    // Y axis slope
    plane_b = (E * I * J + F * F * K - D * I * K + D * H * L - F * (H * J + E * L)) / denom;
    // Z axis intercept
    plane_c = (F * G * J - E * H * J - E * F * K + D * H * K + E * E * L - D * G * L) / denom;
}

least_squares_plane_fit_All_Points()

void least_squares_plane_fit_All_Points	(	const MultidimArray< double > &	Image,
		double &	plane_A,
		double &	plane_B,
		double &	plane_C
	)

Least-squares-fit a plane to an image

Performs the same computation as least_squares_plane_fit but using a complete image instead of only a set of points.

Definition at line 154 of file geometry.cpp.

{
    double  D = 0;
    double  E = 0;
    double  F = 0;
    double  G = 0;
    double  H = 0;
    double  I = 0;
    double  J = 0;
    double  K = 0;
    double  L = 0;
    double  denom = 0;
    int nIterations=0;
    double sumjValues=0.0;
    for (int j=STARTINGX(Image); j<=FINISHINGX(Image); ++j)
    {
        D += j*j;
        sumjValues += j;
        I += 1;
    }
    F = sumjValues;

    double *ref;
    double sumElements;
    for (int i=STARTINGY(Image); i<=FINISHINGY(Image); ++i)
    {
        nIterations++;
        ref = &A2D_ELEM(Image, i, STARTINGX(Image));
        sumElements = 0.0;
        for (int j=STARTINGX(Image); j<=FINISHINGX(Image); ++j)
        {
            J += j*(*ref);
            sumElements += (*ref);
            ref++;
        }
        K += i*sumElements;
        L += sumElements;
        E += sumjValues*i;
        G += i*i*I;
        H += i*I;
    }
    D = D*nIterations;
    F = F*nIterations;
    I = I*nIterations;

    denom = F * F * G - 2 * E * F * H + D * H * H + E * E * I - D * G * I;

    // X axis slope
    plane_a = (H * H * J - G * I * J + E * I * K + F * G * L - H * (F * K + E * L)) / denom;
    // Y axis slope
    plane_b = (E * I * J + F * F * K - D * I * K + D * H * L - F * (H * J + E * L)) / denom;
    // Z axis intercept
    plane_c = (F * G * J - E * H * J - E * F * K + D * H * K + E * E * L - D * G * L) / denom;
}

line_plane_intersection()

int line_plane_intersection	(	const Matrix1D< double >	normal_plane,
		const Matrix1D< double >	vector_line,
		Matrix1D< double > &	intersection_point,
		const Matrix1D< double >	point_line,
		double	point_plane_at_x_y_zero = 0.
	)

Line Plane Intersection

Let ax+by+cz+D=0 be the equation of your plane (if your plane is defined by a normal vector N + one point M, then (a,b,c) are the coordinates of the normal N, and d is calculated by using the coordinates of M in the above equation).

Let your line be defined by one point P(d,e,f) and a vector V(u,v,w), the points on your line are those which verify

x = d + lu y = e + lv z = f + lw

where l takes all real values.

for this point to be on the plane, you have to have

ax + by + cz + D = 0, so,

a(d + lu) + b(e + lv) + c(f + lw) + D = 0

that is

l(au + bv + cw) = -(ad + be + cf + D)

note that, if au + bv + cw = 0, then your line is either in the plane, or parallel to it... otherwise you get the value of l, and the intersection has coordinates:

x = d + lu y = e + lv z = f + lw

where

l = -(ad + be + cf + D) / (au + bv + cw)

a = XX(normal_plane); b = YY(normal_plane); c = ZZ(normal_plane); D = intersection_point;

d = XX(point_line) e = YY(point_line) f = ZZ(point_line)

u = XX(vector_line); v = YY(vector_line); w = ZZ(vector_line);

XX(intersection_point) = x; YY(intersection_point) = y; ZZ(intersection_point) = z;

return 0 if successful return -1 if line parallel to plane return +1 if line in the plane

Definition at line 588 of file geometry.cpp.

{
    double l;
    intersection_point.resize(3);
    // if au+bv+cw=0, then your line is either in the plane, or
    // parallel to it
    if (ABS(dotProduct(normal_plane, vector_line)) < XMIPP_EQUAL_ACCURACY)
    {
        intersection_point = point_line + vector_line;
        if (ABS(dotProduct(intersection_point, normal_plane) +
                point_plane_at_x_y_zero) < XMIPP_EQUAL_ACCURACY)
            return(1);
        else
            return(-1);

    }
    //compute intersection
    l = -1.0 * dotProduct(point_line, normal_plane) +
        point_plane_at_x_y_zero;
    l /=  dotProduct(normal_plane, vector_line);

    intersection_point = point_line + l * vector_line;
    return(0);
}

point_inside_polygon()

bool point_inside_polygon	(	const std::vector< Matrix1D< double > > &	polygon,
		const Matrix1D< double > &	point
	)

Point inside polygon

Given a polygon described by a list of points (the last one and the first one ust be the same), determine whether another point is inside the polygon or not.

Definition at line 418 of file geometry.cpp.

{
    size_t i, j;
    bool retval = false;
    for (i = 0, j = polygon.size() - 1; i < polygon.size(); j = i++)
    {
        if ((((YY(polygon[i]) <= YY(point)) && (YY(point) < YY(polygon[j]))) ||
             ((YY(polygon[j]) <= YY(point)) && (YY(point) < YY(polygon[i])))) &&
            (XX(point) < (XX(polygon[j]) - XX(polygon[i])) *
             (YY(point) - YY(polygon[i])) /
             (YY(polygon[j]) - YY(polygon[i])) + XX(polygon[i])))
            retval = !retval;
    }
    return retval;
}

point_line_distance_3D()

double point_line_distance_3D	(	const Matrix1D< double > &	p,
		const Matrix1D< double > &	a,
		const Matrix1D< double > &	v
	)

Point to line distance in 3D

Let a line in 3-D be specified by the point a and the vector v, this fuction returns the minimum distance of this line to the point p.

Definition at line 85 of file geometry.cpp.

{
    Matrix1D<double> p_a(3);

    V3_MINUS_V3(p_a, p, a);
    return (vectorProduct(p_a, v).module() / v.module());
}

point_plane_distance_3D()

double point_plane_distance_3D ( const Matrix1D< double > &  p, const Matrix1D< double > &  a, const Matrix1D< double > &  v  )

inline

Point to plane distance in 3D

Let a plane in 3-D be specified by the point a and the perpendicular vector v, this fuction returns the minimum distance of this plane to the point p.

Definition at line 171 of file geometry.h.

{
    static Matrix1D< double > p_a(3);

    V3_MINUS_V3(p_a, p, a);
    return (dotProduct(p_a, v) / v.module());
}

rectangle_enclosing()

void rectangle_enclosing	(	const Matrix1D< double > &	v0,
		const Matrix1D< double > &	vF,
		const Matrix2D< double > &	V,
		Matrix1D< double > &	corner1,
		Matrix1D< double > &	corner2
	)

Rectangle which encloses a deformed rectangle

Given a rectangle characterized by the top-left corner and the right-bottom corner, and given a matrix after which the rectangle is deformed. Which is the minimum rectangle which encloses the preceding one? This function is useful for stablishing for loops which will cover for sure the deformed rectangle. All vectors are supposed to be 2x1 and the deformation matrix is 2x2. The corner (x0,y0) goes to V*(x0,y0)' and (xF,yF) to V*(xf,yF)'. After that you can make a loop from corner1 to corner2.

The v0 and vF vectors can be reused as outputs.

Definition at line 328 of file geometry.cpp.

{
    SPEED_UP_temps01;
    Matrix1D<double> v(2);
    corner1.resize(2);
    corner2.resize(2);

    // Store values for reusing input as output vectors
    double XX_v0 = XX(v0);
    double YY_v0 = YY(v0);
    double XX_vF = XX(vF);
    double YY_vF = YY(vF);

    VECTOR_R2(v, XX_v0, YY_v0);
    M2x2_BY_V2x1(v, V, v);
    XX(corner1) = XX(v);
    XX(corner2) = XX(v);
    YY(corner1) = YY(v);
    YY(corner2) = YY(v);

#define DEFORM_AND_CHOOSE_CORNERS2D \
    M2x2_BY_V2x1(v,V,v); \
    XX(corner1)=XMIPP_MIN(XX(corner1),XX(v)); \
    XX(corner2)=XMIPP_MAX(XX(corner2),XX(v)); \
    YY(corner1)=XMIPP_MIN(YY(corner1),YY(v)); \
    YY(corner2)=XMIPP_MAX(YY(corner2),YY(v));

    VECTOR_R2(v, XX_vF, YY_v0);
    DEFORM_AND_CHOOSE_CORNERS2D;
    VECTOR_R2(v, XX_v0, YY_vF);
    DEFORM_AND_CHOOSE_CORNERS2D;
    VECTOR_R2(v, XX_vF, YY_vF);
    DEFORM_AND_CHOOSE_CORNERS2D;
}

spherical_distance()

double spherical_distance	(	const Matrix1D< double > &	r1,
		const Matrix1D< double > &	r2
	)

Spherical distance

This function returns the distance over a sphere, not the straight line but the line which goes from one point to the other going over the surface of a sphere, supposing that both points lie on the same sphere.

Definition at line 73 of file geometry.cpp.

{
    double r1r2 = XX(r1) * XX(r2) + YY(r1) * YY(r2) + ZZ(r1) * ZZ(r2);
    double R1 = sqrt(XX(r1) * XX(r1) + YY(r1) * YY(r1) + ZZ(r1) * ZZ(r1));
    double R2 = sqrt(XX(r2) * XX(r2) + YY(r2) * YY(r2) + ZZ(r2) * ZZ(r2));
    double argument = r1r2 / (R1 * R2);
    double ang = acos(CLIP(argument,-1.,1.));
    double retVal = ang*R1;
    return retVal;
}

triangle_area()

double triangle_area	(	double	x1,
		double	y1,
		double	x2,
		double	y2,
		double	x3,
		double	y3
	)

Area of a triangle Given the coordinates of three points this function calculates the area of the triangle

Definition at line 511 of file geometry.cpp.

{
    double trigarea = ((x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1))/2;
    return (trigarea > 0) ? trigarea : -trigarea;
}

Uproject_to_plane() [1/3]

void Uproject_to_plane	(	const Matrix1D< double > &	point,
		const Matrix1D< double > &	direction,
		double	distance,
		Matrix1D< double > &	result
	)

Project a point to a plane (direction vector, distance)

Given the coordinates for a vector in R3, and a plane (defined by its direction vector and the minimum distance from the plane to the coordinate system origin). This function computes the perpendicular projection from the vector to the plane. This function has tried to be optimized in speed as it is used in core routines within huge loops, this is why the result is given as an argument, and why no check about the dimensionality of the vectors is performed. The routine performs faster if the result vector is already in R3.

The following example projects the point P=(1,1,1) to the XY-plane storing the result in Pp (belonging to R3), the result is obviously Pp=(1,1,0).

Matrix1D< double > Z = vectorR3(0, 0, 1), P = vector_R3(1, 1, 1), Pp(3);
Uproject_to_plane(P,Z,0,Pp);
std::cout << "After projecting: Pp=" << Pp.transpose() << std::endl;

The starting U in the function name stands for the fact that the plane and the point are in the same reference system (called by default Universal).

The result and point vectors can be the same one.

Definition at line 39 of file geometry.cpp.

{

    if (result.size() != 3)
        result.resize(3);
    double xx = distance - (XX(point) * XX(direction) + YY(point) * YY(direction) +
                            ZZ(point) * ZZ(direction));
    XX(result) = XX(point) + xx * XX(direction);
    YY(result) = YY(point) + xx * YY(direction);
    ZZ(result) = ZZ(point) + xx * ZZ(direction);
}

Uproject_to_plane() [2/3]

void Uproject_to_plane	(	const Matrix1D< double > &	r,
		double	rot,
		double	tilt,
		double	psi,
		Matrix1D< double > &	result
	)

Project a vector to a plane (Euler angles)

These planes are restricted to have 0 distance to the universal coordinate system. In this special case, a plane can be defined by 3 Euler angles (this is specially suited for projections, where the projection plane is defined by its 3 Euler angles). Then, the coordinate system associated to the 3 Euler angles (let's call its vectors X',Y',Z') defines a projection plane where Z' is the direction vector, and X'Y' are in-plane vectors. Actually, X' coincides with the X vector in the Matrix2D definition and Y' with the Y vector.

The resulting vector is in R3, and the function has been optimized for speed, so the result is passed as a parameter. This function is based in the one which projects a point given the Euler matrix of the projection plane. If you are to project several points to the same plane, you should better use the project to plane function where you give the Euler matrix.

The following example projects the point P=(1,1,1) to the XY-plane storing the result in Pp (belonging to R3), the result is obviously Pp=(1,1,0).

Matrix1D< double > P = vectorR3(1, 1, 1), Pp(3);
Uproject_to_plane(P, 0, 0, 0, Pp);
std::cout << "After projecting: Pp=" << Pp.transpose() << std::endl;

The starting U in the function name stands for the fact that the plane, and the point are in the same reference system (called by default Universal).

The result and point vectors can be the same one.

Definition at line 54 of file geometry.cpp.

{
    Matrix2D<double> euler(3, 3);
    Euler_angles2matrix(rot, tilt, psi, euler);
    Uproject_to_plane(r, euler, result);
}

Uproject_to_plane() [3/3]

void Uproject_to_plane	(	const Matrix1D< double > &	r,
		const Matrix2D< double > &	euler,
		Matrix1D< double > &	result
	)

Project a vector to a plane (Euler matrix)

These planes are restricted to have 0 distance to the universal coordinate system. In this special case, a plane can be defined by 3 Euler angles (this is specially suited for projections, where the projection plane is defined by its 3 Euler angles). Then, the coordinate system associated to the 3 Euler angles (let's call its vectors X',Y',Z') defines a projection plane where Z' is the direction vector, and X'Y' are in-plane vectors. Actually, X' coincides with the X vector in the Matrix2D definition and Y' with the Y vector.

The resulting vector is in R3, and the function has been optimized for speed, if the result vector is already 3 dimensional when entering the function, no resize is performed; and the result is passed as a parameter. If you are to project several points to the same plane, you should better use the project to plane function where you give the Euler matrix.

The following example projects the point P=(1,1,1) to the XY-plane storing the result in Pp (belonging to R3), the result is obviously Pp=(1,1,0).

Matrix1D< double > P = vectorR3(1, 1, 1), Pp(3);
Matrix2D< double > euler;
Euler_angles2matrix(0, 0, 0, euler);
Uproject_to_plane(P, euler, Pp);
std::cout << "After projecting: Pp=" << Pp.transpose() << std::endl;

The starting U in the function name stands for the fact that the plane, and the point are in the same reference system (called by default Universal).

The result and point vectors can be the same one.

Definition at line 63 of file geometry.cpp.

{
    SPEED_UP_temps012;
    if (VEC_XSIZE(result) != 3)
        result.resize(3);
    M3x3_BY_V3x1(result, euler, r);
}