Polar coordinates

Collaboration diagram for Polar coordinates:

Classes
class	Polar_fftw_plans
class	Polar< T >
class	RotationalCorrelationAux

Functions
void	fourierTransformRings (Polar< double > &in, Polar< std::complex< double > > &out, Polar_fftw_plans &plans, bool conjugated=DONT_CONJUGATE)
void	inverseFourierTransformRings (Polar< std::complex< double > > &in, Polar< double > &out, Polar_fftw_plans &plans, bool conjugated=DONT_CONJUGATE)
void	rotationalCorrelation (const Polar< std::complex< double > > &M1, const Polar< std::complex< double > > &M2, MultidimArray< double > &angles, RotationalCorrelationAux &aux)
template<bool NORMALIZE>
void	polarFourierTransform (const MultidimArray< double > &in, Polar< double > &inAux, Polar< std::complex< double > > &out, bool conjugated, int first_ring, int last_ring, Polar_fftw_plans *&plans, int BsplineOrder)
template<bool NORMALIZE>
void	polarFourierTransform (const MultidimArray< double > &in, Polar< std::complex< double > > &out, bool flag, int first_ring, int last_ring, Polar_fftw_plans *&plans, int BsplineOrder=3)
void	normalizedPolarFourierTransform (Polar< double > &polarIn, Polar< std::complex< double > > &out, bool conjugated, Polar_fftw_plans *&plans)
double	best_rotation (const Polar< std::complex< double > > &I1, const Polar< std::complex< double > > &I2, RotationalCorrelationAux &aux)
void	alignRotationally (MultidimArray< double > &I1, MultidimArray< double > &I2, RotationalCorrelationAux &aux, int splineOrder=1, int wrap=xmipp_transformation::WRAP)
void	image_convertCartesianToPolar (MultidimArray< double > &in, MultidimArray< double > &out, double Rmin, double Rmax, double deltaR, double angMin, double angMax, double deltaAng)
void	image_convertCartesianToPolar_ZoomAtCenter (const MultidimArray< double > &in, MultidimArray< double > &out, Matrix1D< double > &R, double zoomFactor, double Rmin, double Rmax, int NRSteps, float angMin, double angMax, int NAngSteps)
void	volume_convertCartesianToCylindrical (const MultidimArray< double > &in, MultidimArray< double > &out, double Rmin, double Rmax, double deltaR, float angMin, double angMax, float deltaAng, Matrix1D< double > &axis)
void	volume_convertCartesianToSpherical (const MultidimArray< double > &in, MultidimArray< double > &out, double Rmin, double Rmax, double deltaR=1., double deltaRot=2.*PI/360., double deltaTilt=PI/180.)

Detailed Description

Function Documentation

alignRotationally()

void alignRotationally	(	MultidimArray< double > &	I1,
		MultidimArray< double > &	I2,
		RotationalCorrelationAux &	aux,
		int	splineOrder = 1,
		int	wrap = xmipp_transformation::WRAP
	)

Align I2 rotationally to I1

Definition at line 234 of file polar.cpp.

                                                                  {
    I1.setXmippOrigin();
    I2.setXmippOrigin();

    Polar_fftw_plans *plans = nullptr;
    Polar<std::complex<double> > polarFourierI2;
    Polar<std::complex<double> > polarFourierI1;
    polarFourierTransform<true>(I1, polarFourierI1, false, XSIZE(I1) / 5,
            XSIZE(I1) / 2, plans);
    polarFourierTransform<true>(I2, polarFourierI2, true, XSIZE(I2) / 5,
            XSIZE(I2) / 2, plans);

    MultidimArray<double> rotationalCorr;
    rotationalCorr.resize(2 * polarFourierI2.getSampleNoOuterRing() - 1);
    aux.local_transformer.setReal(rotationalCorr);
    double bestRot = best_rotation(polarFourierI1, polarFourierI2, aux);

    MultidimArray<double> tmp = I2;
    rotate(splineOrder, I2, tmp, -bestRot, 'Z', wrap);

}

best_rotation()

double best_rotation	(	const Polar< std::complex< double > > &	I1,
		const Polar< std::complex< double > > &	I2,
		RotationalCorrelationAux &	aux
	)

Best rotation between two normalized polar Fourier transforms.

Definition at line 212 of file polar.cpp.

                                                                             {
    MultidimArray<double> angles;
    rotationalCorrelation(I1, I2, angles, aux);

    // Compute the maximum of correlation (inside local_transformer)
    const MultidimArray<double> &corr = aux.local_transformer.getReal();
    int imax = 0;
    double maxval = DIRECT_MULTIDIM_ELEM(corr,0);
    double* ptr = nullptr;
    unsigned long int n;
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(corr,n,ptr)
        if (*ptr > maxval) {
            maxval = *ptr;
            imax = n;
        }

    // Return the corresponding angle
    return DIRECT_A1D_ELEM(angles,imax);
}

fourierTransformRings()

void fourierTransformRings	(	Polar< double > &	in,
		Polar< std::complex< double > > &	out,
		Polar_fftw_plans &	plans,
		bool	conjugated = DONT_CONJUGATE
	)

Calculate FourierTransform of all rings

This function returns a polar of complex<double> by calculating the FT of all rings

Note that the Polar_fftw_plans may be re-used for polars of the same size They should initially be calculated as in the example below

MultidimArray<double> angles, corr;
Polar_fftw_plans plans;
Polar<std::complex<double> > F1, F2;

M1.calculateFftwPlans(plans); // M1 is a Polar<double>
fourierTransformRings(M1, F1, plans, false);
fourierTransformRings(M1, F2, plans, true);// complex conjugated

Definition at line 34 of file polar.cpp.

                         {
    MultidimArray<std::complex<double> > Fring;
    out.clear();
    for (int iring = 0; iring < in.getRingNo(); iring++) {

        plans.arrays[iring] = in.rings[iring];
        (plans.transformers[iring])->FourierTransform();
        (plans.transformers[iring])->getFourierAlias(Fring);
        if (conjugated) {
            auto *ptrFring_i = (double*) &DIRECT_A1D_ELEM(Fring,0);
            ++ptrFring_i;
            for (size_t i = 0; i < XSIZE(Fring); ++i, ptrFring_i += 2)
                (*ptrFring_i) *= -1;
        }
        out.rings.push_back(Fring);
    }
    out.mode = in.mode;
    out.ring_radius = in.ring_radius;
}

image_convertCartesianToPolar()

void image_convertCartesianToPolar	(	MultidimArray< double > &	in,
		MultidimArray< double > &	out,
		double	Rmin,
		double	Rmax,
		double	deltaR,
		double	angMin,
		double	angMax,
		double	deltaAng
	)

Produce a polar image from a cartesian image. You can give the minimum and maximum radius for the interpolation, the delta radius and the delta angle. It is assumed that the input image has the Xmipp origin. Delta Ang must be in radians.

image_convertCartesianToPolar_ZoomAtCenter()

void image_convertCartesianToPolar_ZoomAtCenter	(	const MultidimArray< double > &	in,
		MultidimArray< double > &	out,
		Matrix1D< double > &	R,
		double	zoomFactor,
		double	Rmin,
		double	Rmax,
		int	NRSteps,
		float	angMin,
		double	angMax,
		int	NAngSteps
	)

Produce a polar image from a cartesian image. Identical to the previous function. The Zoom factor must be >=1 (==1 means no zoom). If R is with a 0 size, then it is assumed that this routine must compute the sampling radii. If it is not 0, then it is assumed that the sampling radii have already been calculated.

Definition at line 276 of file polar.cpp.

                       {
    /* Octave
     * r=0:64;plot(r,d.^(r/d),r,r,d*((r/d).^2.8));
     */
    float deltaAng = (angMax - angMin + 1) / NAngSteps;
    double deltaR = (Rmax - Rmin + 1) / NRSteps;
    out.initZeros(NAngSteps, NRSteps);
    if (VEC_XSIZE(R) == 0) {
        R.initZeros(NRSteps);
        double Rrange = Rmax - Rmin;
        for (int j = 0; j < NRSteps; ++j)
            VEC_ELEM(R,j) = Rmin
                    + Rrange * pow(j * deltaR / Rrange, zoomFactor);
    }
    for (int i = 0; i < NAngSteps; ++i) {
        float angle = angMin + i * deltaAng;
        float s = sin(angle);
        float c = cos(angle);
        for (int j = 0; j < NRSteps; ++j) {
            float Rj = VEC_ELEM(R,j);
            A2D_ELEM(out,i,j) = in.interpolatedElement2D(Rj * c, Rj * s);
        }
    }
}

inverseFourierTransformRings()

void inverseFourierTransformRings	(	Polar< std::complex< double > > &	in,
		Polar< double > &	out,
		Polar_fftw_plans &	plans,
		bool	conjugated = DONT_CONJUGATE
	)

Calculate inverse FourierTransform of all rings

This function returns a Polar<double> by calculating the inverse FT of all rings in a Polar<std::complex<double> >

Note that the Polar_fftw_plans may be re-used for polars of the same size They should initially be calculated as in the example below

MultidimArray<double> angles, corr;
Polar_fftw_plans plans;
Polar<std::complex<double> > F1, F2;

M1.calculateFftwPlans(plans); // M1 is a Polar<double>
fourierTransformRings(M1, F1, plans);
inverseFourierTransformRings(F1, M1, plans);

Definition at line 87 of file polar.cpp.

                                                                      {
    out.clear();
    for (int iring = 0; iring < in.getRingNo(); iring++) {
        (plans.transformers[iring])->setFourier(in.rings[iring]);
        (plans.transformers[iring])->inverseFourierTransform(); // fReal points to plans.arrays[iring]
        out.rings.push_back(plans.arrays[iring]);
    }
    out.mode = in.mode;
    out.ring_radius = in.ring_radius;
}

normalizedPolarFourierTransform()

void normalizedPolarFourierTransform	(	Polar< double > &	polarIn,
		Polar< std::complex< double > > &	out,
		bool	conjugated,
		Polar_fftw_plans *&	plans
	)

Definition at line 197 of file polar.cpp.

                                  {
    double mean;
    double stddev;
    polarIn.computeAverageAndStddev(mean, stddev);
    polarIn.normalize(mean, stddev);
    if (plans == nullptr) {
        plans = new Polar_fftw_plans();
        polarIn.calculateFftwPlans(*plans);
    }
    fourierTransformRings(polarIn, out, *plans, conjugated);
}

polarFourierTransform() [1/2]

template<bool NORMALIZE>

void polarFourierTransform	(	const MultidimArray< double > &	in,
		Polar< double > &	inAux,
		Polar< std::complex< double > > &	out,
		bool	conjugated,
		int	first_ring,
		int	last_ring,
		Polar_fftw_plans *&	plans,
		int	BsplineOrder
	)

Compute a polar Fourier transform (with/out normalization) of the input image. If plans is NULL, they are computed and returned.

Definition at line 153 of file polar.cpp.

                                                                   {
    if (BsplineOrder == 1)
        inAux.getPolarFromCartesianBSpline(in, first_ring, last_ring, 1);
    else {
        MultidimArray<double> Maux;
        produceSplineCoefficients(3, Maux, in);
        inAux.getPolarFromCartesianBSpline(Maux, first_ring, last_ring,
                BsplineOrder);
    }
    if (NORMALIZE) {
        double mean;
        double stddev;
        inAux.computeAverageAndStddev(mean, stddev);
        inAux.normalize(mean, stddev);
    }
    if (plans == nullptr) {
        plans = new Polar_fftw_plans();
        inAux.calculateFftwPlans(*plans);
    }
    fourierTransformRings(inAux, out, *plans, conjugated);
}

polarFourierTransform() [2/2]

template<bool NORMALIZE>

void polarFourierTransform	(	const MultidimArray< double > &	in,
		Polar< std::complex< double > > &	out,
		bool	flag,
		int	first_ring,
		int	last_ring,
		Polar_fftw_plans *&	plans,
		int	BsplineOrder = 3
	)

Definition at line 183 of file polar.cpp.

                                                                   {
    Polar<double> polarIn;
    polarFourierTransform<NORMALIZE>(in, polarIn, out, conjugated, first_ring, last_ring, plans, BsplineOrder);
}

rotationalCorrelation()

void rotationalCorrelation	(	const Polar< std::complex< double > > &	M1,
		const Polar< std::complex< double > > &	M2,
		MultidimArray< double > &	angles,
		RotationalCorrelationAux &	aux
	)

Fourier-space rotational Cross-Correlation Function

This function returns the rotational cross-correlation function of two complex Polars M1 and M2 using the cross-correlation convolution theorem.

Note that the local_transformer should have corr (with the right size) already in its fReal, and a Fourier Transform already calculated

MultidimArray<double> angles, corr;
FourierTransformer local_transformer;
Polar_fftw_plans plans;
Polar<std::complex<double> > F1, F2;

// M1 and M2 are already Polar<double>

M1.calculateFftwPlans(plans);
fourierTransformRings(M1, F1, plans, false);
fourierTransformRings(M1, F2, plans, true);// complex conjugated
corr.resize(M1.getSampleNoOuterRing());
local_transformer.setReal(corr, true);
local_transformer.FourierTransform();
rotationalCorrelation(F1,F2,angles,corr,local_transformer);

Note that this function can be used for real-space and fourier-space correlations!!

Definition at line 99 of file polar.cpp.

                                       {
    int nrings = M1.getRingNo();
    if (nrings != M2.getRingNo()) {
        char errorMsg[256];
        sprintf(
                errorMsg,
                "rotationalCorrelation: polar structures have unequal number of rings:\
        nrings %d and M2.getRingNo %d",
                nrings, M2.getRingNo());
        std::cerr << errorMsg << std::endl;
        REPORT_ERROR(ERR_VALUE_INCORRECT, errorMsg);
    }
    // Fsum should already be set with the right size in the local_transformer
    // (i.e. through a FourierTransform of corr)
    aux.local_transformer.getFourierAlias(aux.Fsum);
    memset((double*) MULTIDIM_ARRAY(aux.Fsum), 0,
            2 * MULTIDIM_SIZE(aux.Fsum) * sizeof(double));

    // Multiply M1 and M2 over all rings and sum
    // Assume M2 is already complex conjugated!
    std::complex<double> auxC;
    for (int iring = 0; iring < nrings; iring++) {
        double w = (2. * PI * M1.ring_radius[iring]);
        int imax = M1.getSampleNo(iring);
        const MultidimArray<std::complex<double> > &M1_iring = M1.rings[iring];
        const MultidimArray<std::complex<double> > &M2_iring = M2.rings[iring];
        auto *ptr1 = (double*) MULTIDIM_ARRAY(M1_iring);
        auto *ptr2 = (double*) MULTIDIM_ARRAY(M2_iring);
        auto *ptrFsum = (double *) MULTIDIM_ARRAY(aux.Fsum);
        for (int i = 0; i < imax; i++) {
            double a = *ptr1++;
            double b = *ptr1++;
            double c = *ptr2++;
            double d = *ptr2++;
            *(ptrFsum++) += w * (a * c - b * d);
            *(ptrFsum++) += w * (b * c + a * d);
        }
    }

    // Inverse FFT to get real-space correlations
    // The local_transformer should already have corr as setReal!!
    aux.local_transformer.inverseFourierTransform();

    angles.resize(XSIZE(aux.local_transformer.getReal()));
    double Kaux = 360. / XSIZE(angles);
    for (size_t i = 0; i < XSIZE(angles); i++)
        DIRECT_A1D_ELEM(angles,i) = (double) i * Kaux;
}

volume_convertCartesianToCylindrical()

void volume_convertCartesianToCylindrical	(	const MultidimArray< double > &	in,
		MultidimArray< double > &	out,
		double	Rmin,
		double	Rmax,
		double	deltaR,
		float	angMin,
		double	angMax,
		float	deltaAng,
		Matrix1D< double > &	axis
	)

Produce a cylindrical volume from a cartesian volume. You can give the minimum and maximum radius for the interpolation, the delta radius and the delta angle. It is assumed that the input image has the Xmipp origin. Delta Ang must be in radians.

Definition at line 305 of file polar.cpp.

                                                                             {
    auto NAngSteps = (size_t)floor((angMax - angMin) / deltaAng);
    auto NRSteps = (size_t)floor((Rmax - Rmin) / deltaR);
    out.initZeros(ZSIZE(in), NAngSteps, NRSteps);
    STARTINGZ(out) = STARTINGZ(in);
    STARTINGY(out) = STARTINGX(out) = 0;

    // Calculate rotation matrix
    Matrix2D<double> M;
    if (XX(axis)!=0.0 || YY(axis)!=0.0)
    {
        // 1) Calculate rotation axis v=(0,0,1) x axis
        Matrix1D<double> v(3);
        XX(v)=-YY(axis);
        YY(v)=XX(axis);

        // 2) Calculate rotation angle (0,0,1).axis=module(axis).cos(ang)
        double rotAng=acos(ZZ(axis)/axis.module());

        // 3) Calculate matrix
        rotation3DMatrix(RAD2DEG(rotAng), v, M, false);
    }

    Matrix1D<double> p(3);
    Matrix1D<double> pc(3);
    SPEED_UP_temps012;
    for (size_t i = 0; i < NAngSteps; ++i) {
        float angle = angMin + i * deltaAng;
        float s = sin(angle);
        float c = cos(angle);
        for (int k = STARTINGZ(in); k <= FINISHINGZ(in); ++k) {
            ZZ(p)=k;
            for (size_t j = 0; j < NRSteps; ++j) {
                float R = Rmin + j * deltaR;
                YY(p) = R * s;
                XX(p) = R * c;
                if (MAT_XSIZE(M)>0)
                {
                    M3x3_BY_V3x1(pc,M,p);
                    A3D_ELEM(out,k,i,j) = in.interpolatedElement3D(XX(pc),YY(pc),ZZ(pc));
                }
                else
                    A3D_ELEM(out,k,i,j) = in.interpolatedElement3D(XX(p),YY(p),ZZ(p));
            }
        }
    }
}

volume_convertCartesianToSpherical()

void volume_convertCartesianToSpherical	(	const MultidimArray< double > &	in,
		MultidimArray< double > &	out,
		double	Rmin,
		double	Rmax,
		double	deltaR = 1.,
		double	deltaRot = 2.*PI/360.,
		double	deltaTilt = PI/180.
	)

Produce a spherical volume from a cartesian voume. You can give the minimum and maximum radii for the interpolation. It is assumed that the input image has the Xmipp origin. Delta Ang must be in radians.

Definition at line 356 of file polar.cpp.

                                           {
    auto NRotSteps = (size_t)floor(2*PI / deltaRot);
    auto NTiltSteps = (size_t)floor(PI / deltaRot);
    auto NRSteps = (size_t)floor((Rmax - Rmin) / deltaR);
    out.initZeros(NRSteps, NRotSteps, NTiltSteps);
    STARTINGZ(out) = STARTINGY(out) = STARTINGX(out) = 0;

    Matrix1D<double> p(3);
    for (size_t i = 0; i < NRotSteps; ++i) {
        float rot = i * deltaRot;
        float c = cos(rot);
        float s = sin(rot);
        for (size_t j = 0; j < NRotSteps; ++j) {
            float tilt = j * deltaTilt;
            float stilt = sin(tilt);
            float ctilt = cos(tilt);
            float sc = stilt * c;
            float ss = stilt * s;
            for (size_t k = 0; k<NRSteps; ++k) {
                float R = Rmin + k * deltaR;
                ZZ(p)=R*ctilt;
                YY(p)=R*ss;
                XX(p)=R*sc;
                A3D_ELEM(out,k,i,j) = in.interpolatedElement3D(XX(p),YY(p),ZZ(p));
            }
        }
    }
}