EMAN::RotateInFSProcessor Class Reference

#include <processor.h>

Inheritance diagram for EMAN::RotateInFSProcessor:

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

Collaboration graph
[legend]
List of all members.

Public Member Functions

virtual void process_inplace (EMData *image)
 To process an image in-place.
virtual EMDataprocess (const EMData *const image)
 To proccess an image out-of-place.
virtual string get_name () const
 Get the processor's name.
string get_desc () const
 Get the descrition of this specific processor.
virtual TypeDict get_param_types () const
 Get processor parameter information in a dictionary.

Static Public Member Functions

static ProcessorNEW ()

Static Public Attributes

static const string NAME = "rotateinfs"

Detailed Description

Definition at line 7060 of file processor.h.


Member Function Documentation

string EMAN::RotateInFSProcessor::get_desc (  )  const [inline, virtual]

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 7075 of file processor.h.

07076                         {
07077                                 return "Rotates a Fourier object using a kernel.";
07078                         }

virtual string EMAN::RotateInFSProcessor::get_name (  )  const [inline, virtual]

Get the processor's name.

Each processor is identified by a unique name.

Returns:
The processor's name.

Implements EMAN::Processor.

Definition at line 7067 of file processor.h.

References NAME.

07068                         {
07069                                 return NAME;
07070                         }

virtual TypeDict EMAN::RotateInFSProcessor::get_param_types (  )  const [inline, virtual]

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 7079 of file processor.h.

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

07080                         {
07081                                 TypeDict d;     
07082                                 d.put("transform", EMObject::TRANSFORM, "transform");
07083                                 d.put("interpCutoff", EMObject::FLOAT, "cutoff for interpolation");
07084 //                              d.put("offset", EMObject::FLOAT, "offset for FT centering");
07085 //                              d.put("angle", EMObject::FLOAT, "angle");
07086                                 return d;
07087                         }

static Processor* EMAN::RotateInFSProcessor::NEW (  )  [inline, static]

Definition at line 7071 of file processor.h.

07072                         {
07073                                 return new RotateInFSProcessor( );
07074                         }

EMData * RotateInFSProcessor::process ( const EMData *const   image  )  [virtual]

To proccess an image out-of-place.

For those processors which can only be processed out-of-place, override this function to give the right behavior.

Parameters:
image The image will be copied, actual process happen on copy of image.
Returns:
the image processing result, may or may not be the same size of the input image

Reimplemented from EMAN::Processor.

Definition at line 9965 of file processor.cpp.

References copy(), and process_inplace().

09966 {
09967 
09968     EMData* imageCp        = image -> copy(); // This is the rotated image
09969     process_inplace(imageCp);
09970     
09971     return imageCp;
09972 }

void RotateInFSProcessor::process_inplace ( EMData image  )  [virtual]

To process an image in-place.

For those processors which can only be processed out-of-place, override this function to just print out some error message to remind user call the out-of-place version.

Parameters:
image The image to be processed.

Implements EMAN::Processor.

Definition at line 9975 of file processor.cpp.

References abs, copy(), EMAN::EMData::get_value_at(), get_value_at(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), ImageDimensionException, EMAN::Transform::inverse(), EMAN::EMData::is_real(), EMAN::Processor::params, pi, EMAN::Dict::set_default(), set_value_at(), EMAN::EMData::set_value_at(), sqrt(), EMAN::EMData::to_zero(), update(), and write_image().

Referenced by process().

09976 {
09977 //      float angle = params["angle"];
09978         
09979 
09980 //      Transform *rotNow  = params.set_default("transform",&Transform());
09981         Transform *rotNow  = params["transform"];
09982         float interpCutoff = params.set_default("interpCutoff",0.8f);   // JFF, can you move this to input parameter? 
09983 //      float interpCutoff = params["interpCutoff"];   // JFF, can you move this to input parameter? 
09984 //      float interpCutoff =.8;   // JFF, can you move this to input parameter? 
09985         
09986         
09987         int debug=0;
09988         
09989 //      error: conversion from ‘EMAN::EMObject’ to non-scalar type ‘EMAN::Transform’ requested
09990         
09991         
09992         // if 2N is size of image, then sizes of FFT are (2N+2,2N,2N) or  (2N+2,2N,1)
09993         // if 2N+1 is size of image, then sizes of FFT are (2N+2,2N+1,2N+1) or  (2N+2,2N+1,1)
09994         
09995         int x_size = image->get_xsize();  //16
09996         int y_size = image->get_ysize();  int y_odd=  (y_size%2);
09997         int z_size = image->get_zsize();
09998         
09999 //      float size_check = abs(y_size-z_size)+abs(x_size-y_size-2+y_odd);
10000 //      if (size_check!=0) throw ImageDimensionException("Only cubic images");
10001         int size_check = abs(x_size-y_size-2+y_odd)+ abs(z_size-1)*abs(z_size-y_size);
10002         int N = x_size/2-1;
10003         int Mid = N+1;
10004         if (size_check!=0) throw ImageDimensionException("Only square or cubic  images for now");
10005         if (image->is_real()) throw ImageDimensionException("Only for Fourier images");
10006 //      if (y_odd==0) throw ImageDimensionException("Only use odd images for now");
10007         
10008         if (debug) printf("Mid=%d, x_size=%d, y_size=%d, N=%d, z_size=%d \n", Mid,x_size,y_size, N, z_size );
10009         
10010         EMData* RotIm        = image -> copy(); // This is the rotated image
10011         EMData* WeightIm     = image -> copy(); WeightIm     ->to_zero();// This is the Weight image for rotated image
10012         EMData* PhaseOrigIm  = new EMData(N+1,2*N+1,z_size) ; PhaseOrigIm  ->to_zero();// This is the Weight image for rotated image
10013         EMData* PhaseFinalIm = PhaseOrigIm -> copy(); PhaseFinalIm ->to_zero();// This is the Weight image for rotated image
10014         EMData* MagFinalIm   = PhaseOrigIm -> copy(); MagFinalIm   ->to_zero();// This is the Weight image for rotated image
10015         
10016 //      float*   data      = image    -> get_data();
10017 //      float* Rotdata     = RotIm    -> get_data();  // This is the data of the rotated image
10018 //      float* WeightData  = WeightIm -> get_data();  //
10019         
10020         float  WeightNowX, WeightNowY, WeightNowZ ;
10021         int    kxMin,kxMax, kyMin, kyMax,  kzMin, kzMax, kxBefore, kyBefore, kzBefore ;
10022         float  kxRT, kyRT, kzRT ;
10023         
10024         Vec3f PosAfter;
10025         Vec3f PosBefore;
10026         Transform invRotNow;
10027         // Fill out real and imaginary full images
10028 
10029         //int kz=0;
10030 
10031         if (debug) {image -> write_image("OrigImageFFT.hdf");   
10032                     printf("Just wrote OrigImageFFT.hdf \n"); }
10033         
10034 
10035         for (kxBefore = 0; kxBefore <=  N      ; ++kxBefore) {  // These are the  kx, ky coordinates of the original image
10036            for (kyBefore = 0; kyBefore < y_size ; ++kyBefore)  {         //  We need to rephase
10037               for (kzBefore = 0; kzBefore < z_size ; ++kzBefore)  {         //  We need to rephase
10038 
10039             // Now we need a 
10040               float CurrentReal = RotIm -> get_value_at(2*kxBefore  ,kyBefore, kzBefore);
10041               float CurrentImag = RotIm -> get_value_at(2*kxBefore+1,kyBefore, kzBefore);
10042               
10043 //         fftOfxPRB3(1+mod(ik-Mid,N))=Grand*exp(-2*pi*1i*Mid*(ik-Mid)/x_size); % Phase to apply to centered version
10044               
10045 //            float Phase    = -2*pi*(kxBefore+kyBefore + kzBefore)*(Mid)/y_size;
10046               float Phase    = -pi*(kxBefore+kyBefore + kzBefore)*x_size/y_size;
10047               // Phase    = 0;
10048               float CosPhase = cos( Phase);
10049               float SinPhase = sin( Phase);
10050              
10051               float NewRealValue = CosPhase*CurrentReal -SinPhase*CurrentImag;
10052               float NewImagValue = SinPhase*CurrentReal +CosPhase*CurrentImag;
10053               
10054               RotIm ->set_value_at(2*kxBefore  ,kyBefore, kzBefore, NewRealValue);
10055               RotIm ->set_value_at(2*kxBefore+1,kyBefore, kzBefore, NewImagValue);
10056         }}}
10057         
10058         if (debug) {RotIm  -> write_image("OrigImageFFTAfterPhaseCorrection.hdf"); 
10059                     printf("  Just wrote OrigImageFFTAfterPhaseCorrection.hdf \n");}
10060         
10061         // RotIm ->set_value_at(2*Mid-1,0, 0, 0);
10062         if (debug) printf("  Just about to start second loop  \n");
10063         
10064         image ->to_zero();
10065         invRotNow = rotNow ->inverse(); //  no match for ‘operator=’ in ‘PosBefore = EMAN::operator*(const EMAN::Transform&, const EMAN::Transform&)((
10066          
10067         
10068         for (int kxAfter = 0; kxAfter <= N  ; ++kxAfter) {  // These are the  kx, ky, kz coordinates of the rotated image
10069           for (int kyAfter = -N; kyAfter < y_size-N     ; ++kyAfter)  {     // referring to a properly centered version
10070              for (int kzAfter = -z_size/2; kzAfter <= z_size/2  ; ++kzAfter)  {   
10071 
10072             // Now we need a 
10073 
10074               PosAfter = Vec3f(kxAfter, kyAfter, kzAfter);
10075               PosBefore = invRotNow*PosAfter;
10076               kxRT = PosBefore[0]; // This will be the off-lattice site, where the point was rotated from
10077               kyRT = PosBefore[1]; //  
10078               kzRT = PosBefore[2]; //  
10079               
10080               
10081               kxMin = ceil( kxRT-interpCutoff); kxMax = floor(kxRT+interpCutoff);
10082               kyMin = ceil( kyRT-interpCutoff); kyMax = floor(kyRT+interpCutoff);
10083               kzMin = ceil( kzRT-interpCutoff); kzMax = floor(kzRT+interpCutoff);
10084 
10085     
10086 //              printf("Block 0,kx=%d, ky=%d,kxMin=%d, kyMin=%d, kxMax=%d, kyMax=%d, kyAfter=%d  \n",kxAfter,kyAfter,kxMin,kyMin, kxMax, kyMax, kyAfter); 
10087               //continue; 
10088               for (int kxI= kxMin; kxI <= kxMax; ++kxI){  // go through this
10089                 for (int kyI= kyMin; kyI <= kyMax; ++kyI){  // and get values to interp
10090                    for (int kzI= kzMin; kzI <= kzMax; ++kzI){ //  
10091 
10092 // 
10093                      if (abs(kxI) >N) continue; // don't go off the lattice
10094                      if (abs(kyI) >N) continue;
10095                      if (abs(kzI) >z_size/2) continue;
10096                      
10097                      float distx= abs(kxI-kxRT);
10098                      float disty= abs(kyI-kyRT);
10099                      float distz= abs(kzI-kzRT);
10100 
10101                      // fold kxI, kyI back into lattice if possible
10102                      int IsComplexConj= 1;
10103                      
10104                      if (kxI<0) IsComplexConj=-1;
10105                      kxBefore= IsComplexConj*kxI; // The Proper coordinates will be
10106                      kyBefore= IsComplexConj*kyI; // where the original data is stored 
10107                      kzBefore= IsComplexConj*kzI; // At this point kxBefore >=0, but not necessarily kyBefore
10108                      
10109                      if ( kyBefore<0 ) kyBefore += y_size; // makes sure kyBefore is also >0
10110                      if ( kzBefore<0 ) kzBefore += y_size; // makes sure kzBefore is also >0
10111                      
10112                      WeightNowX  = (distx ==0)? 1: (sin(pi*distx) /(pi*distx)) ;
10113                      WeightNowY  = (disty ==0)? 1: (sin(pi*disty) /(pi*disty)) ;
10114                      WeightNowZ  = (distz ==0)? 1: (sin(pi*distz) /(pi*distz)) ;
10115                      
10116                      
10117                      float CurrentValue;
10118                      float ToAdd;
10119                      int kyAfterInd = (kyAfter+y_size)%(y_size);
10120                      int kzAfterInd = (kzAfter+z_size)%(z_size);
10121                      
10122                     // if (kxAfter==0) IsComplexConj*=-1;
10123                      
10124 //                   if ((kxI+kyI)%1 ==0)
10125 //                       printf("Block5,kx=%d, ky=%d,kxI=%d, kyI=%d, kxBefore=%d, kyBefore=%d  \n",kxAfter,kyAfter,kxI,kyI, kxBefore,kyBefore);
10126 //                       printf("  %d,     %d,  %d,  %d,        %d,  %d,       %d, %d \n",IsComplexConj,kxAfter,kyAfter, kyAfterInd,kxI,kyI, kxBefore,kyBefore);
10127                      
10128                      CurrentValue =   image -> get_value_at(2*kxAfter,kyAfterInd, kzAfterInd);  // Update real part of Image
10129                      ToAdd =   WeightNowX*WeightNowY*WeightNowZ*(RotIm -> get_value_at(2*kxBefore,kyBefore, kzBefore));
10130                      image -> set_value_at(2*kxAfter  ,kyAfterInd  , kzAfterInd,  ToAdd   + CurrentValue );
10131                      
10132                      
10133                      CurrentValue =   WeightIm -> get_value_at(kxAfter,kyAfterInd, kzAfterInd);    // Update real  part of Weight image
10134                      ToAdd =   WeightNowX*WeightNowY;
10135                      WeightIm -> set_value_at(kxAfter  , kyAfterInd , kzAfterInd,  abs(ToAdd)   + CurrentValue );
10136                      
10137                      CurrentValue = image -> get_value_at(2*kxAfter+1,kyAfterInd);    // Update imaginary   part of image
10138                      ToAdd =   IsComplexConj*WeightNowX*WeightNowY*RotIm -> get_value_at(2*kxBefore+1,kyBefore, kzBefore );
10139                      image -> set_value_at(2*kxAfter+1  , kyAfterInd , kzAfterInd,  ToAdd   + CurrentValue );
10140                      
10141                       
10142                 }}}
10143 
10144             
10145           }}}
10146 
10147 //        Set the image values to the rotated image, because we do it in place
10148 
10149 
10150         if (debug) { image        -> write_image("RotImageBeforeFinalPhaseCorrection.hdf"); 
10151                      printf("  Just wrote RotImageBeforeFinalPhaseCorrection.hdf \n");   }
10152 
10153 
10154         for (kxBefore = 0; kxBefore <= N     ; ++kxBefore) {      // This is  the normalization step
10155           for (kyBefore = 0; kyBefore < y_size   ; ++kyBefore)  {  // These are the  kx, ky, kz coordinates of the original image       
10156               for (kzBefore = 0; kzBefore < z_size   ; ++kzBefore)  {       
10157 
10158               float CurrentReal = image -> get_value_at(2*kxBefore   , kyBefore, kzBefore);
10159               float CurrentImag = image -> get_value_at(2*kxBefore+1 , kyBefore, kzBefore);
10160               
10161               PhaseFinalIm -> set_value_at(kxBefore,kyBefore, kzBefore, atan2(CurrentImag,CurrentReal));
10162               MagFinalIm   -> set_value_at(kxBefore,kyBefore, kzBefore, sqrt(CurrentImag*CurrentImag+CurrentReal*CurrentReal) );
10163               float WeightNow    = WeightIm -> get_value_at(kxBefore,kyBefore, kzBefore);
10164               if (WeightNow>0) {
10165                  float val =  (image->get_value_at(2*kxBefore   , kyBefore, kzBefore))/WeightNow;
10166                  image -> set_value_at(2*kxBefore   , kyBefore, kzBefore, val);
10167                  val      =  (image->get_value_at(2*kxBefore +1  , kyBefore, kzBefore))/WeightNow;
10168                  image -> set_value_at(2*kxBefore +1  , kyBefore, kzBefore,  val);
10169               }
10170         
10171         }}}
10172         
10173         if (debug) { printf("  Just did normalization step \n");}
10174         
10175 
10176         for ( kxBefore = 0; kxBefore < Mid     ; ++kxBefore) {     //  This is the rephase step
10177           for ( kyBefore = 0; kyBefore < y_size   ; ++kyBefore)  {     
10178             for ( kzBefore = 0; kzBefore < z_size   ; ++kzBefore)  {     
10179 
10180               float CurrentReal = image -> get_value_at(2*kxBefore  ,kyBefore, kzBefore);
10181               float CurrentImag = image -> get_value_at(2*kxBefore+1,kyBefore, kzBefore);
10182               
10183 //            float Phase    = +2*pi*(kxBefore+kyBefore+kzBefore)*(Mid)/y_size;
10184               float Phase    =  pi*(kxBefore + kyBefore + kzBefore)*x_size/y_size;
10185               // Phase    = 0; // Offset should be Mid-1
10186               float CosPhase = cos( Phase);
10187               float SinPhase = sin( Phase);
10188              
10189               float NewRealValue = CosPhase*CurrentReal -SinPhase*CurrentImag;
10190               float NewImagValue = SinPhase*CurrentReal +CosPhase*CurrentImag;
10191               
10192               image ->set_value_at(2*kxBefore,  kyBefore,  kzBefore, NewRealValue);
10193               image ->set_value_at(2*kxBefore+1,kyBefore,  kzBefore, NewImagValue);
10194         }}}
10195         
10196         if (debug) {
10197            image        -> write_image("RotatedImageFFT.hdf");
10198            PhaseFinalIm -> write_image("PhaseImInFS.hdf");   // These are the phases,mags of the image when properly centered
10199            MagFinalIm   -> write_image("MagFinalImInFS.hdf"); 
10200            WeightIm     -> write_image("WeightIm.hdf");
10201            printf("  Just wrote RotatedImageFFT.hdf \n");
10202         }
10203 
10204         image -> update();
10205         
10206 }


Member Data Documentation

const string RotateInFSProcessor::NAME = "rotateinfs" [static]

Definition at line 7088 of file processor.h.

Referenced by get_name().


The documentation for this class was generated from the following files:
Generated on Thu May 3 10:11:04 2012 for EMAN2 by  doxygen 1.4.7