#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 7060 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 7075 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 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.
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 }
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 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.
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 }
const string RotateInFSProcessor::NAME = "rotateinfs" [static] |