Polar< T > Class Template Reference

#include <polar.h>

Collaboration diagram for Polar< T >:

Public Member Functions
	Polar ()
	Polar (const Polar &P)
	~Polar ()
Polar &	operator= (const Polar &P)
Polar &	operator- (const T val) const
Polar &	operator+ (const T val)
Polar &	operator* (const T val)
Polar &	operator/ (const T val)
void	operator-= (const T val)
void	operator+= (const T val)
void	operator*= (const T val)
void	operator/= (const T val)
Polar &	operator-= (const Polar< T > in)
Polar &	operator+= (const Polar< T > in)
Polar &	operator*= (const Polar< T > in)
Polar &	operator/= (const Polar< T > in)
void	rename (const FileName &newName)
void	clear ()
const FileName &	name () const
int	getRingNo () const
int	getMode () const
double	getOversample () const
int	getSampleNo (int iring) const
int	getSampleNoOuterRing () const
double	getRadius (int iring) const
MultidimArray< T > &	getRing (int i)
const MultidimArray< T > &	getRing (int i) const
void	setRing (int i, MultidimArray< T > val)
T &	getPixel (int r, int f) const
T &	operator() (int r, int f) const
void	setPixel (int r, int f, T val) const
void	computeAverageAndStddev (double &avg, double &stddev, int mode=FULL_CIRCLES) const
void	normalize (double average, double stddev)
void	getCartesianCoordinates (std::vector< double > &x, std::vector< double > &y, std::vector< T > &data, const int extra_shell=GRIDDING_K/2)
void	getPolarFromCartesianBSpline (const MultidimArray< T > &M1, int first_ring, int last_ring, int BsplineOrder=3, double xoff=0., double yoff=0., double oversample1=1., int mode1=FULL_CIRCLES)
void	calculateFftwPlans (Polar_fftw_plans &out)

Public Attributes
int	mode
double	oversample
std::vector< double >	ring_radius
std::vector< MultidimArray< T > >	rings

Protected Attributes
FileName	fn_pol

Detailed Description

template<typename T>
class Polar< T >

Class for polar coodinates

Definition at line 67 of file polar.h.

Constructor & Destructor Documentation

Polar() [1/2]

template<typename T>

Polar< T >::Polar ( )

inline

Empty constructor

An empty image with size 0x0 is created.

Polar P;

Definition at line 85 of file polar.h.

    {
        fn_pol = "";
        ring_radius.clear();
        rings.clear();
        mode = FULL_CIRCLES;
        oversample = 1.;
    }

Polar() [2/2]

template<typename T>

Polar< T >::Polar ( const Polar< T > &  P )

inline

Copy constructor

Polar P2(P1);

Definition at line 100 of file polar.h.

    {
        rings = P.rings;
        fn_pol = P.fn_pol;
        ring_radius = P.ring_radius;
        mode = P.mode;
        oversample = P.oversample;
    }

~Polar()

template<typename T>

Polar< T >::~Polar ( )

inline

Destructor.

Definition at line 111 of file polar.h.

    {
        rings.clear();
        ring_radius.clear();
    }

Member Function Documentation

calculateFftwPlans()

template<typename T>

void Polar< T >::calculateFftwPlans ( Polar_fftw_plans &  out )

inline

Precalculate a vector with FFTW plans for all rings

Definition at line 708 of file polar.h.

    {
        (out.arrays).resize(rings.size());
        for (size_t iring = 0; iring < rings.size(); iring++)
        {
            (out.arrays)[iring] = rings[iring];
            auto *ptr_transformer = new FourierTransformer();
            ptr_transformer->setReal((out.arrays)[iring]);
            out.transformers.push_back(ptr_transformer);
        }
    }

clear()

template<typename T>

void Polar< T >::clear ( )

inline

Empty polar

This function clears the polar to an empty vector without name.

P.clear();

Definition at line 295 of file polar.h.

    {
        fn_pol = "";
        rings.clear();
        ring_radius.clear();
        mode = FULL_CIRCLES;
        oversample = 1.;
    }

computeAverageAndStddev()

template<typename T>

void Polar< T >::computeAverageAndStddev ( double &  avg, double &  stddev, int  mode = FULL_CIRCLES  ) const

inline

Compute sum or average of all pixels in polar rings.

Definition at line 488 of file polar.h.

    {
        double aux;
        double sum = 0.;
        double sum2=0.;
        double twopi;
        double w;
        double N = 0;

        if (mode == FULL_CIRCLES)
            twopi = 2.*PI;
        else if (mode == HALF_CIRCLES)
            twopi = PI;
        else
            REPORT_ERROR(ERR_VALUE_INCORRECT,"Incorrect mode for computeAverageAndStddev");

        size_t imax=rings.size();
        const double *ptrRing_radius=&ring_radius[0];
        for (size_t i = 0; i < imax; ++i, ++ptrRing_radius)
        {
            // take varying sampling into account
            w = (twopi * (*ptrRing_radius)) / (double) XSIZE(rings[i]);
            const MultidimArray<T> &rings_i=rings[i];
            for (size_t j = 0; j < XSIZE(rings_i); j++)
            {
                aux = DIRECT_A1D_ELEM(rings_i,j);
                double waux=w*aux;
                sum += waux;
                sum2 += waux * aux;
                N += w;
            }
        }
        if (N>0)
        {
            sum2 = sum2 / N;
            avg = sum / N;
            stddev=sqrt(fabs(sum2-avg*avg));
        }
        else if (N != 0.)
        {
            avg = sum;
            stddev=sqrt(fabs(sum2-avg*avg));
        }
        else
            stddev=avg=0;
    }

getCartesianCoordinates()

template<typename T>

void Polar< T >::getCartesianCoordinates ( std::vector< double > &  x, std::vector< double > &  y, std::vector< T > &  data, const int  extra_shell = GRIDDING_K/2  )

inline

Definition at line 559 of file polar.h.

    {
        // Only for full circles for now!
        if (mode != FULL_CIRCLES)
            REPORT_ERROR(ERR_VALUE_INCORRECT,"VoronoiArea only implemented for FULL_CIRCLES mode of Polar");

        // First fill the vector with the originally sampled coordinates
        x.clear();
        y.clear();
        data.clear();
        for (int i = 0; i < rings.size(); i++)
        {
            int nsam = XSIZE(rings[i]);
            float dphi = TWOPI / (float)nsam;
            float radius = ring_radius[i];
            for (int j = 0; j < nsam; j++)
            {
                float tmp = j * dphi;
                x.push_back(radius * sin(tmp));
                y.push_back(radius * cos(tmp));
                data.push_back(rings[i](j));
            }
        }

        // Add additional points on the inside and outside of the rings
        // Add a maximum of "extra_shell" rings
        // Set data to zero here
        auto first_ring  = (int)floor(ring_radius[0]);
        auto last_ring   = (int)ceil(ring_radius[rings.size()-1]);
        int outer       = last_ring + extra_shell;
        int inner       = XMIPP_MAX(0,first_ring - extra_shell);
        for (int iradius = 0; iradius < outer; iradius +=1)
        {
            if ( (iradius >= inner && iradius < first_ring) ||
                 (iradius <= outer && iradius > last_ring) )
            {
                double radius=iradius;
                int nsam = 2 * (int)( 0.5 * oversample * TWOPI * radius );
                nsam = XMIPP_MAX(1, nsam);
                float dphi = TWOPI / (float)nsam;
                for (int j = 0; j < nsam; j++)
                {
                    float tmp = j * dphi;
                    x.push_back(radius * sin(tmp));
                    y.push_back(radius * cos(tmp));
                    data.push_back(0.);
                }
            }
        }
    }

getMode()

template<typename T>

int Polar< T >::getMode ( ) const

inline

Mode access

This function is used to know the "mode" of the polar.

There are two modes: FULL_CIRCLES = 0 (used for asymmetric functions) HALF_CIRCLES = 0 (used for symmetric functions, e.g. Fourier transforms)

std::cout << "Mode: " << P.getMode() << std::endl;

Definition at line 343 of file polar.h.

    {
        return mode;
    }

getOversample()

template<typename T>

double Polar< T >::getOversample ( ) const

inline

Oversample access

This function is used to know the oversampling factor of the polar.

Oversample = 1 means there is no oversampling Oversample > 1 means oversampling Oversampling < 1 means undersampling

std::cout << "Oversample: " << P.getOversample() << std::endl;

Definition at line 360 of file polar.h.

    {
        return oversample;
    }

getPixel()

template<typename T>

T& Polar< T >::getPixel ( int  r, int  f  ) const

inline

Pixel access

This operator is used to access a pixel within the polar. That means, given the ring and the number of the pixel in the rotational direction.

int ring = 1, phi = 0;
std::cout <<" first pixel of second ring= "<<P.getPixel(ring,phi)<<std::endl;

Definition at line 446 of file polar.h.

    {
        return rings[r](f);
    }

getPolarFromCartesianBSpline()

template<typename T>

void Polar< T >::getPolarFromCartesianBSpline ( const MultidimArray< T > &  M1, int  first_ring, int  last_ring, int  BsplineOrder = 3, double  xoff = 0., double  yoff = 0., double  oversample1 = 1., int  mode1 = FULL_CIRCLES  )

inline

Convert cartesian MultidimArray to Polar using B-spline interpolation

The input MultidimArray is assumed to be pre-processed for B-splines

Polar P;
MultidimArray<double> Maux;
img().produceSplineCoefficients(Maux,3);
P.getPolarFromCartesianBSpline(Maux,1,15);

Definition at line 625 of file polar.h.

    {
        double twopi;
        double xp;
        double yp;
        double minxp;
        double maxxp;
        double minyp;
        double maxyp;

        auto noOfRings = getNoOfRings(first_ring, last_ring);
        rings.resize(noOfRings);
        ring_radius.resize(noOfRings);

        mode = mode1;
        oversample = oversample1;

        if (mode == FULL_CIRCLES)
            twopi = TWOPI;
        else if (mode == HALF_CIRCLES)
            twopi = PI;
        else
            REPORT_ERROR(ERR_VALUE_INCORRECT,"Incorrect mode for getPolarFromCartesian");


        // Limits of the matrix (not oversized!)
        minxp = FIRST_XMIPP_INDEX(XSIZE(M1));
        minyp = FIRST_XMIPP_INDEX(YSIZE(M1));
        maxxp = LAST_XMIPP_INDEX(XSIZE(M1));
        maxyp = LAST_XMIPP_INDEX(YSIZE(M1));
        double minxp_e=minxp-XMIPP_EQUAL_ACCURACY;
        double minyp_e=minyp-XMIPP_EQUAL_ACCURACY;
        double maxxp_e=maxxp+XMIPP_EQUAL_ACCURACY;
        double maxyp_e=maxyp+XMIPP_EQUAL_ACCURACY;

        std::lock_guard<std::mutex> lockGuard(m_mutex);

        ensureAngleCache(first_ring, last_ring);

        // Loop over all polar coordinates
        for (int r = 0; r < noOfRings; ++r)
        {
            float radius = (float) r + first_ring;
            ring_radius.at(r) = radius;
            auto &Mring = rings.at(r);
            // Non-constant sampling!! (always even for convenient Half2Whole of FTs)
            int nsam = getNoOfSamples(twopi, radius);
            Mring.resizeNoCopy(nsam);

            auto &s = m_lastSinRadius.at(r);
            auto &c = m_lastCosRadius.at(r);

            for (int sample = 0; sample < nsam; ++sample)
            {
                // from polar to original cartesian coordinates
                xp = s.at(sample);
                yp = c.at(sample);

                // Origin offsets
                xp += xoff;
                yp += yoff;

                // Wrap coordinates
                if (xp < minxp_e || xp > maxxp_e)
                    xp = realWRAP(xp, minxp - 0.5, maxxp + 0.5);
                if (yp < minyp_e || yp > maxyp_e)
                    yp = realWRAP(yp, minyp - 0.5, maxyp + 0.5);

                // Perform the convolution interpolation
                if (BsplineOrder==1)
                    DIRECT_A1D_ELEM(Mring,sample) = M1.interpolatedElement2DOutsideZero(xp,yp);
                else
                    DIRECT_A1D_ELEM(Mring,sample) = M1.interpolatedElementBSpline2D(xp,yp,BsplineOrder);
            }
        }
    }

getRadius()

template<typename T>

double Polar< T >::getRadius ( int  iring ) const

inline

The radius of each ring access

This function is used to know the radius of a given ring.

std::cout << "Radius of second ring: " << P.getRadius(1) << std::endl;

Definition at line 399 of file polar.h.

    {
        return ring_radius[iring];
    }

getRing() [1/2]

template<typename T>

MultidimArray< T >& Polar< T >::getRing ( int  i )

inline

1D Matrix access

This operator can be used to access any ring of the polar as MultidimArray.

MultidimArray<double> secondring = P.getRing(1);

Definition at line 412 of file polar.h.

    {
        return rings[i];
    }

getRing() [2/2]

template<typename T>

const MultidimArray< T >& Polar< T >::getRing ( int  i ) const

inline

Definition at line 416 of file polar.h.

    {
        return rings[i];
    }

getRingNo()

template<typename T>

int Polar< T >::getRingNo ( ) const

inline

Number of rings access

This function is used to know the number of rings in the polar.

std::cout << "Number of rings: " << P.getRingNo() << std::endl;

Definition at line 326 of file polar.h.

    {
        return rings.size();
    }

getSampleNo()

template<typename T>

int Polar< T >::getSampleNo ( int  iring ) const

inline

Number of samples in each ring access

This function is used to know the number of samples in a given ring.

std::cout << "Number of samples in second ring: " << P.getSampleNo(1) << std::endl;

Definition at line 373 of file polar.h.

    {
        return XSIZE(rings[iring]);
    }

getSampleNoOuterRing()

template<typename T>

int Polar< T >::getSampleNoOuterRing ( ) const

inline

Number of samples in outer ring access

This function is used to know the number of samples in the outer ring.

std::cout << "Number of samples in outer ring: " << P.getSampleNoOuterRing() << std::endl;

Definition at line 386 of file polar.h.

    {
        return XSIZE(rings[rings.size()-1]);
    }

name()

template<typename T>

const FileName& Polar< T >::name ( ) const

inline

Name access

This function is used to know the name of the polar. It cannot be used to assign a new one. You may use rename() for that.

std::cout << "Polar name: " << P.name() << std::endl;

Definition at line 313 of file polar.h.

    {
        return fn_pol;
    }

normalize()

template<typename T>

void Polar< T >::normalize ( double  average, double  stddev  )

inline

Definition at line 536 of file polar.h.

    {
        size_t imax=rings.size();
        double istddev=1.0/stddev;
        for (size_t i = 0; i < imax; i++)
        {
            MultidimArray<T> &rings_i=rings[i];
            for (size_t j = 0; j < XSIZE(rings_i); j++)
                DIRECT_A1D_ELEM(rings_i,j) = (DIRECT_A1D_ELEM(rings_i,j)-average)*istddev;
        }
    }

operator()()

template<typename T>

T& Polar< T >::operator() ( int  r, int  f  ) const

inline

Pixel access

This operator is used to access a pixel within the polar. That means, given the ring and the number of the pixel in the rotational direction.

int ring = 1, phi = 0;
std::cout <<" first pixel of second ring= "<<P(ring,phi)<<std::endl;

Definition at line 462 of file polar.h.

    {
        return DIRECT_A1D_ELEM(rings[r],f);
    }

operator*()

template<typename T>

Polar& Polar< T >::operator* ( const T  val )

inline

Multiply by a constant pixel-by-pixel

Definition at line 160 of file polar.h.

    {
        Polar<T> result;
        result = *(this);
        for (size_t i = 0; i < rings.size(); i++)
            for (size_t j = 0; j < XSIZE(rings[i]); j++)
                result(i,j) *= val;

        return result;
    }

operator*=() [1/2]

template<typename T>

void Polar< T >::operator*= ( const T  val )

inline

Multiply by a constant pixel-by-pixel

Definition at line 212 of file polar.h.

    {
        size_t imax=rings.size();
        for (size_t i = 0; i < imax; i++)
        {
            MultidimArray<T> &rings_i=rings[i];
            for (int j = 0; j < XSIZE(rings_i); j++)
                DIRECT_A1D_ELEM(rings_i,j) *= val;
        }
    }

operator*=() [2/2]

template<typename T>

Polar& Polar< T >::operator*= ( const Polar< T >  in )

inline

Multiply two polars pixel-by-pixel

Definition at line 257 of file polar.h.

    {
        for (size_t i = 0; i < rings.size(); i++)
            for (size_t j = 0; j < XSIZE(rings[i]); j++)
                (*(this))(i,j) *= in(i,j);
    }

operator+()

template<typename T>

Polar& Polar< T >::operator+ ( const T  val )

inline

Add a constant pixel-by-pixel

Definition at line 147 of file polar.h.

    {
        Polar<T> result;
        result = *(this);
        for (size_t i = 0; i < rings.size(); i++)
            for (size_t j = 0; j < XSIZE(rings[i]); j++)
                result(i,j) += val;

        return result;
    }

operator+=() [1/2]

template<typename T>

void Polar< T >::operator+= ( const T  val )

inline

Add a constant pixel-by-pixel

Definition at line 199 of file polar.h.

    {
        size_t imax=rings.size();
        for (size_t i = 0; i < imax; i++)
        {
            MultidimArray<T> &rings_i=rings[i];
            for (int j = 0; j < XSIZE(rings_i); j++)
                DIRECT_A1D_ELEM(rings_i,j) += val;
        }
    }

operator+=() [2/2]

template<typename T>

Polar& Polar< T >::operator+= ( const Polar< T >  in )

inline

Add two polars pixel-by-pixel

Definition at line 248 of file polar.h.

    {
        for (size_t i = 0; i < rings.size(); i++)
            for (size_t j = 0; j < XSIZE(rings[i]); j++)
                (*(this))(i,j) += in(i,j);
    }

operator-()

template<typename T>

Polar& Polar< T >::operator- ( const T  val ) const

inline

Subtract a constant pixel-by-pixel

Definition at line 134 of file polar.h.

    {
        Polar<T> result;
        result = *(this);
        for (size_t i = 0; i < rings.size(); i++)
            for (size_t j = 0; j < XSIZE(rings[i]); j++)
                result(i,j) -= val;

        return result;
    }

operator-=() [1/2]

template<typename T>

void Polar< T >::operator-= ( const T  val )

inline

Subtract a constant pixel-by-pixel

Definition at line 186 of file polar.h.

    {
        size_t imax=rings.size();
        for (size_t i = 0; i < imax; i++)
        {
            MultidimArray<T> &rings_i=rings[i];
            for (size_t j = 0; j < XSIZE(rings_i); j++)
                DIRECT_A1D_ELEM(rings_i,j) -= val;
        }
    }

operator-=() [2/2]

template<typename T>

Polar& Polar< T >::operator-= ( const Polar< T >  in )

inline

Subtract two polars pixel-by-pixel

Definition at line 239 of file polar.h.

    {
        for (size_t i = 0; i < rings.size(); i++)
            for (size_t j = 0; j < XSIZE(rings[i]); j++)
                (*(this))(i,j) -= in(i,j);
    }

operator/()

template<typename T>

Polar& Polar< T >::operator/ ( const T  val )

inline

Divide by a constant pixel-by-pixel

Definition at line 173 of file polar.h.

    {
        Polar<T> result;
        result = *(this);
        for (size_t i = 0; i < rings.size(); i++)
            for (size_t j = 0; j < XSIZE(rings[i]); j++)
                result(i,j) /= val;

        return result;
    }

operator/=() [1/2]

template<typename T>

void Polar< T >::operator/= ( const T  val )

inline

Divide by a constant pixel-by-pixel

Definition at line 225 of file polar.h.

    {
        size_t imax=rings.size();
        double ival=1.0/val;
        for (size_t i = 0; i < imax; i++)
        {
            MultidimArray<T> &rings_i=rings[i];
            for (size_t j = 0; j < XSIZE(rings_i); j++)
                DIRECT_A1D_ELEM(rings_i,j) *= ival;
        }
    }

operator/=() [2/2]

template<typename T>

Polar& Polar< T >::operator/= ( const Polar< T >  in )

inline

Divide two polars pixel-by-pixel

Definition at line 266 of file polar.h.

    {
        for (size_t i = 0; i < rings.size(); i++)
            for (size_t j = 0; j < XSIZE(rings[i]); j++)
                (*(this))(i,j) /= in(i,j);
    }

operator=()

template<typename T>

Polar& Polar< T >::operator= ( const Polar< T > &  P )

inline

Assignment

Definition at line 119 of file polar.h.

    {
        if (this != &P)
        {
            fn_pol = P.fn_pol;
            rings = P.rings;
            ring_radius = P.ring_radius;
            mode = P.mode;
            oversample = P.oversample;
        }
        return *this;
    }

rename()

template<typename T>

void Polar< T >::rename ( const FileName &  newName )

inline

Rename polar

Give a new name to the polar.

P.rename("new_name");

Definition at line 282 of file polar.h.

    {
        fn_pol = newName;
    }

setPixel()

template<typename T>

void Polar< T >::setPixel ( int  r, int  f, T  val  ) const

inline

Pixel access

This operator is used to set a pixel within the polar. That means, given the ring and the number of the pixel in the rotational direction.

int ring = 1, phi = 0;
double val = 1.2;
P.setPixel(ring,phi,val);

Definition at line 479 of file polar.h.

    {
        rings[r](f) = val;
    }

setRing()

template<typename T>

void Polar< T >::setRing ( int  i, MultidimArray< T >  val  )

inline

1D Matrix access

This operator can be used to set any ring of the polar with a MultidimArray.

MultidimArray<double> ring;
P.setRing(1,ring);

Definition at line 430 of file polar.h.

    {
        rings[i] = val;
    }

Member Data Documentation

fn_pol

template<typename T>

FileName Polar< T >::fn_pol

protected

Definition at line 70 of file polar.h.

mode

template<typename T>

int Polar< T >::mode

Definition at line 72 of file polar.h.

oversample

template<typename T>

double Polar< T >::oversample

Definition at line 73 of file polar.h.

ring_radius

template<typename T>

std::vector<double> Polar< T >::ring_radius

Definition at line 74 of file polar.h.

rings

template<typename T>

std::vector<MultidimArray<T> > Polar< T >::rings

Definition at line 75 of file polar.h.

---

The documentation for this class was generated from the following files:

• xmipp/libraries/data/polar.h
• xmipp/libraries/data/polar.cpp