Geometrical transformations

Collaboration diagram for Geometrical transformations:

Namespaces
	applyGeometryImpl

Functions
void	geo2TransformationMatrix (const MDRow &imageHeader, Matrix2D< double > &A, bool only_apply_shifts=false)
bool	getLoopRange (double value, double min, double max, double delta, int loopLimit, int &minIter, int &maxIter)
void	string2TransformationMatrix (const String &matrixStr, Matrix2D< double > &matrix, size_t dim=4)
template<typename T >
void	transformationMatrix2Parameters2D (const Matrix2D< T > &A, bool &flip, T &scale, T &shiftX, T &shiftY, T &psi)
void	transformationMatrix2Parameters3D (const Matrix2D< double > &A, bool &flip, double &scale, double &shiftX, double &shiftY, double &shiftZ, double &rot, double &tilt, double &psi)
void	transformationMatrix2Geo (const Matrix2D< double > &A, MDRow &imageGeo)
template<typename T >
void	rotation2DMatrix (T ang, Matrix2D< T > &m, bool homogeneous=true)
template<typename T >
void	translation2DMatrix (const Matrix1D< T > &translation, Matrix2D< T > &resMatrix, bool inverse=false)
void	rotation3DMatrix (double ang, char axis, Matrix2D< double > &m, bool homogeneous=true)
void	rotation3DMatrix (double ang, const Matrix1D< double > &axis, Matrix2D< double > &m, bool homogeneous=true)
void	alignWithZ (const Matrix1D< double > &axis, Matrix2D< double > &m, bool homogeneous=true)
template<typename T >
void	translation3DMatrix (const Matrix1D< T > &translation, Matrix2D< T > &resMatrix, bool inverse=false)
void	scale3DMatrix (const Matrix1D< double > &sc, Matrix2D< double > &m, bool homogeneous=true)
void	rotation3DMatrixFromIcoOrientations (const char *icoFrom, const char *icoTo, Matrix2D< double > &R)
template<typename T1 , typename T , typename T2 >
void	applyGeometry (int SplineDegree, MultidimArray< T > &__restrict__ V2, const MultidimArray< T1 > &__restrict__ V1, const Matrix2D< T2 > &At, bool inv, bool wrap, T outside=T(0), MultidimArray< double > *BcoeffsPtr=NULL)
template<typename T >
void	applyGeometry (int SplineDegree, MultidimArray< T > &V2, const MultidimArrayGeneric &V1, const Matrix2D< double > &A, bool inv, bool wrap, T outside=0)
template<typename T >
void	selfApplyGeometry (int SplineDegree, MultidimArray< T > &V1, const Matrix2D< double > &A, bool inv, bool wrap, T outside=0)
template<>
void	applyGeometry (int SplineDegree, MultidimArray< std::complex< double > > &V2, const MultidimArray< std::complex< double > > &V1, const Matrix2D< double > &A, bool inv, bool wrap, std::complex< double > outside, MultidimArray< double > *BcoeffsPtr)
template<>
void	selfApplyGeometry (int SplineDegree, MultidimArray< std::complex< double > > &V1, const Matrix2D< double > &A, bool inv, bool wrap, std::complex< double > outside)
void	applyGeometry (int SplineDegree, MultidimArrayGeneric &V2, const MultidimArrayGeneric &V1, const Matrix2D< double > &A, bool inv, bool wrap, double outside)
template<typename T >
void	produceSplineCoefficients (int SplineDegree, MultidimArray< double > &coeffs, const MultidimArray< T > &V1)
void	produceSplineCoefficients (int SplineDegree, MultidimArray< double > &coeffs, const MultidimArray< std::complex< double > > &V1)
template<typename T >
void	produceImageFromSplineCoefficients (int SplineDegree, MultidimArray< T > &img, const MultidimArray< double > &coeffs)
template<typename T >
void	rotate (int SplineDegree, MultidimArray< T > &V2, const MultidimArray< T > &V1, double ang, char axis='Z', bool wrap=xmipp_transformation::DONT_WRAP, T outside=0)
template<typename T >
void	selfRotate (int SplineDegree, MultidimArray< T > &V1, double ang, char axis='Z', bool wrap=xmipp_transformation::DONT_WRAP, T outside=0)
template<typename T >
void	translate (int SplineDegree, MultidimArray< T > &V2, const MultidimArray< T > &V1, const Matrix1D< double > &v, bool wrap=xmipp_transformation::WRAP, T outside=0)
template<typename T >
void	selfTranslate (int SplineDegree, MultidimArray< T > &V1, const Matrix1D< double > &v, bool wrap=xmipp_transformation::WRAP, T outside=0)
template<typename T >
void	translateCenterOfMassToCenter (int SplineDegree, MultidimArray< T > &V2, const MultidimArray< T > &V1, bool wrap=xmipp_transformation::WRAP)
template<typename T >
void	selfTranslateCenterOfMassToCenter (int SplineDegree, MultidimArray< T > &V1, bool wrap=xmipp_transformation::WRAP)
template<typename T1 , typename T >
void	scaleToSize (int SplineDegree, MultidimArray< T > &V2, const MultidimArray< T1 > &V1, size_t Xdim, size_t Ydim, size_t Zdim=1)
template<typename T >
void	scaleToSize (int SplineDegree, MultidimArray< T > &V2, const MultidimArrayGeneric &V1, int Xdim, int Ydim, int Zdim=1)
template<typename T >
void	scaleToSize (int SplineDegree, MultidimArrayGeneric &V2, const MultidimArray< T > &V1, int Xdim, int Ydim, int Zdim=1)
void	scaleToSize (int SplineDegree, MultidimArrayGeneric &V2, const MultidimArrayGeneric &V1, int Xdim, int Ydim, int Zdim=1)
template<typename T >
void	selfScaleToSize (int SplineDegree, MultidimArray< T > &V1, int Xdim, int Ydim, int Zdim=1)
void	selfScaleToSize (int SplineDegree, MultidimArrayGeneric &V1, int Xdim, int Ydim, int Zdim=1)
void	scaleToSize (int SplineDegree, MultidimArray< std::complex< double > > &V2, const MultidimArray< std::complex< double > > &V1, int Xdim, int Ydim, int Zdim=1)
void	selfScaleToSize (int SplineDegree, MultidimArray< std::complex< double > > &V1, int Xdim, int Ydim, int Zdim=1)
template<typename T >
void	reduceBSpline (int SplineDegree, MultidimArray< double > &V2, const MultidimArray< T > &V1)
template<typename T >
void	expandBSpline (int SplineDegree, MultidimArray< double > &V2, const MultidimArray< T > &V1)
template<typename T >
void	pyramidReduce (int SplineDegree, MultidimArray< T > &V2, const MultidimArray< T > &V1, int levels=1)
template<typename T >
void	selfPyramidReduce (int SplineDegree, MultidimArray< T > &V1, int levels=1)
void	selfPyramidReduce (int SplineDegree, MultidimArrayGeneric &V1, int levels=1)
template<typename T >
void	pyramidExpand (int SplineDegree, MultidimArray< T > &V2, const MultidimArray< T > &V1, int levels=1)
void	pyramidExpand (int SplineDegree, MultidimArrayGeneric &V2, const MultidimArrayGeneric &V1, int levels=1)
void	pyramidReduce (int SplineDegree, MultidimArrayGeneric &V2, const MultidimArrayGeneric &V1, int levels=1)
template<typename T >
void	selfPyramidExpand (int SplineDegree, MultidimArray< T > &V1, int levels=1)
void	selfPyramidExpand (int SplineDegree, MultidimArrayGeneric &V1, int levels=1)
template<typename T >
void	radialAverage (const MultidimArray< T > &m, Matrix1D< int > &center_of_rot, MultidimArray< T > &radial_mean, MultidimArray< int > &radial_count, const bool &rounding=false)
template<typename T >
void	radialAveragePrecomputeDistance (const MultidimArray< T > &m, Matrix1D< int > &center_of_rot, MultidimArray< int > &distance, int &dim, const bool &rounding=false)
template<typename T >
void	fastRadialAverage (const MultidimArray< T > &m, const MultidimArray< int > &distance, int dim, MultidimArray< T > &radial_mean, MultidimArray< int > &radial_count)
template<typename T >
void	radialAverageAxis (const MultidimArray< T > &in, char axis, MultidimArray< double > &out)
template<typename T >
void	radialAverageNonCubic (const MultidimArray< T > &m, Matrix1D< int > &center_of_rot, MultidimArray< T > &radial_mean, MultidimArray< int > &radial_count, bool rounding=false)
void	radiallySymmetrize (const MultidimArray< double > &img, MultidimArray< double > &radialImg)
double	interpolatedElementBSplineDiffX (MultidimArray< double > &vol, double x, double y, double z, int SplineDegree=3)
double	interpolatedElementBSplineDiffY (MultidimArray< double > &vol, double x, double y, double z, int SplineDegree=3)
double	interpolatedElementBSplineDiffZ (MultidimArray< double > &vol, double x, double y, double z, int SplineDegree=3)

Detailed Description

Function Documentation

alignWithZ()

void alignWithZ	(	const Matrix1D< double > &	axis,
		Matrix2D< double > &	m,
		bool	homogeneous = true
	)

Matrix which transforms the given axis into Z

A geometrical transformation matrix (4x4) is returned such that the given axis is rotated until it is aligned with the Z axis. This is very useful in order to produce rotational matrices, for instance, around any axis.

Matrix2D< double > A = alignWithZ(axis);
return A.transpose() * rotation3DMatrix(ang, 'Z') * A;

The returned matrix is such that A*axis=Z, where Z and axis are column vectors.

Definition at line 471 of file transformations.cpp.

{
    if (axis.size() != 3)
        REPORT_ERROR(ERR_MATRIX_SIZE, "alignWithZ: Axis is not in R3");
    if (homogeneous)
    {
        result.initZeros(4,4);
        dMij(result,3, 3) = 1;
    }
    else
        result.initZeros(3,3);
    Matrix1D<double>  Axis(axis);
    Axis.selfNormalize();

    // Compute length of the projection on YZ plane
    double proj_mod = sqrt(YY(Axis) * YY(Axis) + ZZ(Axis) * ZZ(Axis));
    if (proj_mod > XMIPP_EQUAL_ACCURACY)
    {   // proj_mod!=0
        // Build Matrix result, which makes the turning axis coincident with Z
        dMij(result,0, 0) = proj_mod;
        dMij(result,0, 1) = -XX(Axis) * YY(Axis) / proj_mod;
        dMij(result,0, 2) = -XX(Axis) * ZZ(Axis) / proj_mod;
        dMij(result,1, 0) = 0;
        dMij(result,1, 1) = ZZ(Axis) / proj_mod;
        dMij(result,1, 2) = -YY(Axis) / proj_mod;
        dMij(result,2, 0) = XX(Axis);
        dMij(result,2, 1) = YY(Axis);
        dMij(result,2, 2) = ZZ(Axis);
    }
    else
    {
        // I know that the Axis is the X axis, EITHER POSITIVE OR NEGATIVE!!
        dMij(result,0, 0) = 0;
        dMij(result,0, 1) = 0;
        dMij(result,0, 2) = (XX(Axis) > 0)? -1 : 1;
        dMij(result,1, 0) = 0;
        dMij(result,1, 1) = 1;
        dMij(result,1, 2) = 0;
        dMij(result,2, 0) = (XX(Axis) > 0)? 1 : -1;
        dMij(result,2, 1) = 0;
        dMij(result,2, 2) = 0;
    }
}

applyGeometry() [1/4]

template<typename T1 , typename T , typename T2 >

void applyGeometry	(	int	SplineDegree,
		MultidimArray< T > &__restrict__	V2,
		const MultidimArray< T1 > &__restrict__	V1,
		const Matrix2D< T2 > &	At,
		bool	inv,
		bool	wrap,
		T	outside = T(0),
		MultidimArray< double > *	BcoeffsPtr = NULL
	)

Applies a geometrical transformation.

Any geometrical transformation defined by the matrix A (double (4x4)!! ie, in homogeneous R3 coordinates) is applied to the volume V1. The result is stored in V2 (it cannot be the same as the input volume). An exception is thrown if the transformation matrix is not 4x4.

Structure of the transformation matrix: It should have the following components

r11 r12 r13 x r21 r22 r23 y r31 r32 r33 z 0 0 0 1

where (x,y,z) is the translation desired, and Rij are the components of the rotation matrix R. If you want to apply a scaling factor to the transformation, then multiply r11, r22 and r33 by it.

The result volume (with ndim=1) is resized to the same dimensions as V1 if V2 is empty (0x0) at the beginning, if it is not, ie, if V2 has got some size then only those values in the volume are filled, this is very useful for resizing the volume, then you manually resize the output volume to the desired size and then call this routine.

The relationship between the output coordinates and the input ones are

out = A * in
(x, y, z) = A * (x', y', z')

This function works independently from the logical indexing of each matrix, it sets the logical center and the physical center of the image and work with these 2 coordinate spaces. At the end the original logical indexing of each matrix is kept.

The procedure followed goes from coordinates in the output volume to the ones in the input one, so the inverse of the A matrix is needed. There is a flag telling if the given matrix is already the inverse one or the normal one. If it is the normal one internally the matrix is inversed. If you are to do many "rotations" then some time is spent in inverting the matrix. Normally the matrix is the normal one.

There is something else to tell about the geometrical tranformation. The value of the voxel in the output volume is computed via bilinear interpolation in the input volume. If any of the voxels participating in the interpolation falls outside the input volume, then automatically the corresponding output voxel is set to 0, unless that the wrap flag has been set to 1. In this case if the voxel falls out by the right hand then it is "wrapped" and the corresponding voxel in the left hand is used. The same is appliable to top-bottom. Usually wrap mode is off. Wrap mode is interesting for translations but not for rotations, for example.

The inverse mode and wrapping mode should be taken by default by the routine, g++ seems to have problems with template functions outside a class with default parameters. So, I'm sorry, you will have to put them always. The usual combination is

applyGeometry(..., IS_NOT_INV, DONT_WRAP).

Although you can also use the constants IS_INV, or WRAP.

Matrix2D< double > A(4,4);
A.initIdentity;
applyGeometry(V2, A, V1);

Definition at line 531 of file transformations.h.

{
#ifndef RELEASE_MODE
    if (&V1 == (MultidimArray<T1>*)&V2)
        REPORT_ERROR(ERR_VALUE_INCORRECT,"ApplyGeometry: Input array cannot be the same as output array");

    if ( V1.getDim()==2 && ((MAT_XSIZE(At) != 3) || (MAT_YSIZE(At) != 3)) )
        REPORT_ERROR(ERR_MATRIX_SIZE,"ApplyGeometry: 2D transformation matrix is not 3x3");

    if ( V1.getDim()==3 && ((MAT_XSIZE(At) != 4) || (MAT_YSIZE(At) != 4)) )
        REPORT_ERROR(ERR_MATRIX_SIZE,"ApplyGeometry: 3D transformation matrix is not 4x4");
#endif

    if (At.isIdentity() && ( XSIZE(V2) == 0 || SAME_SHAPE3D(V1,V2) ) )
    {
        typeCast(V1,V2);
        return;
    }

    if (XSIZE(V1) == 0)
    {
        V2.clear();
        return;
    }

    MultidimArray<double> Bcoeffs;
    MultidimArray<double> *BcoeffsToUse=NULL;
    Matrix2D<double> A, Ainv;
    typeCast(At, A);
    const Matrix2D<double> * Aptr=&A;
    if (!inv)
    {
        Ainv = A.inv();
        Aptr=&Ainv;
    }
    const Matrix2D<double> &Aref=*Aptr;

    // For scalings the output matrix is resized outside to the final
    // size instead of being resized inside the routine with the
    // same size as the input matrix
    if (XSIZE(V2) == 0)
        V2.resizeNoCopy(V1);

    if (outside != 0.)
    {
        // Initialize output matrix with value=outside
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(V2)
        DIRECT_MULTIDIM_ELEM(V2, n) = outside;
    }
    else
        V2.initZeros();

    if (V1.getDim() == 2)
    {
        // this version should be slightly more optimized than the general one bellow
        if (1 == SplineDegree) {
            if (wrap) {
                applyGeometryImpl::applyGeometry2DDegree1<T1, T, true>(V2, V1, Aref);
            } else {
                applyGeometryImpl::applyGeometry2DDegree1<T1, T, false>(V2, V1, Aref);
            }
            return;
        }

        // 2D transformation
        double Aref00=MAT_ELEM(Aref,0,0);
        double Aref10=MAT_ELEM(Aref,1,0);

        // Find center and limits of image
        double cen_y  = (int)(YSIZE(V2) / 2);
        double cen_x  = (int)(XSIZE(V2) / 2);
        double cen_yp = (int)(YSIZE(V1) / 2);
        double cen_xp = (int)(XSIZE(V1) / 2);
        double minxp  = -cen_xp;
        double minyp  = -cen_yp;
        double minxpp = minxp-XMIPP_EQUAL_ACCURACY;
        double minypp = minyp-XMIPP_EQUAL_ACCURACY;
        double maxxp  = XSIZE(V1) - cen_xp - 1;
        double maxyp  = YSIZE(V1) - cen_yp - 1;
        double maxxpp = maxxp+XMIPP_EQUAL_ACCURACY;
        double maxypp = maxyp+XMIPP_EQUAL_ACCURACY;
        size_t Xdim   = XSIZE(V1);
        size_t Ydim   = YSIZE(V1);

        if (SplineDegree > 1)
        {
            // Build the B-spline coefficients
            if (BcoeffsPtr!=NULL)
                BcoeffsToUse=BcoeffsPtr;
            else
            {
                produceSplineCoefficients(SplineDegree, Bcoeffs, V1); //Bcoeffs is a single image
                BcoeffsToUse = &Bcoeffs;
            }
            STARTINGX(*BcoeffsToUse) = (int) minxp;
            STARTINGY(*BcoeffsToUse) = (int) minyp;
        }

        // Now we go from the output image to the input image, ie, for any pixel
        // in the output image we calculate which are the corresponding ones in
        // the original image, make an interpolation with them and put this value
        // at the output pixel

#ifdef DEBUG_APPLYGEO
        std::cout << "A\n" << Aref << std::endl
        << "(cen_x ,cen_y )=(" << cen_x  << "," << cen_y  << ")\n"
        << "(cen_xp,cen_yp)=(" << cen_xp << "," << cen_yp << ")\n"
        << "(min_xp,min_yp)=(" << minxp  << "," << minyp  << ")\n"
        << "(max_xp,max_yp)=(" << maxxp  << "," << maxyp  << ")\n";
#endif
        // Calculate position of the beginning of the row in the output image
        double x = -cen_x;
        double y = -cen_y;
        for (size_t i = 0; i < YSIZE(V2); i++)
        {
            // Calculate this position in the input image according to the
            // geometrical transformation
            // they are related by
            // coords_output(=x,y) = A * coords_input (=xp,yp)
            double xp = x * MAT_ELEM(Aref, 0, 0) + y * MAT_ELEM(Aref, 0, 1) + MAT_ELEM(Aref, 0, 2);
            double yp = x * MAT_ELEM(Aref, 1, 0) + y * MAT_ELEM(Aref, 1, 1) + MAT_ELEM(Aref, 1, 2);

            // Inner loop boundaries.
            int globalMin=0, globalMax=XSIZE(V2);

            if (!wrap)
            {
                // First and last iteration with valid values for x and y coordinates.
                int minX, maxX, minY, maxY;

                // Compute valid iterations in x and y coordinates. If one of them is always out
                // of boundaries then the inner loop is not executed this iteration.
                if (!getLoopRange( xp, minxpp, maxxpp, Aref00, XSIZE(V2), minX, maxX) ||
                    !getLoopRange( yp, minypp, maxypp, Aref10, XSIZE(V2), minY, maxY))
                {
                    y++;
                    continue;
                }
                else
                {
                    // Compute initial iteration.
                    globalMin = minX;
                    if (minX < minY)
                    {
                        globalMin = minY;
                    }

                    // Compute last iteration.
                    globalMax = maxX;
                    if (maxX > maxY)
                    {
                        globalMax = maxY;
                    }
                    globalMax++;

                    // Check max iteration is not higher than image.
                    if ((globalMax >= 0) && ((size_t)globalMax > XSIZE(V2)))
                    {
                        globalMax = XSIZE(V2);
                    }

                    xp += globalMin*Aref00;
                    yp += globalMin*Aref10;
                }
            }

            // Loop over j is splitted according to wrap (wrap==true is not
            // vectorizable) and also according to SplineDegree value
            // (I have not fully analyzed vector dependences for 
            // SplineDegree==3 and else branch)
            if (wrap) 
            {
                // This is original implementation
                for (int j=globalMin; j<globalMax ;j++)
                {
#ifdef DEBUG_APPLYGEO

                    std::cout << "Computing (" << i << "," << j << ")\n";
                    std::cout << "   (y, x) =(" << y << "," << x << ")\n"
                    << "   before wrapping (y',x')=(" << yp << "," << xp << ") "
                    << std::endl;
#endif
                    bool x_isOut = XMIPP_RANGE_OUTSIDE_FAST(xp, minxpp, maxxpp);
                    bool y_isOut = XMIPP_RANGE_OUTSIDE_FAST(yp, minypp, maxypp);

                    if (x_isOut)
                    {
                        xp = realWRAP(xp, minxp - 0.5, maxxp + 0.5);
                    }

                    if (y_isOut)
                    {
                        yp = realWRAP(yp, minyp - 0.5, maxyp + 0.5);
                    }

#ifdef DEBUG_APPLYGEO

                    std::cout << "   after wrapping (y',x')=(" << yp << "," << xp << ") "
                    << std::endl;
                    std::cout << "   Interp = " << interp << std::endl;
                    // The following line sounds dangerous...
                    //x++;
#endif

                    if (SplineDegree==1)
                    {
                        // Linear interpolation

                        // Calculate the integer position in input image, be 
                        // careful that it is not the nearest but the one 
                        // at the top left corner of the interpolation square. 
                        // Ie, (0.7,0.7) would give (0,0)
                        // Calculate also weights for point m1+1,n1+1
                        double wx = xp + cen_xp;
                        int m1 = (int) wx;
                        wx = wx - m1;
                        int m2 = m1 + 1;
                        double wy = yp + cen_yp;
                        int n1 = (int) wy;
                        wy = wy - n1;
                        int n2 = n1 + 1;

                        // m2 and n2 can be out by 1 so wrap must be check here
                        if ((m2 >= 0) && ((size_t)m2 >= Xdim))
                            m2 = 0;
                        if ((n2 >=0) && ((size_t)n2 >= Ydim))
                            n2 = 0;

#ifdef DEBUG_APPLYGEO
                        std::cout << "   From (" << n1 << "," << m1 << ") and ("
                                << n2 << "," << m2 << ")\n";
                        std::cout << "   wx= " << wx << " wy= " << wy 
                                << std::endl;
#endif

                        // Perform interpolation
                        // if wx == 0 means that the rightest point is useless 
                        // for this interpolation, and even it might not be 
                        // defined if m1=xdim-1
                        // The same can be said for wy.
                        double wx_1 = (1-wx);
                        double wy_1 = (1-wy);
                        double aux2=wy_1* wx_1 ;
                        double tmp  = aux2 * DIRECT_A2D_ELEM(V1, n1, m1);

                        if ((wx != 0) && ((m2 < 0) || ((size_t)m2 < V1.xdim)))
                            tmp += (wy_1-aux2) * DIRECT_A2D_ELEM(V1, n1, m2);

                        if ((wy != 0) && ((n2 < 0) || ((size_t)n2 < V1.ydim)))
                        {
                            aux2=wy * wx_1;
                            tmp += aux2 * DIRECT_A2D_ELEM(V1, n2, m1);

                            if ((wx != 0) && ((m2 < 0) || ((size_t)m2 < V1.xdim)))
                                tmp += (wy-aux2) * DIRECT_A2D_ELEM(V1, n2, m2);
                        }

                        dAij(V2, i, j) = (T) tmp;
                    }
                    else if (SplineDegree==0)
                    {
                        dAij(V2, i, j) = (T) A2D_ELEM(V1,(int)trunc(yp),(int)trunc(xp));
                    }
                    else if (SplineDegree==3)
                    {
                        // B-spline interpolation
                        dAij(V2, i, j) = (T) BcoeffsToUse->interpolatedElementBSpline2D_Degree3(xp, yp);
                    }
                    else
                    {
                        // B-spline interpolation
                        dAij(V2, i, j) = (T) BcoeffsToUse->interpolatedElementBSpline2D(xp, yp, SplineDegree);
                    }
#ifdef DEBUG_APPYGEO
                    std::cout << "   val= " << dAij(V2, i, j) << std::endl;
#endif

                    // Compute new point inside input image
                    xp += Aref00;
                    yp += Aref10;
                }
            } /* wrap == true */
            else
            {
                if (SplineDegree==1)
                {
                    #pragma simd reduction (+:xp,yp)
                    for (int j=globalMin; j<globalMax ;j++)
                    {
#ifdef DEBUG_APPLYGEO

                        std::cout << "Computing (" << i << "," << j << ")\n";
                        std::cout << "   (y, x) =(" << y << "," << x << ")\n"
                        << "   before wrapping (y',x')=(" << yp << "," << xp << ") "
                        << std::endl;

                        std::cout << "   after wrapping (y',x')=(" << yp << "," << xp << ") "
                        << std::endl;
                        std::cout << "   Interp = " << interp << std::endl;
                        // The following line sounds dangerous...
                        //x++;
#endif

                        // Linear interpolation

                        // Calculate the integer position in input image, be careful
                        // that it is not the nearest but the one at the top left corner
                        // of the interpolation square. Ie, (0.7,0.7) would give (0,0)
                        // Calculate also weights for point m1+1,n1+1
                        double wx = xp + cen_xp;
                        int m1 = (int) wx;
                        wx = wx - m1;
                        int m2 = m1 + 1;
                        double wy = yp + cen_yp;
                        int n1 = (int) wy;
                        wy = wy - n1;
                        int n2 = n1 + 1;

#ifdef DEBUG_APPLYGEO
                        std::cout << "   From (" << n1 << "," << m1 << ") and ("
                            << n2 << "," << m2 << ")\n";
                        std::cout << "   wx= " << wx << " wy= " << wy << std::endl;
#endif

                        // Perform interpolation
                        // if wx == 0 means that the rightest point is useless for this
                        // interpolation, and even it might not be defined if m1=xdim-1
                        // The same can be said for wy.
                        double wx_1 = (1-wx);
                        double wy_1 = (1-wy);
                        double aux2=wy_1* wx_1 ;
                        double tmp  = aux2 * DIRECT_A2D_ELEM(V1, n1, m1);

                        if ((wx != 0) && ((m2 < 0) || ((size_t)m2 < V1.xdim)))
                            tmp += (wy_1-aux2) * DIRECT_A2D_ELEM(V1, n1, m2);

                        if ((wy != 0) && ((n2 < 0) || ((size_t)n2 < V1.ydim)))
                        {
                            aux2=wy * wx_1;
                            tmp += aux2 * DIRECT_A2D_ELEM(V1, n2, m1);
    
                            if ((wx != 0) && ((m2 < 0) || ((size_t)m2 < V1.xdim)))
                                tmp += (wy-aux2) * DIRECT_A2D_ELEM(V1, n2, m2);
                        }

                        dAij(V2, i, j) = (T) tmp;

#ifdef DEBUG_APPYGEO
                        std::cout << "   val= " << dAij(V2, i, j) << std::endl;
#endif

                        // Compute new point inside input image
                        xp += Aref00;
                        yp += Aref10;
                    }
                }
                else if (SplineDegree==0)
                {
                    #pragma simd reduction (+:xp,yp)
                    for (int j=globalMin; j<globalMax ;j++) 
                    {
#ifdef DEBUG_APPLYGEO

                        std::cout << "Computing (" << i << "," << j << ")\n";
                        std::cout << "   (y, x) =(" << y << "," << x << ")\n"
                        << "   before wrapping (y',x')=(" << yp << "," << xp << ") "
                        << std::endl;

                        std::cout << "   after wrapping (y',x')=(" << yp << "," << xp << ") "
                        << std::endl;
                        std::cout << "   Interp = " << interp << std::endl;
                        // The following line sounds dangerous...
                        //x++;
#endif
                        dAij(V2, i, j) = (T) A2D_ELEM(V1,(int)trunc(yp),(int)trunc(xp));

#ifdef DEBUG_APPYGEO
                        std::cout << "   val= " << dAij(V2, i, j) << std::endl;
#endif

                        // Compute new point inside input image
                        xp += Aref00;
                        yp += Aref10;
                    }
                }
                else if (SplineDegree==3)
                {
                    for (int j=globalMin; j<globalMax ;j++)
                    {
#ifdef DEBUG_APPLYGEO

                        std::cout << "Computing (" << i << "," << j << ")\n";
                        std::cout << "   (y, x) =(" << y << "," << x << ")\n"
                        << "   before wrapping (y',x')=(" << yp << "," << xp << ") "
                        << std::endl;

                        std::cout << "   after wrapping (y',x')=(" << yp << "," << xp << ") "
                        << std::endl;
                        std::cout << "   Interp = " << interp << std::endl;
                        // The following line sounds dangerous...
                        //x++;
#endif

                        // B-spline interpolation
                        dAij(V2, i, j) = (T) BcoeffsToUse->interpolatedElementBSpline2D_Degree3(xp, yp);

#ifdef DEBUG_APPYGEO
                        std::cout << "   val= " << dAij(V2, i, j) << std::endl;
#endif

                        // Compute new point inside input image
                        xp += Aref00;
                        yp += Aref10;
                    }
                }
                else
                {
                    for (int j=globalMin; j<globalMax ;j++)
                    {
#ifdef DEBUG_APPLYGEO

                        std::cout << "Computing (" << i << "," << j << ")\n";
                        std::cout << "   (y, x) =(" << y << "," << x << ")\n"
                        << "   before wrapping (y',x')=(" << yp << "," << xp << ") "
                        << std::endl;

                        std::cout << "   after wrapping (y',x')=(" << yp << "," << xp << ") "
                        << std::endl;
                        std::cout << "   Interp = " << interp << std::endl;
                        // The following line sounds dangerous...
                        //x++;
#endif

                        // B-spline interpolation
                        dAij(V2, i, j) = (T) BcoeffsToUse->interpolatedElementBSpline2D(xp, yp, SplineDegree);
#ifdef DEBUG_APPYGEO
                        std::cout << "   val= " << dAij(V2, i, j) << std::endl;
#endif

                        // Compute new point inside input image
                        xp += Aref00;
                        yp += Aref10;
                    }
                }
            } /* wrap == false */

            y++;
        }
    }
    else
    {
        // 3D transformation
        size_t m1, n1, o1, m2, n2, o2;
        double x, y, z, xp, yp, zp;
        double minxp, minyp, maxxp, maxyp, minzp, maxzp;
        double cen_x, cen_y, cen_z, cen_xp, cen_yp, cen_zp;
        double wx, wy, wz;
        double Aref00=MAT_ELEM(Aref,0,0);
        double Aref10=MAT_ELEM(Aref,1,0);
        double Aref20=MAT_ELEM(Aref,2,0);

        // Find center of MultidimArray
        cen_z = (int)(V2.zdim / 2);
        cen_y = (int)(V2.ydim / 2);
        cen_x = (int)(V2.xdim / 2);
        cen_zp = (int)(V1.zdim / 2);
        cen_yp = (int)(V1.ydim / 2);
        cen_xp = (int)(V1.xdim / 2);
        minxp = -cen_xp;
        minyp = -cen_yp;
        minzp = -cen_zp;
        maxxp = V1.xdim - cen_xp - 1;
        maxyp = V1.ydim - cen_yp - 1;
        maxzp = V1.zdim - cen_zp - 1;

#ifdef DEBUG

        std::cout << "Geometry 2 center=("
        << cen_z  << "," << cen_y  << "," << cen_x  << ")\n"
        << "Geometry 1 center=("
        << cen_zp << "," << cen_yp << "," << cen_xp << ")\n"
        << "           min=("
        << minzp  << "," << minyp  << "," << minxp  << ")\n"
        << "           max=("
        << maxzp  << "," << maxyp  << "," << maxxp  << ")\n"
        ;
#endif

        if (SplineDegree > 1)
        {
            // Build the B-spline coefficients
            if (BcoeffsPtr!=NULL)
                BcoeffsToUse=BcoeffsPtr;
            else
            {
                produceSplineCoefficients(SplineDegree, Bcoeffs, V1); //Bcoeffs is a single image
                BcoeffsToUse = &Bcoeffs;
            }
            STARTINGX(*BcoeffsToUse) = (int) minxp;
            STARTINGY(*BcoeffsToUse) = (int) minyp;
            STARTINGZ(*BcoeffsToUse) = (int) minzp;
        }

        // Now we go from the output MultidimArray to the input MultidimArray, ie, for any
        // voxel in the output MultidimArray we calculate which are the corresponding
        // ones in the original MultidimArray, make an interpolation with them and put
        // this value at the output voxel

        // V2 is not initialised to 0 because all its pixels are rewritten
        for (size_t k = 0; k < V2.zdim; k++)
            for (size_t i = 0; i < V2.ydim; i++)
            {
                // Calculate position of the beginning of the row in the output
                // MultidimArray
                x = -cen_x;
                y = i - cen_y;
                z = k - cen_z;

                // Calculate this position in the input image according to the
                // geometrical transformation they are related by
                // coords_output(=x,y) = A * coords_input (=xp,yp)
                xp = x * MAT_ELEM(Aref, 0, 0) + y * MAT_ELEM(Aref, 0, 1) + z * MAT_ELEM(Aref, 0, 2) + MAT_ELEM(Aref, 0, 3);
                yp = x * MAT_ELEM(Aref, 1, 0) + y * MAT_ELEM(Aref, 1, 1) + z * MAT_ELEM(Aref, 1, 2) + MAT_ELEM(Aref, 1, 3);
                zp = x * MAT_ELEM(Aref, 2, 0) + y * MAT_ELEM(Aref, 2, 1) + z * MAT_ELEM(Aref, 2, 2) + MAT_ELEM(Aref, 2, 3);

                for (size_t j = 0; j < V2.xdim; j++)
                {
                    bool interp;
                    double tmp;

#ifdef DEBUG

                    bool show_debug = false;
                    if ((i == 0 && j == 0 && k == 0) ||
                        (i == V2.ydim - 1 && j == V2.xdim - 1 && k == V2.zdim - 1))
                        show_debug = true;

                    if (show_debug)
                        std::cout << "(x,y,z)-->(xp,yp,zp)= "
                        << "(" << x  << "," << y  << "," << z  << ") "
                        << "(" << xp << "," << yp << "," << zp << ")\n";
#endif

                    // If the point is outside the volume, apply a periodic
                    // extension of the volume, what exits by one side enters by
                    // the other
                    interp  = true;
                    bool x_isOut = XMIPP_RANGE_OUTSIDE(xp, minxp, maxxp);
                    bool y_isOut = XMIPP_RANGE_OUTSIDE(yp, minyp, maxyp);
                    bool z_isOut = XMIPP_RANGE_OUTSIDE(zp, minzp, maxzp);

                    if (wrap)
                    {
                        if (x_isOut)
                            xp = realWRAP(xp, minxp - 0.5, maxxp + 0.5);

                        if (y_isOut)
                            yp = realWRAP(yp, minyp - 0.5, maxyp + 0.5);

                        if (z_isOut)
                            zp = realWRAP(zp, minzp - 0.5, maxzp + 0.5);
                    }
                    else if (x_isOut || y_isOut || z_isOut)
                        interp = false;

                    if (interp)
                    {
                        if (SplineDegree == 1)
                        {
                            // Linear interpolation

                            // Calculate the integer position in input volume, be
                            // careful that it is not the nearest but the one at the
                            // top left corner of the interpolation square. Ie,
                            // (0.7,0.7) would give (0,0)
                            // Calculate also weights for point m1+1,n1+1
                            wx = xp + cen_xp;
                            m1 = (int) wx;
                            wx = wx - m1;
                            m2 = m1 + 1;
                            wy = yp + cen_yp;
                            n1 = (int) wy;
                            wy = wy - n1;
                            n2 = n1 + 1;
                            wz = zp + cen_zp;
                            o1 = (int) wz;
                            wz = wz - o1;
                            o2 = o1 + 1;

#ifdef DEBUG

                            if (show_debug)
                            {
                                std::cout << "After wrapping(xp,yp,zp)= "
                                << "(" << xp << "," << yp << "," << zp << ")\n";
                                std::cout << "(m1,n1,o1)-->(m2,n2,o2)="
                                << "(" << m1 << "," << n1 << "," << o1 << ") "
                                << "(" << m2 << "," << n2 << "," << o2 << ")\n";
                                std::cout << "(wx,wy,wz)="
                                << "(" << wx << "," << wy << "," << wz << ")\n";
                            }
#endif

                            // Perform interpolation
                            // if wx == 0 means that the rightest point is useless for
                            // this interpolation, and even it might not be defined if
                            // m1=xdim-1
                            // The same can be said for wy.
                            double wx_1=1-wx;
                            double wy_1=1-wy;
                            double wz_1=1-wz;

                            double aux1=wz_1 * wy_1;
                            double aux2=aux1*wx_1;
                            tmp  =  aux2 * DIRECT_A3D_ELEM(V1, o1, n1, m1);

                            if (wx != 0 && m2 < V1.xdim)
                                tmp += (aux1-aux2)* DIRECT_A3D_ELEM(V1, o1, n1, m2);

                            if (wy != 0 && n2 < V1.ydim)
                            {
                                aux1=wz_1 * wy;
                                aux2=aux1*wx_1;
                                tmp += aux2 * DIRECT_A3D_ELEM(V1, o1, n2, m1);
                                if (wx != 0 && m2 < V1.xdim)
                                    tmp += (aux1-aux2) * DIRECT_A3D_ELEM(V1, o1, n2, m2);
                            }

                            if (wz != 0 && o2 < V1.zdim)
                            {
                                aux1=wz * wy_1;
                                aux2=aux1*wx_1;
                                tmp += aux2 * DIRECT_A3D_ELEM(V1, o2, n1, m1);
                                if (wx != 0 && m2 < V1.xdim)
                                    tmp += (aux1-aux2) * DIRECT_A3D_ELEM(V1, o2, n1, m2);
                                if (wy != 0 && n2 < V1.ydim)
                                {
                                    aux1=wz * wy;
                                    aux2=aux1*wx_1;
                                    tmp += aux2 * DIRECT_A3D_ELEM(V1, o2, n2, m1);
                                    if (wx != 0 && m2 < V1.xdim)
                                        tmp += (aux1-aux2) * DIRECT_A3D_ELEM(V1, o2, n2, m2);
                                }
                            }

#ifdef DEBUG
                            if (show_debug)
                                std::cout <<
                                "tmp1=" << DIRECT_A3D_ELEM(V1, o1, n1, m1) << " "
                                << (T)(wz_1 *wy_1 *wx_1 * DIRECT_A3D_ELEM(V1, o1, n1, m1))
                                << std::endl <<
                                "tmp2=" << DIRECT_A3D_ELEM(V1, o1, n1, m2) << " "
                                << (T)(wz_1 *wy_1 * wx * DIRECT_A3D_ELEM(V1, o1, n1, m2))
                                << std::endl <<
                                "tmp3=" << DIRECT_A3D_ELEM(V1, o1, n2, m1) << " "
                                << (T)(wz_1 * wy *wx_1 * DIRECT_A3D_ELEM(V1, o1, n2, m1))
                                << std::endl <<
                                "tmp4=" << DIRECT_A3D_ELEM(V1, o1, n2, m2) << " "
                                << (T)(wz_1 * wy * wx * DIRECT_A3D_ELEM(V1, o2, n1, m1))
                                << std::endl <<
                                "tmp6=" << DIRECT_A3D_ELEM(V1, o2, n1, m2) << " "
                                << (T)(wz * wy_1 * wx * DIRECT_A3D_ELEM(V1, o2, n1, m2))
                                << std::endl <<
                                "tmp7=" << DIRECT_A3D_ELEM(V1, o2, n2, m1) << " "
                                << (T)(wz * wy *wx_1 * DIRECT_A3D_ELEM(V1, o2, n2, m1))
                                << std::endl <<
                                "tmp8=" << DIRECT_A3D_ELEM(V1, o2, n2, m2) << " "
                                << (T)(wz * wy * wx * DIRECT_A3D_ELEM(V1, o2, n2, m2))
                                << std::endl <<
                                "tmp= " << tmp << std::endl;
#endif

                            dAkij(V2 , k, i, j) = (T)tmp;
                        }
                        else if (SplineDegree==0)
                        {
                            dAkij(V2, k, i, j)=(T)A3D_ELEM(V1,(int)trunc(zp),(int)trunc(yp),(int)trunc(xp));
                        }
                        else
                        {
                            // B-spline interpolation
                            dAkij(V2, k, i, j) =
                                (T) BcoeffsToUse->interpolatedElementBSpline3D(xp, yp, zp, SplineDegree);
                        }
                    }
                    else
                        dAkij(V2, k, i, j) = outside;

                    // Compute new point inside input image
                    xp += Aref00;
                    yp += Aref10;
                    zp += Aref20;
                }
            }
    }
}

applyGeometry() [2/4]

template<typename T >

void applyGeometry	(	int	SplineDegree,
		MultidimArray< T > &	V2,
		const MultidimArrayGeneric &	V1,
		const Matrix2D< double > &	A,
		bool	inv,
		bool	wrap,
		T	outside = 0
	)

Definition at line 1234 of file transformations.h.

{
#define APPLYGEO(type)  applyGeometry(SplineDegree,V2, (*(MultidimArray<type>*)(V1.im)), A, inv, wrap, outside);
    SWITCHDATATYPE(V1.datatype, APPLYGEO)
#undef APPLYGEO
}

applyGeometry() [3/4]

template<>

void applyGeometry	(	int	SplineDegree,
		MultidimArray< std::complex< double > > &	V2,
		const MultidimArray< std::complex< double > > &	V1,
		const Matrix2D< double > &	A,
		bool	inv,
		bool	wrap,
		std::complex< double >	outside,
		MultidimArray< double > *	BcoeffsPtr
	)

Definition at line 668 of file transformations.cpp.

{

    if (SplineDegree > 1)
    {
        MultidimArray<double> re, im, rotre, rotim;
        MultidimArray<std::complex<double> > oneImg;
        double outre, outim;
        re.resize(ZSIZE(V1), YSIZE(V1), XSIZE(V1));
        im.resize(ZSIZE(V1), YSIZE(V1), XSIZE(V1));
        outre = outside.real();
        outim = outside.imag();
        oneImg=V1;
        Complex2RealImag(MULTIDIM_ARRAY(oneImg),
                         MULTIDIM_ARRAY(re), MULTIDIM_ARRAY(im),
                         MULTIDIM_SIZE(oneImg));
        applyGeometry(SplineDegree, rotre, re, A, inv, wrap, outre);
        applyGeometry(SplineDegree, rotim, im, A, inv, wrap, outim);
        V2.resize(oneImg);
        RealImag2Complex(MULTIDIM_ARRAY(rotre), MULTIDIM_ARRAY(rotim),
                         MULTIDIM_ARRAY(V2), MULTIDIM_SIZE(re));
    }
    else
    { //FIXME I do not think you want to recall your self
        REPORT_ERROR(ERR_NOT_IMPLEMENTED,"I do not think you want to recall your self");
       // applyGeometry(SplineDegree, V2, V1, A, inv, wrap, outside); // this was causing crash of the sonarcloud analyzer
    }
}

applyGeometry() [4/4]

void applyGeometry	(	int	SplineDegree,
		MultidimArrayGeneric &	V2,
		const MultidimArrayGeneric &	V1,
		const Matrix2D< double > &	A,
		bool	inv,
		bool	wrap,
		double	outside
	)

Definition at line 713 of file transformations.cpp.

{
#define APPLYGEO(type)  applyGeometry(SplineDegree,(*(MultidimArray<type>*)(V2.im)), \
                        (*(MultidimArray<type>*)(V1.im)), A, inv, wrap, (type) outside);
    SWITCHDATATYPE(V1.datatype, APPLYGEO)
#undef APPLYGEO

}

expandBSpline()

template<typename T >

void expandBSpline	(	int	SplineDegree,
		MultidimArray< double > &	V2,
		const MultidimArray< T > &	V1
	)

Expand a set of B-spline coefficients.

Knowing that V1 is a (single image of) B-spline coefficients, produce the expanded set of B-spline coefficients using the two-scale relationship.

Definition at line 144 of file transformations.cpp.

{
    double g[200]; // Coefficients of the reduce filter
    long ng; // Number of coefficients of the reduce filter
    double h[200]; // Coefficients of the expansion filter
    long nh; // Number of coefficients of the expansion filter
    short IsCentered; // Equal TRUE if the filter is a centered spline, FALSE otherwise */

    // Get the filter
    if (GetPyramidFilter("Centered Spline", SplineDegree, g, &ng, h, &nh,
                         &IsCentered))
        REPORT_ERROR(ERR_UNCLASSIFIED, "Unable to load the filter coefficients");

    MultidimArray< double > aux;
    typeCast(V1, aux);

    if (V1.getDim() == 2)
    {
        V2.initZeros(2 * YSIZE(aux), 2 * XSIZE(aux));
        Expand_2D(MULTIDIM_ARRAY(aux), XSIZE(aux), YSIZE(aux),
                  MULTIDIM_ARRAY(V2), h, nh, IsCentered);
    }
    else if (V1.getDim() == 3)
    {
        V2.initZeros(2 * ZSIZE(aux), 2 * YSIZE(aux), 2 * XSIZE(aux));
        Expand_3D(MULTIDIM_ARRAY(aux), XSIZE(aux), YSIZE(aux), ZSIZE(aux),
                  MULTIDIM_ARRAY(V2), h, nh, IsCentered);
    }
    else
        REPORT_ERROR(ERR_MULTIDIM_DIM,"expandBSpline ERROR: only valid for 2D or 3D arrays");
}

fastRadialAverage()

template<typename T >

void fastRadialAverage	(	const MultidimArray< T > &	m,
		const MultidimArray< int > &	distance,
		int	dim,
		MultidimArray< T > &	radial_mean,
		MultidimArray< int > &	radial_count
	)

Definition at line 1854 of file transformations.h.

{
    // Define the vectors
    radial_mean.initZeros(dim);
    radial_count.initZeros(dim);

    // Perform the radial sum and count pixels that contribute to every
    // distance
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(m)
    {
        int d=DIRECT_MULTIDIM_ELEM(distance,n);
        A1D_ELEM(radial_mean,d) += DIRECT_MULTIDIM_ELEM(m,n);
        ++A1D_ELEM(radial_count,d);
    }

    // Perform the mean
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY1D(radial_mean)
      if (DIRECT_MULTIDIM_ELEM(radial_count,i) > 0)
        DIRECT_MULTIDIM_ELEM(radial_mean,i) /= DIRECT_MULTIDIM_ELEM(radial_count,i);
}

geo2TransformationMatrix()

void geo2TransformationMatrix	(	const MDRow &	imageHeader,
		Matrix2D< double > &	A,
		bool	only_apply_shifts = false
	)

Get geometric transformation matrix from image geometry Now this info is stored in metadata The order is: first flip (if necessary), then rotate, then shift A=Translation*Rotation*(Mirror)

Definition at line 178 of file transformations.cpp.

{
    // This has only been implemented for 2D images...
    double psi = 0, shiftX = 0., shiftY = 0., scale = 1.;
    bool flip = false;

    imageGeo.getValue(MDL_ANGLE_PSI, psi);
    imageGeo.getValue(MDL_SHIFT_X, shiftX);
    imageGeo.getValue(MDL_SHIFT_Y, shiftY);
    imageGeo.getValue(MDL_SCALE, scale);
    imageGeo.getValue(MDL_FLIP, flip);

    psi = realWRAP(psi, 0., 360.);

    int dim = A.Xdim() - 1;
    //This check the case when matrix A is not initialized with correct size
    if (dim < 2 || dim > 3)
    {
        dim = 3;
        A.resizeNoCopy(dim + 1, dim + 1);
    }

    if (only_apply_shifts)
        A.initIdentity();
    else if (dim == 2) //2D geometry
        rotation2DMatrix(psi, A, true);
    else if (dim == 3)//3D geometry
    {
        double rot = 0., tilt = 0., shiftZ = 0.;
        imageGeo.getValue(MDL_ANGLE_ROT, rot);
        imageGeo.getValue(MDL_ANGLE_TILT, tilt);
        imageGeo.getValue(MDL_SHIFT_Z, shiftZ);
        Euler_angles2matrix(rot, tilt, psi, A, true);
        dMij(A, 2, dim) = shiftZ;
    }
    dMij(A, 0, dim) = shiftX;
    dMij(A, 1, dim) = shiftY;

    if (scale != 1.)
    {
        if (scale==0.) // Protection against badly formed metadatas
            scale=1.0;
        if (dim == 2)
        {
            M3x3_BY_CT(A, A, scale);
        }
        else if (dim == 3)
        {
            M4x4_BY_CT(A, A, scale);
        }
        dMij(A, dim, dim) = 1.;
    }

    if (flip)
    {
        dMij(A, 0, 0) *= -1.;
        dMij(A, 0, 1) *= -1.;
        if (dim == 3)
            dMij(A, 0, 2) *= -1.;
    }
}

getLoopRange()

bool getLoopRange	(	double	value,
		double	min,
		double	max,
		double	delta,
		int	loopLimit,
		int &	minIter,
		int &	maxIter
	)

Definition at line 603 of file transformations.cpp.

{
    bool validRange=true;// Return value. TRUE if input value is into range boundaries.

    // Value is under lower boundary.
    if (value < min)
    {
        // If delta is negative -> moving to lower values and never into valid range.
        if (delta <= DELTA_THRESHOLD)
        {
            validRange = false;
        }
        // Compute first and last iterations into valid values.
        else
        {
            minIter = (int) (fabs(min - value) / delta);
            minIter++;
            maxIter = (int) (fabs(max - value) / delta);
        }
    }
    // Value is over upper boundary.
    else if (value >= max)
    {
        // If delta is negative -> moving to lower values and never into valid range.
        if (delta >= -DELTA_THRESHOLD)
        {
            validRange = false;
        }
        // Compute first and last iterations into valid values.
        else
        {
            minIter = (int) (fabs(value - max) / -delta);
            minIter++;
            maxIter = (int) (fabs(value - min) / -delta);
        }
    }
    // First value into valid range.
    else
    {
        // Compute first and last iterations into valid values.
        if (delta > DELTA_THRESHOLD)
        {
            minIter = 0;
            maxIter = (int) (fabs(max - value) / delta);
        }
        // Compute first and last iterations into valid values.
        else if (delta < -DELTA_THRESHOLD)
        {
            minIter = 0;
            maxIter = (int) (fabs(value - min) / -delta);
        }
        // If delta is zero then always in valid range.
        else
        {
            minIter = 0;
            maxIter = loopLimit;
        }
    }

    return(validRange);
}

interpolatedElementBSplineDiffX()

double interpolatedElementBSplineDiffX	(	MultidimArray< double > &	vol,
		double	x,
		double	y,
		double	z,
		int	SplineDegree
	)

Interpolates the value of the 3D matrix M at the point (x,y,z) knowing that this image is a set of B-spline coefficients. And making the diff of x, such-> V=sum(Coef diff(Bx) By Bz) Only for BSplines of degree 3!!

(x,y,z) are in logical coordinates.

Definition at line 853 of file transformations.cpp.

{
    int SplineDegree_1 = SplineDegree - 1;
    double aux;

    // Logical to physical
    z -= STARTINGZ(vol);
    y -= STARTINGY(vol);
    x -= STARTINGX(vol);

    int l1 = CEIL(x - SplineDegree_1);
    int l2 = l1 + SplineDegree;

    int m1 = CEIL(y - SplineDegree_1);
    int m2 = m1 + SplineDegree;

    int n1 = CEIL(z - SplineDegree_1);
    int n2 = n1 + SplineDegree;

    double zyxsum = 0.0;
    int Zdim=(int)ZSIZE(vol);
    int Ydim=(int)YSIZE(vol);
    int Xdim=(int)XSIZE(vol);
    for (int n = n1; n <= n2; n++)
    {
        int equivalent_n=n;
        if      (n<0)
            equivalent_n=-n-1;
        else if (n>=Zdim)
            equivalent_n=2*Zdim-n-1;
        double yxsum = 0.0;
        for (int m = m1; m <= m2; m++)
        {
            int equivalent_m=m;
            if      (m<0)
                equivalent_m=-m-1;
            else if (m>=Ydim)
                equivalent_m=2*Ydim-m-1;
            double xsum = 0.0;
            for (int l = l1; l <= l2; l++)
            {
                double xminusl = x - (double) l;
                int equivalent_l=l;
                if      (l<0)
                    equivalent_l=-l-1;
                else if (l>=Xdim)
                    equivalent_l=2*Xdim-l-1;
                double Coeff = (double) DIRECT_A3D_ELEM(vol, equivalent_n,
                                                             equivalent_m,
                                                             equivalent_l );
                switch (SplineDegree)
                {
                case 2:
                    xsum += Coeff * Bspline02(xminusl);
                    break;
                case 3:
                    BSPLINE03DIFF1(aux,xminusl);
                    xsum += Coeff * aux;
                    break;
                case 4:
                    xsum += Coeff * Bspline04(xminusl);
                    break;
                case 5:
                    xsum += Coeff * Bspline05(xminusl);
                    break;
                case 6:
                    xsum += Coeff * Bspline06(xminusl);
                    break;
                case 7:
                    xsum += Coeff * Bspline07(xminusl);
                    break;
                case 8:
                    xsum += Coeff * Bspline08(xminusl);
                    break;
                case 9:
                    xsum += Coeff * Bspline09(xminusl);
                    break;
                }
            }

            double yminusm = y - (double) m;
            switch (SplineDegree)
            {
            case 2:
                yxsum += xsum * Bspline02(yminusm);
                break;
            case 3:
                BSPLINE03(aux,yminusm);
                yxsum += xsum * aux;
                break;
            case 4:
                yxsum += xsum * Bspline04(yminusm);
                break;
            case 5:
                yxsum += xsum * Bspline05(yminusm);
                break;
            case 6:
                yxsum += xsum * Bspline06(yminusm);
                break;
            case 7:
                yxsum += xsum * Bspline07(yminusm);
                break;
            case 8:
                yxsum += xsum * Bspline08(yminusm);
                break;
            case 9:
                yxsum += xsum * Bspline09(yminusm);
                break;
            }
        }

        double zminusn = z - (double) n;
        switch (SplineDegree)
        {
        case 2:
            zyxsum += yxsum * Bspline02(zminusn);
            break;
        case 3:
            BSPLINE03(aux,zminusn);
            zyxsum += yxsum * aux;
            break;
        case 4:
            zyxsum += yxsum * Bspline04(zminusn);
            break;
        case 5:
            zyxsum += yxsum * Bspline05(zminusn);
            break;
        case 6:
            zyxsum += yxsum * Bspline06(zminusn);
            break;
        case 7:
            zyxsum += yxsum * Bspline07(zminusn);
            break;
        case 8:
            zyxsum += yxsum * Bspline08(zminusn);
            break;
        case 9:
            zyxsum += yxsum * Bspline09(zminusn);
            break;
        }
    }

    return zyxsum;
}

interpolatedElementBSplineDiffY()

double interpolatedElementBSplineDiffY	(	MultidimArray< double > &	vol,
		double	x,
		double	y,
		double	z,
		int	SplineDegree
	)

Interpolates the value of the 3D matrix M at the point (x,y,z) knowing that this image is a set of B-spline coefficients. And making the diff of y, such-> V=sum(Coef Bx diff(By) Bz) Only for BSplines of degree 3!!

(x,y,z) are in logical coordinates.

Definition at line 1008 of file transformations.cpp.

{
    int SplineDegree_1 = SplineDegree - 1;
    double aux;

    // Logical to physical
    z -= STARTINGZ(vol);
    y -= STARTINGY(vol);
    x -= STARTINGX(vol);

    int l1 = CEIL(x - SplineDegree_1);
    int l2 = l1 + SplineDegree;

    int m1 = CEIL(y - SplineDegree_1);
    int m2 = m1 + SplineDegree;

    int n1 = CEIL(z - SplineDegree_1);
    int n2 = n1 + SplineDegree;

    double zyxsum = 0.0;
    int Zdim=(int)ZSIZE(vol);
    int Ydim=(int)YSIZE(vol);
    int Xdim=(int)XSIZE(vol);
    for (int n = n1; n <= n2; n++)
    {
        int equivalent_n=n;
        if      (n<0)
            equivalent_n=-n-1;
        else if (n>=Zdim)
            equivalent_n=2*Zdim-n-1;
        double yxsum = 0.0;
        for (int m = m1; m <= m2; m++)
        {
            int equivalent_m=m;
            if      (m<0)
                equivalent_m=-m-1;
            else if (m>=Ydim)
                equivalent_m=2*Ydim-m-1;
            double xsum = 0.0;
            for (int l = l1; l <= l2; l++)
            {
                double xminusl = x - (double) l;
                int equivalent_l=l;
                if      (l<0)
                    equivalent_l=-l-1;
                else if (l>=Xdim)
                    equivalent_l=2*Xdim-l-1;
                double Coeff = (double) DIRECT_A3D_ELEM(vol, equivalent_n,
                                                             equivalent_m,
                                                             equivalent_l );
                double aux;
                switch (SplineDegree)
                {
                case 2:
                    xsum += Coeff * Bspline02(xminusl);
                    break;
                case 3:
                    BSPLINE03(aux,xminusl);
                    xsum += Coeff * aux;
                    break;
                case 4:
                    xsum += Coeff * Bspline04(xminusl);
                    break;
                case 5:
                    xsum += Coeff * Bspline05(xminusl);
                    break;
                case 6:
                    xsum += Coeff * Bspline06(xminusl);
                    break;
                case 7:
                    xsum += Coeff * Bspline07(xminusl);
                    break;
                case 8:
                    xsum += Coeff * Bspline08(xminusl);
                    break;
                case 9:
                    xsum += Coeff * Bspline09(xminusl);
                    break;
                }
            }

            double yminusm = y - (double) m;
            switch (SplineDegree)
            {
            case 2:
                yxsum += xsum * Bspline02(yminusm);
                break;
            case 3:
                BSPLINE03DIFF1(aux,yminusm);
                yxsum += xsum * aux;
                break;
            case 4:
                yxsum += xsum * Bspline04(yminusm);
                break;
            case 5:
                yxsum += xsum * Bspline05(yminusm);
                break;
            case 6:
                yxsum += xsum * Bspline06(yminusm);
                break;
            case 7:
                yxsum += xsum * Bspline07(yminusm);
                break;
            case 8:
                yxsum += xsum * Bspline08(yminusm);
                break;
            case 9:
                yxsum += xsum * Bspline09(yminusm);
                break;
            }
        }

        double zminusn = z - (double) n;
        switch (SplineDegree)
        {
        case 2:
            zyxsum += yxsum * Bspline02(zminusn);
            break;
        case 3:
            BSPLINE03(aux,zminusn);
            zyxsum += yxsum * aux;
            break;
        case 4:
            zyxsum += yxsum * Bspline04(zminusn);
            break;
        case 5:
            zyxsum += yxsum * Bspline05(zminusn);
            break;
        case 6:
            zyxsum += yxsum * Bspline06(zminusn);
            break;
        case 7:
            zyxsum += yxsum * Bspline07(zminusn);
            break;
        case 8:
            zyxsum += yxsum * Bspline08(zminusn);
            break;
        case 9:
            zyxsum += yxsum * Bspline09(zminusn);
            break;
        }
    }

    return zyxsum;
}

interpolatedElementBSplineDiffZ()

double interpolatedElementBSplineDiffZ	(	MultidimArray< double > &	vol,
		double	x,
		double	y,
		double	z,
		int	SplineDegree
	)

Interpolates the value of the 3D matrix M at the point (x,y,z) knowing that this image is a set of B-spline coefficients. And making the diff of z, such-> V=sum(Coef Bx By diff(Bz)) Only for BSplines of degree 3!!

(x,y,z) are in logical coordinates.

Definition at line 1164 of file transformations.cpp.

{
    int SplineDegree_1 = SplineDegree - 1;
    double aux;

    // Logical to physical
    z -= STARTINGZ(vol);
    y -= STARTINGY(vol);
    x -= STARTINGX(vol);

    int l1 = CEIL(x - SplineDegree_1);
    int l2 = l1 + SplineDegree;

    int m1 = CEIL(y - SplineDegree_1);
    int m2 = m1 + SplineDegree;

    int n1 = CEIL(z - SplineDegree_1);
    int n2 = n1 + SplineDegree;

    double zyxsum = 0.0;
    int Zdim=(int)ZSIZE(vol);
    int Ydim=(int)YSIZE(vol);
    int Xdim=(int)XSIZE(vol);
    for (int n = n1; n <= n2; n++)
    {
        int equivalent_n=n;
        if      (n<0)
            equivalent_n=-n-1;
        else if (n>=Zdim)
            equivalent_n=2*Zdim-n-1;
        double yxsum = 0.0;
        for (int m = m1; m <= m2; m++)
        {
            int equivalent_m=m;
            if      (m<0)
                equivalent_m=-m-1;
            else if (m>=Ydim)
                equivalent_m=2*Ydim-m-1;
            double xsum = 0.0;
            for (int l = l1; l <= l2; l++)
            {
                double xminusl = x - (double) l;
                int equivalent_l=l;
                if      (l<0)
                    equivalent_l=-l-1;
                else if (l>=Xdim)
                    equivalent_l=2*Xdim-l-1;
                double Coeff = (double) DIRECT_A3D_ELEM(vol, equivalent_n,
                                                             equivalent_m,
                                                             equivalent_l );
                double aux;
                switch (SplineDegree)
                {
                case 2:
                    xsum += Coeff * Bspline02(xminusl);
                    break;
                case 3:
                    BSPLINE03(aux,xminusl);
                    xsum += Coeff * aux;
                    break;
                case 4:
                    xsum += Coeff * Bspline04(xminusl);
                    break;
                case 5:
                    xsum += Coeff * Bspline05(xminusl);
                    break;
                case 6:
                    xsum += Coeff * Bspline06(xminusl);
                    break;
                case 7:
                    xsum += Coeff * Bspline07(xminusl);
                    break;
                case 8:
                    xsum += Coeff * Bspline08(xminusl);
                    break;
                case 9:
                    xsum += Coeff * Bspline09(xminusl);
                    break;
                }
            }

            double yminusm = y - (double) m;
            switch (SplineDegree)
            {
            case 2:
                yxsum += xsum * Bspline02(yminusm);
                break;
            case 3:
                BSPLINE03(aux,yminusm);
                yxsum += xsum * aux;
                break;
            case 4:
                yxsum += xsum * Bspline04(yminusm);
                break;
            case 5:
                yxsum += xsum * Bspline05(yminusm);
                break;
            case 6:
                yxsum += xsum * Bspline06(yminusm);
                break;
            case 7:
                yxsum += xsum * Bspline07(yminusm);
                break;
            case 8:
                yxsum += xsum * Bspline08(yminusm);
                break;
            case 9:
                yxsum += xsum * Bspline09(yminusm);
                break;
            }
        }

        double zminusn = z - (double) n;
        switch (SplineDegree)
        {
        case 2:
            zyxsum += yxsum * Bspline02(zminusn);
            break;
        case 3:
            BSPLINE03DIFF1(aux,zminusn);
            zyxsum += yxsum * aux;
            break;
        case 4:
            zyxsum += yxsum * Bspline04(zminusn);
            break;
        case 5:
            zyxsum += yxsum * Bspline05(zminusn);
            break;
        case 6:
            zyxsum += yxsum * Bspline06(zminusn);
            break;
        case 7:
            zyxsum += yxsum * Bspline07(zminusn);
            break;
        case 8:
            zyxsum += yxsum * Bspline08(zminusn);
            break;
        case 9:
            zyxsum += yxsum * Bspline09(zminusn);
            break;
        }
    }

    return zyxsum;
}

produceImageFromSplineCoefficients()

template<typename T >

void produceImageFromSplineCoefficients	(	int	SplineDegree,
		MultidimArray< T > &	img,
		const MultidimArray< double > &	coeffs
	)

Produce image from B-spline coefficients.

Note that the coeffs and img are only single images!

Definition at line 64 of file transformations.cpp.

{
    MultidimArray< double > imgD;
    imgD.initZeros(ZSIZE(coeffs), YSIZE(coeffs), XSIZE(coeffs));
    STARTINGX(img) = STARTINGX(coeffs);
    STARTINGY(img) = STARTINGY(coeffs);
    STARTINGZ(img) = STARTINGZ(coeffs);

    int Status;
    MultidimArray< double > aux(coeffs);

    ChangeBasisVolume(MULTIDIM_ARRAY(aux), MULTIDIM_ARRAY(imgD),
                      XSIZE(coeffs), YSIZE(coeffs), ZSIZE(coeffs),
                      BasicSpline, CardinalSpline, SplineDegree,
                      MirrorOnBounds, DBL_EPSILON, &Status);
    if (Status)
        REPORT_ERROR(ERR_UNCLASSIFIED, "Error in ImageFromSplineCoefficients...");
    typeCast(imgD, img);
}

produceSplineCoefficients() [1/2]

template<typename T >

void produceSplineCoefficients	(	int	SplineDegree,
		MultidimArray< double > &	coeffs,
		const MultidimArray< T > &	V1
	)

Produce spline coefficients.

Create a single image with spline coefficients for the nth image

Definition at line 41 of file transformations.cpp.

{

    coeffs.initZeros(ZSIZE(V1), YSIZE(V1), XSIZE(V1));
    STARTINGX(coeffs) = STARTINGX(V1);
    STARTINGY(coeffs) = STARTINGY(V1);
    STARTINGZ(coeffs) = STARTINGZ(V1);

    int Status;
    MultidimArray< double > aux;
    typeCast(V1, aux); // This will create a single volume!

    ChangeBasisVolume(MULTIDIM_ARRAY(aux), MULTIDIM_ARRAY(coeffs),
                      XSIZE(V1), YSIZE(V1), ZSIZE(V1),
                      CardinalSpline, BasicSpline, SplineDegree,
                      MirrorOffBounds, DBL_EPSILON, &Status);
    if (Status)
        REPORT_ERROR(ERR_UNCLASSIFIED, "Error in produceSplineCoefficients...");
}

produceSplineCoefficients() [2/2]

void produceSplineCoefficients	(	int	SplineDegree,
		MultidimArray< double > &	coeffs,
		const MultidimArray< std::complex< double > > &	V1
	)

Definition at line 727 of file transformations.cpp.

{
    // TODO Implement
    REPORT_ERROR(ERR_NOT_IMPLEMENTED,"Spline coefficients of a complex matrix is not implemented.");
}

pyramidExpand() [1/2]

template<typename T >

void pyramidExpand	(	int	SplineDegree,
		MultidimArray< T > &	V2,
		const MultidimArray< T > &	V1,
		int	levels = 1
	)

Expand the nth volume by 2 using a BSpline pyramid.

Definition at line 1640 of file transformations.h.

{
    MultidimArray< double > coeffs, coeffs2;

    produceSplineCoefficients(SplineDegree, coeffs, V1);

    for (int i = 0; i < levels; i++)
    {
        expandBSpline(SplineDegree, coeffs2, coeffs);
        coeffs = coeffs2;
    }

    produceImageFromSplineCoefficients(SplineDegree, V2, coeffs);

}

pyramidExpand() [2/2]

void pyramidExpand	(	int	SplineDegree,
		MultidimArrayGeneric &	V2,
		const MultidimArrayGeneric &	V1,
		int	levels = 1
	)

Definition at line 819 of file transformations.cpp.

{
  if (V1.datatype != V2.datatype)
    REPORT_ERROR(ERR_PARAM_INCORRECT, "pyramidExpand: MultidimArrayGeneric requires same datatype");
#define PYRAMIDEXPAND(type) pyramidExpand(SplineDegree,MULTIDIM_ARRAY_TYPE(V2,type),\
                                              MULTIDIM_ARRAY_TYPE(V1,type), levels);
    SWITCHDATATYPE(V1.datatype, PYRAMIDEXPAND)
#undef PYRAMIDEXPAND
}

pyramidReduce() [1/2]

template<typename T >

void pyramidReduce	(	int	SplineDegree,
		MultidimArray< T > &	V2,
		const MultidimArray< T > &	V1,
		int	levels = 1
	)

Reduce the nth volume by 2 using a BSpline pyramid.

Definition at line 1597 of file transformations.h.

{
    MultidimArray< double > coeffs, coeffs2;

    produceSplineCoefficients(SplineDegree, coeffs, V1);

    for (int i = 0; i < levels; i++)
    {
        reduceBSpline(SplineDegree, coeffs2, coeffs);
        coeffs = coeffs2;
    }

    produceImageFromSplineCoefficients(SplineDegree, V2, coeffs);
}

pyramidReduce() [2/2]

void pyramidReduce	(	int	SplineDegree,
		MultidimArrayGeneric &	V2,
		const MultidimArrayGeneric &	V1,
		int	levels = 1
	)

Definition at line 832 of file transformations.cpp.

{
  if (V1.datatype != V2.datatype)
    REPORT_ERROR(ERR_PARAM_INCORRECT, "pyramidReduce: MultidimArrayGeneric requires same datatype");
#define PYRAMIDREDUCE(type) pyramidReduce(SplineDegree,MULTIDIM_ARRAY_TYPE(V2,type),\
                                              MULTIDIM_ARRAY_TYPE(V1,type), levels);
    SWITCHDATATYPE(V1.datatype, PYRAMIDREDUCE)
#undef PYRAMIDREDUCE
}

radialAverage()

template<typename T >

void radialAverage	(	const MultidimArray< T > &	m,
		Matrix1D< int > &	center_of_rot,
		MultidimArray< T > &	radial_mean,
		MultidimArray< int > &	radial_count,
		const bool &	rounding = false
	)

Does a radial average of a 2D/3D image, around the voxel where is the origin.

A vector radial_mean is returned where:

• the first element is the mean of the voxels whose distance to the origin is (0-1),
• the second element is the mean of the voxels whose distance to the origin is (1-2)
• and so on.

A second vector radial_count is returned containing the number of voxels over which each radial average was calculated.

Sjors nov2003: if rounding=true, element=round(distance);

• so the first element is the mean of the voxels whose distance to the origin is (0.5-1.5),
• the second element is the mean of the voxels whose distance to the origin is (1.5-2.5)
• and so on.

Definition at line 1711 of file transformations.h.

{
    Matrix1D< double > idx(3);

    // If center_of_rot was written for 2D image
    if (center_of_rot.size() < 3)
        center_of_rot.resize(3);

    // First determine the maximum distance that one should expect, to set the
    // dimension of the radial average vector
    MultidimArray< int > distances(8);

    double z = STARTINGZ(m) - ZZ(center_of_rot);
    double y = STARTINGY(m) - YY(center_of_rot);
    double x = STARTINGX(m) - XX(center_of_rot);

    distances(0) = (int) floor(sqrt(x * x + y * y + z * z));
    x = FINISHINGX(m) - XX(center_of_rot);

    distances(1) = (int) floor(sqrt(x * x + y * y + z * z));
    y = FINISHINGY(m) - YY(center_of_rot);

    distances(2) = (int) floor(sqrt(x * x + y * y + z * z));
    x = STARTINGX(m) - XX(center_of_rot);

    distances(3) = (int) floor(sqrt(x * x + y * y + z * z));
    z = FINISHINGZ(m) - ZZ(center_of_rot);

    distances(4) = (int) floor(sqrt(x * x + y * y + z * z));
    x = FINISHINGX(m) - XX(center_of_rot);

    distances(5) = (int) floor(sqrt(x * x + y * y + z * z));
    y = STARTINGY(m) - YY(center_of_rot);

    distances(6) = (int) floor(sqrt(x * x + y * y + z * z));
    x = STARTINGX(m) - XX(center_of_rot);

    distances(7) = (int) floor(sqrt(x * x + y * y + z * z));

    int dim = (int) CEIL(distances.computeMax()) + 1;
    if (rounding)
        dim++;

    // Define the vectors
    radial_mean.initZeros(dim);
    radial_count.initZeros(dim);

    // Perform the radial sum and count pixels that contribute to every
    // distance
    FOR_ALL_ELEMENTS_IN_ARRAY3D(m)
    {
        ZZ(idx) = k - ZZ(center_of_rot);
        YY(idx) = i - YY(center_of_rot);
        XX(idx) = j - XX(center_of_rot);

        // Determine distance to the center
        ;
        double module = sqrt(ZZ(idx)*ZZ(idx)+YY(idx)*YY(idx)+XX(idx)*XX(idx));
        int distance = (rounding) ? (int) round(module) : (int) floor(module);

        // Sum the value to the pixels with the same distance
        DIRECT_MULTIDIM_ELEM(radial_mean,distance) += A3D_ELEM(m, k, i, j);

        // Count the pixel
        DIRECT_MULTIDIM_ELEM(radial_count,distance)++;
    }

    // Perform the mean
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY1D(radial_mean)
      if (DIRECT_MULTIDIM_ELEM(radial_count,i) > 0)
        DIRECT_MULTIDIM_ELEM(radial_mean,i) /= DIRECT_MULTIDIM_ELEM(radial_count,i);
}

radialAverageAxis()

template<typename T >

void radialAverageAxis	(	const MultidimArray< T > &	in,
		char	axis,
		MultidimArray< double > &	out
	)

Definition at line 1880 of file transformations.h.

{
    MultidimArray<double> inCentered;
    inCentered.alias(in);
    inCentered.setXmippOrigin();
    if (axis=='z')
    {
        out.initZeros(ZSIZE(in),XSIZE(in));
        out.setXmippOrigin();
        for (int i=STARTINGY(inCentered); i<=FINISHINGY(inCentered); ++i)
        {
            double z=i;
            for (int j=0; j<XSIZE(out)/2; ++j)
            {
                for (double ang=0; ang<TWOPI; ang+=TWOPI/72)
                {
                    double x=j*cos(ang);
                    double y=j*sin(ang);
                    double val=inCentered.interpolatedElement3D(x,y,z);
                    A2D_ELEM(out,i,j)+=val;
                }
                A2D_ELEM(out,i,j)/=73;
                A2D_ELEM(out,i,-j)=A2D_ELEM(out,i,j);
            }
        }
    }

    else
        REPORT_ERROR(ERR_ARG_INCORRECT,"Not implemented yet");
}

radialAverageNonCubic()

template<typename T >

void radialAverageNonCubic	(	const MultidimArray< T > &	m,
		Matrix1D< int > &	center_of_rot,
		MultidimArray< T > &	radial_mean,
		MultidimArray< int > &	radial_count,
		bool	rounding = false
	)

Definition at line 1375 of file transformations.cpp.

{
    Matrix1D< double > idx(3);

    size_t sizemax = std::max({XSIZE(m), YSIZE(m), ZSIZE(m)});
    double scalex = double(XSIZE(m)/sizemax);
    double scaley = double(YSIZE(m)/sizemax);
    double scalez = double(ZSIZE(m)/sizemax);

    // If center_of_rot was written for 2D image
    if (center_of_rot.size() < 3)
        center_of_rot.resize(3);

    // First determine the maximum distance that one should expect, to set the
    // dimension of the radial average vector
    MultidimArray< int > distances(8);

    const double z0 = STARTINGZ(m) - ZZ(center_of_rot);
    const double y0 = STARTINGY(m) - YY(center_of_rot);
    const double x0 = STARTINGX(m) - XX(center_of_rot);

    const double xf = FINISHINGX(m) - XX(center_of_rot);
    const double yf = FINISHINGY(m) - YY(center_of_rot);
    const double zf = FINISHINGZ(m) - ZZ(center_of_rot);

    distances(0) = (int) floor(sqrt(x0 * x0 + y0 * y0 + z0 * z0));
    distances(1) = (int) floor(sqrt(xf * xf + y0 * y0 + z0 * z0));
    distances(2) = (int) floor(sqrt(xf * xf + yf * yf + z0 * z0));
    distances(3) = (int) floor(sqrt(x0 * x0 + yf * yf + z0 * z0));
    distances(4) = (int) floor(sqrt(x0 * x0 + yf * yf + zf * zf));
    distances(5) = (int) floor(sqrt(xf * xf + yf * yf + zf * zf));
    distances(6) = (int) floor(sqrt(xf * xf + y0 * y0 + zf * zf));
    distances(7) = (int) floor(sqrt(x0 * x0 + y0 * y0 + zf * zf));

    int dim = CEIL(distances.computeMax()) + 1;
    if (rounding)
        dim++;

    // Define the vectors
    radial_mean.initZeros(dim);
    radial_count.initZeros(dim);

    // Perform the radial sum and count pixels that contribute to every
    // distance
    FOR_ALL_ELEMENTS_IN_ARRAY3D(m)
    {
        ZZ(idx) = scalez * (k - ZZ(center_of_rot));
        YY(idx) = scaley * (i - YY(center_of_rot));
        XX(idx) = scalex * (j - XX(center_of_rot));

        // Determine distance to the center
        double mod = sqrt(ZZ(idx)*ZZ(idx)+YY(idx)*YY(idx)+XX(idx)*XX(idx));
        int distance = rounding ? (int) round(mod) : (int) floor(mod);

        // Sum the value to the pixels with the same distance
        DIRECT_MULTIDIM_ELEM(radial_mean,distance) += A3D_ELEM(m, k, i, j);

        // Count the pixel
        DIRECT_MULTIDIM_ELEM(radial_count,distance)++;
    }

    // Perform the mean
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY1D(radial_mean)
    if (DIRECT_MULTIDIM_ELEM(radial_count,i) > 0)
        DIRECT_MULTIDIM_ELEM(radial_mean,i) /= DIRECT_MULTIDIM_ELEM(radial_count,i);
}

radialAveragePrecomputeDistance()

template<typename T >

void radialAveragePrecomputeDistance	(	const MultidimArray< T > &	m,
		Matrix1D< int > &	center_of_rot,
		MultidimArray< int > &	distance,
		int &	dim,
		const bool &	rounding = false
	)

Definition at line 1789 of file transformations.h.

{
    Matrix1D< double > idx(3);

    // If center_of_rot was written for 2D image
    if (center_of_rot.size() < 3)
        center_of_rot.resize(3);

    // First determine the maximum distance that one should expect, to set the
    // dimension of the radial average vector
    MultidimArray< int > distances(8);

    double z = STARTINGZ(m) - ZZ(center_of_rot);
    double y = STARTINGY(m) - YY(center_of_rot);
    double x = STARTINGX(m) - XX(center_of_rot);

    distances(0) = (int) floor(sqrt(x * x + y * y + z * z));
    x = FINISHINGX(m) - XX(center_of_rot);

    distances(1) = (int) floor(sqrt(x * x + y * y + z * z));
    y = FINISHINGY(m) - YY(center_of_rot);

    distances(2) = (int) floor(sqrt(x * x + y * y + z * z));
    x = STARTINGX(m) - XX(center_of_rot);

    distances(3) = (int) floor(sqrt(x * x + y * y + z * z));
    z = FINISHINGZ(m) - ZZ(center_of_rot);

    distances(4) = (int) floor(sqrt(x * x + y * y + z * z));
    x = FINISHINGX(m) - XX(center_of_rot);

    distances(5) = (int) floor(sqrt(x * x + y * y + z * z));
    y = STARTINGY(m) - YY(center_of_rot);

    distances(6) = (int) floor(sqrt(x * x + y * y + z * z));
    x = STARTINGX(m) - XX(center_of_rot);

    distances(7) = (int) floor(sqrt(x * x + y * y + z * z));

    dim = (int) CEIL(distances.computeMax()) + 1;
    if (rounding)
        dim++;

    // Perform the radial sum and count pixels that contribute to every
    // distance
    distance.initZeros(m);
    FOR_ALL_ELEMENTS_IN_ARRAY3D(m)
    {
        ZZ(idx) = k - ZZ(center_of_rot);
        YY(idx) = i - YY(center_of_rot);
        XX(idx) = j - XX(center_of_rot);

        // Determine distance to the center
        if (rounding)
            A3D_ELEM(distance,k,i,j) = (int) round(idx.module());
        else
            A3D_ELEM(distance,k,i,j) = (int) floor(idx.module());
    }
}

radiallySymmetrize()

void radiallySymmetrize	(	const MultidimArray< double > &	img,
		MultidimArray< double > &	radialImg
	)

Definition at line 1312 of file transformations.cpp.

{
    Matrix1D<int> center(2);
    center.initZeros();
    MultidimArray<int> distance, radial_count;
    MultidimArray<double> radial_mean;
    int dim;
    radialAveragePrecomputeDistance(img, center, distance, dim);
    fastRadialAverage(img, distance, dim, radial_mean, radial_count);

    radialImg.initZeros(img);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(radialImg)
    {
        int d=DIRECT_MULTIDIM_ELEM(distance,n);
        DIRECT_MULTIDIM_ELEM(radialImg,n)=A1D_ELEM(radial_mean,d);
    }
}

reduceBSpline()

template<typename T >

void reduceBSpline	(	int	SplineDegree,
		MultidimArray< double > &	V2,
		const MultidimArray< T > &	V1
	)

Reduce a set of B-spline coefficients.

Knowing that V1 is a (single image of) B-spline coefficients, produce the reduced set of B-spline coefficients using the two-scale relationship.

Definition at line 87 of file transformations.cpp.

{
    double g[200]; // Coefficients of the reduce filter
    long ng; // Number of coefficients of the reduce filter
    double h[200]; // Coefficients of the expansion filter
    long nh; // Number of coefficients of the expansion filter
    short IsCentered; // Equal TRUE if the filter is a centered spline

    // Get the filter
    const char *splineType="Centered Spline";
    if (GetPyramidFilter(splineType, SplineDegree,
                         g, &ng, h, &nh, &IsCentered))
        REPORT_ERROR(ERR_UNCLASSIFIED, "Unable to load the filter coefficients");

    MultidimArray< double>  aux;
    typeCast(V1, aux);
    if (V1.getDim() == 2)
    {
        if (XSIZE(aux) % 2 != 0 && YSIZE(aux) % 2 != 0)
            aux.resize(YSIZE(aux) - 1, XSIZE(aux) - 1);
        else if (YSIZE(aux) % 2 != 0)
            aux.resize(YSIZE(aux) - 1, XSIZE(aux));
        else if (XSIZE(aux) % 2 != 0)
            aux.resize(YSIZE(aux), XSIZE(aux) - 1);

        V2.initZeros(YSIZE(aux) / 2, XSIZE(aux) / 2);
        Reduce_2D(MULTIDIM_ARRAY(aux), XSIZE(aux), YSIZE(aux),
                  MULTIDIM_ARRAY(V2), g, ng, IsCentered);
    }
    else if (V1.getDim() == 3)
    {
        if (XSIZE(aux) % 2 != 0 && YSIZE(aux) % 2 != 0 && ZSIZE(aux) % 2 != 0)
            aux.resize(ZSIZE(aux - 1), YSIZE(aux) - 1, XSIZE(aux) - 1);
        else if (XSIZE(aux) % 2 != 0 && YSIZE(aux) % 2 != 0 && ZSIZE(aux) % 2 == 0)
            aux.resize(ZSIZE(aux), YSIZE(aux) - 1, XSIZE(aux) - 1);
        else if (XSIZE(aux) % 2 != 0 && YSIZE(aux) % 2 == 0 && ZSIZE(aux) % 2 != 0)
            aux.resize(ZSIZE(aux) - 1, YSIZE(aux), XSIZE(aux) - 1);
        else if (XSIZE(aux) % 2 != 0 && YSIZE(aux) % 2 == 0 && ZSIZE(aux) % 2 == 0)
            aux.resize(ZSIZE(aux), YSIZE(aux), XSIZE(aux) - 1);
        else if (XSIZE(aux) % 2 == 0 && YSIZE(aux) % 2 != 0 && ZSIZE(aux) % 2 != 0)
            aux.resize(ZSIZE(aux) - 1, YSIZE(aux) - 1, XSIZE(aux));
        else if (XSIZE(aux) % 2 == 0 && YSIZE(aux) % 2 != 0 && ZSIZE(aux) % 2 == 0)
            aux.resize(ZSIZE(aux), YSIZE(aux) - 1, XSIZE(aux));
        else if (XSIZE(aux) % 2 == 0 && YSIZE(aux) % 2 == 0 && ZSIZE(aux) % 2 != 0)
            aux.resize(ZSIZE(aux) - 1, YSIZE(aux), XSIZE(aux));

        V2.initZeros(ZSIZE(aux) / 2, YSIZE(aux) / 2, XSIZE(aux) / 2);
        Reduce_3D(MULTIDIM_ARRAY(aux), XSIZE(aux), YSIZE(aux), ZSIZE(aux),
                  MULTIDIM_ARRAY(V2), g, ng, IsCentered);
    }
    else
        REPORT_ERROR(ERR_MULTIDIM_DIM,"reduceBSpline ERROR: only valid for 2D or 3D arrays");
}

rotate()

template<typename T >

void rotate	(	int	SplineDegree,
		MultidimArray< T > &	V2,
		const MultidimArray< T > &	V1,
		double	ang,
		char	axis = 'Z',
		bool	wrap = xmipp_transformation::DONT_WRAP,
		T	outside = 0
	)

Rotate an array around a given system axis.

The rotation angle is in degrees, and the rotational axis is either 'X', 'Y' or 'Z' for 3D arrays and only 'Z' for 2D arrays. An exception is thrown if the axis given is not one of these.

V2 = V1.rotate(60);

Definition at line 1323 of file transformations.h.

{
    Matrix2D< double > tmp;
    if (V1.getDim()==2)
    {
        rotation2DMatrix(ang,tmp);
    }
    else if (V1.getDim()==3)
    {
        rotation3DMatrix(ang, axis, tmp);
    }
    else
        REPORT_ERROR(ERR_MULTIDIM_DIM,"rotate ERROR: rotate only valid for 2D or 3D arrays");

    applyGeometry(SplineDegree, V2, V1, tmp, xmipp_transformation::IS_NOT_INV, wrap, outside);
}

rotation2DMatrix()

template<typename T >

void rotation2DMatrix	(	T	ang,
		Matrix2D< T > &	m,
		bool	homogeneous = true
	)

Creates a rotational matrix (3x3) for images

The rotation angle is in degrees. m must have been already resized to 3x3

rotation2DMatrix(60,m);

Definition at line 361 of file transformations.cpp.

{
    // FIXME DS this might not be true, but just to make sure
    static_assert(std::is_floating_point<T>::value,
            "Only float and double are allowed as template parameters");

    T rad = DEG2RAD(ang);
    T cosine = cos(rad);
    T sine = sin(rad);

    if (homogeneous)
    {
        result.resizeNoCopy(3,3); // sizes will be tested inside
        // now we have 3x3 matrix row wise matrix
        result.mdata[0] = cosine;
        result.mdata[1] = sine;
        result.mdata[2] = 0;

        result.mdata[3] = -sine;
        result.mdata[4] = cosine;
        result.mdata[5] = 0;

        result.mdata[6] = 0;
        result.mdata[7] = 0;
        result.mdata[8] = 1;
    } else {
        result.resizeNoCopy(2,2); // sizes will be tested inside
        // now we have 2x2 matrix row wise matrix
        result.mdata[0] = cosine;
        result.mdata[1] = sine;

        result.mdata[2] = -sine;
        result.mdata[3] = cosine;
    }
}

rotation3DMatrix() [1/2]

void rotation3DMatrix	(	double	ang,
		char	axis,
		Matrix2D< double > &	m,
		bool	homogeneous = true
	)

Creates a rotational matrix (4x4) for volumes around system axis

The rotation angle is in degrees, and the rotational axis is either 'X', 'Y' or 'Z'. An exception is thrown if the axis given is not one of these.

The returned matrices are respectively alpha degrees around Z

[ cos(A) -sin(A)     0   ]
[ sin(A)  cos(A)     0   ]
[   0       0        1   ]

alpha degrees around Y

[ cos(A)    0    -sin(A) ]
[   0       1       0    ]
[ sin(A)    0     cos(A) ]

alpha degrees around X

[   1       0       0    ]
[   0     cos(A) -sin(A) ]
[   0     sin(A)  cos(A) ]

m = rotation3DMatrix(60, 'X');

Definition at line 426 of file transformations.cpp.

{
    if (homogeneous)
    {
        result.initZeros(4,4);
        dMij(result,3, 3) = 1;
    }
    else
        result.initZeros(3,3);

    double cosine, sine;
    ang = DEG2RAD(ang);
    cosine = cos(ang);
    sine = sin(ang);

    switch (axis)
    {
    case 'Z':
        dMij(result,0, 0) = cosine;
        dMij(result,0, 1) = sine;
        dMij(result,1, 0) = -sine;
        dMij(result,1, 1) = cosine;
        dMij(result,2, 2) = 1;
        break;
    case 'Y':
        dMij(result,0, 0) = cosine;
        dMij(result,0, 2) = sine;
        dMij(result,2, 0) = -sine;
        dMij(result,2, 2) = cosine;
        dMij(result,1, 1) = 1;
        break;
    case 'X':
        dMij(result,1, 1) = cosine;
        dMij(result,1, 2) = sine;
        dMij(result,2, 1) = -sine;
        dMij(result,2, 2) = cosine;
        dMij(result,0, 0) = 1;
        break;
    default:
        REPORT_ERROR(ERR_VALUE_INCORRECT, "rotation3DMatrix: Unknown axis");
    }
}

rotation3DMatrix() [2/2]

void rotation3DMatrix	(	double	ang,
		const Matrix1D< double > &	axis,
		Matrix2D< double > &	m,
		bool	homogeneous = true
	)

Creates a rotational matrix (4x4) for volumes around any axis

The rotation angle is in degrees, and the rotational axis is given as a R3 vector. An exception is thrown if the axis is not a R3 vector. The axis needs not to be unitary.

m = rotation3DMatrix(60, vectorR3(1, 1, 1));

Definition at line 517 of file transformations.cpp.

{
#ifdef NEVERDEFINED
    // Compute a matrix which makes the turning axis coincident with Z
    // And turn around this axis
    Matrix2D<double> A, R;
    alignWithZ(axis, A, homogeneous);
    rotation3DMatrix(ang, 'Z', R, homogeneous);
    result = A.transpose() * R * A;
#else
    // http://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
    if (homogeneous)
        result.initIdentity(4);
    else
        result.initIdentity(3);
    double s,c;
    //sincos(-DEG2RAD(ang),&s,&c);
    s = sin(-DEG2RAD(ang));
    c = cos(-DEG2RAD(ang));
    double c1=1-c;
    double x=XX(axis);
    double y=YY(axis);
    double z=ZZ(axis);
    double xy=x*y;
    double xz=x*z;
    double yz=y*z;
    double x2=x*x;
    double y2=y*y;
    double z2=z*z;
    dMij(result,0,0)=c+x2*c1;
    dMij(result,0,1)=xy*c1-z*s;
    dMij(result,0,2)=xz*c1+y*s;
    dMij(result,1,0)=xy*c1+z*s;
    dMij(result,1,1)=c+y2*c1;
    dMij(result,1,2)=yz*c1-x*s;
    dMij(result,2,0)=xz*c1-y*s;
    dMij(result,2,1)=yz*c1+x*s;
    dMij(result,2,2)=c+z2*c1;
#endif
}

rotation3DMatrixFromIcoOrientations()

void rotation3DMatrixFromIcoOrientations	(	const char *	icoFrom,
		const char *	icoTo,
		Matrix2D< double > &	R
	)

Creates a rotation matrix (R) to change the orientation from one standard orientation (icoFrom) to another one (icoTo)

Definition at line 1332 of file transformations.cpp.

{
    std::vector<char> symLabel;
    symLabel.push_back((int)icoFrom[1]);
    symLabel.push_back((int)icoTo[1]);
    Matrix1D<double> xyz(3); 
    Matrix2D<double> Rfrom, Rto;

    for(int i=0; i<2; i++)
    {
        switch (symLabel[i])
        {
            case '1':
                xyz = vectorR3(0.0, 90.0, 0.0);
                break;
            case '2':
                xyz = vectorR3(0.0, 0.0, 0.0);
                break;
            case '3':
                xyz = vectorR3(0.0, 31.7175, 0.0);
                break;
            case '4':
                xyz = vectorR3(0.0, -31.7175, 0.0);
                break;
/*            case '5':
                xyz = vectorR3(31.7175, 90.0, 0.0); // not aviable yet
                break;
            case '6':
                xyz = vectorR3(-31.7175, 90.0, 0.0); // not aviable yet
                break;*/
            default:
                REPORT_ERROR(ERR_PARAM_INCORRECT, "Incorrect standard icosahedral orientation");
        }
        if (i==0)
            Euler_angles2matrix(XX(xyz), YY(xyz), ZZ(xyz), Rfrom, true);
        else
            Euler_angles2matrix(XX(xyz), YY(xyz), ZZ(xyz), Rto, true);
    }
    R = Rto * Rfrom.transpose();
}

scale3DMatrix()

void scale3DMatrix	(	const Matrix1D< double > &	sc,
		Matrix2D< double > &	m,
		bool	homogeneous = true
	)

Creates a scaling matrix (4x4) for volumes

The scaling factors for the different axis must be given as a vector. So that, XX(sc)=scale for X axis, YY(sc)=...

Definition at line 584 of file transformations.cpp.

{
    if (VEC_XSIZE(sc) != 3)
        REPORT_ERROR(ERR_MATRIX_SIZE, "Scale3D_matrix: vector is not in R3");

    if (homogeneous)
    {
        result.initZeros(4,4);
        dMij(result,3, 3) = 1;
    }
    else
        result.initZeros(3,3);
    dMij(result,0, 0) = XX(sc);
    dMij(result,1, 1) = YY(sc);
    dMij(result,2, 2) = ZZ(sc);
}

scaleToSize() [1/5]

template<typename T1 , typename T >

void scaleToSize	(	int	SplineDegree,
		MultidimArray< T > &	V2,
		const MultidimArray< T1 > &	V1,
		size_t	Xdim,
		size_t	Ydim,
		size_t	Zdim = 1
	)

Scales to a new size.

The volume is scaled (resampled) to fill a new size. It is not the same as "selfWindow" in this same class. The size can be larger or smaller than the actual one.

V2 = V1.scaleToSize(128, 128, 128);

Definition at line 1451 of file transformations.h.

{
    if (Xdim == XSIZE(V1) && Ydim == YSIZE(V1) && \
        (!(V1.getDim()==3) || (Zdim == ZSIZE(V1))) )
    {
        typeCast(V1,V2);
        return;
    }

    Matrix2D< double > tmp;
    if (V1.getDim()==2)
    {
        tmp.initIdentity(3);
        tmp(0, 0) = (double) Xdim / (double) XSIZE(V1);
        tmp(1, 1) = (double) Ydim / (double) YSIZE(V1);
        V2.initZeros(1, 1, Ydim, Xdim);
    }
    else if (V1.getDim()==3)
    {
        tmp.initIdentity(4);
        tmp(0, 0) = (double) Xdim / (double) XSIZE(V1);
        tmp(1, 1) = (double) Ydim / (double) YSIZE(V1);
        tmp(2, 2) = (double) Zdim / (double) ZSIZE(V1);
        V2.initZeros(1, Zdim, Ydim, Xdim);
    }
    else
        REPORT_ERROR(ERR_MULTIDIM_DIM,"scaleToSize ERROR: scaleToSize only valid for 2D or 3D arrays");

    V2.setXmippOrigin();
    applyGeometry(SplineDegree, V2, V1, tmp, xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP, (T)0);
}

scaleToSize() [2/5]

template<typename T >

void scaleToSize ( int  SplineDegree, MultidimArray< T > &  V2, const MultidimArrayGeneric &  V1, int  Xdim, int  Ydim, int  Zdim = 1  )

inline

Scales to a new size.

The volume is scaled (resampled) to fill a new size. Same as previous, but in this case the input is a MultidimArrayGeneric.

scaleToSize(VolumeOut, VolumeGenericInput, 128, 128, 128);

Definition at line 1497 of file transformations.h.

{
#define SCALETOSIZE(type) scaleToSize(SplineDegree,V2,*((MultidimArray<type>*)(V1.im)),Xdim,Ydim,Zdim);
    SWITCHDATATYPE(V1.datatype,SCALETOSIZE)
#undef SCALETOSIZE
}

scaleToSize() [3/5]

template<typename T >

void scaleToSize ( int  SplineDegree, MultidimArrayGeneric &  V2, const MultidimArray< T > &  V1, int  Xdim, int  Ydim, int  Zdim = 1  )

inline

Scales to a new size.

The volume is scaled (resampled) to fill a new size. Same as previous, but in this case the input is a MultidimArrayGeneric.

scaleToSize(VolumeOut, VolumeGenericInput, 128, 128, 128);

Definition at line 1517 of file transformations.h.

{
#define SCALETOSIZE(type) scaleToSize(SplineDegree,MULTIDIM_ARRAY_TYPE(V2,type),V1,Xdim,Ydim,Zdim);
    SWITCHDATATYPE(V1.datatype,SCALETOSIZE)
#undef SCALETOSIZE
}

scaleToSize() [4/5]

void scaleToSize	(	int	SplineDegree,
		MultidimArrayGeneric &	V2,
		const MultidimArrayGeneric &	V1,
		int	Xdim,
		int	Ydim,
		int	Zdim = 1
	)

Scales to a new size.

The volume is scaled (resampled) to fill a new size. Same as previous, but in this case the input is a MultidimArrayGeneric.

scaleToSize(VolumeOut, VolumeGenericInput, 128, 128, 128);

Definition at line 746 of file transformations.cpp.

{
  if (V1.datatype != V2.datatype)
    REPORT_ERROR(ERR_PARAM_INCORRECT, "scaleToSize: MultidimArrayGeneric requires same datatype");
#define SCALETOSIZE(type) scaleToSize(SplineDegree,MULTIDIM_ARRAY_TYPE(V2,type), \
                          MULTIDIM_ARRAY_TYPE(V1,type),Xdim,Ydim,Zdim);
    SWITCHDATATYPE(V1.datatype, SCALETOSIZE)
#undef SCALETOSIZE
}

scaleToSize() [5/5]

void scaleToSize	(	int	SplineDegree,
		MultidimArray< std::complex< double > > &	V2,
		const MultidimArray< std::complex< double > > &	V1,
		int	Xdim,
		int	Ydim,
		int	Zdim = 1
	)

Definition at line 759 of file transformations.cpp.

{
    if (SplineDegree > 1)
    {
        MultidimArray< double > re, im, aux;
        MultidimArray<std::complex<double> > oneImg;

        re.resize(ZSIZE(V1), YSIZE(V1), XSIZE(V1));
        im.resize(ZSIZE(V1), YSIZE(V1), XSIZE(V1));

        oneImg=V1;
        Complex2RealImag(MULTIDIM_ARRAY(oneImg),
                         MULTIDIM_ARRAY(re), MULTIDIM_ARRAY(im),
                         MULTIDIM_SIZE(oneImg));
        aux = re;
        scaleToSize(SplineDegree, re, aux, Ydim, Xdim, Zdim);
        aux = im;
        scaleToSize(SplineDegree, im, aux, Ydim, Xdim, Zdim);
        RealImag2Complex(MULTIDIM_ARRAY(re), MULTIDIM_ARRAY(im),
                         MULTIDIM_ARRAY(V2), MULTIDIM_SIZE(re));
    }
    else
        scaleToSize(SplineDegree, V2, V1, Xdim, Ydim, Zdim);

}

selfApplyGeometry() [1/2]

template<typename T >

void selfApplyGeometry	(	int	SplineDegree,
		MultidimArray< T > &	V1,
		const Matrix2D< double > &	A,
		bool	inv,
		bool	wrap,
		T	outside = 0
	)

Applies a geometrical transformation and overwrites the input matrix.

The same as the previous function, but input array is overwritten

Definition at line 1253 of file transformations.h.

{
    MultidimArray<T> aux = V1;
    V1.initZeros();
    applyGeometry(SplineDegree, V1, aux, A, inv, wrap, outside);
}

selfApplyGeometry() [2/2]

template<>

void selfApplyGeometry	(	int	SplineDegree,
		MultidimArray< std::complex< double > > &	V1,
		const Matrix2D< double > &	A,
		bool	inv,
		bool	wrap,
		std::complex< double >	outside
	)

Definition at line 704 of file transformations.cpp.

{
    MultidimArray<std::complex<double> > aux = V1;
    applyGeometry(Splinedegree, V1, aux, A, inv, wrap, outside);
}

selfPyramidExpand() [1/2]

template<typename T >

void selfPyramidExpand	(	int	SplineDegree,
		MultidimArray< T > &	V1,
		int	levels = 1
	)

Expand the nth volume by 2 using a BSpline pyramid.

The same as the previous function, but input array is overwritten

Definition at line 1675 of file transformations.h.

{
    MultidimArray<T> aux = V1;
    V1.initZeros();
    pyramidExpand(SplineDegree, V1, aux, levels);
}

selfPyramidExpand() [2/2]

void selfPyramidExpand	(	int	SplineDegree,
		MultidimArrayGeneric &	V1,
		int	levels
	)

Same as previous but for MultidimArrayGeneric

Definition at line 809 of file transformations.cpp.

{
#define SELFPYRAMIDEXPAND(type) selfPyramidExpand(SplineDegree, \
                                      *((MultidimArray<type>*)(V1.im)), levels);
    SWITCHDATATYPE(V1.datatype,SELFPYRAMIDEXPAND);
#undef SELFPYRAMIDEXPAND
}

selfPyramidReduce() [1/2]

template<typename T >

void selfPyramidReduce	(	int	SplineDegree,
		MultidimArray< T > &	V1,
		int	levels = 1
	)

Reduce the nth volume by 2 using a BSpline pyramid.

The same as the previous function, but input array is overwritten

Definition at line 1621 of file transformations.h.

{
    MultidimArray<T> aux = V1;
    V1.initZeros();
    pyramidReduce(SplineDegree, V1, aux, levels);
}

selfPyramidReduce() [2/2]

void selfPyramidReduce	(	int	SplineDegree,
		MultidimArrayGeneric &	V1,
		int	levels
	)

Same as previous but for MultidimArrayGeneric

Same as template version but for MultidimArrayGeneric

Definition at line 798 of file transformations.cpp.

{
#define SELFPYRAMIDREDUCE(type) selfPyramidReduce(SplineDegree, \
                                      *((MultidimArray<type>*)(V1.im)), levels);
    SWITCHDATATYPE(V1.datatype,SELFPYRAMIDREDUCE);
#undef SELFPYRAMIDREDUCE
}

selfRotate()

template<typename T >

void selfRotate	(	int	SplineDegree,
		MultidimArray< T > &	V1,
		double	ang,
		char	axis = 'Z',
		bool	wrap = xmipp_transformation::DONT_WRAP,
		T	outside = 0
	)

Rotate an array around a given system axis.

The same as the previous function, but input array is overwritten

Definition at line 1350 of file transformations.h.

{
    MultidimArray<T> aux = V1;
    V1.initZeros();
    rotate(SplineDegree, V1, aux, ang, axis, wrap, outside);
}

selfScaleToSize() [1/3]

template<typename T >

void selfScaleToSize	(	int	SplineDegree,
		MultidimArray< T > &	V1,
		int	Xdim,
		int	Ydim,
		int	Zdim = 1
	)

Scales to a new size.

The same as the previous function, but input array is overwritten

Definition at line 1546 of file transformations.h.

{
    MultidimArray<T> aux = V1;
    scaleToSize(SplineDegree, V1, aux, Xdim, Ydim, Zdim);
}

selfScaleToSize() [2/3]

void selfScaleToSize	(	int	SplineDegree,
		MultidimArrayGeneric &	V1,
		int	Xdim,
		int	Ydim,
		int	Zdim = 1
	)

Definition at line 736 of file transformations.cpp.

{
#define SELFSCALETOSIZE(type) selfScaleToSize(SplineDegree,MULTIDIM_ARRAY_TYPE(V1,type), \
                                                                Xdim,Ydim,Zdim);
    SWITCHDATATYPE(V1.datatype,SELFSCALETOSIZE)
#undef SELFSCALETOSIZE
}

selfScaleToSize() [3/3]

void selfScaleToSize	(	int	SplineDegree,
		MultidimArray< std::complex< double > > &	V1,
		int	Xdim,
		int	Ydim,
		int	Zdim = 1
	)

Definition at line 789 of file transformations.cpp.

{
    MultidimArray<std::complex<double> > aux;
    scaleToSize(SplineDegree, V1, aux, Xdim, Ydim, Zdim);
}

selfTranslate()

template<typename T >

void selfTranslate	(	int	SplineDegree,
		MultidimArray< T > &	V1,
		const Matrix1D< double > &	v,
		bool	wrap = xmipp_transformation::WRAP,
		T	outside = 0
	)

Translate an array.

The same as the previous function, but input array is overwritten

Definition at line 1394 of file transformations.h.

{
    MultidimArray<T> aux = V1;
    V1.initZeros();
    translate(SplineDegree, V1, aux, v, wrap, outside);
}

selfTranslateCenterOfMassToCenter()

template<typename T >

void selfTranslateCenterOfMassToCenter	(	int	SplineDegree,
		MultidimArray< T > &	V1,
		bool	wrap = xmipp_transformation::WRAP
	)

Translate center of mass to center

The same as the previous function, but input array is overwritten

Definition at line 1430 of file transformations.h.

{
    MultidimArray<T> aux = V1;
    V1.initZeros();
    translateCenterOfMassToCenter(SplineDegree, V1, aux, wrap);
}

string2TransformationMatrix()

void string2TransformationMatrix	(	const String &	matrixStr,
		Matrix2D< double > &	matrix,
		size_t	dim = 4
	)

Retrieve the matrix from an string representation Valid formats are: [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]] or 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 or 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1

Definition at line 241 of file transformations.cpp.

{
  matrix.resizeNoCopy(dim, dim);

  String matrixStrCopy(matrixStr);
  char c;

  for (size_t i = 0; i < matrixStr.size(); i++)
  {
    c = matrixStr[i];
    if (c == '[' or c == ']' or c == ',')
      matrixStrCopy[i] = ' ';
  }

  size_t n = 4; // EMX matrix are always 4x4
  std::stringstream ss(matrixStrCopy);
  size_t d_1 = dim - 1;

  for (size_t i = 0; i < n; ++i)
    for (size_t j = 0; j < n; ++j)
    {
        //TODO validate that M(dim, dim) is 1
        //if (i == dim-1 && j == dim-1)

      ss >> dMij(matrix, i < dim ? i : i-1, j < dim ? j : j-1);
    }
  dMij(matrix, d_1, d_1) = 1.;

}

transformationMatrix2Geo()

void transformationMatrix2Geo	(	const Matrix2D< double > &	A,
		MDRow &	imageGeo
	)

Retrieve the geometry transfomations from matrix

Definition at line 330 of file transformations.cpp.

{
    bool flip = false;
    double scale = 1, shiftX = 0, shiftY = 0, psi = 0, shiftZ = 0, rot = 0, tilt = 0;

    int dim = A.Xdim() - 1;

    if (dim == 2)
        transformationMatrix2Parameters2D(A, flip, scale, shiftX, shiftY, psi);
    else if (dim == 3)
    {
        transformationMatrix2Parameters3D(A, flip, scale, shiftX, shiftY, shiftZ,
                                          rot, tilt, psi);
        ADD_IF_EXIST_NONZERO(MDL_ANGLE_ROT, rot);
        ADD_IF_EXIST_NONZERO(MDL_ANGLE_TILT, tilt);
        ADD_IF_EXIST_NONZERO(MDL_SHIFT_Z, shiftZ);
    }

    ADD_IF_EXIST_NONZERO(MDL_ANGLE_PSI, psi);
    ADD_IF_EXIST_NONZERO(MDL_SHIFT_X, shiftX);
    ADD_IF_EXIST_NONZERO(MDL_SHIFT_Y, shiftY);

    if (imageGeo.containsLabel(MDL_SCALE) || !XMIPP_EQUAL_REAL(scale, 1.))
        imageGeo.setValue(MDL_SCALE, scale);

    if (imageGeo.containsLabel(MDL_FLIP) || flip)
        imageGeo.setValue(MDL_FLIP, flip);
}

transformationMatrix2Parameters2D()

template<typename T >

void transformationMatrix2Parameters2D	(	const Matrix2D< T > &	A,
		bool &	flip,
		T &	scale,
		T &	shiftX,
		T &	shiftY,
		T &	psi
	)

Retrieve the geometry transformations from matrix for 2D.

Definition at line 272 of file transformations.cpp.

{
    // FIXME DS this might not be true, but just to make sure
    static_assert(std::is_floating_point<T>::value,
            "Only float and double are allowed as template parameters");
    //Calculate determinant for getting flip
    flip = ((dMij(A, 0, 0) * dMij(A, 1, 1) - dMij(A, 0, 1) * dMij(A, 1, 0) ) < 0);
    int sgn = flip ? -1 : 1;
    T cosine = sgn * dMij(A, 0, 0);
    T sine = sgn * dMij(A, 0, 1);
    T scale2 = cosine * cosine +  sine * sine;
    scale = sqrt(scale2);
    T invScale = 1 / scale;
    shiftX = dMij(A, 0, 2) * invScale;
    shiftY = dMij(A, 1, 2) * invScale;
    psi = RAD2DEG(atan2(sine, cosine));
}

transformationMatrix2Parameters3D()

void transformationMatrix2Parameters3D	(	const Matrix2D< double > &	A,
		bool &	flip,
		double &	scale,
		double &	shiftX,
		double &	shiftY,
		double &	shiftZ,
		double &	rot,
		double &	tilt,
		double &	psi
	)

Retrieve the geometry transformations from matrix for D.

Definition at line 295 of file transformations.cpp.

{
    scale =  sqrt(dMij(A,2,0)*dMij(A,2,0) \
                  + dMij(A,2,1)*dMij(A,2,1)\
                  + dMij(A,2,2)*dMij(A,2,2) );
    double invScale = 1./ scale;

    Matrix2D<double> tmpMatrix(4,4);
    M4x4_BY_CT(tmpMatrix, A, invScale);

    Matrix2D<double> eulerMatrix(3,3);

    FOR_ALL_ELEMENTS_IN_MATRIX2D(eulerMatrix)
    dMij(eulerMatrix,i,j) = dMij(tmpMatrix, i, j);
    //check determinant if -1 then flip = true
    flip = tmpMatrix.det3x3() < 0;
    if (flip)
    {
        dMij(eulerMatrix, 0, 0) *= -1.;
        dMij(eulerMatrix, 0, 1) *= -1.;
        dMij(eulerMatrix, 0, 2) *= -1.;
    }
    Euler_matrix2angles(eulerMatrix, rot, tilt, psi);

    shiftX = dMij(tmpMatrix,0,3);
    shiftY = dMij(tmpMatrix,1,3);
    shiftZ = dMij(tmpMatrix,2,3);
}

translate()

template<typename T >

void translate	(	int	SplineDegree,
		MultidimArray< T > &	V2,
		const MultidimArray< T > &	V1,
		const Matrix1D< double > &	v,
		bool	wrap = xmipp_transformation::WRAP,
		T	outside = 0
	)

Translate a array.

The shift is given as a R2 or R3 vector (shift_X, shift_Y, shift_Z) for 2D and 3D arrays, respectively. An exception is thrown if the displacement is not a R3 vector.

// Displacement of 2 pixels down
V2 = V1.translate(vectorR3(0, 0, 2));

Definition at line 1372 of file transformations.h.

{
    Matrix2D< double > tmp;
    if (V1.getDim()==2)
        translation2DMatrix(v, tmp,true);
    else if (V1.getDim()==3)
        translation3DMatrix(v, tmp,true);
    else
        REPORT_ERROR(ERR_MULTIDIM_DIM,"translate ERROR: translate only valid for 2D or 3D arrays");
    applyGeometry(SplineDegree, V2, V1, tmp, xmipp_transformation::IS_INV, wrap, outside);
}

translateCenterOfMassToCenter()

template<typename T >

void translateCenterOfMassToCenter	(	int	SplineDegree,
		MultidimArray< T > &	V2,
		const MultidimArray< T > &	V1,
		bool	wrap = xmipp_transformation::WRAP
	)

Translate center of mass to center

If the input has very high values, it is better to rescale it to be between 0 and 1.

Definition at line 1411 of file transformations.h.

{
    V2 = V1;
    V2.setXmippOrigin();
    Matrix1D< double > center;
    V2.centerOfMass(center);
    center *= -1;
    translate(SplineDegree, V2, V1, center, wrap, 0);
}

translation2DMatrix()

template<typename T >

void translation2DMatrix	(	const Matrix1D< T > &	translation,
		Matrix2D< T > &	resMatrix,
		bool	inverse = false
	)

Creates a translational matrix (3x3) for images

The shift is given as a R2 vector (shift_X, shift_Y). An exception is thrown if the displacement is not a R2 vector.

// Displacement of 1 pixel to the right
 resMatrix = translation2DMatrix(vectorR2(1, 0));

Definition at line 405 of file transformations.cpp.

{
    if (VEC_XSIZE(translation) != 2)
        REPORT_ERROR(ERR_MATRIX_SIZE, "Translation2D_matrix: vector is not in R2");

    resMatrix.initIdentity(3);
    if (inverse)
    {
        dMij(resMatrix,0, 2) = -XX(translation);
        dMij(resMatrix,1, 2) = -YY(translation);
    }
    else
    {
        dMij(resMatrix,0, 2) = XX(translation);
        dMij(resMatrix,1, 2) = YY(translation);
    }
}

translation3DMatrix()

template<typename T >

void translation3DMatrix	(	const Matrix1D< T > &	translation,
		Matrix2D< T > &	resMatrix,
		bool	inverse = false
	)

Creates a translational matrix (4x4) for volumes

The shift is given as a R3 vector (shift_X, shift_Y, shift_Z). An exception is thrown if the displacement is not a R3 vector.

// Displacement of 2 pixels down
resMatrix = translation3DMatrix(vectorR3(0, 0, 2));

Definition at line 563 of file transformations.cpp.

{
    if (VEC_XSIZE(translation) != 3)
        REPORT_ERROR(ERR_MATRIX_SIZE, "Translation3D_matrix: vector is not in R3");

    resMatrix.initIdentity(4);
    if (inverse)
    {
        dMij(resMatrix,0, 3) = -XX(translation);
        dMij(resMatrix,1, 3) = -YY(translation);
        dMij(resMatrix,2, 3) = -ZZ(translation);
    }
    else
    {
        dMij(resMatrix,0, 3) = XX(translation);
        dMij(resMatrix,1, 3) = YY(translation);
        dMij(resMatrix,2, 3) = ZZ(translation);
    }
}