#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 | |
Processor * | NEW () |
Static Public Attributes | |
const string | NAME = "rotateinfs" |
|
Get the descrition of this specific processor. This function must be overwritten by a subclass.
Implements EMAN::Processor. Definition at line 7077 of file processor.h. 07078 { 07079 return "Rotates a Fourier object using a kernel."; 07080 }
|
|
Get the processor's name. Each processor is identified by a unique name.
Implements EMAN::Processor. Definition at line 7069 of file processor.h. 07070 {
07071 return NAME;
07072 }
|
|
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 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 }
|
|
Definition at line 7073 of file processor.h. 07074 { 07075 return new RotateInFSProcessor( ); 07076 }
|
|
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.
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 }
|
|
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.
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 }
|
|
Definition at line 245 of file processor.cpp. |