EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Member Functions | List of all members
EMAN::TransformProcessor Class Reference
a function or class that is CUDA enabled

Transform the image using a Transform object. More...

#include <processor.h>

Inheritance diagram for EMAN::TransformProcessor:
Inheritance graph
[legend]
Collaboration diagram for EMAN::TransformProcessor:
Collaboration graph
[legend]

Public Member Functions

virtual string get_name () const
 Get the processor's name. More...
 
virtual void process_inplace (EMData *image)
 
virtual EMData * process (const EMData *const image)
 
virtual TypeDict get_param_types () const
 Get processor parameter information in a dictionary. More...
 
virtual string get_desc () const
 Get the descrition of this specific processor. More...
 
float * transform (const EMData *const image, const Transform &t) const
 
- Public Member Functions inherited from EMAN::Processor
virtual ~Processor ()
 
virtual void process_list_inplace (vector< EMData * > &images)
 To process multiple images using the same algorithm. More...
 
virtual Dict get_params () const
 Get the processor parameters in a key/value dictionary. More...
 
virtual void set_params (const Dict &new_params)
 Set the processor parameters using a key/value dictionary. More...
 

Static Public Member Functions

static Processor * NEW ()
 
- Static Public Member Functions inherited from EMAN::Processor
static string get_group_desc ()
 Get the description of this group of processors. More...
 
static void EMFourierFilterInPlace (EMData *fimage, Dict params)
 Compute a Fourier-filter processed image in place. More...
 
static EMData * EMFourierFilter (EMData *fimage, Dict params)
 Compute a Fourier-processor processed image without altering the original image. More...
 

Static Public Attributes

static const string NAME = "xform"
 

Private Member Functions

void assert_valid_aspect (const EMData *const image) const
 

Additional Inherited Members

- Public Types inherited from EMAN::Processor
enum  fourier_filter_types {
  TOP_HAT_LOW_PASS , TOP_HAT_HIGH_PASS , TOP_HAT_BAND_PASS , TOP_HOMOMORPHIC ,
  GAUSS_LOW_PASS , GAUSS_HIGH_PASS , GAUSS_BAND_PASS , GAUSS_INVERSE ,
  GAUSS_HOMOMORPHIC , BUTTERWORTH_LOW_PASS , BUTTERWORTH_HIGH_PASS , BUTTERWORTH_HOMOMORPHIC ,
  KAISER_I0 , KAISER_SINH , KAISER_I0_INVERSE , KAISER_SINH_INVERSE ,
  SHIFT , TANH_LOW_PASS , TANH_HIGH_PASS , TANH_HOMOMORPHIC ,
  TANH_BAND_PASS , RADIAL_TABLE , CTF_
}
 Fourier filter Processor type enum. More...
 
- Protected Attributes inherited from EMAN::Processor
Dict params
 

Detailed Description

Transform the image using a Transform object.

Author
David Woolford
Date
September 2008
Parameters
transformThe Transform object that will be applied to the image

Definition at line 2825 of file processor.h.

Member Function Documentation

◆ assert_valid_aspect()

void TransformProcessor::assert_valid_aspect ( const EMData *const  image) const
private

Definition at line 11620 of file processor.cpp.

11620 {
11621 int ndim = image->get_ndim();
11622 if (ndim != 2 && ndim != 3) throw ImageDimensionException("Transforming an EMData only works if it's 2D or 3D");
11623
11624// if (! params.has_key("transform") ) throw InvalidParameterException("You must specify a Transform in order to perform this operation");
11625}
ImageDimensionException
#define ImageDimensionException(desc)
Definition: exception.h:166

References ImageDimensionException.

Referenced by process(), and process_inplace().

◆ get_desc()

virtual string EMAN::TransformProcessor::get_desc ( ) const
inlinevirtual

Get the descrition of this specific processor.

This function must be overwritten by a subclass.

Returns
The description of this processor.

Implements EMAN::Processor.

Definition at line 2864 of file processor.h.

2865 {
2866 return "The image is transformed using transform parameter, or alternatively specific numbers alpha,tx,ty or az,alt,phi,tx,ty,tz (shortcut for convenience)";
2867 }

◆ get_name()

virtual string EMAN::TransformProcessor::get_name ( ) const
inlinevirtual

Get the processor's name.

Each processor is identified by a unique name.

Returns
The processor's name.

Implements EMAN::Processor.

Definition at line 2828 of file processor.h.

2829 {
2830 return NAME;
2831 }
EMAN::TransformProcessor::NAME
static const string NAME
Definition: processor.h:2869

References NAME.

◆ get_param_types()

virtual TypeDict EMAN::TransformProcessor::get_param_types ( ) const
inlinevirtual

Get processor parameter information in a dictionary.

Each parameter has one record in the dictionary. Each record contains its name, data-type, and description.

Returns
A dictionary containing the parameter info.

Reimplemented from EMAN::Processor.

Definition at line 2849 of file processor.h.

2850 {
2851 TypeDict d;
2852 d.put("transform", EMObject::TRANSFORM, "The Transform object that will be applied to the image" );
2853 d.put("alpha", EMObject::FLOAT, "2-D alpha angle" );
2854 d.put("az", EMObject::FLOAT, "3-D Azimuth" );
2855 d.put("alt", EMObject::FLOAT, "3-D Altitude" );
2856 d.put("phi", EMObject::FLOAT, "3-D Phi" );
2857 d.put("tx", EMObject::FLOAT, "x translation" );
2858 d.put("ty", EMObject::FLOAT, "y translation" );
2859 d.put("tz", EMObject::FLOAT, "y translation" );
2860 d.put("zerocorners",EMObject::INT,"If set, corners (anything beyond radius/2-1) may be zeroed out in real or Fourier space. This will produce a considerable speedup in Fourier rotations. ");
2861 return d;
2862 }
EMAN::TypeDict
TypeDict is a dictionary to store <string, EMObject::ObjectType> pair.
Definition: emobject.h:305
EMAN::TypeDict::put
void put(const string &key, EMObject::ObjectType o, const string &desc="")
Definition: emobject.h:330

References EMAN::EMObject::FLOAT, EMAN::EMObject::INT, EMAN::TypeDict::put(), and EMAN::EMObject::TRANSFORM.

◆ NEW()

static Processor * EMAN::TransformProcessor::NEW ( )
inlinestatic

Definition at line 2832 of file processor.h.

2833 {
2834 return new TransformProcessor();
2835 }
EMAN::TransformProcessor
Transform the image using a Transform object.
Definition: processor.h:2826

◆ process()

EMData * TransformProcessor::process ( const EMData *const  image)
virtual
Exceptions
ImageDimensionExceptionif the image is not 2D or 3D
InvalidParameterExceptionif the Transform parameter is not specified

Reimplemented from EMAN::Processor.

Definition at line 11659 of file processor.cpp.

11659 {
11660 ENTERFUNC;
11661
11662 assert_valid_aspect(image);
11663
11664 Transform t;
11665 if (params.has_key("transform")) t=*(Transform *)params["transform"];
11666 else {
11667 if (params.has_key("alpha")) params["phi"]=params["alpha"];
11668 float az=params.set_default("az",0.0f);
11669 float alt=params.set_default("alt",0.0f);
11670 float phi=params.set_default("phi",0.0f);
11671 float tx=params.set_default("tx",0.0f);
11672 float ty=params.set_default("ty",0.0f);
11673 float tz=params.set_default("tz",0.0f);
11674
11675 t.set_rotation(Dict("type","eman","az",az,"alt",alt,"phi",phi));
11676 t.set_trans(tx,ty,tz);
11677 }
11678
11679 EMData* p = 0;
11680#ifdef EMAN2_USING_CUDA
11681 if(EMData::usecuda == 1 && image->isrodataongpu()){
11682 //cout << "using CUDA xform" << endl;
11683 p = new EMData(0,0,image->get_xsize(),image->get_ysize(),image->get_zsize(),image->get_attr_dict());
11684 float * m = new float[12];
11685 Transform inv = t.inverse();
11686 inv.copy_matrix_into_array(m);
11687 image->bindcudaarrayA(true);
11688 p->runcuda(emdata_transform_cuda(m,image->get_xsize(),image->get_ysize(),image->get_zsize()));
11689 image->unbindcudaarryA();
11690 delete [] m;
11691 p->update();
11692 }
11693#endif
11694
11695 if ( p == 0 ) {
11696 float* des_data = transform(image,t);
11697 p = new EMData(des_data,image->get_xsize(),image->get_ysize(),image->get_zsize(),image->get_attr_dict());
11698 }
11699
11700 // all_translation += transform.get_trans();
11701
11702 float scale = t.get_scale();
11703 if (scale != 1.0) {
11704 p->scale_pixel(1.0f/scale);
11705// update_emdata_attributes(p,image->get_attr_dict(),scale);
11706 }
11707
11708 EXITFUNC;
11709 return p;
11710}
EMAN::Dict
Dict is a dictionary to store <string, EMObject> pair.
Definition: emobject.h:385
EMAN::Dict::set_default
type set_default(const string &key, type val)
Default setting behavior This can be achieved using a template - d.woolford Jan 2008 (before there wa...
Definition: emobject.h:569
EMAN::Dict::has_key
bool has_key(const string &key) const
Ask the Dictionary if it as a particular key.
Definition: emobject.h:511
EMAN::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
EMAN::TransformProcessor::transform
float * transform(const EMData *const image, const Transform &t) const
Definition: processor.cpp:10696
EMAN::TransformProcessor::assert_valid_aspect
void assert_valid_aspect(const EMData *const image) const
Definition: processor.cpp:11620
EMAN::Transform
A Transform object is a somewhat specialized object designed specifically for EMAN2/Sparx storage of ...
Definition: transform.h:75
EMAN::Transform::copy_matrix_into_array
void copy_matrix_into_array(float *const) const
Definition: transform.cpp:158
EMAN::Transform::set_rotation
void set_rotation(const Dict &rotation)
Set a rotation using a specific Euler type and the dictionary interface Works for all Euler types.
Definition: transform.cpp:519
EMAN::Transform::get_scale
float get_scale() const
Get the scale that was applied.
Definition: transform.cpp:1145
EMAN::Transform::inverse
Transform inverse() const
Get the inverse of this transformation matrix.
Definition: transform.cpp:1327
EMAN::Transform::set_trans
void set_trans(const float &x, const float &y, const float &z=0)
Set the post translation component.
Definition: transform.cpp:1036
emdata_transform_cuda
float * emdata_transform_cuda(const float *const m, const int nx, const int ny, const int nz)
ENTERFUNC
#define ENTERFUNC
Definition: log.h:48
EXITFUNC
#define EXITFUNC
Definition: log.h:49

References assert_valid_aspect(), EMAN::Transform::copy_matrix_into_array(), emdata_transform_cuda(), ENTERFUNC, EXITFUNC, EMAN::Transform::get_scale(), EMAN::Dict::has_key(), EMAN::Transform::inverse(), EMAN::Processor::params, EMAN::Dict::set_default(), EMAN::Transform::set_rotation(), EMAN::Transform::set_trans(), and transform().

◆ process_inplace()

void TransformProcessor::process_inplace ( EMData *  image)
virtual
Exceptions
ImageDimensionExceptionif the image is not 2D or 3D
InvalidParameterExceptionif the Transform parameter is not specified

Implements EMAN::Processor.

Definition at line 11712 of file processor.cpp.

11712 {
11713 ENTERFUNC;
11714
11715 assert_valid_aspect(image);
11716
11717 Transform t;
11718 if (params.has_key("transform")) t=*(Transform *)params["transform"];
11719 else {
11720 if (params.has_key("alpha")) params["phi"]=params["alpha"];
11721 float az=params.set_default("az",0.0f);
11722 float alt=params.set_default("alt",0.0f);
11723 float phi=params.set_default("phi",0.0f);
11724 float tx=params.set_default("tx",0.0f);
11725 float ty=params.set_default("ty",0.0f);
11726 float tz=params.set_default("tz",0.0f);
11727
11728 t.set_rotation(Dict("type","eman","az",az,"alt",alt,"phi",phi));
11729 t.set_trans(tx,ty,tz);
11730 }
11731
11732 // all_translation += transform.get_trans();
11733 bool use_cpu = true;
11734
11735#ifdef EMAN2_USING_CUDA
11736 if(EMData::usecuda == 1 && image->isrodataongpu()){
11737 //cout << "CUDA xform inplace" << endl;
11738 image->bindcudaarrayA(false);
11739 float * m = new float[12];
11740 Transform inv = t.inverse();
11741 inv.copy_matrix_into_array(m);
11742 image->runcuda(emdata_transform_cuda(m,image->get_xsize(),image->get_ysize(),image->get_zsize()));
11743 image->unbindcudaarryA();
11744 delete [] m;
11745 use_cpu = false;
11746 image->update();
11747 }
11748#endif
11749 if ( use_cpu ) {
11750 float* des_data = transform(image,t);
11751 image->set_data(des_data,image->get_xsize(),image->get_ysize(),image->get_zsize());
11752 image->update();
11753 }
11754 float scale = t.get_scale();
11755 if (scale != 1.0f) {
11756 image->scale_pixel(1.0f/scale);
11757// update_emdata_attributes(image,image->get_attr_dict(),scale);
11758 }
11759
11760 EXITFUNC;
11761}

References assert_valid_aspect(), EMAN::Transform::copy_matrix_into_array(), emdata_transform_cuda(), ENTERFUNC, EXITFUNC, EMAN::Transform::get_scale(), EMAN::Dict::has_key(), EMAN::Transform::inverse(), EMAN::Processor::params, EMAN::Dict::set_default(), EMAN::Transform::set_rotation(), EMAN::Transform::set_trans(), and transform().

◆ transform()

float * TransformProcessor::transform ( const EMData *const  image,
const Transform &  t 
) const

Definition at line 10696 of file processor.cpp.

10696 {
10697
10698 ENTERFUNC;
10699
10700 Transform inv = t.inverse();
10701 int nx = image->get_xsize();
10702 int ny = image->get_ysize();
10703 int nz = image->get_zsize();
10704 int nxy = nx*ny;
10705 int N = ny;
10706
10707 int zerocorners = params.set_default("zerocorners",0);
10708
10709
10710 const float * const src_data = image->get_const_data();
10711 float *des_data = (float *) EMUtil::em_calloc(sizeof(float)*nx,ny*nz);
10712
10713 if ((nz == 1)&&(image -> is_real())) {
10714 Vec2f offset(nx/2,ny/2);
10715 for (int j = 0; j < ny; j++) {
10716 for (int i = 0; i < nx; i++) {
10717 Vec2f coord(i-nx/2,j-ny/2);
10718 Vec2f soln = inv*coord;
10719 soln += offset;
10720
10721 float x2 = soln[0];
10722 float y2 = soln[1];
10723
10724 if (x2 < 0 || x2 >= nx || y2 < 0 || y2 >= ny ) {
10725 des_data[i + j * nx] = 0; // It may be tempting to set this value to the
10726 // mean but in fact this is not a good thing to do. Talk to S.Ludtke about it.
10727 }
10728 else {
10729 int ii = Util::fast_floor(x2);
10730 int jj = Util::fast_floor(y2);
10731 int k0 = ii + jj * nx;
10732 int k1 = k0 + 1;
10733 int k2 = k0 + nx;
10734 int k3 = k0 + nx + 1;
10735
10736 if (ii == nx - 1) {
10737 k1--;
10738 k3--;
10739 }
10740 if (jj == ny - 1) {
10741 k2 -= nx;
10742 k3 -= nx;
10743 }
10744
10745 float t = x2 - ii;
10746 float u = y2 - jj;
10747
10748 des_data[i + j * nx] = Util::bilinear_interpolate(src_data[k0],src_data[k1], src_data[k2], src_data[k3],t,u);
10749 }
10750 }
10751 }
10752 }
10753 if ((nz == 1)&&(image -> is_complex())&&(nx%2==0)&&((2*(nx-ny)-3)*(2*(nx-ny)-3)==1)&&(zerocorners==0) ) {
10754 //printf("Hello 2-d complex TransformProcessor \n");
10755 // make sure there was a realImage.process('xform.phaseorigin.tocorner')
10756 // before transformation to Fourier Space
10757// This rotates a complex image that is a FT of a real space image: it has Friedel symmetries
10758// An arbitrary complex image, F, can be decomposed into two Friedel images
10759// G(k) = (F(k) + F*(-k))/2 , H(k) = (-i(F(k) - F*(k)))/2;
10760// via F(k) = G(k) + i H(k); notice G and H are Friedel symmetric (not necessarily real)
10761// But g,h are real, using f=g+ih.
10762// First make sure that image has proper size;
10763// if 2N is side of real image, then sizes of FFT are (2N+2,2N)
10764// if 2N+1 is side of real image, then sizes of FFT are (2N+2,2N+1)
10765// so we need nx =ny+2, and ny even
10766// or nx = ny+1 and ny odd
10767// so nx even, and 2 *(nx -ny) -3= +- 1; So abs(2*(nx-ny)-3) == 1
10768 float theta = t.get_rotation("eman").get("phi"); theta=theta*pi/180;
10769 Vec3f transNow= t.get_trans();
10770 float xshift= transNow[0]; float yshift= transNow[1];
10771 float tempR; float tempI;float tempW;
10772// int kNy= ny; //size of the real space image
10773// int kNx= nx/2; //
10774 Vec2f offset(nx/2,ny/2);
10775 float Mid =(N+1.0)/2.0; // Check this
10776 for (int kyN = 0; kyN < ny; kyN++) {
10777 int kyNew = kyN;
10778 if (kyN>=nx/2) kyNew=kyN-ny; // Step 0 Unalias
10779 float kxOldPre= - sin(theta)* kyNew;
10780 float kyOldPre= cos(theta)* kyNew;
10781 float phaseY = -2*pi*kyNew*yshift/ny;
10782
10783 for (int kxN = 0; kxN < (nx/2); kxN++) {
10784 int kxNew=kxN;
10785 if (kxN >= nx/2) kxNew=kxN-ny;// Step 0 Unalias
10786 // Step 1, Do rotation and find 4 nn
10787 float kxOld= kxOldPre + cos(theta)* kxNew ;
10788 float kyOld= kyOldPre + sin(theta)* kxNew ;
10789 //
10790 int kxLower= floor(kxOld); int kxUpper= kxLower+1; // All these values
10791 int kyLower= floor(kyOld); int kyUpper= kyLower+1; // may be neg
10792 float dkxLower= (kxUpper-kxOld); float dkxUpper= (kxOld-kxLower);
10793 float dkyLower= (kyUpper-kyOld); float dkyUpper= (kyOld-kyLower);
10794//
10795 int kxL= kxLower; int kyL=kyLower;
10796 float dataLL_R= 0; float dataLL_I=0; int flag=1;
10797 if ((abs(kxL)<Mid) && (abs(kyL)<Mid)) { // Step 2 Make sure to be in First BZ
10798 kxL = (N+kxL)%N; kyL = (N+kyL)%N;
10799 if (kxL> N/2){ kxL=(N-kxL)%N; kyL=(N-kyL)%N ;flag=-1;} // Step 3: if nec, use Friedel paired
10800 dataLL_R= src_data[2*kxL+ kyL*nx]; // image -> get_value_at(2*kxL,kyL);
10801 dataLL_I= flag*src_data[2*kxL+1+ kyL*nx]; // image -> get_value_at(2*kxL+1,kyL);
10802 }
10803
10804 kxL=kxLower; int kyU=kyUpper;
10805 float dataLU_R= 0; float dataLU_I=0; flag=1;
10806 if ((abs(kxL)<Mid) && (abs(kyU)<Mid)){ // Step 2 Make sure to be in First BZ
10807 kxL = (N+kxL)%N; kyU = (N+kyU)%N;
10808 if (kxL> N/2){ kxL=(N-kxL)%N; kyU=(N-kyU)%N;flag=-1;} // Step 3
10809 dataLU_R= src_data[2*kxL+ kyU*nx];// image -> get_value_at(2*kxL,kyU);
10810 dataLU_I= flag*src_data[2*kxL+1+ kyU*nx];// flag*image -> get_value_at(2*kxL+1,kyU);
10811 }
10812
10813 int kxU= kxUpper; kyL=kyLower;
10814 float dataUL_R= 0; float dataUL_I=0; flag=1;
10815 if ((abs(kxU)<Mid) && (abs(kyL)<Mid)) { // Step 2
10816 kxU = (N+kxU)%N; kyL = (N+kyL)%N;
10817 if (kxU> N/2) { kxU=(N-kxU)%N; kyL=(N-kyL)%N;flag=-1;} // Step 3
10818 dataUL_R= src_data[2*kxU + kyL*nx]; // image -> get_value_at(2*kxU,kyL);
10819 dataUL_I= flag*src_data[2*kxU+1 + kyL*nx]; // flag*image -> get_value_at(2*kxU+1,kyL);
10820 }
10821
10822 kxU= kxUpper; kyU=kyUpper;
10823 float dataUU_R= 0; float dataUU_I=0; flag=1;
10824 if ((abs(kxU)<Mid) & (abs(kyU)<Mid)){ // Step 2
10825 kxU = (N+kxU)%N; kyU = (N+kyU)%N;
10826 if (kxU> N/2) { kxU=(N-kxU)%N; kyU=(N-kyU)%N;flag=-1;} // Step 3
10827 dataUU_R= src_data[2*kxU + kyU*nx]; // image -> get_value_at(2*kxU,kyU);
10828 dataUU_I= flag*src_data[2*kxU+1 + kyU*nx]; //flag*image -> get_value_at(2*kxU+1,kyU);
10829 }
10830 // Step 4 Assign Weights
10831 float WLL = dkxLower*dkyLower ;
10832 float WLU = dkxLower*dkyUpper ;// make more intricated weightings here
10833 float WUL = dkxUpper*dkyLower ;// WLL(dkxLower,dkyLower)
10834 float WUU = dkxUpper*dkyUpper ;// etc
10835 tempW = WLL + WLU + WUL + WUU ;
10836
10837 // Step 5 Assign Real, then Imaginar Values
10838 tempR = WLL*dataLL_R + WLU*dataLU_R + WUL* dataUL_R + WUU * dataUU_R ;
10839 tempI = WLL*dataLL_I + WLU*dataLU_I + WUL* dataUL_I + WUU * dataUU_I ;
10840 //
10841 float phase = phaseY -2*pi*kxNew*xshift/ny;
10842 float tempRb=tempR*cos(phase) - tempI*sin(phase);
10843 float tempIb=tempR*sin(phase) + tempI*cos(phase);
10844 //
10845 des_data[2*kxN + nx* kyN] = tempRb/tempW;
10846 des_data[2*kxN+1 + nx* kyN] = tempIb/tempW;
10847 //printf(" kxNew = %d, kyNew = %d,kxOld = %3.2f, kyOld = %3.2f, xl = %d,xU = %d,yl = %d,yu = %d, tempR = %3.2f, tempI=%3.2f, \n",
10848 // kxNew,kyNew, kxOld, kyOld, kxLower,kxUpper,kyLower,kyUpper, tempR, tempI);
10849 }
10850 }
10851 }
10852
10853 if ((nz == 1)&&(image -> is_complex())&&(nx%2==0)&&((2*(nx-ny)-3)*(2*(nx-ny)-3)==1)&&(zerocorners==1) ) {
10854 //printf("Hello 2-d complex TransformProcessor \n");
10855 // make sure there was a realImage.process('xform.phaseorigin.tocorner')
10856 // before transformation to Fourier Space
10857// This rotates a complex image that is a FT of a real space image: it has Friedel symmetries
10858// An arbitrary complex image, F, can be decomposed into two Friedel images
10859// G(k) = (F(k) + F*(-k))/2 , H(k) = (-i(F(k) - F*(k)))/2;
10860// via F(k) = G(k) + i H(k); notice G and H are Friedel symmetric (not necessarily real)
10861// But g,h are real, using f=g+ih.
10862// First make sure that image has proper size;
10863// if 2N is side of real image, then sizes of FFT are (2N+2,2N)
10864// if 2N+1 is side of real image, then sizes of FFT are (2N+2,2N+1)
10865// so we need nx =ny+2, and ny even
10866// or nx = ny+1 and ny odd
10867// so nx even, and 2 *(nx -ny) -3= +- 1; So abs(2*(nx-ny)-3) == 1
10868 float theta = t.get_rotation("2d").get("alpha"); theta=theta*pi/180;
10869 float Ctheta= cos(theta);
10870 float Stheta= sin(theta);
10871 int mirror= t.get_mirror()?(-1.0f):1.0f; // added by steve on 11/7/16
10872 Vec3f transNow= t.get_trans();
10873 float xshift= transNow[0]; float yshift= transNow[1];
10874// printf("%f\t%f\t%f\t%d\n",theta,xshift,yshift,mirror);
10875 float tempR; float tempI; // float tempW;
10876 float Mid =N*N/4; // steve changed to N^2
10877// int kNy= ny; //size of the real space image
10878// int kNx= nx/2; //
10879 Vec2f offset(nx/2,ny/2);
10880 float phaseConstx = -2*pi*xshift/ny ;
10881 float k1= cos(phaseConstx); float k2= sin(phaseConstx);
10882 float k3= 1.0/k1; float k4= k2/k1; // that is 1/cos and tan()
10883 int nxy=nx*ny;
10884
10885 for (int kyN = 0; kyN < ny; kyN++) {
10886 int kyNew = kyN;
10887 if (kyN>=nx/2) kyNew=kyN-ny; // Step 0 Unalias
10888 float kxOld= - Stheta* kyNew - Ctheta;
10889 float kyOld= Ctheta* kyNew - Stheta;
10890 float phase = -2*pi*kyNew*yshift/ny - phaseConstx ;
10891 float Cphase = cos(phase);
10892 float Sphase = sin(phase);
10893 int kx,ky;
10894 int IndexOut;
10895 if (mirror==-1.0) {
10896 if (kyN>0) IndexOut = -2 + nx * (ny-kyN); // Added by steve on 11/7/16
10897 else IndexOut = -2+ nx* kyN;
10898 }
10899 else IndexOut = -2+ nx* kyN;
10900 for (int kxN = 0; kxN < (nx/2); kxN++) {
10901 IndexOut += 2;
10902 // int kxNew=kxN;
10903 // if (kxN >= nx/2) kxNew=kxN-ny;// Step 0 Unalias, but this won't happen
10904 // Step 1, Do rotation and find 4 nn
10905 kxOld += Ctheta ;// was using kxNew instead of kxN
10906 kyOld += Stheta ;
10907 phase += phaseConstx;
10908 Cphase = Cphase*k1 -Sphase*k2; //update using trig addition; this is cos = cos cos -sin sin
10909 Sphase = Sphase*k3+ Cphase*k4; // and sin = sin (1/ cos) + cos * tan;
10910
10911// if ((abs(kxOld)>=Mid) || (abs(kyOld)>=Mid)) { // out of bounds
10912 if (Util::square_sum(kxOld,kyOld)>=Mid) {
10913 des_data[IndexOut] = 0;
10914 des_data[IndexOut+1] = 0;
10915 continue;}
10916
10917 //
10918 int kxLower= floor(kxOld); // All these values
10919 int kyLower= floor(kyOld); // may be neg
10920 //float dkxLower= (kxLower+1-kxOld);
10921 float dkxUpper= (kxOld-kxLower);
10922 //float dkyLower= (kyLower+1-kyOld);
10923 float dkyUpper= (kyOld-kyLower);
10924//
10925 kx= kxLower; ky =kyLower;
10926 int flagLL=1;
10927 if (kx<0) kx += N;
10928 if (ky<0) ky += N;
10929 if (kx> N/2){ kx=(N-kx) ; ky=(N-ky) ;flagLL=-1;} // Step 3: if nec, use Friedel paired
10930 int kLL =2*kx+ky*nx;
10931
10932 kx = kxLower; ky = kyLower+1;
10933 int flagLU =1;
10934 if (kx<0) kx += N;
10935 if (ky<0) ky += N;
10936 if (kx> N/2){ kx=(N-kx) ; ky=(N-ky) ;flagLU=-1;} // Step 3: if nec, use Friedel paired
10937 int kLU =2*kx+ky*nx;
10938
10939 kx = kxLower+1; ky=kyLower;
10940 int flagUL=1;
10941 if (kx<0) kx += N;
10942 if (ky<0) ky += N;
10943 if (kx> N/2){ kx=(N-kx) ; ky=(N-ky) ;flagUL=-1;} // Step 3: if nec, use Friedel paired
10944 int kUL =2*kx+ky*nx;
10945
10946 kx = kxLower+1; ky = kyLower+1;
10947 int flagUU =1;
10948 if (kx<0) kx += N;
10949 if (ky<0) ky += N;
10950 if (kx> N/2){ kx=(N-kx) ; ky=(N-ky) ;flagUU=-1;} // Step 3: if nec, use Friedel paired
10951 int kUU =2*kx+ky*nx;
10952
10953 // Step 4 Assign Weights
10954// float WLL = dkxLower*dkyLower ;
10955// float WLU = dkxLower*dkyUpper ;// make more intricated weightings here
10956// float WUL = dkxUpper*dkyLower ;// WLL(dkxLower,dkyLower)
10957// float WUU = dkxUpper*dkyUpper ;// etc
10958// tempW = WLL + WLU + WUL + WUU ;
10959
10960 // Added by steve. MAJOR limit checking issues here! Hope this doesn't cause any other problems
10961 if (kLL<0||kUL<0||kLU<0||kUU<0||kLL>=nxy||kUL>=nxy||kLU>=nxy||kUU>=nxy) continue;
10962 //printf("%d %d %d %d\n",kLL,kUL,kLU,kUU);
10963
10964 // Step 5 Assign Real, then Imaginary Values
10965 tempR = Util::bilinear_interpolate(src_data[kLL],src_data[kUL], src_data[kLU], src_data[kUU],dkxUpper,dkyUpper);//
10966 //WLL*dataLL_R + WLU*dataLU_R + WUL* dataUL_R + WUU * dataUU_R ;
10967 tempI = Util::bilinear_interpolate(flagLL* src_data[kLL+1],flagUL*src_data[kUL+1],
10968 flagLU*src_data[kLU+1], flagUU*src_data[kUU+1],dkxUpper,dkyUpper);
10969 //WLL*dataLL_I + WLU*dataLU_I + WUL* dataUL_I + WUU * dataUU_I ;
10970 //
10971
10972 float tempRb=tempR*Cphase - tempI*Sphase;
10973 float tempIb=tempR*Sphase + tempI*Cphase;
10974 //
10975 des_data[IndexOut] = tempRb;
10976// des_data[IndexOut+1] = tempIb;
10977 des_data[IndexOut+1] = tempIb*mirror;
10978 //printf(" kxNew = %d, kyNew = %d,kxOld = %3.2f, kyOld = %3.2f, xl = %d,xU = %d,yl = %d,yu = %d, tempR = %3.2f, tempI=%3.2f, \n",
10979 // kxNew,kyNew, kxOld, kyOld, kxLower,kxUpper,kyLower,kyUpper, tempR, tempI);
10980 }
10981 }
10982 }
10983
10984 if ((nz > 1)&&(image -> is_complex())&&(zerocorners==3)) {
10985 //printf("Hello 3-d complex TransformProcessor \n");
10986 float phi = t.get_rotation("eman").get("phi"); phi=pi*phi/180;
10987 float alt = t.get_rotation("eman").get("alt"); alt=pi*alt/180;
10988 float az = t.get_rotation("eman").get("az"); az=pi*az /180;
10989 Vec3f transNow= t.get_trans();
10990 float xshift= transNow[0]; float yshift= transNow[1];float zshift= transNow[2];
10991
10992 float phaseConstx = -2*pi*xshift/ny ;
10993 float k1= cos(phaseConstx); float k2= sin(phaseConstx);
10994 float k3= 1.0/k1; float k4= k2/k1; // that is 1/cos and tan()
10995
10996 float MatXX = (cos(az)*cos(phi) - sin(az)*cos(alt)*sin(phi) );
10997 float MatXY = (- cos(az)*sin(phi) - sin(az)*cos(alt)*cos(phi) ) ;
10998 float MatXZ = sin(az)*sin(alt) ;
10999 float MatYX = (sin(az)*cos(phi) + cos(az)*cos(alt)*sin(phi) );
11000 float MatYY = (- sin(az)*sin(phi) + cos(az)*cos(alt)*cos(phi) ) ;
11001 float MatYZ = - cos(az)*sin(alt) ;
11002 float MatZX = sin(alt)*sin(phi);
11003 float MatZY = sin(alt)*cos(phi);
11004 float MatZZ = cos(alt) ;
11005 float tempR; float tempI; float tempW;
11006 float Mid =(N+1.0)/2.0; // Check this
11007 int nxny = nx*ny;
11008
11009 for (int kzN = 0; kzN < ny; kzN++) {
11010 int kzNew=kzN;
11011 if (kzN >= nx/2) kzNew=kzN-N; // Step 0 Unalias new coords;
11012 for (int kyN = 0; kyN < ny; kyN++) {// moves them to lesser mag
11013 int kyNew=kyN;
11014 if (kyN>=nx/2) kyNew=kyN-ny; // Step 0 Unalias
11015 float kxPre = MatXY * kyNew + MatXZ *kzNew;
11016 float kyPre = MatYY * kyNew + MatYZ*kzNew;
11017 float kzPre = MatZY * kyNew + MatZZ*kzNew;
11018 float phase = -2*pi*kzNew*zshift/ny-2*pi*kyNew*yshift/ny - phaseConstx ;
11019 float Cphase = cos(phase);
11020 float Sphase = sin(phase);
11021
11022 float OutBounds2= (2*Mid*Mid- (kyNew*kyNew+kzNew*kzNew)) ;
11023 int kxNewMax= nx/2;
11024 if (OutBounds2< 0) kxNewMax=0;
11025 else if (OutBounds2<(nx*nx/4)) kxNewMax=sqrt(OutBounds2);
11026 for (int kxN = kxNewMax; kxN < nx/2 ; kxN++ ) {
11027 des_data[2*kxN + nx* kyN +nxny*kzN] = 0;
11028 des_data[2*kxN + nx* kyN +nxny*kzN+1] = 0;
11029 }
11030
11031
11032 for (int kxN = 0; kxN < kxNewMax; kxN++ ) {
11033 Cphase = Cphase*k1 -Sphase*k2; //update using trig addition; this is cos = cos cos -sin sin
11034 Sphase = Sphase*k3+ Cphase*k4; // and sin = sin (1/ cos) + cos * tan;
11035 // Step 1: Do inverse Rotation to get former values, and alias Step1
11036 float kxOld= MatXX * kxN + kxPre;
11037 float kyOld= MatYX * kxN + kyPre;
11038 float kzOld= MatZX * kxN + kzPre;
11039 //
11040 if ((abs(kxOld)>=Mid) || (abs(kyOld)>=Mid) || (abs(kzOld)>=Mid) ) { // out of bounds
11041 des_data[2*kxN + nx* kyN +nxny*kzN] = 0;
11042 des_data[2*kxN + nx* kyN +nxny*kzN+1] = 0;
11043 continue;}
11044 //
11045 int kxLower= floor(kxOld); int kxUpper= kxLower+1;
11046 int kyLower= floor(kyOld); int kyUpper= kyLower+1;
11047 int kzLower= floor(kzOld); int kzUpper= kzLower+1;
11048 //
11049 float dkxLower= (kxUpper-kxOld); float dkxUpper= (kxOld-kxLower);
11050 float dkyLower= (kyUpper-kyOld); float dkyUpper= (kyOld-kyLower);
11051 float dkzLower= (kzUpper-kzOld); float dkzUpper= (kzOld-kzLower);
11052 // LLL 1
11053 int kxL= kxLower; int kyL=kyLower; int kzL=kzLower;
11054 float dataLLL_R= 0; float dataLLL_I=0; int flag=1;
11055 if ( (abs(kxL)<Mid) && (abs(kyL)<Mid) && (abs(kzL)<Mid) ) { // Step 2 Make sure to be in First BZ
11056 kxL = (N+kxL)%N; kyL = (N+kyL)%N; kzL = (N+kzL)%N;
11057 if (kxL> N/2){kxL=(N-kxL)%N; kyL=(N-kyL)%N ; kzL=(N-kzL)%N ;flag=-1;} // Step 3: use Friedel paired
11058 dataLLL_R= src_data[ 2*kxL + nx*kyL+ nxny*kzL ]; //get_value_at(2*kxL,kyL,kzL);
11059 dataLLL_I=flag*src_data[ 2*kxL+1 + nx*kyL+ nxny*kzL ];//get_value_at(2*kxL+1,kyL,kzL);
11060 }
11061 // LLU 2
11062 kxL= kxLower; kyL=kyLower; int kzU=kzUpper;
11063 float dataLLU_R= 0; float dataLLU_I=0; flag=1;
11064 if ( (abs(kxL)<Mid) && (abs(kyL)<Mid) && (abs(kzU)<Mid) ) {// Step 2 Make sure to be in First BZ
11065 kxL = (N+kxL)%N; kyL = (N+kyL)%N; kzU = (N+kzU)%N;
11066 if (kxL> N/2){kxL=(N-kxL)%N; kyL=(N-kyL)%N ; kzU=(N-kzU)%N ;flag=-1;} // Step 3: use Friedel paired
11067 dataLLU_R= src_data[ 2*kxL + nx*kyL+ nxny*kzU ]; // image -> get_value_at(2*kxL ,kyL,kzU);
11068 dataLLU_I= flag*src_data[ 2*kxL+1 + nx*kyL+ nxny*kzU ]; // image -> get_value_at(2*kxL+1,kyL,kzU);
11069 }
11070 // LUL 3
11071 kxL= kxLower; int kyU=kyUpper; kzL=kzLower;
11072 float dataLUL_R= 0; float dataLUL_I=0; flag=1;
11073 if ( (abs(kxL)<Mid) && (abs(kyU)<Mid)&& (abs(kzL)<Mid) ) {// Step 2 Make sure to be in First BZ
11074 kxL = (N+kxL)%N; kyU = (N+kyU)%N; kzL = (N+kzL)%N;
11075 if (kxL> N/2){ kxL=(N-kxL)%N; kyU=(N-kyU)%N; kzL=(N-kzL)%N ;flag=-1;}// Step 3
11076 dataLUL_R= src_data[ 2*kxL + nx*kyU+ nxny*kzL ]; // image -> get_value_at(2*kxL ,kyU,kzL);
11077 dataLUL_I=flag*src_data[ 2*kxL+1 + nx*kyU+ nxny*kzL ]; // image -> get_value_at(2*kxL+1,kyU,kzL);
11078 }
11079 // LUU 4
11080 kxL= kxLower; kyU=kyUpper; kzL=kzUpper;
11081 float dataLUU_R= 0; float dataLUU_I=0; flag=1;
11082 if ( (abs(kxL)<Mid) && (abs(kyU)<Mid)&& (abs(kzU)<Mid)) {// Step 2 Make sure to be in First BZ
11083 kxL = (N+kxL)%N; kyU = (N+kyU)%N; kzU = (N+kzU)%N;
11084 if (kxL> N/2){kxL=(N-kxL)%N; kyU=(N-kyU)%N; kzL=(N-kzL)%N ;flag=-1;} // Step 3
11085 dataLUU_R= src_data[ 2*kxL + nx*kyU+ nxny*kzU ]; // image -> get_value_at(2*kxL ,kyU,kzU);
11086 dataLUU_I=flag*src_data[ 2*kxL+1 + nx*kyU+ nxny*kzU ]; // image -> get_value_at(2*kxL+1,kyU,kzU);
11087 }
11088 // ULL 5
11089 int kxU= kxUpper; kyL=kyLower; kzL=kzLower;
11090 float dataULL_R= 0; float dataULL_I=0; flag=1;
11091 if ( (abs(kxU)<Mid) && (abs(kyL)<Mid) && (abs(kzL)<Mid) ) {// Step 2
11092 kxU = (N+kxU)%N; kyL = (N+kyL)%N; kzL = (N+kzL)%N;
11093 if (kxU> N/2){kxU=(N-kxU)%N; kyL=(N-kyL)%N; kzL=(N-kzL)%N ;flag=-1;} // Step 3
11094 dataULL_R= src_data[ 2*kxU + nx*kyL+ nxny*kzL ]; // image -> get_value_at(2*kxU ,kyL,kzL);
11095 dataULL_I=flag*src_data[ 2*kxU+1 + nx*kyL+ nxny*kzL ]; // image -> get_value_at(2*kxU+1,kyL,kzL);
11096 }
11097 // ULU 6
11098 kxU= kxUpper; kyL=kyLower; kzU=kzUpper;
11099 float dataULU_R= 0; float dataULU_I=0; flag=1;
11100 if ( (abs(kxU)<Mid) && (abs(kyL)<Mid)&& (abs(kzU)<Mid) ) {// Step 2
11101 kxU = (N+kxU)%N; kyL = (N+kyL)%N; kzU = (N+kzU)%N;
11102 if (kxU> N/2){kxU=(N-kxU)%N; kyL=(N-kyL)%N; kzU=(N-kzU)%N ;flag=-1;} // Step 3
11103 dataULU_R= src_data[ 2*kxU + nx*kyL+ nxny*kzU ]; // image -> get_value_at(2*kxU ,kyL,kzU);
11104 dataULU_I=flag*src_data[ 2*kxU+1 + nx*kyL+ nxny*kzU ]; // image -> get_value_at(2*kxU+1,kyL,kzU);
11105 }
11106 // UUL 7
11107 kxU= kxUpper; kyU=kyUpper; kzL=kzLower;
11108 float dataUUL_R= 0; float dataUUL_I=0; flag=1;
11109 if ( (abs(kxU)<Mid) && (abs(kyU)<Mid) && (abs(kzL)<Mid) ) {// Step 2
11110 kxU = (N+kxU)%N; kyU = (N+kyU)%N; kzL = (N+kzL)%N;
11111 if (kxU> N/2){kxU=(N-kxU)%N; kyU=(N-kyU)%N; kzL=(N-kzL)%N ;flag=-1;} // Step 3
11112 dataUUL_R= src_data[ 2*kxU + nx*kyU+ nxny*kzL ]; // image -> get_value_at(2*kxU ,kyU,kzL);
11113 dataUUL_I=flag*src_data[ 2*kxU+1 + nx*kyU+ nxny*kzL ]; // image -> get_value_at(2*kxU+1,kyU,kzL);
11114 }
11115 // UUU 8
11116 kxU= kxUpper; kyU=kyUpper; kzU=kzUpper;
11117 float dataUUU_R= 0; float dataUUU_I=0; flag=1;
11118 if ( (abs(kxU)<Mid) && (abs(kyU)<Mid) && (abs(kzU)<Mid) ) { // Step 2
11119 kxU = (N+kxU)%N; kyU = (N+kyU)%N; kzU = (N+kzU)%N;
11120 if (kxU> N/2) {kxU=(N-kxU)%N; kyU=(N-kyU)%N; kzU=(N-kzU)%N ;flag=-1;} // Step 3
11121 dataUUU_R= src_data[ 2*kxU + nx*kyU+ nxny*kzU ]; // image -> get_value_at(2*kxU ,kyU,kzU);
11122 dataUUU_I=flag*src_data[ 2*kxU+1 + nx*kyU+ nxny*kzU ]; // image -> get_value_at(2*kxU+1,kyU,kzU);
11123 }
11124 // Step 4 Assign Weights
11125 float WLLL = dkxLower*dkyLower*dkzLower ;
11126 // WLLL = sqrt(pow(dkxUpper,2)+ pow(dkyUpper,2) + pow(dkzUpper,2)) ;
11127 // WLLL = sin(pi* WLLL+.00000001)/(pi* WLLL+.00000001);
11128 float WLLU = dkxLower*dkyLower*dkzUpper ;
11129 float WLUL = dkxLower*dkyUpper*dkzLower ;
11130 float WLUU = dkxLower*dkyUpper*dkzUpper ;
11131 float WULL = dkxUpper*dkyLower*dkzLower ;
11132 float WULU = dkxUpper*dkyLower*dkzUpper ;
11133 float WUUL = dkxUpper*dkyUpper*dkzLower;
11134 float WUUU = dkxUpper*dkyUpper*dkzUpper;
11135 tempW = WLLL + WLLU + WLUL + WLUU + WULL + WULU + WUUL + WUUU ;
11136 // Step 5 Assign Real, then Imaginary Values
11137 tempR = WLLL*dataLLL_R + WLLU*dataLLU_R + WLUL*dataLUL_R + WLUU*dataLUU_R ;
11138 tempR += WULL*dataULL_R + WULU*dataULU_R + WUUL*dataUUL_R + WUUU*dataUUU_R ;
11139 //
11140 tempI = WLLL*dataLLL_I + WLLU*dataLLU_I + WLUL*dataLUL_I + WLUU*dataLUU_I ;
11141 tempI += WULL*dataULL_I + WULU*dataULU_I + WUUL*dataUUL_I + WUUU*dataUUU_I ;
11142 //
11143// float phase = -2*pi*(kxNew*xshift+kyNew*yshift+kzNew*zshift)/ny;
11144 float tempRb=tempR*Cphase - tempI*Sphase;
11145 float tempIb=tempR*Sphase + tempI*Cphase;
11146 //
11147 des_data[2*kxN + nx* kyN +nxny*kzN] = tempRb/tempW;
11148 des_data[2*kxN+1 + nx* kyN +nxny*kzN] = tempIb/tempW;
11149 }}} // end z, y, x loops through new coordinates
11150 } // end rotations in Fourier Space 3D
11151 // Steve trying for more optimization
11152 if ((nz > 1)&&(image -> is_complex())&&(zerocorners==2)) {
11153 //printf("Hello 3-d complex TransformProcessor \n");
11154 float phi = t.get_rotation("eman").get("phi"); phi=pi*phi/180;
11155 float alt = t.get_rotation("eman").get("alt"); alt=pi*alt/180;
11156 float az = t.get_rotation("eman").get("az"); az=pi*az /180;
11157 Vec3f transNow= t.get_trans();
11158 float xshift= transNow[0]; float yshift= transNow[1];float zshift= transNow[2];
11159
11160 float phaseConstx = -2*pi*xshift/ny ;
11161 float k1= cos(phaseConstx); float k2= sin(phaseConstx);
11162 float k3= 1.0/k1; float k4= k2/k1; // that is 1/cos and tan()
11163
11164 float MatXX = (cos(az)*cos(phi) - sin(az)*cos(alt)*sin(phi) );
11165 float MatXY = (- cos(az)*sin(phi) - sin(az)*cos(alt)*cos(phi) ) ;
11166 float MatXZ = sin(az)*sin(alt) ;
11167 float MatYX = (sin(az)*cos(phi) + cos(az)*cos(alt)*sin(phi) );
11168 float MatYY = (- sin(az)*sin(phi) + cos(az)*cos(alt)*cos(phi) ) ;
11169 float MatYZ = - cos(az)*sin(alt) ;
11170 float MatZX = sin(alt)*sin(phi);
11171 float MatZY = sin(alt)*cos(phi);
11172 float MatZZ = cos(alt) ;
11173 //float tempR; float tempI; float tempW; # commented because variables were unused and produced warnings during compilation.
11174 float Mid =(N+1.0)/2.0; // Check this
11175 int lim=(N/2)*(N/2); // this is NOT N*N/4 !
11176 int nxny = nx*ny;
11177
11178 for (int kzN = 0; kzN < ny; kzN++) {
11179 int kzNew=kzN;
11180 if (kzN >= nx/2) kzNew=kzN-N; // Step 0 Unalias new coords;
11181 for (int kyN = 0; kyN < ny; kyN++) {// moves them to lesser mag
11182 int kyNew=kyN;
11183 if (kyN>=nx/2) kyNew=kyN-ny; // Step 0 Unalias
11184
11185 // Establish limits for this row and skip if necessary
11186 int kyz2=kyNew*kyNew+kzNew*kzNew;
11187 if (kyz2>lim) continue; // Whole y,z row is 'missing'
11188 int kxNewMax=(int)floor(sqrt((float)(lim-kyz2)));
11189 float kxPre = MatXY * kyNew + MatXZ *kzNew;
11190 float kyPre = MatYY * kyNew + MatYZ*kzNew;
11191 float kzPre = MatZY * kyNew + MatZZ*kzNew;
11192 float phase = -2*pi*kzNew*zshift/ny-2*pi*kyNew*yshift/ny - phaseConstx ;
11193 float Cphase = cos(phase);
11194 float Sphase = sin(phase);
11195
11196// for (int kxN = kxNewMax; kxN < nx/2 ; kxN++ ) {
11197// des_data[2*kxN + nx* kyN +nxny*kzN] = 0;
11198// des_data[2*kxN + nx* kyN +nxny*kzN+1] = 0;
11199// }
11200
11201
11202 for (int kxN = 0; kxN < kxNewMax; kxN++ ) {
11203 Cphase = Cphase*k1 -Sphase*k2; //update using trig addition; this is cos = cos cos -sin sin
11204 Sphase = Sphase*k3+ Cphase*k4; // and sin = sin (1/ cos) + cos * tan;
11205 // Step 1: Do inverse Rotation to get former values, and alias Step1
11206 float kxOld= MatXX * kxN + kxPre;
11207 float kyOld= MatYX * kxN + kyPre;
11208 float kzOld= MatZX * kxN + kzPre;
11209 //
11210 // This is impossible
11211// if ((abs(kxOld)>=Mid) || (abs(kyOld)>=Mid) || (abs(kzOld)>=Mid) ) { // out of bounds
11212// des_data[2*kxN + nx* kyN +nxny*kzN] = 0;
11213// des_data[2*kxN + nx* kyN +nxny*kzN+1] = 0;
11214// continue;}
11215 //
11216 int kx0= floor(kxOld);
11217 int ky0= floor(kyOld);
11218 int kz0= floor(kzOld);
11219 //
11220 float dkx0= (kxOld-kx0);
11221 float dky0= (kyOld-ky0);
11222 float dkz0= (kzOld-kz0);
11223
11224 std::complex<float> c0 = image->get_complex_at(kx0 ,ky0 ,kz0 );
11225 std::complex<float> c1 = image->get_complex_at(kx0+1,ky0 ,kz0 );
11226 std::complex<float> c2 = image->get_complex_at(kx0 ,ky0+1,kz0 );
11227 std::complex<float> c3 = image->get_complex_at(kx0+1,ky0+1,kz0 );
11228 std::complex<float> c4 = image->get_complex_at(kx0 ,ky0 ,kz0+1);
11229 std::complex<float> c5 = image->get_complex_at(kx0+1,ky0 ,kz0+1);
11230 std::complex<float> c6 = image->get_complex_at(kx0 ,ky0+1,kz0+1);
11231 std::complex<float> c7 = image->get_complex_at(kx0+1,ky0+1,kz0+1);
11232
11233 std::complex<float> nwv = Util::trilinear_interpolate_complex(c0,c1,c2,c3,c4,c5,c6,c7,dkx0,dky0,dkz0);
11234
11235 des_data[2*kxN + nx* kyN +nxny*kzN] = nwv.real()*Cphase - nwv.imag()*Sphase;
11236 des_data[2*kxN+1 + nx* kyN +nxny*kzN] = nwv.real()*Sphase + nwv.imag()*Cphase;
11237
11238// // LLL 1
11239// int kxL= kxLower; int kyL=kyLower; int kzL=kzLower;
11240// float dataLLL_R= 0; float dataLLL_I=0; int flag=1;
11241// if ( (abs(kxL)<Mid) && (abs(kyL)<Mid) && (abs(kzL)<Mid) ) { // Step 2 Make sure to be in First BZ
11242// kxL = (N+kxL)%N; kyL = (N+kyL)%N; kzL = (N+kzL)%N;
11243// if (kxL> N/2){kxL=(N-kxL)%N; kyL=(N-kyL)%N ; kzL=(N-kzL)%N ;flag=-1;} // Step 3: use Friedel paired
11244// dataLLL_R= src_data[ 2*kxL + nx*kyL+ nxny*kzL ]; //get_value_at(2*kxL,kyL,kzL);
11245// dataLLL_I=flag*src_data[ 2*kxL+1 + nx*kyL+ nxny*kzL ];//get_value_at(2*kxL+1,kyL,kzL);
11246// }
11247// // LLU 2
11248// kxL= kxLower; kyL=kyLower; int kzU=kzUpper;
11249// float dataLLU_R= 0; float dataLLU_I=0; flag=1;
11250// if ( (abs(kxL)<Mid) && (abs(kyL)<Mid) && (abs(kzU)<Mid) ) {// Step 2 Make sure to be in First BZ
11251// kxL = (N+kxL)%N; kyL = (N+kyL)%N; kzU = (N+kzU)%N;
11252// if (kxL> N/2){kxL=(N-kxL)%N; kyL=(N-kyL)%N ; kzU=(N-kzU)%N ;flag=-1;} // Step 3: use Friedel paired
11253// dataLLU_R= src_data[ 2*kxL + nx*kyL+ nxny*kzU ]; // image -> get_value_at(2*kxL ,kyL,kzU);
11254// dataLLU_I= flag*src_data[ 2*kxL+1 + nx*kyL+ nxny*kzU ]; // image -> get_value_at(2*kxL+1,kyL,kzU);
11255// }
11256// // LUL 3
11257// kxL= kxLower; int kyU=kyUpper; kzL=kzLower;
11258// float dataLUL_R= 0; float dataLUL_I=0; flag=1;
11259// if ( (abs(kxL)<Mid) && (abs(kyU)<Mid)&& (abs(kzL)<Mid) ) {// Step 2 Make sure to be in First BZ
11260// kxL = (N+kxL)%N; kyU = (N+kyU)%N; kzL = (N+kzL)%N;
11261// if (kxL> N/2){ kxL=(N-kxL)%N; kyU=(N-kyU)%N; kzL=(N-kzL)%N ;flag=-1;}// Step 3
11262// dataLUL_R= src_data[ 2*kxL + nx*kyU+ nxny*kzL ]; // image -> get_value_at(2*kxL ,kyU,kzL);
11263// dataLUL_I=flag*src_data[ 2*kxL+1 + nx*kyU+ nxny*kzL ]; // image -> get_value_at(2*kxL+1,kyU,kzL);
11264// }
11265// // LUU 4
11266// kxL= kxLower; kyU=kyUpper; kzL=kzUpper;
11267// float dataLUU_R= 0; float dataLUU_I=0; flag=1;
11268// if ( (abs(kxL)<Mid) && (abs(kyU)<Mid)&& (abs(kzU)<Mid)) {// Step 2 Make sure to be in First BZ
11269// kxL = (N+kxL)%N; kyU = (N+kyU)%N; kzU = (N+kzU)%N;
11270// if (kxL> N/2){kxL=(N-kxL)%N; kyU=(N-kyU)%N; kzL=(N-kzL)%N ;flag=-1;} // Step 3
11271// dataLUU_R= src_data[ 2*kxL + nx*kyU+ nxny*kzU ]; // image -> get_value_at(2*kxL ,kyU,kzU);
11272// dataLUU_I=flag*src_data[ 2*kxL+1 + nx*kyU+ nxny*kzU ]; // image -> get_value_at(2*kxL+1,kyU,kzU);
11273// }
11274// // ULL 5
11275// int kxU= kxUpper; kyL=kyLower; kzL=kzLower;
11276// float dataULL_R= 0; float dataULL_I=0; flag=1;
11277// if ( (abs(kxU)<Mid) && (abs(kyL)<Mid) && (abs(kzL)<Mid) ) {// Step 2
11278// kxU = (N+kxU)%N; kyL = (N+kyL)%N; kzL = (N+kzL)%N;
11279// if (kxU> N/2){kxU=(N-kxU)%N; kyL=(N-kyL)%N; kzL=(N-kzL)%N ;flag=-1;} // Step 3
11280// dataULL_R= src_data[ 2*kxU + nx*kyL+ nxny*kzL ]; // image -> get_value_at(2*kxU ,kyL,kzL);
11281// dataULL_I=flag*src_data[ 2*kxU+1 + nx*kyL+ nxny*kzL ]; // image -> get_value_at(2*kxU+1,kyL,kzL);
11282// }
11283// // ULU 6
11284// kxU= kxUpper; kyL=kyLower; kzU=kzUpper;
11285// float dataULU_R= 0; float dataULU_I=0; flag=1;
11286// if ( (abs(kxU)<Mid) && (abs(kyL)<Mid)&& (abs(kzU)<Mid) ) {// Step 2
11287// kxU = (N+kxU)%N; kyL = (N+kyL)%N; kzU = (N+kzU)%N;
11288// if (kxU> N/2){kxU=(N-kxU)%N; kyL=(N-kyL)%N; kzU=(N-kzU)%N ;flag=-1;} // Step 3
11289// dataULU_R= src_data[ 2*kxU + nx*kyL+ nxny*kzU ]; // image -> get_value_at(2*kxU ,kyL,kzU);
11290// dataULU_I=flag*src_data[ 2*kxU+1 + nx*kyL+ nxny*kzU ]; // image -> get_value_at(2*kxU+1,kyL,kzU);
11291// }
11292// // UUL 7
11293// kxU= kxUpper; kyU=kyUpper; kzL=kzLower;
11294// float dataUUL_R= 0; float dataUUL_I=0; flag=1;
11295// if ( (abs(kxU)<Mid) && (abs(kyU)<Mid) && (abs(kzL)<Mid) ) {// Step 2
11296// kxU = (N+kxU)%N; kyU = (N+kyU)%N; kzL = (N+kzL)%N;
11297// if (kxU> N/2){kxU=(N-kxU)%N; kyU=(N-kyU)%N; kzL=(N-kzL)%N ;flag=-1;} // Step 3
11298// dataUUL_R= src_data[ 2*kxU + nx*kyU+ nxny*kzL ]; // image -> get_value_at(2*kxU ,kyU,kzL);
11299// dataUUL_I=flag*src_data[ 2*kxU+1 + nx*kyU+ nxny*kzL ]; // image -> get_value_at(2*kxU+1,kyU,kzL);
11300// }
11301// // UUU 8
11302// kxU= kxUpper; kyU=kyUpper; kzU=kzUpper;
11303// float dataUUU_R= 0; float dataUUU_I=0; flag=1;
11304// if ( (abs(kxU)<Mid) && (abs(kyU)<Mid) && (abs(kzU)<Mid) ) { // Step 2
11305// kxU = (N+kxU)%N; kyU = (N+kyU)%N; kzU = (N+kzU)%N;
11306// if (kxU> N/2) {kxU=(N-kxU)%N; kyU=(N-kyU)%N; kzU=(N-kzU)%N ;flag=-1;} // Step 3
11307// dataUUU_R= src_data[ 2*kxU + nx*kyU+ nxny*kzU ]; // image -> get_value_at(2*kxU ,kyU,kzU);
11308// dataUUU_I=flag*src_data[ 2*kxU+1 + nx*kyU+ nxny*kzU ]; // image -> get_value_at(2*kxU+1,kyU,kzU);
11309// }
11310// // Step 4 Assign Weights
11311// float WLLL = dkxLower*dkyLower*dkzLower ;
11312// // WLLL = sqrt(pow(dkxUpper,2)+ pow(dkyUpper,2) + pow(dkzUpper,2)) ;
11313// // WLLL = sin(pi* WLLL+.00000001)/(pi* WLLL+.00000001);
11314// float WLLU = dkxLower*dkyLower*dkzUpper ;
11315// float WLUL = dkxLower*dkyUpper*dkzLower ;
11316// float WLUU = dkxLower*dkyUpper*dkzUpper ;
11317// float WULL = dkxUpper*dkyLower*dkzLower ;
11318// float WULU = dkxUpper*dkyLower*dkzUpper ;
11319// float WUUL = dkxUpper*dkyUpper*dkzLower;
11320// float WUUU = dkxUpper*dkyUpper*dkzUpper;
11321// tempW = WLLL + WLLU + WLUL + WLUU + WULL + WULU + WUUL + WUUU ;
11322// // Step 5 Assign Real, then Imaginary Values
11323// tempR = WLLL*dataLLL_R + WLLU*dataLLU_R + WLUL*dataLUL_R + WLUU*dataLUU_R ;
11324// tempR += WULL*dataULL_R + WULU*dataULU_R + WUUL*dataUUL_R + WUUU*dataUUU_R ;
11325// //
11326// tempI = WLLL*dataLLL_I + WLLU*dataLLU_I + WLUL*dataLUL_I + WLUU*dataLUU_I ;
11327// tempI += WULL*dataULL_I + WULU*dataULU_I + WUUL*dataUUL_I + WUUU*dataUUU_I ;
11328// //
11329// // float phase = -2*pi*(kxNew*xshift+kyNew*yshift+kzNew*zshift)/ny;
11330// float tempRb=tempR*Cphase - tempI*Sphase;
11331// float tempIb=tempR*Sphase + tempI*Cphase;
11332// //
11333// des_data[2*kxN + nx* kyN +nxny*kzN] = tempRb/tempW;
11334// des_data[2*kxN+1 + nx* kyN +nxny*kzN] = tempIb/tempW;
11335 }}} // end z, y, x loops through new coordinates
11336 } // end rotations in Fourier Space 3D
11337 if ((nz > 1)&&(image -> is_complex())&&(zerocorners<=1)) {
11338 //printf("Hello 3-d complex TransformProcessor \n");
11339 float phi = t.get_rotation("eman").get("phi"); phi=pi*phi/180;
11340 float alt = t.get_rotation("eman").get("alt"); alt=pi*alt/180;
11341 float az = t.get_rotation("eman").get("az"); az=pi*az /180;
11342 Vec3f transNow= t.get_trans();
11343 float xshift= transNow[0]; float yshift= transNow[1];float zshift= transNow[2];
11344
11345 float MatXX = (cos(az)*cos(phi) - sin(az)*cos(alt)*sin(phi) );
11346 float MatXY = (- cos(az)*sin(phi) - sin(az)*cos(alt)*cos(phi) ) ;
11347 float MatXZ = sin(az)*sin(alt) ;
11348 float MatYX = (sin(az)*cos(phi) + cos(az)*cos(alt)*sin(phi) );
11349 float MatYY = (- sin(az)*sin(phi) + cos(az)*cos(alt)*cos(phi) ) ;
11350 float MatYZ = - cos(az)*sin(alt) ;
11351 float MatZX = sin(alt)*sin(phi);
11352 float MatZY = sin(alt)*cos(phi);
11353 float MatZZ = cos(alt) ;
11354 float tempR=0; float tempI=0;
11355 float Mid =(N+1.0)/2.0; // Check this
11356 float phaseConstx = -2*pi*xshift/ny ;
11357 float k1= cos(phaseConstx); float k2= sin(phaseConstx);
11358 float k3= 1.0/k1; float k4= k2/k1; // that is 1/cos and tan()
11359// float dataRLLL, dataRLLU, dataRLUL, dataRLUU;
11360// float dataRULL, dataRULU, dataRUUL, dataRUUU;
11361// float dataILLL, dataILLU, dataILUL, dataILUU;
11362// float dataIULL, dataIULU, dataIUUL, dataIUUU;
11363 int nxny = nx*ny;
11364
11365 for (int kzN = 0; kzN < ny; kzN++) {
11366 int kzNew=kzN;
11367 if (kzN >= nx/2) kzNew=kzN-N; // Step 0 Unalias new coords;
11368 for (int kyN = 0; kyN < ny; kyN++) {// moves them to lesser mag
11369 int kyNew=kyN;
11370 if (kyN>=nx/2) kyNew=kyN-ny; // Step 0 Unalias
11371 // Step 1: Do inverse Rotation to get former values, and alias Step1
11372 float kxOld = MatXY * kyNew + MatXZ *kzNew - MatXX;
11373 float kyOld = MatYY * kyNew + MatYZ*kzNew - MatYX;
11374 float kzOld = MatZY * kyNew + MatZZ*kzNew - MatZX;
11375 float phase = -2*pi*kzNew*zshift/ny-2*pi*kyNew*yshift/ny - phaseConstx ;
11376 float Cphase = cos(phase);
11377 float Sphase = sin(phase);
11378 int kx, ky,kz, II;
11379
11380 int IndexOut= -2+ nx* kyN +nxny*kzN;
11381 float OutBounds2 = (Mid*Mid- (kyOld*kyOld+kzOld*kzOld)) ;
11382
11383 int kxNewMax= nx/2;
11384 if (OutBounds2< 0) kxNewMax=0;
11385 else if (OutBounds2<(nx*nx/4)) kxNewMax=sqrt(OutBounds2);
11386 for (int kxN = kxNewMax; kxN < nx/2 ; kxN++ ) {
11387 des_data[2*kxN + nx* kyN +nxny*kzN] = 0;
11388 des_data[2*kxN + nx* kyN +nxny*kzN+1] = 0;
11389 }
11390
11391
11392 for (int kxNew = 0; kxNew < kxNewMax; kxNew++ ) {
11393/* printf(" kxNew = %d, kyNew = %d, kzNew = %d,kxOld = %3.2f, kyOld = %3.2f, kzOld = %3.2f \n",
11394 kxNew,kyNew,kzNew, kxOld, kyOld, kzOld);*/
11395
11396 IndexOut +=2;
11397
11398 kxOld += MatXX ;
11399 kyOld += MatYX ;
11400 kzOld += MatZX ;
11401 phase += phaseConstx;
11402 Cphase = Cphase*k1 -Sphase*k2; //update using trig addition; this is cos = cos cos -sin sin
11403 Sphase = Sphase*k3+ Cphase*k4; // and sin = sin (1/ cos) + cos * tan;
11404
11405 //
11406 if ((abs(kxOld)>=Mid) || (abs(kyOld)>=Mid) || (abs(kzOld)>=Mid) ) { // out of bounds
11407 des_data[IndexOut] = 0;
11408 des_data[IndexOut+1] = 0;
11409 continue;}
11410/* if (kxOld*kxOld+OutBounds> 2*Mid*Mid){ // out of bounds
11411 des_data[IndexOut] = 0;
11412 des_data[IndexOut+1] = 0;
11413 continue;} */
11414 //
11415 int kxLower= floor(kxOld);
11416 int kyLower= floor(kyOld);
11417 int kzLower= floor(kzOld);
11418 //
11419 float dkxUpper= (kxOld-kxLower);
11420 float dkyUpper= (kyOld-kyLower);
11421 float dkzUpper= (kzOld-kzLower);
11422
11423 // LLL 1
11424 kx= kxLower; ky =kyLower; kz=kzLower;
11425// dataRLLL=0; dataILLL=0;
11426// if ( (abs(kx)<Mid) && (abs(ky)<Mid) && (abs(kz)<Mid) ) { // Step 2 Make sure to be in First BZ
11427 int flagLLL=1;
11428 if (kx<0) kx += N; if (ky<0) ky += N; if (kz<0) kz += N;
11429 if (kx> N/2){kx=N-kx;ky=(N-ky)%N;kz=(N-kz)%N; flagLLL=-1;} // Step 3: if nec, use Friedel paired
11430 int kLLL =2*kx + nx*ky+ nxny*kz ;
11431// dataRLLL = src_data[ kLLL ];
11432// dataILLL = flagLLL*src_data[ 1 + kLLL ]; //
11433// }
11434
11435 // LLU 2
11436 kx= kxLower; ky =kyLower; kz=kzLower+1;
11437// dataRLLU=0; dataILLU=0;
11438// if ( (abs(kx)<Mid) && (abs(ky)<Mid) && (abs(kz)<Mid) ) { // Step 2 Make sure to be in First BZ
11439 int flagLLU=1;
11440 if (kx<0) kx += N; if (ky<0) ky += N; if (kz<0) kz += N;
11441 if (kx> N/2){kx=N-kx;ky=(N-ky)%N;kz=(N-kz)%N; flagLLU=-1;} // Step 3: if nec, use Friedel paired
11442 int kLLU =2*kx + nx*ky+ nxny*kz ;
11443// dataRLLU = src_data[ kLLU];
11444// dataILLU = flagLLU*src_data[1+kLLU];
11445// }
11446
11447 // LUL 3
11448 kx= kxLower; ky =kyLower+1; kz=kzLower;
11449// dataRLUL=0; dataILUL=0;
11450// if ( (abs(kx)<Mid) && (abs(ky)<Mid) && (abs(kz)<Mid) ) { // Step 2 Make sure to be in First BZ
11451 int flagLUL=1;
11452 if (kx<0) kx += N; if (ky<0) ky += N; if (kz<0) kz += N;
11453 if (kx> N/2){kx=N-kx;ky=(N-ky)%N;kz=(N-kz)%N; flagLUL=-1;} // Step 3: if nec, use Friedel paired
11454 int kLUL =2*kx + nx*ky+ nxny*kz ;
11455// dataRLUL = src_data[ kLUL];
11456// dataILUL = flagLUL*src_data[1+kLUL];
11457// }
11458
11459 // LUU 4
11460 kx= kxLower; ky =kyLower+1; kz=kzLower+1;
11461// dataRLUU=0; dataILUU=0;
11462// if ( (abs(kx)<Mid) && (abs(ky)<Mid) && (abs(kz)<Mid) ) { // Step 2 Make sure to be in First BZ
11463 int flagLUU=1;
11464 if (kx<0) kx += N; if (ky<0) ky += N; if (kz<0) kz += N;
11465 if (kx> N/2){kx=N-kx;ky=(N-ky)%N;kz=(N-kz)%N; flagLUU=-1;} // Step 3: if nec, use Friedel paired
11466 int kLUU =2*kx + nx*ky+ nxny*kz ;
11467// dataRLUU = src_data[ kLUU];
11468// dataILUU = flagLUU*src_data[1+kLUU];
11469// }
11470
11471 // ULL 5
11472 kx= kxLower+1; ky =kyLower; kz=kzLower;
11473// dataRULL=0; dataIULL=0;
11474// if ( (abs(kx)<Mid) && (abs(ky)<Mid) && (abs(kz)<Mid) ) { // Step 2 Make sure to be in First BZ
11475 int flagULL=1;
11476 if (kx<0) kx += N; if (ky<0) ky += N; if (kz<0) kz += N;
11477 if (kx> N/2){kx=N-kx;ky=(N-ky)%N;kz=(N-kz)%N; flagULL=-1;} // Step 3: if nec, use Friedel paired
11478 int kULL =2*kx + nx*ky+ nxny*kz ;
11479// dataRULL = src_data[ kULL];
11480// dataIULL = flagULL*src_data[1+kULL];
11481// }
11482
11483 // ULU 6
11484 kx= kxLower+1; ky =kyLower; kz=kzLower+1;
11485// dataRULU=0; dataIULU=0;
11486// if ( (abs(kx)<Mid) && (abs(ky)<Mid) && (abs(kz)<Mid) ) { // Step 2 Make sure to be in First BZ
11487 int flagULU=1;
11488 if (kx<0) kx += N; if (ky<0) ky += N;if (kz<0) kz += N;
11489 if (kx> N/2){kx=N-kx;ky=(N-ky)%N;kz=(N-kz)%N; flagULU=-1;} // Step 3: if nec, use Friedel paired
11490 int kULU =2*kx + nx*ky+ nxny*kz ;
11491// dataRULU = src_data[ kULU];
11492// dataIULU = flagULU*src_data[1+kULU];
11493// }
11494
11495 // UUL 7
11496 kx= kxLower+1; ky =kyLower+1; kz=kzLower;
11497// dataRUUL=0; dataIUUL=0;
11498// if ( (abs(kx)<Mid) && (abs(ky)<Mid) && (abs(kz)<Mid) ) { // Step 2 Make sure to be in First BZ
11499 int flagUUL=1;
11500 if (kx<0) kx += N; if (ky<0) ky += N;if (kz<0) kz += N;
11501 if (kx> N/2){kx=N-kx;ky=(N-ky)%N;kz=(N-kz)%N; flagUUL=-1;} // Step 3: if nec, use Friedel paired
11502 int kUUL =2*kx + nx*ky+ nxny*kz ;
11503// dataRUUL = src_data[ kUUL];
11504// dataIUUL = flagUUL*src_data[1+kUUL];
11505// }
11506
11507 // UUU 8
11508 kx= kxLower+1; ky =kyLower+1; kz=kzLower+1;
11509// dataRUUU=0; dataIUUU=0;
11510// if ( (abs(kx)<Mid) && (abs(ky)<Mid) && (abs(kz)<Mid) ) { // Step 2 Make sure to be in First BZ
11511 int flagUUU=1;
11512 if (kx<0) kx += N; if (ky<0) ky += N;if (kz<0) kz += N;
11513 if (kx> N/2){kx=N-kx;ky=(N-ky)%N;kz=(N-kz)%N; flagUUU=-1;} // Step 3: if nec, use Friedel paired
11514 int kUUU =2*kx + nx*ky+ nxny*kz ;
11515// dataRUUU = src_data[ kUUU];
11516// dataIUUU = flagUUU*src_data[1+kUUU];
11517// }
11518
11519 // Step 4 Assign Real, then Imaginary Values
11520/* printf(" kxNew = %d, kyNew = %d, kzNew = %d,kxOld = %3.2f, kyOld = %3.2f, kzOld = %3.2f \n",
11521 kxNew,kyNew,kzNew, kxOld, kyOld, kzOld);
11522 printf(" dataRLLL= %f, dataRULL = %f, dataRLUL = %f,dataRUUL = %f, dataRLLU = %f, dataRULU = %f, dataRLUU = %f, dataRUUU = %f \n",
11523 dataRLLL,dataRULL,dataRLUL, dataRUUL,dataRLLU, dataRULU,dataRLUU, dataRUUU);
11524 printf(" dataILLL= %f, dataIULL = %f, dataILUL = %f,dataIUUL = %f, dataILLU = %f, dataIULU = %f, dataILUU = %f, dataIUUU = %f \n",
11525 dataILLL,dataIULL,dataILUL, dataIUUL,dataILLU, dataIULU,dataILUU, dataIUUU);*/
11526/* printf(" kLLL= %f, kULL = %d, kLUL = %d,kUUL = %d, kLLU = %d, kULU = %d, kLUU = %d, kUUU = %d \n",
11527 kLLL,kULL,kLUL, kUUL,kLLU, kULU,kLUU, kUUU);*/
11528
11529 tempR = Util::trilinear_interpolate(
11530 src_data[kLLL] , src_data[kULL],
11531 src_data[kLUL] , src_data[kUUL],
11532 src_data[kLLU] , src_data[kULU],
11533 src_data[kLUU] , src_data[kUUU],
11534 dkxUpper,dkyUpper,dkzUpper);
11535
11536 tempI = Util::trilinear_interpolate(
11537 flagLLL*src_data[kLLL+1], flagULL*src_data[kULL+1],
11538 flagLUL*src_data[kLUL+1], flagUUL*src_data[kUUL+1],
11539 flagLLU*src_data[kLLU+1], flagULU*src_data[kULU+1],
11540 flagLUU*src_data[kLUU+1], flagUUU*src_data[kUUU+1],
11541 dkxUpper,dkyUpper,dkzUpper);
11542
11543 // Step 5 Apply translation
11544 float tempRb=tempR*Cphase - tempI*Sphase;
11545 float tempIb=tempR*Sphase + tempI*Cphase;
11546 //
11547 des_data[IndexOut] = tempRb;
11548 des_data[IndexOut+1] = tempIb;
11549 }}} // end z, y, x loops through new coordinates
11550 } // end rotations in Fourier Space 3D
11551 if ((nz > 1)&&(image -> is_real())) {
11552 size_t l=0, ii, k0, k1, k2, k3, k4, k5, k6, k7;
11553 Vec3f offset(nx/2,ny/2,nz/2);
11554 float x2, y2, z2, tuvx, tuvy, tuvz;
11555 int ix, iy, iz;
11556 for (int k = 0; k < nz; ++k) {
11557 for (int j = 0; j < ny; ++j) {
11558 for (int i = 0; i < nx; ++i,++l) {
11559 Vec3f coord(i-nx/2,j-ny/2,k-nz/2);
11560 Vec3f soln = inv*coord;
11561 soln += offset;
11562
11563 x2 = soln[0];
11564 y2 = soln[1];
11565 z2 = soln[2];
11566
11567 if (x2 < 0 || y2 < 0 || z2 < 0 || x2 >= nx || y2 >= ny || z2>= nz ) {
11568 des_data[l] = 0;
11569 }
11570 else {
11571 ix = Util::fast_floor(x2);
11572 iy = Util::fast_floor(y2);
11573 iz = Util::fast_floor(z2);
11574 tuvx = x2-ix;
11575 tuvy = y2-iy;
11576 tuvz = z2-iz;
11577 ii = ix + iy * nx + iz * nxy;
11578
11579 k0 = ii;
11580 k1 = k0 + 1;
11581 k2 = k0 + nx;
11582 k3 = k0 + nx+1;
11583 k4 = k0 + nxy;
11584 k5 = k1 + nxy;
11585 k6 = k2 + nxy;
11586 k7 = k3 + nxy;
11587
11588 if (ix == nx - 1) {
11589 k1--;
11590 k3--;
11591 k5--;
11592 k7--;
11593 }
11594 if (iy == ny - 1) {
11595 k2 -= nx;
11596 k3 -= nx;
11597 k6 -= nx;
11598 k7 -= nx;
11599 }
11600 if (iz == nz - 1) {
11601 k4 -= nxy;
11602 k5 -= nxy;
11603 k6 -= nxy;
11604 k7 -= nxy;
11605 }
11606
11607 des_data[l] = Util::trilinear_interpolate(src_data[k0],
11608 src_data[k1], src_data[k2], src_data[k3], src_data[k4],
11609 src_data[k5], src_data[k6], src_data[k7], tuvx, tuvy, tuvz);
11610 }
11611 }
11612 }
11613 }
11614 }
11615
11616 EXITFUNC;
11617 return des_data;
11618}
EMAN::Dict::get
EMObject get(const string &key) const
Get the EMObject corresponding to the particular key Probably better to just use operator[].
Definition: emobject.h:527
EMAN::EMUtil::em_calloc
static void * em_calloc(const size_t nmemb, const size_t size)
Definition: emutil.h:370
EMAN::Transform::get_mirror
bool get_mirror() const
Query whether x_mirroring is occurring.
Definition: transform.cpp:1250
EMAN::Transform::get_rotation
Dict get_rotation(const string &euler_type="eman") const
Get a rotation in any Euler format.
Definition: transform.cpp:829
EMAN::Transform::get_trans
Vec3f get_trans() const
Get the post trans as a vec3f.
Definition: transform.cpp:1046
EMAN::Util::trilinear_interpolate_complex
static std::complex< float > trilinear_interpolate_complex(std::complex< float > p1, std::complex< float > p2, std::complex< float > p3, std::complex< float > p4, std::complex< float > p5, std::complex< float > p6, std::complex< float > p7, std::complex< float > p8, float t, float u, float v)
Calculate trilinear interpolation.
Definition: util.h:645
EMAN::Util::fast_floor
static int fast_floor(float x)
A fast way to calculate a floor, which is largest integral value not greater than argument.
Definition: util.h:874
EMAN::Util::trilinear_interpolate
static float trilinear_interpolate(float p1, float p2, float p3, float p4, float p5, float p6, float p7, float p8, float t, float u, float v)
Calculate trilinear interpolation.
Definition: util.h:619
EMAN::Util::bilinear_interpolate
static float bilinear_interpolate(float p1, float p2, float p3, float p4, float t, float u)
Calculate bilinear interpolation.
Definition: util.h:543
EMAN::Util::square_sum
static float square_sum(float x, float y)
Calcuate (x*x + y*y).
Definition: util.h:764
EMAN::Vec2
The Vec2 is precisely the same as Vec3 except it works exclusively in 2D Note there are convenient ty...
Definition: vec3.h:710
sqrt
EMData * sqrt() const
return square root of current image
Definition: emdata_core.cpp:1084
phase
EMData * phase() const
return phase part of a complex image as a real image format
Definition: emdata_core.cpp:1384
is_real
bool is_real() const
Is this a real image?
Definition: emdata_metadata.h:667
is_complex
bool is_complex() const
Is this a complex image?
Definition: emdata_metadata.h:648

References EMAN::Util::bilinear_interpolate(), EMAN::EMUtil::em_calloc(), ENTERFUNC, EXITFUNC, EMAN::Util::fast_floor(), EMAN::Dict::get(), EMAN::Transform::get_mirror(), EMAN::Transform::get_rotation(), EMAN::Transform::get_trans(), EMAN::Transform::inverse(), is_complex(), is_real(), EMAN::Processor::params, phase(), EMAN::Dict::set_default(), sqrt(), EMAN::Util::square_sum(), EMAN::Util::trilinear_interpolate(), and EMAN::Util::trilinear_interpolate_complex().

Referenced by process(), and process_inplace().

Member Data Documentation

◆ NAME

const string TransformProcessor::NAME = "xform"
static

Definition at line 2869 of file processor.h.

Referenced by get_name().

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.2