#include <rotational_spectrum.h>

Collaboration diagram for Cylindrical_Wave_Decomposition:

Public Member Functions
double	interpolate (MultidimArray< double > &img, double y, double x)
	Interpolate image value (bilinear) More...

void	compute_cwd (MultidimArray< double > &img)
	Compute the Cylindrical Wave decomposition of an image. More...

Public Attributes
int	numin
	Minimum harmonics. More...

int	numax
	Maximum harmonics. More...

double	x0
	Center of symmetry (x) More...

double	y0
	Center of symmetry (x) More...

double	r1
	Minimum integration radius. More...

double	r2
	Maximum integration radius. More...

double	r3
	Integration increment. More...

int	ir
	Ir. More...

MultidimArray< double >	out_ampcos
	Ampcos. More...

MultidimArray< double >	out_ampsin
	Ampsin. More...

Friends
std::ostream &	operator<< (std::ostream &_out, const Cylindrical_Wave_Decomposition &_cwd)
	Show this object. More...

Detailed Description

Cylindrical_Wave_Decomposition class.

Definition at line 39 of file rotational_spectrum.h.

Member Function Documentation

◆ compute_cwd()

void Cylindrical_Wave_Decomposition::compute_cwd ( MultidimArray< double > & img )

Compute the Cylindrical Wave decomposition of an image.

Definition at line 65 of file rotational_spectrum.cpp.

 {
     std::array<double,1281> coseno = {};
     std::array<double,5191> ampcos = {};
     std::array<double,5191> ampsen = {};
 
     double rh;
     rh = XMIPP_MIN(r2, x0 - 4.);
     rh = XMIPP_MIN(rh, y0 - 4.);
     rh = XMIPP_MIN(rh, YSIZE(img) - x0 - 3.);
     rh = XMIPP_MIN(rh, XSIZE(img) - y0 - 3.);
     int ir_internal;
     ir_internal = (int)((rh - r1) / r3 + 1);
     int ind;
     ind = 0;
     int numin1;
     int numax1;
     numin1 = numin + 1;
     numax1 = numax + 1;
     int kk;
     for (kk = numin1;kk <= numax1;kk++)
     {
         int k;
         k = kk - 1;
         double coefca;
         double coefcb;
         double coefsa;
         double coefsb;
         double h;
         int ntot;
         int my;
         int my2;
         int my4;
         int my5;
         int j;
         if (k != 0)
         {
             my = (int)(1 + PI * rh / 2. / k);
             my2 = 2 * my;
             my4 = my * k;
             my5 = my4 - 1;
             ntot = 4 * my4;
             h = 2. * PI / ntot;
             double hdpi;
             hdpi = h / PI;
             double th;
             th = k * h;
             double ys;
             ys = sin(th);
             double zs;
             zs = cos(th);
             double ys2;
             ys2 = sin(2. * th);
             double b1;
             b1 = 2. / (th * th) * (1. + zs * zs - ys2 / th);
             double g1;
             g1 = 4. / (th * th) * (ys / th - zs);
             double d1;
             d1 = 2. * th / 45.;
             double e1;
             e1 = d1 * ys * 2.;
             d1 *= ys2;
             coefca = (b1 + e1) * hdpi;
             coefcb = (g1 - d1) * hdpi;
             coefsa = (b1 - e1) * hdpi;
             coefsb = (g1 + d1) * hdpi;
         }
         else
         {
             my = (int)(1 + PI * rh / 2.);
             my2 = 2 * my;
             my4 = my;
             my5 = my4 - 1;
             ntot = 4 * my4;
             h = 2. * PI / ntot;
             coefca = h / PI / 2.;
         }
         int i;
         for (i = 1;i <= my5;i++)
         {
             double fi;
             fi = i * h;
             coseno[i] = sin(fi);
         }
         coseno[my4] = 1.;
         int my3;
         my3 = 2 * my4;
         for (i = 1;i <= my5;i++)
         {
             coseno[my3-i] = coseno[i];
             coseno[my3+i] = -coseno[i];
             coseno[ntot-i] = -coseno[i];
             coseno[ntot+i] = coseno[i];
         }
         coseno[my3] = 0.;
         coseno[my3+my4] = -1.;
         coseno[ntot] = 0.;
         coseno[ntot+my4] = 1.;
         double r11;
         r11 = r1 - r3;
         double ac;
         double r;
         for (int jr = 1;jr <= ir_internal;jr++)
         {
             ind++;
             r = r11 + r3 * jr;
             ac = 0.;
             int i1c;
             i1c = my4;
             int i1s;
             i1s = 0;
             double x;
             double y;
             double z;
             if (k != 0)
             {
                 double as;
                 double bc;
                 double bs;
                 as = bc = bs = 0.;
                 for (i = 1;i <= k;i++)
                 {
                     int i2c;
                     i2c = my4;
                     int i2s;
                     i2s = 0;
                     for (j = 1;j <= my2;j++)
                     {
                         i1c++;
                         i1s++;
                         i2c += k;
                         i2s += k;
                         x = x0 + r * coseno[i1c];
                         y = y0 + r * coseno[i1s];
                         z = interpolate(img, y, x);
                         bc += z * coseno[i2c];
                         bs += z * coseno[i2s];
                         i1c++;
                         i1s++;
                         i2c += k;
                         i2s += k;
                         x = x0 + r * coseno[i1c];
                         y = y0 + r * coseno[i1s];
                         z = interpolate(img, y, x);
                         ac += z * coseno[i2c];
                         as += z * coseno[i2s];
                     }
                 }
                 ampcos[ind] = coefca * ac + coefcb * bc;
                 ampsen[ind] = coefsa * as + coefsb * bs;
             }
             else
                 for (j = 1;j <= ntot;j++)
                 {
                     i1c++;
                     i1s++;
                     x = x0 + r * coseno[i1c];
                     y = y0 + r * coseno[i1s];
                     z = interpolate(img, y, x);
                     ac += z;
                     ampcos[ind] = coefca * ac;
                     ampsen[ind] = 0.;
                 }
         }
         if (k == 0)
         {
             ac = 0.;
             for (j = 1;j <= ir_internal;j++)
             {
                 r = r11 + j * r3;
                 ac += ampcos[j] * 2. * PI * r;
             }
             ac /= (PI * (r * r - r1 * r1));
             for (j = 1;j <= ir_internal;j++)
                 ampcos[j] -= ac;
         }
     }
 
     out_ampcos.initZeros((numax - numin + 1)*ir_internal);
     out_ampsin = out_ampcos;
     FOR_ALL_ELEMENTS_IN_ARRAY1D(out_ampcos)
     {
         A1D_ELEM(out_ampcos, i) = ampcos[i+1];
         A1D_ELEM(out_ampsin, i) = ampsen[i+1];
     }
 }

◆ interpolate()

double Cylindrical_Wave_Decomposition::interpolate	(	MultidimArray< double > &	img,
		double	y,
		double	x
	)

Interpolate image value (bilinear)

Definition at line 49 of file rotational_spectrum.cpp.

 {
     y = y - 1; // Since Fortran starts indexing at 1
     x = x - 1;
     auto iy = (int)y; /* Trunc the x, y coordinates to int */
     auto ix = (int)x;
     double scale =  y - iy;
     double *ptr_yx = &DIRECT_A2D_ELEM(img, iy, ix);
     double *ptr_y1x = ptr_yx + XSIZE(img);
     double introw1 = *ptr_yx  + scale * (*(ptr_yx + 1)  - *ptr_yx);
     double introw2 = *ptr_y1x + scale * (*(ptr_y1x + 1) - *ptr_y1x);
     return introw1 + (x - ix)*(introw2 - introw1);
 }

Friends And Related Function Documentation

◆ operator<<

std::ostream& operator<<	(	std::ostream &	_out,
		const Cylindrical_Wave_Decomposition &	_cwd
	)

friend

Show this object.

Definition at line 31 of file rotational_spectrum.cpp.

 {
     _out << "ir=" << _cwd.ir << std::endl
     << "numin=" << _cwd.numin << std::endl
     << "numax=" << _cwd.numax << std::endl
     << "x0=" << _cwd.x0 << std::endl
     << "y0=" << _cwd.y0 << std::endl
     << "r1=" << _cwd.r1 << std::endl
     << "r2=" << _cwd.r2 << std::endl
     << "r3=" << _cwd.r3 << std::endl;
     FOR_ALL_ELEMENTS_IN_ARRAY1D(_cwd.out_ampcos)
     _out << _cwd.out_ampcos(i) << " "
     << _cwd.out_ampsin(i) << std::endl;
     return _out;
 }

Member Data Documentation

◆ ir

int Cylindrical_Wave_Decomposition::ir

Ir.

Definition at line 64 of file rotational_spectrum.h.

◆ numax

int Cylindrical_Wave_Decomposition::numax

Maximum harmonics.

Definition at line 46 of file rotational_spectrum.h.

◆ numin

int Cylindrical_Wave_Decomposition::numin

Minimum harmonics.

Definition at line 43 of file rotational_spectrum.h.

◆ out_ampcos

MultidimArray< double > Cylindrical_Wave_Decomposition::out_ampcos

Ampcos.

Definition at line 67 of file rotational_spectrum.h.

◆ out_ampsin

MultidimArray< double > Cylindrical_Wave_Decomposition::out_ampsin

Ampsin.

Definition at line 70 of file rotational_spectrum.h.

◆ r1

double Cylindrical_Wave_Decomposition::r1

Minimum integration radius.

Definition at line 55 of file rotational_spectrum.h.

◆ r2

double Cylindrical_Wave_Decomposition::r2

Maximum integration radius.

Definition at line 58 of file rotational_spectrum.h.

◆ r3

double Cylindrical_Wave_Decomposition::r3

Integration increment.

Definition at line 61 of file rotational_spectrum.h.

◆ x0

double Cylindrical_Wave_Decomposition::x0

Center of symmetry (x)

Definition at line 49 of file rotational_spectrum.h.

◆ y0

double Cylindrical_Wave_Decomposition::y0