Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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

ProcessorNEW ()

Static Public Attributes

const string NAME = "rotateinfs"

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

07078                         {
07079                                 return "Rotates a Fourier object using a kernel.";
07080                         }

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

07070                         {
07071                                 return NAME;
07072                         }

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

References EMAN::TypeDict::put().

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

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

Definition at line 7073 of file processor.h.

07074                         {
07075                                 return new RotateInFSProcessor( );
07076                         }

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 9980 of file processor.cpp.

References copy(), and process_inplace().

09981 {
09982 
09983     EMData* imageCp        = image -> copy(); // This is the rotated image
09984     process_inplace(imageCp);
09985     
09986     return imageCp;
09987 }

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 9990 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(), pi, EMAN::Dict::set_default(), set_value_at(), EMAN::EMData::set_value_at(), sqrt(), EMAN::EMData::to_zero(), update(), EMAN::Vec3f, and write_image().

Referenced by process().

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


Member Data Documentation

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

Definition at line 245 of file processor.cpp.


The documentation for this class was generated from the following files:
Generated on Fri Aug 10 16:38:20 2012 for EMAN2 by  doxygen 1.3.9.1