#include <processor.h>
Inheritance diagram for EMAN::RotateInFSProcessor:
Public Member Functions | |
virtual void | process_inplace (EMData *image) |
To process an image in-place. | |
virtual EMData * | process (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 Processor * | NEW () |
Static Public Attributes | |
static const string | NAME = "rotateinfs" |
Definition at line 7101 of file processor.h.
string EMAN::RotateInFSProcessor::get_desc | ( | ) | const [inline, virtual] |
Get the descrition of this specific processor.
This function must be overwritten by a subclass.
Implements EMAN::Processor.
Definition at line 7116 of file processor.h.
virtual string EMAN::RotateInFSProcessor::get_name | ( | ) | const [inline, virtual] |
Get the processor's name.
Each processor is identified by a unique name.
Implements EMAN::Processor.
Definition at line 7108 of file processor.h.
References NAME.
07109 { 07110 return NAME; 07111 }
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.
Reimplemented from EMAN::Processor.
Definition at line 7120 of file processor.h.
References EMAN::EMObject::FLOAT, EMAN::TypeDict::put(), and EMAN::EMObject::TRANSFORM.
07121 { 07122 TypeDict d; 07123 d.put("transform", EMObject::TRANSFORM, "transform"); 07124 d.put("interpCutoff", EMObject::FLOAT, "cutoff for interpolation"); 07125 // d.put("offset", EMObject::FLOAT, "offset for FT centering"); 07126 // d.put("angle", EMObject::FLOAT, "angle"); 07127 return d; 07128 }
static Processor* EMAN::RotateInFSProcessor::NEW | ( | ) | [inline, static] |
Definition at line 7112 of file processor.h.
07113 { 07114 return new RotateInFSProcessor( ); 07115 }
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.
image | The image will be copied, actual process happen on copy of image. |
Reimplemented from EMAN::Processor.
Definition at line 10043 of file processor.cpp.
References copy(), and process_inplace().
10044 { 10045 10046 EMData* imageCp = image -> copy(); // This is the rotated image 10047 process_inplace(imageCp); 10048 10049 return imageCp; 10050 }
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.
image | The image to be processed. |
Implements EMAN::Processor.
Definition at line 10053 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().
10054 { 10055 // float angle = params["angle"]; 10056 10057 10058 // Transform *rotNow = params.set_default("transform",&Transform()); 10059 Transform *rotNow = params["transform"]; 10060 float interpCutoff = params.set_default("interpCutoff",0.8f); // JFF, can you move this to input parameter? 10061 // float interpCutoff = params["interpCutoff"]; // JFF, can you move this to input parameter? 10062 // float interpCutoff =.8; // JFF, can you move this to input parameter? 10063 10064 10065 int debug=0; 10066 10067 // error: conversion from ‘EMAN::EMObject’ to non-scalar type ‘EMAN::Transform’ requested 10068 10069 10070 // if 2N is size of image, then sizes of FFT are (2N+2,2N,2N) or (2N+2,2N,1) 10071 // if 2N+1 is size of image, then sizes of FFT are (2N+2,2N+1,2N+1) or (2N+2,2N+1,1) 10072 10073 int x_size = image->get_xsize(); //16 10074 int y_size = image->get_ysize(); int y_odd= (y_size%2); 10075 int z_size = image->get_zsize(); 10076 10077 // float size_check = abs(y_size-z_size)+abs(x_size-y_size-2+y_odd); 10078 // if (size_check!=0) throw ImageDimensionException("Only cubic images"); 10079 int size_check = abs(x_size-y_size-2+y_odd)+ abs(z_size-1)*abs(z_size-y_size); 10080 int N = x_size/2-1; 10081 int Mid = N+1; 10082 if (size_check!=0) throw ImageDimensionException("Only square or cubic images for now"); 10083 if (image->is_real()) throw ImageDimensionException("Only for Fourier images"); 10084 // if (y_odd==0) throw ImageDimensionException("Only use odd images for now"); 10085 10086 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 ); 10087 10088 EMData* RotIm = image -> copy(); // This is the rotated image 10089 EMData* WeightIm = image -> copy(); WeightIm ->to_zero();// This is the Weight image for rotated image 10090 EMData* PhaseOrigIm = new EMData(N+1,2*N+1,z_size) ; PhaseOrigIm ->to_zero();// This is the Weight image for rotated image 10091 EMData* PhaseFinalIm = PhaseOrigIm -> copy(); PhaseFinalIm ->to_zero();// This is the Weight image for rotated image 10092 EMData* MagFinalIm = PhaseOrigIm -> copy(); MagFinalIm ->to_zero();// This is the Weight image for rotated image 10093 10094 // float* data = image -> get_data(); 10095 // float* Rotdata = RotIm -> get_data(); // This is the data of the rotated image 10096 // float* WeightData = WeightIm -> get_data(); // 10097 10098 float WeightNowX, WeightNowY, WeightNowZ ; 10099 int kxMin,kxMax, kyMin, kyMax, kzMin, kzMax, kxBefore, kyBefore, kzBefore ; 10100 float kxRT, kyRT, kzRT ; 10101 10102 Vec3f PosAfter; 10103 Vec3f PosBefore; 10104 Transform invRotNow; 10105 // Fill out real and imaginary full images 10106 10107 //int kz=0; 10108 10109 if (debug) {image -> write_image("OrigImageFFT.hdf"); 10110 printf("Just wrote OrigImageFFT.hdf \n"); } 10111 10112 10113 for (kxBefore = 0; kxBefore <= N ; ++kxBefore) { // These are the kx, ky coordinates of the original image 10114 for (kyBefore = 0; kyBefore < y_size ; ++kyBefore) { // We need to rephase 10115 for (kzBefore = 0; kzBefore < z_size ; ++kzBefore) { // We need to rephase 10116 10117 // Now we need a 10118 float CurrentReal = RotIm -> get_value_at(2*kxBefore ,kyBefore, kzBefore); 10119 float CurrentImag = RotIm -> get_value_at(2*kxBefore+1,kyBefore, kzBefore); 10120 10121 // fftOfxPRB3(1+mod(ik-Mid,N))=Grand*exp(-2*pi*1i*Mid*(ik-Mid)/x_size); % Phase to apply to centered version 10122 10123 // float Phase = -2*pi*(kxBefore+kyBefore + kzBefore)*(Mid)/y_size; 10124 float Phase = -pi*(kxBefore+kyBefore + kzBefore)*x_size/y_size; 10125 // Phase = 0; 10126 float CosPhase = cos( Phase); 10127 float SinPhase = sin( Phase); 10128 10129 float NewRealValue = CosPhase*CurrentReal -SinPhase*CurrentImag; 10130 float NewImagValue = SinPhase*CurrentReal +CosPhase*CurrentImag; 10131 10132 RotIm ->set_value_at(2*kxBefore ,kyBefore, kzBefore, NewRealValue); 10133 RotIm ->set_value_at(2*kxBefore+1,kyBefore, kzBefore, NewImagValue); 10134 }}} 10135 10136 if (debug) {RotIm -> write_image("OrigImageFFTAfterPhaseCorrection.hdf"); 10137 printf(" Just wrote OrigImageFFTAfterPhaseCorrection.hdf \n");} 10138 10139 // RotIm ->set_value_at(2*Mid-1,0, 0, 0); 10140 if (debug) printf(" Just about to start second loop \n"); 10141 10142 image ->to_zero(); 10143 invRotNow = rotNow ->inverse(); // no match for ‘operator=’ in ‘PosBefore = EMAN::operator*(const EMAN::Transform&, const EMAN::Transform&)(( 10144 10145 10146 for (int kxAfter = 0; kxAfter <= N ; ++kxAfter) { // These are the kx, ky, kz coordinates of the rotated image 10147 for (int kyAfter = -N; kyAfter < y_size-N ; ++kyAfter) { // referring to a properly centered version 10148 for (int kzAfter = -z_size/2; kzAfter <= z_size/2 ; ++kzAfter) { 10149 10150 // Now we need a 10151 10152 PosAfter = Vec3f(kxAfter, kyAfter, kzAfter); 10153 PosBefore = invRotNow*PosAfter; 10154 kxRT = PosBefore[0]; // This will be the off-lattice site, where the point was rotated from 10155 kyRT = PosBefore[1]; // 10156 kzRT = PosBefore[2]; // 10157 10158 10159 kxMin = ceil( kxRT-interpCutoff); kxMax = floor(kxRT+interpCutoff); 10160 kyMin = ceil( kyRT-interpCutoff); kyMax = floor(kyRT+interpCutoff); 10161 kzMin = ceil( kzRT-interpCutoff); kzMax = floor(kzRT+interpCutoff); 10162 10163 10164 // 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); 10165 //continue; 10166 for (int kxI= kxMin; kxI <= kxMax; ++kxI){ // go through this 10167 for (int kyI= kyMin; kyI <= kyMax; ++kyI){ // and get values to interp 10168 for (int kzI= kzMin; kzI <= kzMax; ++kzI){ // 10169 10170 // 10171 if (abs(kxI) >N) continue; // don't go off the lattice 10172 if (abs(kyI) >N) continue; 10173 if (abs(kzI) >z_size/2) continue; 10174 10175 float distx= abs(kxI-kxRT); 10176 float disty= abs(kyI-kyRT); 10177 float distz= abs(kzI-kzRT); 10178 10179 // fold kxI, kyI back into lattice if possible 10180 int IsComplexConj= 1; 10181 10182 if (kxI<0) IsComplexConj=-1; 10183 kxBefore= IsComplexConj*kxI; // The Proper coordinates will be 10184 kyBefore= IsComplexConj*kyI; // where the original data is stored 10185 kzBefore= IsComplexConj*kzI; // At this point kxBefore >=0, but not necessarily kyBefore 10186 10187 if ( kyBefore<0 ) kyBefore += y_size; // makes sure kyBefore is also >0 10188 if ( kzBefore<0 ) kzBefore += y_size; // makes sure kzBefore is also >0 10189 10190 WeightNowX = (distx ==0)? 1: (sin(pi*distx) /(pi*distx)) ; 10191 WeightNowY = (disty ==0)? 1: (sin(pi*disty) /(pi*disty)) ; 10192 WeightNowZ = (distz ==0)? 1: (sin(pi*distz) /(pi*distz)) ; 10193 10194 10195 float CurrentValue; 10196 float ToAdd; 10197 int kyAfterInd = (kyAfter+y_size)%(y_size); 10198 int kzAfterInd = (kzAfter+z_size)%(z_size); 10199 10200 // if (kxAfter==0) IsComplexConj*=-1; 10201 10202 // if ((kxI+kyI)%1 ==0) 10203 // printf("Block5,kx=%d, ky=%d,kxI=%d, kyI=%d, kxBefore=%d, kyBefore=%d \n",kxAfter,kyAfter,kxI,kyI, kxBefore,kyBefore); 10204 // printf(" %d, %d, %d, %d, %d, %d, %d, %d \n",IsComplexConj,kxAfter,kyAfter, kyAfterInd,kxI,kyI, kxBefore,kyBefore); 10205 10206 CurrentValue = image -> get_value_at(2*kxAfter,kyAfterInd, kzAfterInd); // Update real part of Image 10207 ToAdd = WeightNowX*WeightNowY*WeightNowZ*(RotIm -> get_value_at(2*kxBefore,kyBefore, kzBefore)); 10208 image -> set_value_at(2*kxAfter ,kyAfterInd , kzAfterInd, ToAdd + CurrentValue ); 10209 10210 10211 CurrentValue = WeightIm -> get_value_at(kxAfter,kyAfterInd, kzAfterInd); // Update real part of Weight image 10212 ToAdd = WeightNowX*WeightNowY; 10213 WeightIm -> set_value_at(kxAfter , kyAfterInd , kzAfterInd, abs(ToAdd) + CurrentValue ); 10214 10215 CurrentValue = image -> get_value_at(2*kxAfter+1,kyAfterInd); // Update imaginary part of image 10216 ToAdd = IsComplexConj*WeightNowX*WeightNowY*RotIm -> get_value_at(2*kxBefore+1,kyBefore, kzBefore ); 10217 image -> set_value_at(2*kxAfter+1 , kyAfterInd , kzAfterInd, ToAdd + CurrentValue ); 10218 10219 10220 }}} 10221 10222 10223 }}} 10224 10225 // Set the image values to the rotated image, because we do it in place 10226 10227 10228 if (debug) { image -> write_image("RotImageBeforeFinalPhaseCorrection.hdf"); 10229 printf(" Just wrote RotImageBeforeFinalPhaseCorrection.hdf \n"); } 10230 10231 10232 for (kxBefore = 0; kxBefore <= N ; ++kxBefore) { // This is the normalization step 10233 for (kyBefore = 0; kyBefore < y_size ; ++kyBefore) { // These are the kx, ky, kz coordinates of the original image 10234 for (kzBefore = 0; kzBefore < z_size ; ++kzBefore) { 10235 10236 float CurrentReal = image -> get_value_at(2*kxBefore , kyBefore, kzBefore); 10237 float CurrentImag = image -> get_value_at(2*kxBefore+1 , kyBefore, kzBefore); 10238 10239 PhaseFinalIm -> set_value_at(kxBefore,kyBefore, kzBefore, atan2(CurrentImag,CurrentReal)); 10240 MagFinalIm -> set_value_at(kxBefore,kyBefore, kzBefore, sqrt(CurrentImag*CurrentImag+CurrentReal*CurrentReal) ); 10241 float WeightNow = WeightIm -> get_value_at(kxBefore,kyBefore, kzBefore); 10242 if (WeightNow>0) { 10243 float val = (image->get_value_at(2*kxBefore , kyBefore, kzBefore))/WeightNow; 10244 image -> set_value_at(2*kxBefore , kyBefore, kzBefore, val); 10245 val = (image->get_value_at(2*kxBefore +1 , kyBefore, kzBefore))/WeightNow; 10246 image -> set_value_at(2*kxBefore +1 , kyBefore, kzBefore, val); 10247 } 10248 10249 }}} 10250 10251 if (debug) { printf(" Just did normalization step \n");} 10252 10253 10254 for ( kxBefore = 0; kxBefore < Mid ; ++kxBefore) { // This is the rephase step 10255 for ( kyBefore = 0; kyBefore < y_size ; ++kyBefore) { 10256 for ( kzBefore = 0; kzBefore < z_size ; ++kzBefore) { 10257 10258 float CurrentReal = image -> get_value_at(2*kxBefore ,kyBefore, kzBefore); 10259 float CurrentImag = image -> get_value_at(2*kxBefore+1,kyBefore, kzBefore); 10260 10261 // float Phase = +2*pi*(kxBefore+kyBefore+kzBefore)*(Mid)/y_size; 10262 float Phase = pi*(kxBefore + kyBefore + kzBefore)*x_size/y_size; 10263 // Phase = 0; // Offset should be Mid-1 10264 float CosPhase = cos( Phase); 10265 float SinPhase = sin( Phase); 10266 10267 float NewRealValue = CosPhase*CurrentReal -SinPhase*CurrentImag; 10268 float NewImagValue = SinPhase*CurrentReal +CosPhase*CurrentImag; 10269 10270 image ->set_value_at(2*kxBefore, kyBefore, kzBefore, NewRealValue); 10271 image ->set_value_at(2*kxBefore+1,kyBefore, kzBefore, NewImagValue); 10272 }}} 10273 10274 if (debug) { 10275 image -> write_image("RotatedImageFFT.hdf"); 10276 PhaseFinalIm -> write_image("PhaseImInFS.hdf"); // These are the phases,mags of the image when properly centered 10277 MagFinalIm -> write_image("MagFinalImInFS.hdf"); 10278 WeightIm -> write_image("WeightIm.hdf"); 10279 printf(" Just wrote RotatedImageFFT.hdf \n"); 10280 } 10281 10282 image -> update(); 10283 10284 }
const string RotateInFSProcessor::NAME = "rotateinfs" [static] |