00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include <stack>
00037 #include "ctf.h"
00038 #include "emdata.h"
00039 #include <iostream>
00040 #include <cmath>
00041 #include <cstring>
00042
00043 #include <gsl/gsl_sf_bessel.h>
00044 #include <gsl/gsl_errno.h>
00045 #include <vector>
00046 #include <boost/shared_ptr.hpp>
00047
00048 using boost::shared_ptr;
00049 using std::vector;
00050 using std::cout;
00051 using namespace EMAN;
00052 using namespace std;
00053
00054 EMData *EMData::real2FH(float OverSamplekB)
00055 {
00056 int nx = get_xsize();
00057 int ny = get_ysize();
00058 int nz = get_zsize();
00059 int Center = (int) floor( (nx+1.0)/2.0 +.01);
00060 #ifdef DEBUG
00061 printf("nx=%d, ny=%d, nz=%d Center=%d\n", nx,ny,nz, Center);
00062 #endif //DEBUG
00063 float ScalFactor=4.1f;
00064 gsl_set_error_handler_off();
00065
00066 if ( (nz==1) && (nx==ny) && (!is_complex()) && (Center*2)==(nx+1)){
00067 #ifdef DEBUG
00068 printf("entered if \n");fflush(stdout);
00069 #endif //DEBUG
00070
00071 EMData* ImBW = this ;
00072 int Size=nx;
00073 int iMax = (int) floor( (Size-1.0)/2 +.01);
00074 int CountMax = (iMax+2)*(iMax+1)/2;
00075 int *PermMatTr = new int[CountMax];
00076 float *RValsSorted = new float[CountMax];
00077 float *weightofkValsSorted = new float[CountMax];
00078 int *SizeReturned = new int[1];
00079 Util::Radialize(PermMatTr, RValsSorted,weightofkValsSorted,Size, SizeReturned);
00080 int RIntMax= SizeReturned[0];
00081
00082 int mMax = (int) floor( ScalFactor*RValsSorted[RIntMax-1]+10.0);
00083
00084 int kIntMax=2+ (int) floor( RValsSorted[RIntMax-1]*OverSamplekB);
00085 float *kVec2Use= new float[kIntMax];
00086 for (int kk=0; kk<kIntMax; kk++){
00087 kVec2Use[kk]= ((float) kk)/OverSamplekB;}
00088
00089 float *krVec= new float[kIntMax*RIntMax];
00090 int Count=0;
00091 for (int jk=0; jk<kIntMax; jk++ ){
00092 for (int jR=0; jR<RIntMax; jR++ ){
00093 krVec[Count]=2.0f*M_PI*RValsSorted[jR]
00094 *kVec2Use[jk]/( (float) Size);
00095 Count++;
00096
00097 }}
00098 float krVecMin= kVec2Use[1]*RValsSorted[1];
00099 float krVecMax = krVec[kIntMax*RIntMax-1]+krVecMin;
00100 int Number2Use = (int) floor(OverSamplekB*krVecMax+1.0);
00101 float *krVec2Use = new float[Number2Use+1];
00102 float *sampledBesselJ = new float[Number2Use+1];
00103 #ifdef DEBUG
00104 printf("Size=%d, iMax=%d, SizeReturned=%d, RIntMax=%d, \n"
00105 "mMax=%d, kIntMax=%d, krVecMin=%f, krVecMax=%f, Number2Use=%d \n\n",
00106 Size, iMax, SizeReturned[0], RIntMax, mMax, kIntMax,
00107 krVecMin,krVecMax,Number2Use);fflush(stdout);
00108 #endif //DEBUG
00109 for (int jkr=0; jkr<= Number2Use; jkr++) {
00110 krVec2Use[jkr] =((float)jkr)*krVecMax/
00111 ((float)Number2Use);
00112
00113 }
00114
00115
00116 EMData* rhoOfkmB = copy();
00117
00118 rhoOfkmB->set_size(2*(mMax+1),kIntMax);
00119 rhoOfkmB->to_zero();
00120
00121
00122 int CenterM= Center-1;
00123 std::complex <float> *rhoOfRandmTemp = new std::complex <float>[RIntMax];
00124 std::complex <float> rhoTemp;
00125
00126 int PCount=0;
00127
00128 for (int m=0; m <=mMax; m++){
00129
00130 std::complex <float> tempF(0.0f,-1.0f);
00131 std::complex <float> overallFactor = pow(tempF,m);
00132 std::complex <float> mI(0.0f,static_cast<float>(m));
00133 for (int ii=0; ii< RIntMax; ii++){ rhoOfRandmTemp[ii]=0;}
00134 for (int jx=0; jx <Center ; jx++) {
00135 for (int jy=0; jy <=jx; jy++){
00136 float fjx=float(jx);
00137 float fjy= float(jy);
00138 Count = (jx*jx+jx)/2 +1 +jy;
00139 PCount = PermMatTr[Count-1];
00140
00141 rhoTemp = std::complex <float> ((*ImBW)(CenterM+jx,CenterM+jy)) *exp(mI* std::complex <float> (atan2(+fjy,+fjx)))
00142 + std::complex <float> ((*ImBW)(CenterM+jx,CenterM-jy)) * exp(mI*std::complex <float>(atan2(-fjy,+fjx)))
00143 + std::complex <float> ((*ImBW)(CenterM-jx,CenterM+jy)) * exp(mI*std::complex <float>(atan2(+fjy,-fjx)))
00144 + std::complex <float> ((*ImBW)(CenterM-jx,CenterM-jy)) * exp(mI*std::complex <float>(atan2(-fjy,-fjx)))
00145 + std::complex <float> ((*ImBW)(CenterM+jy,CenterM+jx)) * exp(mI*std::complex <float>(atan2(+fjx,+fjy)))
00146 + std::complex <float> ((*ImBW)(CenterM+jy,CenterM-jx)) * exp(mI*std::complex <float>(atan2(-fjx,+fjy)))
00147 + std::complex <float> ((*ImBW)(CenterM-jy,CenterM+jx)) * exp(mI*std::complex <float>(atan2(+fjx,-fjy)))
00148 + std::complex <float> ((*ImBW)(CenterM-jy,CenterM-jx)) * exp(mI*std::complex <float>(atan2(-fjx,-fjy)));
00149 if (((jx+jy)==0)&&(m>0) ){
00150 rhoTemp=0;}
00151
00152
00153
00154
00155 rhoOfRandmTemp[PCount-1] +=
00156 rhoTemp/((float)pow(2.,(int)( (jx==0) +(jy==0)+ (jy==jx))));
00157
00158 }}
00159
00160
00161
00162
00163
00164
00165 float tempp;
00166
00167 for (int st=0; st<= Number2Use; st++){
00168 tempp=krVec2Use[st];
00169 sampledBesselJ[st] = static_cast<float>(gsl_sf_bessel_Jn(m,tempp));
00170
00171 }
00172
00173
00174 float *tempMB = new float [kIntMax*RIntMax];
00175 Util::spline_mat(krVec2Use, sampledBesselJ, Number2Use+1,krVec,tempMB,kIntMax*RIntMax );
00176
00177 std::complex <float> *rowV = new std::complex <float> [kIntMax];
00178
00179
00180
00181
00182
00183 for (int st=0; st < kIntMax; st++) {
00184 rowV[st]=0;
00185 for (int sv=0; sv < RIntMax; sv++) {
00186 rowV[st]+= rhoOfRandmTemp[sv] *tempMB[sv+st*RIntMax];
00187 }
00188 rowV[st] *= overallFactor;
00189
00190 }
00191 for (int st=0; st < kIntMax; st++) {
00192 (*rhoOfkmB)(2*m ,st) = rowV[st].real();
00193 (*rhoOfkmB)(2*m+1,st) = rowV[st].imag();
00194 }
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204 }
00205
00206 update();
00207 rhoOfkmB-> update();
00208 rhoOfkmB->set_complex(true);
00209 if(rhoOfkmB->get_ysize()==1 && rhoOfkmB->get_zsize()==1) {
00210 rhoOfkmB->set_complex_x(true);
00211 }
00212 rhoOfkmB->set_ri(true);
00213 rhoOfkmB->set_FH(true);
00214 rhoOfkmB->set_fftodd(true);
00215 return rhoOfkmB;
00216 } else {
00217 LOGERR("2D real square odd image expected.");
00218 throw ImageFormatException("2D real square odd image expected.");
00219 }
00220 }
00221
00222 EMData *EMData::copy_empty_head() const
00223 {
00224 ENTERFUNC;
00225 EMData *ret = new EMData();
00226 ret->attr_dict = attr_dict;
00227 ret->flags = flags;
00228 ret->all_translation = all_translation;
00229 ret->path = path;
00230 ret->pathnum = pathnum;
00231
00232
00233
00234
00235
00236
00237
00238 ret->update();
00239
00240 EXITFUNC;
00241 return ret;
00242 }
00243
00244
00245 EMData *EMData::FH2F(int Size, float OverSamplekB, int IntensityFlag)
00246 {
00247 int nx=get_xsize();
00248 int ny=get_ysize();
00249 int nz=get_zsize();
00250 float ScalFactor=4.1f;
00251 int Center = (int) floor((Size+1.0)/2.0 +.1);
00252 int CenterM= Center-1;
00253 int CountMax = (Center+1)*Center/2;
00254
00255 int *PermMatTr = new int[CountMax];
00256 float *RValsSorted = new float[CountMax];
00257 float *weightofkValsSorted = new float[CountMax];
00258 int *SizeReturned = new int[1];
00259 Util::Radialize(PermMatTr, RValsSorted,weightofkValsSorted,Size, SizeReturned);
00260 int RIntMax= SizeReturned[0];
00261
00262
00263 int mMax = (int) floor( ScalFactor*RValsSorted[RIntMax-1]+10.0);
00264
00265 int kIntMax = 2+ (int) floor( RValsSorted[RIntMax-1]*OverSamplekB);
00266 float *kVec2Use = new float[kIntMax];
00267 for (int kk=0; kk<kIntMax; kk++){
00268 kVec2Use[kk]= ((float) kk)/OverSamplekB;}
00269
00270
00271
00272 #ifdef DEBUG
00273 printf("nx=%d, ny=%d, nz=%d Center=%d mMax=%d CountMax=%d kIntMax=%d Centerm1=%d Size=%d\n\n",
00274 nx,ny,nz, Center, mMax, CountMax, kIntMax, CenterM, Size);
00275 #endif
00276
00277 EMData * rhoOfkmB = this;
00278
00279
00280
00281
00282 if ( (nx==2*(mMax+1)) && (ny==kIntMax) &&(nz==1) ) {
00283
00284 EMData *rhoOfkandm = copy();
00285 rhoOfkandm ->set_size(2*(mMax+1),RIntMax);
00286 rhoOfkandm ->to_zero();
00287
00288
00289 for (int mr=0; mr <2*(mMax+1); mr++){
00290 float *Row= new float[kIntMax];
00291 float *RowOut= new float[RIntMax];
00292 for (int ii=0; ii<kIntMax; ii++){ Row[ii]=(*rhoOfkmB)(mr,ii);}
00293 Util::spline_mat(kVec2Use, Row, kIntMax, RValsSorted, RowOut, RIntMax );
00294 for (int ii=0; ii<RIntMax; ii++){
00295 (*rhoOfkandm)(mr,ii) = RowOut[ii];
00296
00297 }
00298
00299
00300 }
00301 rhoOfkandm ->update();
00302
00303
00304
00305 EMData* outCopy = rhoOfkandm ->copy();
00306 outCopy->set_size(2*Size,Size,1);
00307 outCopy->to_zero();
00308
00309
00310 int Count =0, kInt, kIntm1;
00311 std::complex <float> ImfTemp;
00312 float kValue, thetak;
00313
00314 for (int jkx=0; jkx <Center; jkx++) {
00315 for (int jky=0; jky<=jkx; jky++){
00316 kInt = PermMatTr[Count];
00317 kIntm1= kInt-1;
00318 Count++;
00319 float fjkx = float(jkx);
00320 float fjky = float(jky);
00321
00322 kValue = std::sqrt(fjkx*fjkx + fjky*fjky ) ;
00323
00324
00325
00326
00327 thetak = atan2(fjky,fjkx);
00328 ImfTemp = (*rhoOfkandm)(0, kIntm1) ;
00329 for (int mm= 1; mm <mMax;mm++) {
00330 std::complex <float> fact(0,-mm*thetak);
00331 std::complex <float> expfact= exp(fact);
00332 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1),(*rhoOfkandm)(2*mm+1,kIntm1));
00333 float mmFac = float(1-2*(mm%2));
00334 if (IntensityFlag==1){ mmFac=1;}
00335 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00336 }
00337 (*outCopy)(2*(CenterM+jkx),CenterM+jky) = ImfTemp.real();
00338 (*outCopy)(2*(CenterM+jkx)+1,CenterM+jky) = ImfTemp.imag();
00339
00340
00341 if (jky>0) {
00342 thetak = atan2(-fjky,fjkx);
00343 ImfTemp = (*rhoOfkandm)(0,kIntm1);
00344 for (int mm= 1; mm<mMax; mm++) {
00345 std::complex <float> fact(0,-mm*thetak);
00346 std::complex <float> expfact= exp(fact);
00347 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1), (*rhoOfkandm)(2*mm+1,kIntm1));
00348 float mmFac = float(1-2*(mm%2));
00349 if (IntensityFlag==1){ mmFac=1;}
00350 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00351 }
00352 (*outCopy)(2*(CenterM+jkx),CenterM-jky) = ImfTemp.real();
00353
00354 (*outCopy)(2*(CenterM+jkx)+1,CenterM-jky) = ImfTemp.imag();
00355 }
00356
00357 if (jkx>0) {
00358 thetak = atan2(fjky,-fjkx);
00359 ImfTemp = (*rhoOfkandm)(0,kIntm1);
00360 for (int mm= 1; mm<mMax; mm++) {
00361 std::complex <float> fact(0,-mm*thetak);
00362 std::complex <float> expfact= exp(fact);
00363 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1), (*rhoOfkandm)(2*mm+1,kIntm1));
00364 float mmFac = float(1-2*(mm%2));
00365 if (IntensityFlag==1){ mmFac=1;}
00366 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00367 }
00368 (*outCopy)(2*(CenterM-jkx) ,CenterM+jky) = ImfTemp.real();
00369 (*outCopy)(2*(CenterM-jkx)+1,CenterM+jky) = ImfTemp.imag();
00370 }
00371
00372 if (jkx>0 && jky>0) {
00373 thetak = atan2(-fjky,-fjkx);
00374 ImfTemp = (*rhoOfkandm)(0 , kIntm1);
00375 for (int mm= 1; mm<mMax; mm++) {
00376 std::complex <float> fact(0,-mm*thetak);
00377 std::complex <float> expfact= exp(fact);
00378 std::complex <float> tempRho( (*rhoOfkandm)(2*mm,kIntm1),(*rhoOfkandm)(2*mm+1,kIntm1) );
00379 float mmFac = float(1-2*(mm%2));
00380 if (IntensityFlag==1){ mmFac=1;}
00381 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00382 }
00383 (*outCopy)(2*(CenterM-jkx) ,CenterM-jky) = ImfTemp.real();
00384 (*outCopy)(2*(CenterM-jkx)+1,CenterM-jky) = ImfTemp.imag();
00385 }
00386
00387 if (jky< jkx) {
00388 thetak = atan2(fjkx,fjky);
00389 ImfTemp = (*rhoOfkandm)(0,kIntm1);
00390 for (int mm= 1; mm<mMax; mm++){
00391 std::complex <float> fact(0,-mm*thetak);
00392 std::complex <float> expfact= exp(fact);
00393 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1),(*rhoOfkandm)(2*mm+1,kIntm1));
00394 float mmFac = float(1-2*(mm%2));
00395 if (IntensityFlag==1){ mmFac=1;}
00396 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00397 }
00398 (*outCopy)(2*(CenterM+jky) ,CenterM+jkx) = ImfTemp.real();
00399 (*outCopy)(2*(CenterM+jky)+1,CenterM+jkx) = ImfTemp.imag();
00400
00401 if (jky>0){
00402 thetak = atan2(fjkx,-fjky);
00403 ImfTemp = (*rhoOfkandm)(0, kIntm1);
00404 for (int mm= 1; mm <mMax; mm++) {
00405 std::complex <float> fact(0,-mm*thetak);
00406 std::complex <float> expfact= exp(fact);
00407 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1),(*rhoOfkandm)(2*mm+1,kIntm1));
00408 float mmFac = float(1-2*(mm%2));
00409 if (IntensityFlag==1){ mmFac=1;}
00410 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00411 }
00412 (*outCopy)(2*(CenterM-jky) ,CenterM+jkx) = ImfTemp.real();
00413 (*outCopy)(2*(CenterM-jky)+1,CenterM+jkx) = ImfTemp.imag();
00414 }
00415
00416 if (jkx>0) {
00417 thetak = atan2(-fjkx,fjky);
00418 ImfTemp = (*rhoOfkandm)(0,kIntm1);
00419 for (int mm= 1; mm <mMax; mm++) {
00420 std::complex <float> fact(0,-mm*thetak);
00421 std::complex <float> expfact= exp(fact);
00422 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1),(*rhoOfkandm)(2*mm+1,kIntm1));
00423 float mmFac = float(1-2*(mm%2));
00424 if (IntensityFlag==1){ mmFac=1;}
00425 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00426 }
00427 (*outCopy)(2*(CenterM+jky) ,CenterM-jkx) = ImfTemp.real();
00428 (*outCopy)(2*(CenterM+jky)+1,CenterM-jkx) = ImfTemp.imag();
00429 }
00430
00431 if (jkx>0 && jky>0) {
00432 thetak = atan2(-fjkx,-fjky);
00433 ImfTemp = (*rhoOfkandm)(0,kIntm1) ;
00434 for (int mm= 1; mm <mMax; mm++) {
00435 std::complex <float> fact(0,-mm*thetak);
00436 std::complex <float> expfact= exp(fact);
00437 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1) ,(*rhoOfkandm)(2*mm+1,kIntm1) );
00438 float mmFac = float(1-2*(mm%2));
00439 if (IntensityFlag==1){ mmFac=1;}
00440 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00441 }
00442 (*outCopy)(2*(CenterM-jky) ,CenterM-jkx) = ImfTemp.real();
00443 (*outCopy)(2*(CenterM-jky)+1,CenterM-jkx) = ImfTemp.imag();
00444 }
00445 }
00446
00447
00448 }
00449 }
00450 outCopy->update();
00451 outCopy->set_complex(true);
00452 if(outCopy->get_ysize()==1 && outCopy->get_zsize()==1) outCopy->set_complex_x(true);
00453 outCopy->set_ri(true);
00454 outCopy->set_FH(false);
00455 outCopy->set_fftodd(true);
00456 outCopy->set_shuffled(true);
00457 return outCopy;
00458 } else {
00459 LOGERR("can't be an FH image not this size");
00460 throw ImageFormatException("something strange about this image: not a FH");
00461
00462 }
00463 }
00464
00465
00466 EMData *EMData::FH2Real(int Size, float OverSamplekB, int)
00467 {
00468 EMData* FFT= FH2F(Size,OverSamplekB,0);
00469 FFT->process_inplace("xform.fourierorigin.tocorner");
00470 EMData* eguess= FFT ->do_ift();
00471 return eguess;
00472 }
00473
00474 float dist(int lnlen, const float* line_1, const float* line_2)
00475 {
00476 float dis2=0.0;
00477 for( int i=0; i < lnlen; ++i) {
00478 float tmp = line_1[i] - line_2[i];
00479 dis2 += tmp*tmp;
00480 }
00481
00482 return dis2;
00483 }
00484
00485 float dist_r(int lnlen, const float* line_1, const float* line_2)
00486 {
00487 double dis2 = 0.0;
00488 for( int i=0; i < lnlen; ++i ) {
00489 float tmp = line_1[lnlen-1-i] - line_2[i];
00490 dis2 += tmp*tmp;
00491 }
00492 return static_cast<float>(std::sqrt(dis2));
00493 }
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525 float EMData::cm_euc(EMData* sinoj, int n1, int n2)
00526 {
00527 int lnlen = get_xsize();
00528 float* line_1 = get_data() + n1 * lnlen;
00529 float* line_2 = sinoj->get_data() + n2 * lnlen;
00530 return dist(lnlen, line_1, line_2);
00531 }
00532
00533 EMData* EMData::rotavg() {
00534
00535 int rmax;
00536
00537 ENTERFUNC;
00538
00539
00540 if (ny<2 && nz <2) {
00541 LOGERR("No 1D images.");
00542 throw ImageDimensionException("No 1D images!");
00543 }
00544
00545 float apix[3];
00546 apix[0] = get_attr_default("apix_x",1.0);
00547 apix[1] = get_attr_default("apix_y",1.0);
00548 apix[2] = get_attr_default("apix_z",1.0);
00549 float min_apix = *std::min_element(&apix[0],&apix[3]);
00550
00551 float apix_x = apix[0]/min_apix;
00552 float apix_y = apix[1]/min_apix;
00553 float apix_z = 1.0;
00554 if( nz > 1)
00555 apix_z=apix[2]/min_apix;
00556 float apix_x2 = apix_x*apix_x;
00557 float apix_y2 = apix_y*apix_y;
00558 float apix_z2 = apix_z*apix_z;
00559
00560 vector<int> saved_offsets = get_array_offsets();
00561 set_array_offsets(-nx/2,-ny/2,-nz/2);
00562
00563
00564 #ifdef _WIN32
00565
00566 if ( nz == 1 ) {
00567 rmax = _cpp_min( nx/2 + nx%2, ny/2 + ny%2);
00568 } else {
00569 rmax = _cpp_min(nx/2 + nx%2, _cpp_min(ny/2 + ny%2, nz/2 + nz%2));
00570 }
00571 #else
00572
00573 if ( nz == 1 ) {
00574 rmax = std::min(nx/2 + nx%2, ny/2 + ny%2);
00575 } else {
00576 rmax = std::min(nx/2 + nx%2, std::min(ny/2 + ny%2, nz/2 + nz%2));
00577 }
00578 #endif //_WIN32
00579
00580 float rmax_ratio = 0.0;
00581 if (rmax == nx/2 + nx%2 )
00582 rmax_ratio = apix_x;
00583 else if (rmax == ny/2 + ny%2)
00584 rmax_ratio = apix_y;
00585 else
00586 rmax_ratio = apix_z;
00587
00588 EMData* ret = new EMData();
00589 ret->set_size(rmax+1, 1, 1);
00590 ret->to_zero();
00591 vector<float> count(rmax+1);
00592 for (int k = -nz/2; k < nz/2 + nz%2; k++) {
00593 if (abs( k*apix_z) > rmax*rmax_ratio ) continue;
00594 for (int j = -ny/2; j < ny/2 + ny%2; j++) {
00595 if (abs( j*apix_y ) > rmax*rmax_ratio) continue;
00596 for (int i = -nx/2; i < nx/2 + nx%2; i++) {
00597 float r = std::sqrt(float(k*k*apix_z2) + float(j*j*apix_y2) + float(i*i*apix_x2))/rmax_ratio;
00598 int ir = int(r);
00599 if (ir >= rmax) continue;
00600 float frac = r - float(ir);
00601 (*ret)(ir) += (*this)(i,j,k)*(1.0f - frac);
00602 (*ret)(ir+1) += (*this)(i,j,k)*frac;
00603 count[ir] += 1.0f - frac;
00604 count[ir+1] += frac;
00605 }
00606 }
00607 }
00608 for (int ir = 0; ir <= rmax; ir++) {
00609 #ifdef _WIN32
00610 (*ret)(ir) /= _cpp_max(count[ir],1.0f);
00611 #else
00612 (*ret)(ir) /= std::max(count[ir],1.0f);
00613 #endif //_WIN32
00614 }
00615
00616 set_array_offsets(saved_offsets);
00617 ret->update();
00618 EXITFUNC;
00619 return ret;
00620 }
00621
00622 EMData* EMData::rotavg_i() {
00623
00624 int rmax;
00625 ENTERFUNC;
00626 if ( ny == 1 && nz == 1 ) {
00627 LOGERR("Input image must be 2-D or 3-D!");
00628 throw ImageDimensionException("Input image must be 2-D or 3-D!");
00629 }
00630
00631 EMData* avg1D = new EMData();
00632 EMData* result = new EMData();
00633
00634 result->set_size(nx,ny,nz);
00635 result->to_zero();
00636 result->set_array_offsets(-nx/2, -ny/2, -nz/2);
00637
00638 if ( nz == 1 ) {
00639 #ifdef _WIN32
00640 rmax = _cpp_min(nx/2 + nx%2, ny/2 + ny%2);
00641 } else {
00642 rmax = _cpp_min(nx/2 + nx%2, _cpp_min(ny/2 + ny%2, nz/2 + nz%2));
00643 #else
00644 rmax = std::min(nx/2 + nx%2, ny/2 + ny%2);
00645 } else {
00646 rmax = std::min(nx/2 + nx%2, std::min(ny/2 + ny%2, nz/2 + nz%2));
00647 #endif //_WIN32
00648 }
00649
00650 avg1D = rotavg();
00651 float padded_value = 0.0, r;
00652 int i, j, k, ir;
00653 size_t number_of_pixels = 0;
00654 for ( k = -nz/2; k < nz/2 + nz%2; k++) {
00655 if (abs(k) > rmax) continue;
00656 for ( j = -ny/2; j < ny/2 + ny%2; j++) {
00657 if (abs(j) > rmax) continue;
00658 for (i = -nx/2; i < nx/2 + nx%2; i++) {
00659 r = std::sqrt(float(k*k) + float(j*j) + float(i*i));
00660 ir = int(r);
00661 if (ir > rmax || ir < rmax-2 ) continue ;
00662 else {
00663 padded_value += (*avg1D)(ir) ;
00664 number_of_pixels++ ;
00665 }
00666 }
00667 }
00668 }
00669 padded_value /= number_of_pixels;
00670 for ( k = -nz/2; k < nz/2 + nz%2; k++) {
00671 for ( j = -ny/2; j < ny/2 + ny%2; j++) {
00672 for ( i = -nx/2; i < nx/2 + nx%2; i++) {
00673 r = std::sqrt(float(k*k) + float(j*j) + float(i*i));
00674 ir = int(r);
00675 if (ir >= rmax) (*result)(i,j,k) = padded_value ;
00676 else (*result)(i,j,k) = (*avg1D)(ir)+((*avg1D)(ir+1)-(*avg1D)(ir))*(r - float(ir));
00677
00678 }
00679 }
00680 }
00681 result->update();
00682 result->set_array_offsets(0,0,0);
00683 EXITFUNC;
00684 return result;
00685 }
00686
00687
00688 EMData* EMData::mult_radial(EMData* radial) {
00689
00690 ENTERFUNC;
00691 if ( ny == 1 && nz == 1 ) {
00692 LOGERR("Input image must be 2-D or 3-D!");
00693 throw ImageDimensionException("Input image must be 2-D or 3-D!");
00694 }
00695
00696 EMData* result = this->copy_head();
00697
00698 result->to_zero();
00699 result->set_array_offsets(-nx/2, -ny/2, -nz/2);
00700 this->set_array_offsets(-nx/2, -ny/2, -nz/2);
00701 int rmax = radial->get_xsize();
00702 int i, j, k, ir;
00703 float r;
00704 for ( k = -nz/2; k < nz/2+nz%2; k++) {
00705 for ( j = -ny/2; j < ny/2+ny%2; j++) {
00706 for ( i = -nx/2; i < nx/2+nx%2; i++) {
00707 r = std::sqrt(float(k*k) + float(j*j) + float(i*i));
00708 ir = int(r);
00709 if(ir < rmax-1) (*result)(i,j,k) = (*this)(i,j,k) * ((*radial)(ir)+((*radial)(ir+1)-(*radial)(ir))*(r - float(ir)));
00710 }
00711 }
00712 }
00713 result->update();
00714 result->set_array_offsets(0,0,0);
00715 this->set_array_offsets(0,0,0);
00716 EXITFUNC;
00717 return result;
00718 }
00719
00720 #define rdata(i,j,k) rdata[(i-1)+((j-1)+(k-1)*ny)*(size_t)nx]
00721 #define square(x) ((x)*(x))
00722 vector<float> EMData::cog() {
00723
00724 vector<float> cntog;
00725 int ndim = get_ndim();
00726 int i=1,j=1,k=1;
00727 float val,sum1=0.f,MX=0.f,RG=0.f,MY=0.f,MZ=0.f,r=0.f;
00728
00729 if (ndim == 1) {
00730 for ( i = 1;i <= nx; i++) {
00731 val = rdata(i,j,k);
00732 sum1 += val;
00733 MX += ((i-1)*val);
00734 }
00735 MX=(MX/sum1);
00736 for ( i = 1;i <= nx; i++) {
00737 val = rdata(i,j,k);
00738 sum1 += val;
00739 RG += val*(square(MX - (i-1)));
00740 }
00741 RG=std::sqrt(RG/sum1);
00742 MX=MX-(nx/2);
00743 cntog.push_back(MX);
00744 cntog.push_back(RG);
00745 #ifdef _WIN32
00746 cntog.push_back((float)Util::round(MX));
00747 #else
00748 cntog.push_back(round(MX));
00749 #endif //_WIN32
00750 } else if (ndim == 2) {
00751 for (j=1;j<=ny;j++) {
00752 for (i=1;i<=nx;i++) {
00753 val = rdata(i,j,k);
00754 sum1 += val;
00755 MX += ((i-1)*val);
00756 MY += ((j-1)*val);
00757 }
00758 }
00759 MX=(MX/sum1);
00760 MY=(MY/sum1);
00761 sum1=0.f;
00762 RG=0.f;
00763 for (j=1;j<=ny;j++) {
00764 r = (square(MY-(j-1)));
00765 for (i=1;i<=nx;i++) {
00766 val = rdata(i,j,k);
00767 sum1 += val;
00768 RG += val*(square(MX - (i-1)) + r);
00769 }
00770 }
00771 RG = std::sqrt(RG/sum1);
00772 MX = MX - nx/2;
00773 MY = MY - ny/2;
00774 cntog.push_back(MX);
00775 cntog.push_back(MY);
00776 cntog.push_back(RG);
00777 #ifdef _WIN32
00778 cntog.push_back((float)Util::round(MX));cntog.push_back((float)Util::round(MY));
00779 #else
00780 cntog.push_back(round(MX));cntog.push_back(round(MY));
00781 #endif //_WIN32
00782 } else {
00783 for (k = 1;k <= nz;k++) {
00784 for (j=1;j<=ny;j++) {
00785 for (i=1;i<=nx;i++) {
00786 val = rdata(i,j,k);
00787 sum1 += val;
00788 MX += ((i-1)*val);
00789 MY += ((j-1)*val);
00790 MZ += ((k-1)*val);
00791 }
00792 }
00793 }
00794 MX = MX/sum1;
00795 MY = MY/sum1;
00796 MZ = MZ/sum1;
00797 sum1=0.f;
00798 RG=0.f;
00799 for (k = 1;k <= nz;k++) {
00800 for (j=1;j<=ny;j++) {
00801 float r = (square(MY-(j-1)) + square(MZ - (k-1)));
00802 for (i=1;i<=nx;i++) {
00803 val = rdata(i,j,k);
00804 sum1 += val;
00805 RG += val*(square(MX - (i-1)) + r);
00806 }
00807 }
00808 }
00809 RG = std::sqrt(RG/sum1);
00810 MX = MX - nx/2;
00811 MY = MY - ny/2;
00812 MZ = MZ - nz/2;
00813 cntog.push_back(MX);
00814 cntog.push_back(MY);
00815 cntog.push_back(MZ);
00816 cntog.push_back(RG);
00817 #ifdef _WIN32
00818 cntog.push_back((float)Util::round(MX));cntog.push_back((float)Util::round(MY));cntog.push_back((float)Util::round(MZ));
00819 #else
00820 cntog.push_back(round(MX));cntog.push_back(round(MY));cntog.push_back(round(MZ));
00821 #endif //_WIN32
00822 }
00823 return cntog;
00824 }
00825 #undef square
00826 #undef rdata
00827
00828
00829
00830
00831
00832 vector < float >EMData::calc_fourier_shell_correlation(EMData * with, float w)
00833 {
00834 ENTERFUNC;
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859 int needfree=0, kz, ky, ii;
00860 float argx, argy, argz;
00861
00862 if (!with) {
00863 throw NullPointerException("NULL input image");
00864 }
00865
00866
00867 EMData *f = this;
00868 EMData *g = with;
00869
00870 int nx = f->get_xsize();
00871 int ny = f->get_ysize();
00872 int nz = f->get_zsize();
00873
00874 if (ny==0 && nz==0) {
00875 throw ImageFormatException( "Cannot calculate FSC for 1D images");
00876 }
00877
00878 if (f->is_complex()) nx = (nx - 2 + f->is_fftodd());
00879 int lsd2 = (nx + 2 - nx%2) ;
00880
00881
00882 EMData* fpimage = NULL;
00883 if (f->is_complex()) fpimage = f;
00884 else {
00885 fpimage= f->norm_pad(false, 1);
00886 fpimage->do_fft_inplace();
00887 needfree|=1;
00888 }
00889
00890
00891 EMData* gpimage = NULL;
00892 if (g->is_complex()) gpimage = g;
00893 else {
00894 gpimage= g->norm_pad(false, 1);
00895 gpimage->do_fft_inplace();
00896 needfree|=2;
00897 }
00898
00899 float *d1 = fpimage->get_data();
00900 float *d2 = gpimage->get_data();
00901
00902 int nx2 = nx/2;
00903 int ny2 = ny/2;
00904 int nz2 = nz/2;
00905
00906 float dx2 = 1.0f/float(nx2)/float(nx2);
00907 float dy2 = 1.0f/float(ny2)/float(ny2);
00908
00909 #ifdef _WIN32
00910 float dz2 = 1.0f / _cpp_max(float(nz2),1.0f)/_cpp_max(float(nz2),1.0f);
00911 int inc = Util::round(float( _cpp_max( _cpp_max(nx2,ny2),nz2) )/w );
00912 #else
00913 float dz2 = 1.0f/std::max(float(nz2),1.0f)/std::max(float(nz2),1.0f);
00914 int inc = Util::round(float(std::max(std::max(nx2,ny2),nz2))/w);
00915 #endif //_WIN32
00916
00917 double* ret = new double[inc+1];
00918 double* n1 = new double[inc+1];
00919 double* n2 = new double[inc+1];
00920 float* lr = new float[inc+1];
00921 for (int i = 0; i <= inc; i++) {
00922 ret[i] = 0; n1[i] = 0; n2[i] = 0; lr[i]=0;
00923 }
00924
00925 for (int iz = 0; iz <= nz-1; iz++) {
00926 if(iz>nz2) kz=iz-nz; else kz=iz; argz = float(kz*kz)*dz2;
00927 for (int iy = 0; iy <= ny-1; iy++) {
00928 if(iy>ny2) ky=iy-ny; else ky=iy; argy = argz + float(ky*ky)*dy2;
00929 for (int ix = 0; ix <= lsd2-1; ix+=2) {
00930
00931 if (ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00932 argx = 0.5f*std::sqrt(argy + float(ix*ix)*0.25f*dx2);
00933 int r = Util::round(inc*2*argx);
00934 if(r <= inc) {
00935 ii = ix + (iy + iz * ny)* lsd2;
00936 ret[r] += d1[ii] * double(d2[ii]) + d1[ii + 1] * double(d2[ii + 1]);
00937 n1[r] += d1[ii] * double(d1[ii]) + d1[ii + 1] * double(d1[ii + 1]);
00938 n2[r] += d2[ii] * double(d2[ii]) + d2[ii + 1] * double(d2[ii + 1]);
00939 lr[r] += 2;
00940 }
00941 }
00942 }
00943 }
00944 }
00945
00946 int linc = 0;
00947 for (int i = 0; i <= inc; i++) if(lr[i]>0) linc++;
00948
00949 vector<float> result(linc*3);
00950
00951 ii = -1;
00952 for (int i = 0; i <= inc; i++) {
00953 if(lr[i]>0) {
00954 ii++;
00955 result[ii] = float(i)/float(2*inc);
00956 result[ii+linc] = float(ret[i] / (std::sqrt(n1[i] * n2[i])));
00957 result[ii+2*linc] = lr[i] ;
00958 }
00959
00960
00961
00962
00963 }
00964
00965 if (needfree&1) {
00966 if (fpimage) {
00967 delete fpimage;
00968 fpimage = 0;
00969 }
00970 }
00971 if (needfree&2) {
00972 if (gpimage) {
00973 delete gpimage;
00974 gpimage = 0;
00975 }
00976 }
00977 delete[] ret; delete[] n1; delete[] n2; delete[] lr;
00978
00979 EXITFUNC;
00980 return result;
00981 }
00982
00983
00984 EMData* EMData::symvol(string symString) {
00985 ENTERFUNC;
00986 int nsym = Transform::get_nsym(symString);
00987 Transform sym;
00988
00989 EMData *svol = new EMData;
00990 svol->set_size(nx, ny, nz);
00991 svol->to_zero();
00992
00993
00994
00995
00996
00997 for (int isym = 0; isym < nsym; isym++) {
00998 Transform rm = sym.get_sym(symString, isym);
00999 EMData* symcopy = this -> rot_scale_trans(rm);
01000 *svol += (*symcopy);
01001 delete symcopy;
01002
01003 }
01004
01005 *svol /= ((float) nsym);
01006 svol->update();
01007 EXITFUNC;
01008 return svol;
01009 }
01010
01011 #define proj(ix,iy,iz) proj[ix-1+(iy-1+(iz-1)*ny)*(size_t)nx]
01012 #define pnewimg(ix,iy,iz) pnewimg[ix-1+(iy-1+(iz-1)*ny)*(size_t)nx]
01013 EMData* EMData::average_circ_sub() const
01014 {
01015
01016
01017 ENTERFUNC;
01018 const EMData* const image = this;
01019 EMData* newimg = copy_head();
01020 float *proj = image->get_data();
01021 float *pnewimg = newimg->get_data();
01022
01023 float r2 = static_cast<float>( (nx/2)*(nx/2) );
01024 float qs=0.0f;
01025 int m=0;
01026 int ncz = nz/2 + 1;
01027 int ncy = ny/2 + 1;
01028 int ncx = nx/2 + 1;
01029 for (int iz = 1; iz <= nz; iz++) {
01030 float yy = static_cast<float>( (iz-ncz)*(iz-ncz) );
01031 for (int iy = 1; iy <=ny; iy++) { float xx = yy + (iy-ncy)*(iy-ncy);
01032 for (int ix = 1; ix <= nx; ix++) {
01033 if ( xx+float((ix-ncx)*(ix-ncx)) > r2 ) {
01034 qs += proj(ix,iy,iz);
01035 m++;
01036 }
01037 }
01038 }
01039 }
01040
01041
01042 if( m > 0 ) qs /= m;
01043
01044 for (int iz = 1; iz <= nz; iz++)
01045 for (int iy = 1; iy <= ny; iy++)
01046 for (int ix = 1; ix <= nx; ix++)
01047 pnewimg(ix,iy,iz) = proj(ix,iy,iz) - qs;
01048 newimg->update();
01049 return newimg;
01050 EXITFUNC;
01051 }
01052
01053
01054
01055
01056
01057 void EMData::onelinenn(int j, int n, int n2, EMData* wptr, EMData* bi, const Transform& tf)
01058 {
01059
01060 int jp = (j >= 0) ? j+1 : n+j+1;
01061
01062
01063 for (int i = 0; i <= n2; i++) {
01064 if (((i*i+j*j) < n*n/4) && !((0 == i) && (j < 0))) {
01065
01066 float xnew = i*tf[0][0] + j*tf[1][0];
01067 float ynew = i*tf[0][1] + j*tf[1][1];
01068 float znew = i*tf[0][2] + j*tf[1][2];
01069 std::complex<float> btq;
01070 if (xnew < 0.) {
01071 xnew = -xnew;
01072 ynew = -ynew;
01073 znew = -znew;
01074 btq = conj(bi->cmplx(i,jp));
01075 } else {
01076 btq = bi->cmplx(i,jp);
01077 }
01078 int ixn = int(xnew + 0.5 + n) - n;
01079 int iyn = int(ynew + 0.5 + n) - n;
01080 int izn = int(znew + 0.5 + n) - n;
01081
01082 int iza, iya;
01083 if (izn >= 0) iza = izn + 1;
01084 else iza = n + izn + 1;
01085
01086 if (iyn >= 0) iya = iyn + 1;
01087 else iya = n + iyn + 1;
01088
01089 cmplx(ixn,iya,iza) += btq;
01090
01091 (*wptr)(ixn,iya,iza)++;
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118 }
01119 }
01120 }
01121
01122
01123 void EMData::onelinenn_mult(int j, int n, int n2, EMData* wptr, EMData* bi, const Transform& tf, int mult)
01124 {
01125
01126 int jp = (j >= 0) ? j+1 : n+j+1;
01127
01128
01129 for (int i = 0; i <= n2; i++) {
01130 if (((i*i+j*j) < n*n/4) && !((0 == i) && (j < 0))) {
01131
01132 float xnew = i*tf[0][0] + j*tf[1][0];
01133 float ynew = i*tf[0][1] + j*tf[1][1];
01134 float znew = i*tf[0][2] + j*tf[1][2];
01135 std::complex<float> btq;
01136 if (xnew < 0.) {
01137 xnew = -xnew;
01138 ynew = -ynew;
01139 znew = -znew;
01140 btq = conj(bi->cmplx(i,jp));
01141 } else {
01142 btq = bi->cmplx(i,jp);
01143 }
01144 int ixn = int(xnew + 0.5 + n) - n;
01145 int iyn = int(ynew + 0.5 + n) - n;
01146 int izn = int(znew + 0.5 + n) - n;
01147
01148
01149 int iza, iya;
01150 if (izn >= 0) iza = izn + 1;
01151 else iza = n + izn + 1;
01152
01153 if (iyn >= 0) iya = iyn + 1;
01154 else iya = n + iyn + 1;
01155
01156 cmplx(ixn,iya,iza) += btq*float(mult);
01157 (*wptr)(ixn,iya,iza)+=float(mult);
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182 }
01183 }
01184 }
01185
01186 void EMData::nn(EMData* wptr, EMData* myfft, const Transform& tf, int mult)
01187 {
01188 ENTERFUNC;
01189 int nxc = attr_dict["nxc"];
01190
01191 vector<int> saved_offsets = get_array_offsets();
01192 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01193 set_array_offsets(0,1,1);
01194 myfft->set_array_offsets(0,1);
01195
01196 if( mult == 1 ) {
01197 for (int iy = -ny/2 + 1; iy <= ny/2; iy++) onelinenn(iy, ny, nxc, wptr, myfft, tf);
01198 } else {
01199 for (int iy = -ny/2 + 1; iy <= ny/2; iy++) onelinenn_mult(iy, ny, nxc, wptr, myfft, tf, mult);
01200 }
01201
01202 set_array_offsets(saved_offsets);
01203 myfft->set_array_offsets(myfft_saved_offsets);
01204 EXITFUNC;
01205 }
01206
01207
01208 void EMData::insert_rect_slice(EMData* w, EMData* myfft, const Transform& trans, int sizeofprojection, float xratio, float yratio, int npad, int mult)
01209 {
01210 ENTERFUNC;
01211 vector<int> saved_offsets = get_array_offsets();
01212 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01213 set_array_offsets(0,1,1);
01214 myfft->set_array_offsets(0,1);
01215
01216
01217
01218 Vec2f coordinate_2d_square;
01219 Vec3f coordinate_3dnew;
01220 Vec3f axis_newx;
01221 Vec3f axis_newy;
01222 Vec3f tempv;
01223
01224
01225
01226 axis_newx[0] = xratio*0.5f*(sizeofprojection*npad)*trans[0][0];
01227 axis_newx[1] = yratio*0.5f*(sizeofprojection*npad)*trans[0][1];
01228 axis_newx[2] = 0.5f*(sizeofprojection*npad)*trans[0][2];
01229
01230 float ellipse_length_x = std::sqrt(axis_newx[0]*axis_newx[0]+axis_newx[1]*axis_newx[1]+axis_newx[2]*axis_newx[2]);
01231
01232 int ellipse_length_x_int = int(ellipse_length_x);
01233 float ellipse_step_x = 0.5f*(sizeofprojection*npad)/float(ellipse_length_x_int);
01234 float xscale = ellipse_step_x;
01235
01236 axis_newy[0] = xratio*0.5f*(sizeofprojection*npad)*trans[1][0];
01237 axis_newy[1] = yratio*0.5f*(sizeofprojection*npad)*trans[1][1];
01238 axis_newy[2] = 0.5f*(sizeofprojection*npad)*trans[1][2];
01239
01240
01241
01242 float ellipse_length_y = std::sqrt(axis_newy[0]*axis_newy[0]+axis_newy[1]*axis_newy[1]+axis_newy[2]*axis_newy[2]);
01243 int ellipse_length_y_int = int(ellipse_length_y);
01244 float ellipse_step_y = 0.5f*(sizeofprojection*npad)/float(ellipse_length_y_int);
01245 float yscale = ellipse_step_y;
01246
01247 std::complex<float> c1;
01248 nz = get_zsize();
01249
01250 float r2=0.25f*sizeofprojection*npad*sizeofprojection*npad;
01251 float r2_at_point;
01252
01253 for(int i=0;i<ellipse_length_x_int;i++) {
01254 for(int j=-1*ellipse_length_y_int+1; j<=ellipse_length_y_int; j++) {
01255
01256 r2_at_point=i*xscale*i*xscale+j*yscale*j*yscale;
01257 if(r2_at_point<=r2 ) {
01258
01259
01260 coordinate_2d_square[0] = xscale*float(i);
01261 coordinate_2d_square[1] = yscale*float(j);
01262 float xnew = coordinate_2d_square[0]*trans[0][0] + coordinate_2d_square[1]*trans[1][0];
01263 float ynew = coordinate_2d_square[0]*trans[0][1] + coordinate_2d_square[1]*trans[1][1];
01264 float znew = coordinate_2d_square[0]*trans[0][2] + coordinate_2d_square[1]*trans[1][2];
01265 coordinate_3dnew[0] =xnew*xratio;
01266 coordinate_3dnew[1] = ynew*yratio;
01267 coordinate_3dnew[2] = znew;
01268
01269
01270 float xp = coordinate_2d_square[0];
01271 float yp = ( coordinate_2d_square[1] >= 0) ? coordinate_2d_square[1]+1 : nz+coordinate_2d_square[1]+1;
01272 std::complex<float> lin_interpolated(0,0);
01273 int xlow=int(xp),xhigh=int(xp)+1;
01274 int ylow=int(yp),yhigh=int(yp)+1;
01275 float tx=xp-xlow,ty=yp-ylow;
01276
01277
01278 if(j == -1) {
01279
01280 if(ylow<yp)
01281 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)*(1-ty) + myfft->cmplx(xlow,yhigh)*(1-tx)*ty
01282 + myfft->cmplx(xhigh,ylow)*tx*(1-ty) + myfft->cmplx(xhigh,yhigh)*tx*ty;
01283 else
01284 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)
01285 + myfft->cmplx(xhigh,ylow)*tx;
01286
01287 }
01288 else {
01289 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)*(1-ty) + myfft->cmplx(xlow,yhigh)*(1-tx)*ty
01290 + myfft->cmplx(xhigh,ylow)*tx*(1-ty)+ myfft->cmplx(xhigh,yhigh)*tx*ty;
01291
01292 }
01293
01294 c1 = lin_interpolated;
01295
01296
01297
01298 std::complex<float> btq;
01299 if ( coordinate_3dnew[0] < 0.) {
01300 coordinate_3dnew[0] = -coordinate_3dnew[0];
01301 coordinate_3dnew[1] = -coordinate_3dnew[1];
01302 coordinate_3dnew[2] = -coordinate_3dnew[2];
01303 btq = conj(c1);
01304 } else {
01305 btq = c1;
01306 }
01307 int ixn = int(coordinate_3dnew[0] + 0.5 + nx) - nx;
01308 int iyn = int(coordinate_3dnew[1] + 0.5 + ny) - ny;
01309 int izn = int(coordinate_3dnew[2] + 0.5 + nz) - nz;
01310
01311 int iza, iya;
01312 if (izn >= 0) iza = izn + 1;
01313 else iza = nz + izn + 1;
01314
01315 if (iyn >= 0) iya = iyn + 1;
01316 else iya = ny + iyn + 1;
01317
01318 cmplx(ixn,iya,iza) += btq*float(mult);
01319 (*w)(ixn,iya,iza) += mult;
01320
01321 }
01322 }
01323
01324 }
01325
01326
01327
01328
01329 set_array_offsets(saved_offsets);
01330 myfft->set_array_offsets(myfft_saved_offsets);
01331 EXITFUNC;
01332
01333 }
01334
01335
01336
01337 void EMData::nn_SSNR(EMData* wptr, EMData* wptr2, EMData* myfft, const Transform& tf, int)
01338 {
01339 ENTERFUNC;
01340 int nxc = attr_dict["nxc"];
01341
01342 vector<int> saved_offsets = get_array_offsets();
01343 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01344
01345 set_array_offsets(0,1,1);
01346 myfft->set_array_offsets(0,1);
01347
01348 int iymin = is_fftodd() ? -ny/2 : -ny/2 + 1 ;
01349 int iymax = ny/2;
01350 int izmin = is_fftodd() ? -nz/2 : -nz/2 + 1 ;
01351 int izmax = nz/2;
01352
01353 for (int iy = iymin; iy <= iymax; iy++) {
01354 int jp = iy >= 0 ? iy+1 : ny+iy+1;
01355 for (int ix = 0; ix <= nxc; ix++) {
01356 if (( 4*(ix*ix+iy*iy) < ny*ny ) && !( ix == 0 && iy < 0 ) ) {
01357 float xnew = ix*tf[0][0] + iy*tf[1][0];
01358 float ynew = ix*tf[0][1] + iy*tf[1][1];
01359 float znew = ix*tf[0][2] + iy*tf[1][2];
01360 std::complex<float> btq;
01361 if (xnew < 0.0) {
01362 xnew = -xnew;
01363 ynew = -ynew;
01364 znew = -znew;
01365 btq = conj(myfft->cmplx(ix,jp));
01366 } else {
01367 btq = myfft->cmplx(ix,jp);
01368 }
01369 int ixn = int(xnew + 0.5 + nx) - nx;
01370 int iyn = int(ynew + 0.5 + ny) - ny;
01371 int izn = int(znew + 0.5 + nz) - nz;
01372 if ((ixn <= nxc) && (iyn >= iymin) && (iyn <= iymax) && (izn >= izmin) && (izn <= izmax)) {
01373 if (ixn >= 0) {
01374 int iza, iya;
01375 if (izn >= 0) iza = izn + 1;
01376 else iza = nz + izn + 1;
01377
01378 if (iyn >= 0) iya = iyn + 1;
01379 else iya = ny + iyn + 1;
01380
01381 cmplx(ixn,iya,iza) += btq;
01382 (*wptr)(ixn,iya,iza)++;
01383 (*wptr2)(ixn,iya,iza) += norm(btq);
01384 } else {
01385 int izt, iyt;
01386 if (izn > 0) izt = nz - izn + 1;
01387 else izt = -izn + 1;
01388
01389 if (iyn > 0) iyt = ny - iyn + 1;
01390 else iyt = -iyn + 1;
01391
01392 cmplx(-ixn,iyt,izt) += conj(btq);
01393 (*wptr)(-ixn,iyt,izt)++;
01394 (*wptr2)(-ixn,iyt,izt) += norm(btq);
01395 }
01396 }
01397 }
01398 }
01399 }
01400 set_array_offsets(saved_offsets);
01401 myfft->set_array_offsets(myfft_saved_offsets);
01402 EXITFUNC;
01403 }
01404
01405
01406
01407 void EMData::symplane0(EMData* wptr) {
01408 ENTERFUNC;
01409 int nxc = attr_dict["nxc"];
01410 int n = nxc*2;
01411
01412 vector<int> saved_offsets = get_array_offsets();
01413 set_array_offsets(0,1,1);
01414 for (int iza = 2; iza <= nxc; iza++) {
01415 for (int iya = 2; iya <= nxc; iya++) {
01416 cmplx(0,iya,iza) += conj(cmplx(0,n-iya+2,n-iza+2));
01417 (*wptr)(0,iya,iza) += (*wptr)(0,n-iya+2,n-iza+2);
01418 cmplx(0,n-iya+2,n-iza+2) = conj(cmplx(0,iya,iza));
01419 (*wptr)(0,n-iya+2,n-iza+2) = (*wptr)(0,iya,iza);
01420 cmplx(0,n-iya+2,iza) += conj(cmplx(0,iya,n-iza+2));
01421 (*wptr)(0,n-iya+2,iza) += (*wptr)(0,iya,n-iza+2);
01422 cmplx(0,iya,n-iza+2) = conj(cmplx(0,n-iya+2,iza));
01423 (*wptr)(0,iya,n-iza+2) = (*wptr)(0,n-iya+2,iza);
01424 }
01425 }
01426 for (int iya = 2; iya <= nxc; iya++) {
01427 cmplx(0,iya,1) += conj(cmplx(0,n-iya+2,1));
01428 (*wptr)(0,iya,1) += (*wptr)(0,n-iya+2,1);
01429 cmplx(0,n-iya+2,1) = conj(cmplx(0,iya,1));
01430 (*wptr)(0,n-iya+2,1) = (*wptr)(0,iya,1);
01431 }
01432 for (int iza = 2; iza <= nxc; iza++) {
01433 cmplx(0,1,iza) += conj(cmplx(0,1,n-iza+2));
01434 (*wptr)(0,1,iza) += (*wptr)(0,1,n-iza+2);
01435 cmplx(0,1,n-iza+2) = conj(cmplx(0,1,iza));
01436 (*wptr)(0,1,n-iza+2) = (*wptr)(0,1,iza);
01437 }
01438 EXITFUNC;
01439 }
01440
01441 void EMData::symplane1(EMData* wptr, EMData* wptr2) {
01442 ENTERFUNC;
01443 int nxc = attr_dict["nxc"];
01444 int n = nxc*2;
01445 vector<int> saved_offsets = get_array_offsets();
01446 set_array_offsets(0,1,1);
01447 for (int iza = 2; iza <= nxc; iza++) {
01448 for (int iya = 2; iya <= nxc; iya++) {
01449 cmplx(0,iya,iza) += conj(cmplx(0,n-iya+2,n-iza+2));
01450 (*wptr)(0,iya,iza) += (*wptr)(0,n-iya+2,n-iza+2);
01451 (*wptr2)(0,iya,iza) += (*wptr2)(0,n-iya+2,n-iza+2);
01452 cmplx(0,n-iya+2,n-iza+2) = conj(cmplx(0,iya,iza));
01453 (*wptr)(0,n-iya+2,n-iza+2) = (*wptr)(0,iya,iza);
01454 (*wptr2)(0,n-iya+2,n-iza+2) = (*wptr2)(0,iya,iza);
01455 cmplx(0,n-iya+2,iza) += conj(cmplx(0,iya,n-iza+2));
01456 (*wptr)(0,n-iya+2,iza) += (*wptr)(0,iya,n-iza+2);
01457 (*wptr2)(0,n-iya+2,iza) += (*wptr2)(0,iya,n-iza+2);
01458 cmplx(0,iya,n-iza+2) = conj(cmplx(0,n-iya+2,iza));
01459 (*wptr)(0,iya,n-iza+2) = (*wptr)(0,n-iya+2,iza);
01460 (*wptr2)(0,iya,n-iza+2) = (*wptr2)(0,n-iya+2,iza);
01461 }
01462 }
01463 for (int iya = 2; iya <= nxc; iya++) {
01464 cmplx(0,iya,1) += conj(cmplx(0,n-iya+2,1));
01465 (*wptr)(0,iya,1) += (*wptr)(0,n-iya+2,1);
01466 (*wptr2)(0,iya,1) += (*wptr2)(0,n-iya+2,1);
01467 cmplx(0,n-iya+2,1) = conj(cmplx(0,iya,1));
01468 (*wptr)(0,n-iya+2,1) = (*wptr)(0,iya,1);
01469 (*wptr2)(0,n-iya+2,1) = (*wptr2)(0,iya,1);
01470 }
01471 for (int iza = 2; iza <= nxc; iza++) {
01472 cmplx(0,1,iza) += conj(cmplx(0,1,n-iza+2));
01473 (*wptr)(0,1,iza) += (*wptr)(0,1,n-iza+2);
01474 (*wptr2)(0,1,iza) += (*wptr2)(0,1,n-iza+2);
01475 cmplx(0,1,n-iza+2) = conj(cmplx(0,1,iza));
01476 (*wptr)(0,1,n-iza+2) = (*wptr)(0,1,iza);
01477 (*wptr2)(0,1,n-iza+2) = (*wptr2)(0,1,iza);
01478 }
01479 EXITFUNC;
01480 }
01481
01482 void EMData::symplane2(EMData* wptr, EMData* wptr2, EMData* wptr3) {
01483 ENTERFUNC;
01484 int nxc = attr_dict["nxc"];
01485 int n = nxc*2;
01486 vector<int> saved_offsets = get_array_offsets();
01487 set_array_offsets(0,1,1);
01488 for (int iza = 2; iza <= nxc; iza++) {
01489 for (int iya = 2; iya <= nxc; iya++) {
01490 cmplx(0,iya,iza) += conj(cmplx(0,n-iya+2,n-iza+2));
01491 (*wptr)(0,iya,iza) += (*wptr)(0,n-iya+2,n-iza+2);
01492 (*wptr2)(0,iya,iza) += (*wptr2)(0,n-iya+2,n-iza+2);
01493 (*wptr3)(0,iya,iza) += (*wptr3)(0,n-iya+2,n-iza+2);
01494
01495 cmplx(0,n-iya+2,n-iza+2) = conj(cmplx(0,iya,iza));
01496 (*wptr)(0,n-iya+2,n-iza+2) = (*wptr)(0,iya,iza);
01497 (*wptr2)(0,n-iya+2,n-iza+2) = (*wptr2)(0,iya,iza);
01498 (*wptr3)(0,n-iya+2,n-iza+2) = (*wptr3)(0,iya,iza);
01499
01500 cmplx(0,n-iya+2,iza) += conj(cmplx(0,iya,n-iza+2));
01501 (*wptr)(0,n-iya+2,iza) += (*wptr)(0,iya,n-iza+2);
01502 (*wptr2)(0,n-iya+2,iza) += (*wptr2)(0,iya,n-iza+2);
01503 (*wptr3)(0,n-iya+2,iza) += (*wptr3)(0,iya,n-iza+2);
01504
01505 cmplx(0,iya,n-iza+2) = conj(cmplx(0,n-iya+2,iza));
01506 (*wptr)(0,iya,n-iza+2) = (*wptr)(0,n-iya+2,iza);
01507 (*wptr2)(0,iya,n-iza+2) = (*wptr2)(0,n-iya+2,iza);
01508 (*wptr3)(0,iya,n-iza+2) = (*wptr3)(0,n-iya+2,iza);
01509 }
01510 }
01511 for (int iya = 2; iya <= nxc; iya++) {
01512 cmplx(0,iya,1) += conj(cmplx(0,n-iya+2,1));
01513 (*wptr)(0,iya,1) += (*wptr)(0,n-iya+2,1);
01514 (*wptr2)(0,iya,1) += (*wptr2)(0,n-iya+2,1);
01515 (*wptr3)(0,iya,1) += (*wptr3)(0,n-iya+2,1);
01516
01517 cmplx(0,n-iya+2,1) = conj(cmplx(0,iya,1));
01518 (*wptr)(0,n-iya+2,1) = (*wptr)(0,iya,1);
01519 (*wptr2)(0,n-iya+2,1) = (*wptr2)(0,iya,1);
01520 (*wptr3)(0,n-iya+2,1) = (*wptr3)(0,iya,1);
01521 }
01522 for (int iza = 2; iza <= nxc; iza++) {
01523 cmplx(0,1,iza) += conj(cmplx(0,1,n-iza+2));
01524 (*wptr)(0,1,iza) += (*wptr)(0,1,n-iza+2);
01525 (*wptr2)(0,1,iza) += (*wptr2)(0,1,n-iza+2);
01526 (*wptr3)(0,1,iza) += (*wptr3)(0,1,n-iza+2);
01527
01528 cmplx(0,1,n-iza+2) = conj(cmplx(0,1,iza));
01529 (*wptr)(0,1,n-iza+2) = (*wptr)(0,1,iza);
01530 (*wptr2)(0,1,n-iza+2) = (*wptr2)(0,1,iza);
01531 (*wptr3)(0,1,n-iza+2) = (*wptr3)(0,1,iza);
01532 }
01533 EXITFUNC;
01534 }
01535
01536
01537 class ctf_store
01538 {
01539 public:
01540
01541 static void init( int winsize, const Ctf* ctf )
01542 {
01543 Dict params = ctf->to_dict();
01544
01545 m_winsize = winsize;
01546
01547 m_voltage = params["voltage"];
01548 m_pixel = params["apix"];
01549 m_cs = params["cs"];
01550 m_ampcont = params["ampcont"];
01551 m_bfactor = params["bfactor"];
01552 m_defocus = params["defocus"];
01553 m_winsize2= m_winsize*m_winsize;
01554 m_vecsize = m_winsize2/4;
01555 }
01556
01557 static float get_ctf( int r2 )
01558 {
01559 float ak = std::sqrt( r2/float(m_winsize2) )/m_pixel;
01560 return Util::tf( m_defocus, ak, m_voltage, m_cs, m_ampcont, m_bfactor, 1 );
01561 }
01562
01563 private:
01564
01565 static int m_winsize, m_winsize2, m_vecsize;
01566 static float m_cs;
01567 static float m_voltage;
01568 static float m_pixel;
01569 static float m_ampcont;
01570 static float m_bfactor;
01571 static float m_defocus;
01572 };
01573
01574
01575 int ctf_store::m_winsize, ctf_store::m_winsize2, ctf_store::m_vecsize;
01576
01577 float ctf_store::m_cs, ctf_store::m_voltage, ctf_store::m_pixel;
01578 float ctf_store::m_ampcont, ctf_store::m_bfactor;
01579 float ctf_store::m_defocus;
01580
01581
01582 class ctf_store_new
01583 {
01584 public:
01585
01586 static void init( int winsize, const Ctf* ctf )
01587 {
01588 Dict params = ctf->to_dict();
01589
01590 m_winsize = winsize;
01591
01592 m_voltage = params["voltage"];
01593 m_pixel = params["apix"];
01594 m_cs = params["cs"];
01595 m_ampcont = params["ampcont"];
01596 m_bfactor = params["bfactor"];
01597 m_defocus = params["defocus"];
01598 m_winsize2= m_winsize*m_winsize;
01599 m_vecsize = m_winsize2/4;
01600 }
01601
01602 static float get_ctf( float r2 )
01603 {
01604 float ak = std::sqrt( r2/float(m_winsize2) )/m_pixel;
01605 return Util::tf( m_defocus, ak, m_voltage, m_cs, m_ampcont, m_bfactor, 1 );
01606 }
01607
01608 private:
01609
01610 static int m_winsize, m_winsize2, m_vecsize;
01611 static float m_cs;
01612 static float m_voltage;
01613 static float m_pixel;
01614 static float m_ampcont;
01615 static float m_bfactor;
01616 static float m_defocus;
01617 };
01618
01619
01620 int ctf_store_new::m_winsize, ctf_store_new::m_winsize2, ctf_store_new::m_vecsize;
01621
01622 float ctf_store_new::m_cs, ctf_store_new::m_voltage, ctf_store_new::m_pixel;
01623 float ctf_store_new::m_ampcont, ctf_store_new::m_bfactor;
01624 float ctf_store_new::m_defocus;
01625
01626
01627
01628
01629 void EMData::onelinenn_ctf(int j, int n, int n2,
01630 EMData* w, EMData* bi, const Transform& tf, int mult) {
01631
01632 int remove = bi->get_attr_default( "remove", 0 );
01633
01634 int jp = (j >= 0) ? j+1 : n+j+1;
01635
01636 for (int i = 0; i <= n2; i++) {
01637 int r2 = i*i+j*j;
01638 if ( (r2<n*n/4) && !((0==i) && (j<0)) ) {
01639 float ctf = ctf_store::get_ctf( r2 );
01640 float xnew = i*tf[0][0] + j*tf[1][0];
01641 float ynew = i*tf[0][1] + j*tf[1][1];
01642 float znew = i*tf[0][2] + j*tf[1][2];
01643 std::complex<float> btq;
01644 if (xnew < 0.) {
01645 xnew = -xnew;
01646 ynew = -ynew;
01647 znew = -znew;
01648 btq = conj(bi->cmplx(i,jp));
01649 } else btq = bi->cmplx(i,jp);
01650 int ixn = int(xnew + 0.5 + n) - n;
01651 int iyn = int(ynew + 0.5 + n) - n;
01652 int izn = int(znew + 0.5 + n) - n;
01653
01654 int iza, iya;
01655 if (izn >= 0) iza = izn + 1;
01656 else iza = n + izn + 1;
01657
01658 if (iyn >= 0) iya = iyn + 1;
01659 else iya = n + iyn + 1;
01660
01661 if(remove > 0 ) {
01662 cmplx(ixn,iya,iza) -= btq*ctf*float(mult);
01663 (*w)(ixn,iya,iza) -= ctf*ctf*mult;
01664 } else {
01665 cmplx(ixn,iya,iza) += btq*ctf*float(mult);
01666 (*w)(ixn,iya,iza) += ctf*ctf*mult;
01667 }
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708 }
01709 }
01710 }
01711
01712 void EMData::onelinenn_ctf_applied(int j, int n, int n2,
01713 EMData* w, EMData* bi, const Transform& tf, int mult) {
01714
01715 int remove = bi->get_attr_default( "remove", 0 );
01716
01717 int jp = (j >= 0) ? j+1 : n+j+1;
01718
01719 for (int i = 0; i <= n2; i++) {
01720 int r2 = i*i + j*j;
01721 if ( (r2< n*n/4) && !((0==i) && (j< 0)) ) {
01722 float ctf = ctf_store::get_ctf(r2);
01723
01724
01725 float xnew = i*tf[0][0] + j*tf[1][0];
01726 float ynew = i*tf[0][1] + j*tf[1][1];
01727 float znew = i*tf[0][2] + j*tf[1][2];
01728 std::complex<float> btq;
01729 if (xnew < 0.) {
01730 xnew = -xnew;
01731 ynew = -ynew;
01732 znew = -znew;
01733 btq = conj(bi->cmplx(i,jp));
01734 } else btq = bi->cmplx(i,jp);
01735 int ixn = int(xnew + 0.5 + n) - n;
01736 int iyn = int(ynew + 0.5 + n) - n;
01737 int izn = int(znew + 0.5 + n) - n;
01738
01739 int iza, iya;
01740 if (izn >= 0) iza = izn + 1;
01741 else iza = n + izn + 1;
01742
01743 if (iyn >= 0) iya = iyn + 1;
01744 else iya = n + iyn + 1;
01745
01746 if( remove > 0 ) {
01747 cmplx(ixn,iya,iza) -= btq*float(mult);
01748 (*w)(ixn,iya,iza) -= mult*ctf*ctf;
01749 } else {
01750 cmplx(ixn,iya,iza) += btq*float(mult);
01751 (*w)(ixn,iya,iza) += mult*ctf*ctf;
01752 }
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791 }
01792 }
01793 }
01794
01795 void
01796 EMData::nn_ctf(EMData* w, EMData* myfft, const Transform& tf, int mult) {
01797 ENTERFUNC;
01798 int nxc = attr_dict["nxc"];
01799
01800 vector<int> saved_offsets = get_array_offsets();
01801 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01802 set_array_offsets(0,1,1);
01803 myfft->set_array_offsets(0,1);
01804
01805 Ctf* ctf = myfft->get_attr("ctf");
01806 ctf_store::init( ny, ctf );
01807 if(ctf) {delete ctf; ctf=0;}
01808
01809
01810 for (int iy = -ny/2 + 1; iy <= ny/2; iy++) onelinenn_ctf(iy, ny, nxc, w, myfft, tf, mult);
01811 set_array_offsets(saved_offsets);
01812 myfft->set_array_offsets(myfft_saved_offsets);
01813 EXITFUNC;
01814 }
01815
01816 void
01817 EMData::nn_ctf_applied(EMData* w, EMData* myfft, const Transform& tf, int mult) {
01818 ENTERFUNC;
01819 int nxc = attr_dict["nxc"];
01820
01821 vector<int> saved_offsets = get_array_offsets();
01822 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01823 set_array_offsets(0,1,1);
01824 myfft->set_array_offsets(0,1);
01825
01826 Ctf* ctf = myfft->get_attr( "ctf" );
01827 ctf_store::init( ny, ctf );
01828 if(ctf) {delete ctf; ctf=0;}
01829
01830
01831
01832 for (int iy = -ny/2 + 1; iy <= ny/2; iy++) onelinenn_ctf_applied(iy, ny, nxc, w, myfft, tf, mult);
01833 set_array_offsets(saved_offsets);
01834 myfft->set_array_offsets(myfft_saved_offsets);
01835 EXITFUNC;
01836 }
01837
01838
01839 void EMData::insert_rect_slice_ctf(EMData* w, EMData* myfft, const Transform& trans, int sizeofprojection, float xratio, float yratio, int npad, int mult)
01840 {
01841 ENTERFUNC;
01842 vector<int> saved_offsets = get_array_offsets();
01843 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01844 set_array_offsets(0,1,1);
01845 myfft->set_array_offsets(0,1);
01846
01847
01848
01849 Vec2f coordinate_2d_square;
01850 Vec3f coordinate_3dnew;
01851 Vec3f axis_newx;
01852 Vec3f axis_newy;
01853 Vec3f tempv;
01854
01855
01856
01857 axis_newx[0] = xratio*0.5f*(sizeofprojection*npad)*trans[0][0];
01858 axis_newx[1] = yratio*0.5f*(sizeofprojection*npad)*trans[0][1];
01859 axis_newx[2] = 0.5f*(sizeofprojection*npad)*trans[0][2];
01860
01861 float ellipse_length_x = std::sqrt(axis_newx[0]*axis_newx[0]+axis_newx[1]*axis_newx[1]+axis_newx[2]*axis_newx[2]);
01862
01863 int ellipse_length_x_int = int(ellipse_length_x);
01864 float ellipse_step_x = 0.5f*(sizeofprojection*npad)/float(ellipse_length_x_int);
01865 float xscale = ellipse_step_x;
01866
01867 axis_newy[0] = xratio*0.5f*(sizeofprojection*npad)*trans[1][0];
01868 axis_newy[1] = yratio*0.5f*(sizeofprojection*npad)*trans[1][1];
01869 axis_newy[2] = 0.5f*(sizeofprojection*npad)*trans[1][2];
01870
01871
01872
01873 float ellipse_length_y = std::sqrt(axis_newy[0]*axis_newy[0]+axis_newy[1]*axis_newy[1]+axis_newy[2]*axis_newy[2]);
01874 int ellipse_length_y_int = int(ellipse_length_y);
01875 float ellipse_step_y = 0.5f*(sizeofprojection*npad)/float(ellipse_length_y_int);
01876 float yscale = ellipse_step_y;
01877
01878 std::complex<float> c1;
01879 nz = get_zsize();
01880 Ctf* ctf = myfft->get_attr( "ctf" );
01881 ctf_store_new::init( nz, ctf );
01882 if(ctf) {delete ctf; ctf=0;}
01883 int remove = myfft->get_attr_default( "remove", 0 );
01884
01885 float r2=0.25f*sizeofprojection*npad*sizeofprojection*npad;
01886 float r2_at_point;
01887
01888 for(int i=0;i<ellipse_length_x_int;i++) {
01889 for(int j=-1*ellipse_length_y_int+1; j<=ellipse_length_y_int; j++) {
01890
01891 r2_at_point=i*xscale*i*xscale+j*yscale*j*yscale;
01892 if(r2_at_point<=r2 && ! ((0==i) && (j<0))) {
01893
01894 float ctf_value = ctf_store_new::get_ctf( r2_at_point );
01895 coordinate_2d_square[0] = xscale*float(i);
01896 coordinate_2d_square[1] = yscale*float(j);
01897 float xnew = coordinate_2d_square[0]*trans[0][0] + coordinate_2d_square[1]*trans[1][0];
01898 float ynew = coordinate_2d_square[0]*trans[0][1] + coordinate_2d_square[1]*trans[1][1];
01899 float znew = coordinate_2d_square[0]*trans[0][2] + coordinate_2d_square[1]*trans[1][2];
01900 coordinate_3dnew[0] =xnew*xratio;
01901 coordinate_3dnew[1] = ynew*yratio;
01902 coordinate_3dnew[2] = znew;
01903
01904
01905 float xp = coordinate_2d_square[0];
01906 float yp = ( coordinate_2d_square[1] >= 0) ? coordinate_2d_square[1]+1 : nz+coordinate_2d_square[1]+1;
01907 std::complex<float> lin_interpolated(0,0);
01908 int xlow=int(xp),xhigh=int(xp)+1;
01909 int ylow=int(yp),yhigh=int(yp)+1;
01910 float tx=xp-xlow,ty=yp-ylow;
01911
01912
01913 if(j == -1) {
01914
01915 if(ylow<yp)
01916 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)*(1-ty) + myfft->cmplx(xlow,yhigh)*(1-tx)*ty
01917 + myfft->cmplx(xhigh,ylow)*tx*(1-ty) + myfft->cmplx(xhigh,yhigh)*tx*ty;
01918 else
01919 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)
01920 + myfft->cmplx(xhigh,ylow)*tx;
01921
01922 }
01923 else {
01924 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)*(1-ty) + myfft->cmplx(xlow,yhigh)*(1-tx)*ty
01925 + myfft->cmplx(xhigh,ylow)*tx*(1-ty)+ myfft->cmplx(xhigh,yhigh)*tx*ty;
01926
01927 }
01928
01929 c1 = lin_interpolated;
01930
01931
01932
01933 std::complex<float> btq;
01934 if ( coordinate_3dnew[0] < 0.) {
01935 coordinate_3dnew[0] = -coordinate_3dnew[0];
01936 coordinate_3dnew[1] = -coordinate_3dnew[1];
01937 coordinate_3dnew[2] = -coordinate_3dnew[2];
01938 btq = conj(c1);
01939 } else {
01940 btq = c1;
01941 }
01942 int ixn = int(coordinate_3dnew[0] + 0.5 + nx) - nx;
01943 int iyn = int(coordinate_3dnew[1] + 0.5 + ny) - ny;
01944 int izn = int(coordinate_3dnew[2] + 0.5 + nz) - nz;
01945
01946 int iza, iya;
01947 if (izn >= 0) iza = izn + 1;
01948 else iza = nz + izn + 1;
01949
01950 if (iyn >= 0) iya = iyn + 1;
01951 else iya = ny + iyn + 1;
01952
01953 if(remove > 0 ) {
01954 cmplx(ixn,iya,iza) -= btq*ctf_value*float(mult);
01955 (*w)(ixn,iya,iza) -= ctf_value*ctf_value*mult;
01956 } else {
01957 cmplx(ixn,iya,iza) += btq*ctf_value*float(mult);
01958 (*w)(ixn,iya,iza) += ctf_value*ctf_value*mult;
01959 }
01960
01961 }
01962 }
01963
01964 }
01965
01966
01967
01968
01969 set_array_offsets(saved_offsets);
01970 myfft->set_array_offsets(myfft_saved_offsets);
01971 EXITFUNC;
01972
01973 }
01974
01975
01976 void EMData::insert_rect_slice_ctf_applied(EMData* w, EMData* myfft,const Transform& trans,int sizeofprojection,float xratio,float yratio,int npad,int mult)
01977 {
01978 ENTERFUNC;
01979 vector<int> saved_offsets = get_array_offsets();
01980 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01981 set_array_offsets(0,1,1);
01982 myfft->set_array_offsets(0,1);
01983
01984
01985
01986 Vec2f coordinate_2d_square;
01987 Vec3f coordinate_3dnew;
01988 Vec3f axis_newx;
01989 Vec3f axis_newy;
01990 Vec3f tempv;
01991
01992
01993
01994 axis_newx[0] = xratio*0.5f*(sizeofprojection*npad)*trans[0][0];
01995 axis_newx[1] = yratio*0.5f*(sizeofprojection*npad)*trans[0][1];
01996 axis_newx[2] = 0.5f*(sizeofprojection*npad)*trans[0][2];
01997
01998 float ellipse_length_x = std::sqrt(axis_newx[0]*axis_newx[0]+axis_newx[1]*axis_newx[1]+axis_newx[2]*axis_newx[2]);
01999
02000 int ellipse_length_x_int = int(ellipse_length_x);
02001 float ellipse_step_x = 0.5f*(sizeofprojection*npad)/float(ellipse_length_x_int);
02002 float xscale = ellipse_step_x;
02003
02004 axis_newy[0] = xratio*0.5f*(sizeofprojection*npad)*trans[1][0];
02005 axis_newy[1] = yratio*0.5f*(sizeofprojection*npad)*trans[1][1];
02006 axis_newy[2] = 0.5f*(sizeofprojection*npad)*trans[1][2];
02007
02008
02009
02010 float ellipse_length_y = std::sqrt(axis_newy[0]*axis_newy[0]+axis_newy[1]*axis_newy[1]+axis_newy[2]*axis_newy[2]);
02011 int ellipse_length_y_int = int(ellipse_length_y);
02012 float ellipse_step_y = 0.5f*(sizeofprojection*npad)/float(ellipse_length_y_int);
02013 float yscale = ellipse_step_y;
02014
02015 std::complex<float> c1;
02016 nz = get_zsize();
02017 Ctf* ctf = myfft->get_attr( "ctf" );
02018 ctf_store_new::init( nz, ctf );
02019 if(ctf) {delete ctf; ctf=0;}
02020 int remove = myfft->get_attr_default( "remove", 0 );
02021
02022 float r2=0.25f*sizeofprojection*npad*sizeofprojection*npad;
02023 float r2_at_point;
02024
02025 for(int i=0;i<ellipse_length_x_int;i++) {
02026 for(int j=-1*ellipse_length_y_int+1; j<=ellipse_length_y_int; j++) {
02027
02028 r2_at_point=i*xscale*i*xscale+j*yscale*j*yscale;
02029 if(r2_at_point<=r2 && ! ((0==i) && (j<0))) {
02030
02031 float ctf_value = ctf_store_new::get_ctf( r2_at_point );
02032 coordinate_2d_square[0] = xscale*float(i);
02033 coordinate_2d_square[1] = yscale*float(j);
02034 float xnew = coordinate_2d_square[0]*trans[0][0] + coordinate_2d_square[1]*trans[1][0];
02035 float ynew = coordinate_2d_square[0]*trans[0][1] + coordinate_2d_square[1]*trans[1][1];
02036 float znew = coordinate_2d_square[0]*trans[0][2] + coordinate_2d_square[1]*trans[1][2];
02037 coordinate_3dnew[0] =xnew*xratio;
02038 coordinate_3dnew[1] = ynew*yratio;
02039 coordinate_3dnew[2] = znew;
02040
02041
02042 float xp = coordinate_2d_square[0];
02043 float yp = ( coordinate_2d_square[1] >= 0) ? coordinate_2d_square[1]+1 : nz+coordinate_2d_square[1]+1;
02044 std::complex<float> lin_interpolated(0,0);
02045 int xlow=int(xp),xhigh=int(xp)+1;
02046 int ylow=int(yp),yhigh=int(yp)+1;
02047 float tx=xp-xlow,ty=yp-ylow;
02048
02049
02050 if(j == -1) {
02051
02052 if(ylow<yp)
02053 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)*(1-ty) + myfft->cmplx(xlow,yhigh)*(1-tx)*ty
02054 + myfft->cmplx(xhigh,ylow)*tx*(1-ty) + myfft->cmplx(xhigh,yhigh)*tx*ty;
02055 else
02056 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)
02057 + myfft->cmplx(xhigh,ylow)*tx;
02058
02059 }
02060 else {
02061 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)*(1-ty) + myfft->cmplx(xlow,yhigh)*(1-tx)*ty
02062 + myfft->cmplx(xhigh,ylow)*tx*(1-ty)+ myfft->cmplx(xhigh,yhigh)*tx*ty;
02063
02064 }
02065
02066 c1 = lin_interpolated;
02067
02068
02069
02070 std::complex<float> btq;
02071 if ( coordinate_3dnew[0] < 0.) {
02072 coordinate_3dnew[0] = -coordinate_3dnew[0];
02073 coordinate_3dnew[1] = -coordinate_3dnew[1];
02074 coordinate_3dnew[2] = -coordinate_3dnew[2];
02075 btq = conj(c1);
02076 } else {
02077 btq = c1;
02078 }
02079 int ixn = int(coordinate_3dnew[0] + 0.5 + nx) - nx;
02080 int iyn = int(coordinate_3dnew[1] + 0.5 + ny) - ny;
02081 int izn = int(coordinate_3dnew[2] + 0.5 + nz) - nz;
02082
02083 int iza, iya;
02084 if (izn >= 0) iza = izn + 1;
02085 else iza = nz + izn + 1;
02086
02087 if (iyn >= 0) iya = iyn + 1;
02088 else iya = ny + iyn + 1;
02089
02090 if(remove > 0 ) {
02091 cmplx(ixn,iya,iza) -= btq*float(mult);
02092 (*w)(ixn,iya,iza) -= ctf_value*ctf_value*mult;
02093 } else {
02094 cmplx(ixn,iya,iza) += btq*float(mult);
02095 (*w)(ixn,iya,iza) += ctf_value*ctf_value*mult;
02096 }
02097
02098 }
02099 }
02100
02101 }
02102
02103
02104
02105
02106 set_array_offsets(saved_offsets);
02107 myfft->set_array_offsets(myfft_saved_offsets);
02108 EXITFUNC;
02109
02110 }
02111
02112
02113
02114
02115
02116
02117
02118
02119
02120
02121
02122
02123
02124
02125
02126
02127
02128
02129
02130
02131
02132
02133
02134
02135
02136
02137
02138
02139
02140
02141
02142
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152
02153
02154
02155
02156
02157
02158
02159
02160
02161
02162
02163
02164
02165
02166
02167
02168
02169
02170
02171
02172
02173
02174
02175
02176
02177
02178
02179
02180
02181 void EMData::nn_SSNR_ctf(EMData* wptr, EMData* wptr2, EMData* wptr3, EMData* myfft, const Transform& tf, int)
02182 {
02183
02184
02185
02186
02187
02188
02189
02190
02191
02192 ENTERFUNC;
02193 int nxc = attr_dict["nxc"];
02194 vector<int> saved_offsets = get_array_offsets();
02195 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
02196 set_array_offsets(0,1,1);
02197 myfft->set_array_offsets(0,1);
02198
02199 Ctf* ctf = myfft->get_attr("ctf");
02200 ctf_store::init( ny, ctf );
02201 int iymin = is_fftodd() ? -ny/2 : -ny/2 + 1;
02202 int iymax = ny/2;
02203 int izmin = is_fftodd() ? -nz/2 : -nz/2 + 1;
02204 int izmax = nz/2;
02205
02206 for (int iy = iymin; iy <= iymax; iy++) {
02207 int jp = iy >= 0 ? iy+1 : ny+iy+1;
02208 for (int ix = 0; ix <= nxc; ix++) {
02209 int r2 = ix*ix+iy*iy;
02210 if (( 4*r2 < ny*ny ) && !( ix == 0 && iy < 0 ) ) {
02211 float ctf = ctf_store::get_ctf( r2 )*10.f;
02212 float xnew = ix*tf[0][0] + iy*tf[1][0];
02213 float ynew = ix*tf[0][1] + iy*tf[1][1];
02214 float znew = ix*tf[0][2] + iy*tf[1][2];
02215 std::complex<float> btq;
02216 if (xnew < 0.0) {
02217 xnew = -xnew;
02218 ynew = -ynew;
02219 znew = -znew;
02220 btq = conj(myfft->cmplx(ix,jp));
02221 } else {
02222 btq = myfft->cmplx(ix,jp);
02223 }
02224 int ixn = int(xnew + 0.5 + nx) - nx;
02225 int iyn = int(ynew + 0.5 + ny) - ny;
02226 int izn = int(znew + 0.5 + nz) - nz;
02227 if ((ixn <= nxc) && (iyn >= iymin) && (iyn <= iymax) && (izn >= izmin) && (izn <= izmax)) {
02228 if (ixn >= 0) {
02229 int iza, iya;
02230 if (izn >= 0) iza = izn + 1;
02231 else iza = nz + izn + 1;
02232
02233 if (iyn >= 0) iya = iyn + 1;
02234 else iya = ny + iyn + 1;
02235
02236 cmplx(ixn,iya,iza) += btq*ctf;
02237 (*wptr)(ixn,iya,iza) += ctf*ctf;
02238 (*wptr2)(ixn,iya,iza) += std::norm(btq);
02239 (*wptr3)(ixn,iya,iza) += 1;
02240 } else {
02241 int izt, iyt;
02242 if (izn > 0) izt = nz - izn + 1;
02243 else izt = -izn + 1;
02244
02245 if (iyn > 0) iyt = ny - iyn + 1;
02246 else iyt = -iyn + 1;
02247
02248 cmplx(-ixn,iyt,izt) += std::conj(btq)*ctf;
02249 (*wptr) (-ixn,iyt,izt) += ctf*ctf;
02250 (*wptr2)(-ixn,iyt,izt) += std::norm(btq);
02251 (*wptr3)(-ixn,iyt,izt) += 1;
02252 }
02253 }
02254 }
02255 }
02256 }
02257 set_array_offsets(saved_offsets);
02258 myfft->set_array_offsets(myfft_saved_offsets);
02259 if(ctf) {delete ctf; ctf=0;}
02260 EXITFUNC;
02261 }
02262
02263
02264
02265
02266
02267
02268
02269
02270
02271
02272
02273
02274
02275
02276
02277
02278
02279
02280
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314
02315
02316
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339
02340
02341
02342
02343
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353
02354
02355
02356
02357
02358
02359
02360
02361 void EMData::symplane0_ctf(EMData* w) {
02362 ENTERFUNC;
02363 int nxc = attr_dict["nxc"];
02364 int n = nxc*2;
02365
02366 vector<int> saved_offsets = get_array_offsets();
02367 set_array_offsets(0,1,1);
02368 for (int iza = 2; iza <= nxc; iza++) {
02369 for (int iya = 2; iya <= nxc; iya++) {
02370 cmplx(0,iya,iza) += conj(cmplx(0,n-iya+2,n-iza+2));
02371 (*w)(0,iya,iza) += (*w)(0,n-iya+2,n-iza+2);
02372 cmplx(0,n-iya+2,n-iza+2) = conj(cmplx(0,iya,iza));
02373 (*w)(0,n-iya+2,n-iza+2) = (*w)(0,iya,iza);
02374 cmplx(0,n-iya+2,iza) += conj(cmplx(0,iya,n-iza+2));
02375 (*w)(0,n-iya+2,iza) += (*w)(0,iya,n-iza+2);
02376 cmplx(0,iya,n-iza+2) = conj(cmplx(0,n-iya+2,iza));
02377 (*w)(0,iya,n-iza+2) = (*w)(0,n-iya+2,iza);
02378 }
02379 }
02380 for (int iya = 2; iya <= nxc; iya++) {
02381 cmplx(0,iya,1) += conj(cmplx(0,n-iya+2,1));
02382 (*w)(0,iya,1) += (*w)(0,n-iya+2,1);
02383 cmplx(0,n-iya+2,1) = conj(cmplx(0,iya,1));
02384 (*w)(0,n-iya+2,1) = (*w)(0,iya,1);
02385 }
02386 for (int iza = 2; iza <= nxc; iza++) {
02387 cmplx(0,1,iza) += conj(cmplx(0,1,n-iza+2));
02388 (*w)(0,1,iza) += (*w)(0,1,n-iza+2);
02389 cmplx(0,1,n-iza+2) = conj(cmplx(0,1,iza));
02390 (*w)(0,1,n-iza+2) = (*w)(0,1,iza);
02391 }
02392 EXITFUNC;
02393 }
02394
02395 void EMData::symplane0_rect(EMData* w) {
02396 ENTERFUNC;
02397 nx=get_xsize();
02398 ny=get_ysize();
02399 nz=get_zsize();
02400 int nzc=nz/2;
02401 int nyc=ny/2;
02402
02403
02404
02405 vector<int> saved_offsets = get_array_offsets();
02406 set_array_offsets(0,1,1);
02407 for (int iza = 2; iza <= nzc; iza++) {
02408 for (int iya = 2; iya <= nyc; iya++) {
02409 cmplx(0,iya,iza) += conj(cmplx(0,ny-iya+2,nz-iza+2));
02410 (*w)(0,iya,iza) += (*w)(0,ny-iya+2,nz-iza+2);
02411 cmplx(0,ny-iya+2,nz-iza+2) = conj(cmplx(0,iya,iza));
02412 (*w)(0,ny-iya+2,nz-iza+2) = (*w)(0,iya,iza);
02413 cmplx(0,ny-iya+2,iza) += conj(cmplx(0,iya,nz-iza+2));
02414 (*w)(0,ny-iya+2,iza) += (*w)(0,iya,nz-iza+2);
02415 cmplx(0,iya,nz-iza+2) = conj(cmplx(0,ny-iya+2,iza));
02416 (*w)(0,iya,nz-iza+2) = (*w)(0,ny-iya+2,iza);
02417 }
02418 }
02419 for (int iya = 2; iya <= nyc; iya++) {
02420 cmplx(0,iya,1) += conj(cmplx(0,ny-iya+2,1));
02421 (*w)(0,iya,1) += (*w)(0,ny-iya+2,1);
02422 cmplx(0,ny-iya+2,1) = conj(cmplx(0,iya,1));
02423 (*w)(0,ny-iya+2,1) = (*w)(0,iya,1);
02424 }
02425 for (int iza = 2; iza <= nzc; iza++) {
02426 cmplx(0,1,iza) += conj(cmplx(0,1,nz-iza+2));
02427 (*w)(0,1,iza) += (*w)(0,1,nz-iza+2);
02428 cmplx(0,1,nz-iza+2) = conj(cmplx(0,1,iza));
02429 (*w)(0,1,nz-iza+2) = (*w)(0,1,iza);
02430 }
02431 EXITFUNC;
02432 }
02433
02434 EMData* EMData::rot_scale_trans2D(float angDeg, float delx, float dely, float scale) {
02435 float ang=angDeg*M_PI/180.0f;
02436 if (1 >= ny)
02437 throw ImageDimensionException("Can't rotate 1D image");
02438 if (nz<2) {
02439 vector<int> saved_offsets = get_array_offsets();
02440 set_array_offsets(0,0,0);
02441 if (0.0f == scale) scale = 1.0f;
02442 EMData* ret = copy_head();
02443 delx = restrict2(delx, nx);
02444 dely = restrict2(dely, ny);
02445
02446 int xc = nx/2;
02447 int yc = ny/2;
02448
02449 float shiftxc = xc + delx;
02450 float shiftyc = yc + dely;
02451
02452 float cang = cos(ang);
02453 float sang = sin(ang);
02454 for (int iy = 0; iy < ny; iy++) {
02455 float y = float(iy) - shiftyc;
02456 float ycang = y*cang/scale + yc;
02457 float ysang = -y*sang/scale + xc;
02458 for (int ix = 0; ix < nx; ix++) {
02459 float x = float(ix) - shiftxc;
02460 float xold = x*cang/scale + ysang ;
02461 float yold = x*sang/scale + ycang ;
02462
02463 (*ret)(ix,iy) = Util::quadri(xold+1.0f, yold+1.0f, nx, ny, get_data());
02464
02465 }
02466 }
02467 set_array_offsets(saved_offsets);
02468 return ret;
02469 } else {
02470 throw ImageDimensionException("Volume not currently supported");
02471 }
02472 }
02473
02474 EMData* EMData::rot_scale_trans2D_background(float angDeg, float delx, float dely, float scale) {
02475 float ang=angDeg*M_PI/180.0f;
02476 if (1 >= ny)
02477 throw ImageDimensionException("Can't rotate 1D image");
02478 if (nz<2) {
02479 vector<int> saved_offsets = get_array_offsets();
02480 set_array_offsets(0,0,0);
02481 if (0.0f == scale) scale = 1.0f;
02482 EMData* ret = copy_head();
02483 delx = restrict2(delx, nx);
02484 dely = restrict2(dely, ny);
02485
02486 int xc = nx/2;
02487 int yc = ny/2;
02488
02489 float shiftxc = xc + delx;
02490 float shiftyc = yc + dely;
02491
02492 float cang = cos(ang);
02493 float sang = sin(ang);
02494 for (int iy = 0; iy < ny; iy++) {
02495 float y = float(iy) - shiftyc;
02496 float ycang = y*cang/scale + yc;
02497 float ysang = -y*sang/scale + xc;
02498 for (int ix = 0; ix < nx; ix++) {
02499 float x = float(ix) - shiftxc;
02500 float xold = x*cang/scale + ysang ;
02501 float yold = x*sang/scale + ycang ;
02502
02503 (*ret)(ix,iy) = Util::quadri_background(xold+1.0f, yold+1.0f, nx, ny, get_data(),ix+1,iy+1);
02504
02505 }
02506 }
02507 set_array_offsets(saved_offsets);
02508 return ret;
02509 } else {
02510 throw ImageDimensionException("Volume not currently supported");
02511 }
02512 }
02513
02514 #define in(i,j,k) in[i+(j+(k*ny))*(size_t)nx]
02515 EMData*
02516 EMData::rot_scale_trans(const Transform &RA) {
02517
02518 EMData* ret = copy_head();
02519 float *in = this->get_data();
02520 vector<int> saved_offsets = get_array_offsets();
02521 set_array_offsets(0,0,0);
02522 Vec3f translations = RA.get_trans();
02523 Transform RAinv = RA.inverse();
02524
02525 if (1 >= ny) throw ImageDimensionException("Can't rotate 1D image");
02526 if (nz < 2) {
02527 float p1, p2, p3, p4;
02528 float delx = translations.at(0);
02529 float dely = translations.at(1);
02530 delx = restrict2(delx, nx);
02531 dely = restrict2(dely, ny);
02532 int xc = nx/2;
02533 int yc = ny/2;
02534
02535 float shiftxc = xc + delx;
02536 float shiftyc = yc + dely;
02537 for (int iy = 0; iy < ny; iy++) {
02538 float y = float(iy) - shiftyc;
02539 float ysang = y*RAinv[0][1]+xc;
02540 float ycang = y*RAinv[1][1]+yc;
02541 for (int ix = 0; ix < nx; ix++) {
02542 float x = float(ix) - shiftxc;
02543 float xold = x*RAinv[0][0] + ysang;
02544 float yold = x*RAinv[1][0] + ycang;
02545
02546 xold = restrict1(xold, nx);
02547 yold = restrict1(yold, ny);
02548
02549 int xfloor = int(xold);
02550 int yfloor = int(yold);
02551 float t = xold-xfloor;
02552 float u = yold-yfloor;
02553 if(xfloor == nx -1 && yfloor == ny -1) {
02554
02555 p1 =in[xfloor + yfloor*ny];
02556 p2 =in[ yfloor*ny];
02557 p3 =in[0];
02558 p4 =in[xfloor];
02559 } else if(xfloor == nx - 1) {
02560
02561 p1 =in[xfloor + yfloor*ny];
02562 p2 =in[ yfloor*ny];
02563 p3 =in[ (yfloor+1)*ny];
02564 p4 =in[xfloor + (yfloor+1)*ny];
02565 } else if(yfloor == ny - 1) {
02566
02567 p1 =in[xfloor + yfloor*ny];
02568 p2 =in[xfloor+1 + yfloor*ny];
02569 p3 =in[xfloor+1 ];
02570 p4 =in[xfloor ];
02571 } else {
02572 p1 =in[xfloor + yfloor*ny];
02573 p2 =in[xfloor+1 + yfloor*ny];
02574 p3 =in[xfloor+1 + (yfloor+1)*ny];
02575 p4 =in[xfloor + (yfloor+1)*ny];
02576 }
02577 (*ret)(ix,iy) = p1 + u * ( p4 - p1) + t * ( p2 - p1 + u *(p3-p2-p4+p1));
02578 }
02579 }
02580 set_array_offsets(saved_offsets);
02581 return ret;
02582 } else {
02583
02584
02585 float delx = translations.at(0);
02586 float dely = translations.at(1);
02587 float delz = translations.at(2);
02588 delx = restrict2(delx, nx);
02589 dely = restrict2(dely, ny);
02590 delz = restrict2(delz, nz);
02591 int xc = nx/2;
02592 int yc = ny/2;
02593 int zc = nz/2;
02594
02595 float shiftxc = xc + delx;
02596 float shiftyc = yc + dely;
02597 float shiftzc = zc + delz;
02598
02599 for (int iz = 0; iz < nz; iz++) {
02600 float z = float(iz) - shiftzc;
02601 float xoldz = z*RAinv[0][2]+xc;
02602 float yoldz = z*RAinv[1][2]+yc;
02603 float zoldz = z*RAinv[2][2]+zc;
02604 for (int iy = 0; iy < ny; iy++) {
02605 float y = float(iy) - shiftyc;
02606 float xoldzy = xoldz + y*RAinv[0][1] ;
02607 float yoldzy = yoldz + y*RAinv[1][1] ;
02608 float zoldzy = zoldz + y*RAinv[2][1] ;
02609 for (int ix = 0; ix < nx; ix++) {
02610 float x = float(ix) - shiftxc;
02611 float xold = xoldzy + x*RAinv[0][0] ;
02612 float yold = yoldzy + x*RAinv[1][0] ;
02613 float zold = zoldzy + x*RAinv[2][0] ;
02614
02615 xold = restrict1(xold, nx);
02616 yold = restrict1(yold, ny);
02617 zold = restrict1(zold, nz);
02618
02619
02620 int IOX = int(xold);
02621 int IOY = int(yold);
02622 int IOZ = int(zold);
02623
02624 #ifdef _WIN32
02625 int IOXp1 = _cpp_min( nx-1 ,IOX+1);
02626 #else
02627 int IOXp1 = std::min( nx-1 ,IOX+1);
02628 #endif //_WIN32
02629
02630 #ifdef _WIN32
02631 int IOYp1 = _cpp_min( ny-1 ,IOY+1);
02632 #else
02633 int IOYp1 = std::min( ny-1 ,IOY+1);
02634 #endif //_WIN32
02635
02636 #ifdef _WIN32
02637 int IOZp1 = _cpp_min( nz-1 ,IOZ+1);
02638 #else
02639 int IOZp1 = std::min( nz-1 ,IOZ+1);
02640 #endif //_WIN32
02641
02642 float dx = xold-IOX;
02643 float dy = yold-IOY;
02644 float dz = zold-IOZ;
02645
02646 float a1 = in(IOX,IOY,IOZ);
02647 float a2 = in(IOXp1,IOY,IOZ) - in(IOX,IOY,IOZ);
02648 float a3 = in(IOX,IOYp1,IOZ) - in(IOX,IOY,IOZ);
02649 float a4 = in(IOX,IOY,IOZp1) - in(IOX,IOY,IOZ);
02650 float a5 = in(IOX,IOY,IOZ) - in(IOXp1,IOY,IOZ) - in(IOX,IOYp1,IOZ) + in(IOXp1,IOYp1,IOZ);
02651 float a6 = in(IOX,IOY,IOZ) - in(IOXp1,IOY,IOZ) - in(IOX,IOY,IOZp1) + in(IOXp1,IOY,IOZp1);
02652 float a7 = in(IOX,IOY,IOZ) - in(IOX,IOYp1,IOZ) - in(IOX,IOY,IOZp1) + in(IOX,IOYp1,IOZp1);
02653 float a8 = in(IOXp1,IOY,IOZ) + in(IOX,IOYp1,IOZ)+ in(IOX,IOY,IOZp1)
02654 - in(IOX,IOY,IOZ)- in(IOXp1,IOYp1,IOZ) - in(IOXp1,IOY,IOZp1)
02655 - in(IOX,IOYp1,IOZp1) + in(IOXp1,IOYp1,IOZp1);
02656 (*ret)(ix,iy,iz) = a1 + dz*(a4 + a6*dx + (a7 + a8*dx)*dy) + a3*dy + dx*(a2 + a5*dy);
02657 }
02658 }
02659 }
02660
02661 set_array_offsets(saved_offsets);
02662 return ret;
02663
02664
02665
02666
02667
02668
02669
02670
02671
02672
02673
02674
02675
02676
02677
02678
02679
02680
02681
02682
02683
02684
02685
02686
02687
02688
02689
02690
02691
02692
02693
02694
02695
02696
02697
02698
02699
02700
02701
02702
02703
02704
02705
02706
02707
02708
02709
02710
02711
02712
02713
02714
02715
02716
02717
02718
02719
02720
02721
02722
02723
02724
02725
02726
02727
02728
02729
02730
02731
02732
02733
02734
02735
02736
02737
02738
02739
02740
02741
02742
02743
02744
02745
02746
02747
02748
02749
02750
02751
02752
02753
02754
02755
02756
02757
02758
02759
02760
02761
02762
02763
02764
02765
02766
02767
02768
02769
02770
02771
02772
02773
02774
02775
02776
02777
02778 }
02779 }
02780 #undef in
02781
02782
02783 #define in(i,j,k) in[i+(j+(k*ny))*(size_t)nx]
02784 EMData*
02785 EMData::rot_scale_trans_background(const Transform &RA) {
02786 EMData* ret = copy_head();
02787 float *in = this->get_data();
02788 vector<int> saved_offsets = get_array_offsets();
02789 set_array_offsets(0,0,0);
02790 Vec3f translations = RA.get_trans();
02791 Transform RAinv = RA.inverse();
02792
02793 if (1 >= ny) throw ImageDimensionException("Can't rotate 1D image");
02794 if (nz < 2) {
02795 float p1, p2, p3, p4;
02796 float delx = translations.at(0);
02797 float dely = translations.at(1);
02798 delx = restrict2(delx, nx);
02799 dely = restrict2(dely, ny);
02800 int xc = nx/2;
02801 int yc = ny/2;
02802
02803 float shiftxc = xc + delx;
02804 float shiftyc = yc + dely;
02805 for (int iy = 0; iy < ny; iy++) {
02806 float y = float(iy) - shiftyc;
02807 float ysang = y*RAinv[0][1]+xc;
02808 float ycang = y*RAinv[1][1]+yc;
02809 for (int ix = 0; ix < nx; ix++) {
02810 float x = float(ix) - shiftxc;
02811 float xold = x*RAinv[0][0] + ysang;
02812 float yold = x*RAinv[1][0] + ycang;
02813
02814
02815
02816 if ( (xold < 0.0f) || (xold >= (float)(nx)) || (yold < 0.0f) || (yold >= (float)(ny)) ){
02817 xold = (float)ix;
02818 yold = (float)iy;
02819 }
02820
02821 int xfloor = int(xold);
02822 int yfloor = int(yold);
02823 float t = xold-xfloor;
02824 float u = yold-yfloor;
02825 if(xfloor == nx -1 && yfloor == ny -1) {
02826
02827 p1 =in[xfloor + yfloor*ny];
02828 p2 =in[ yfloor*ny];
02829 p3 =in[0];
02830 p4 =in[xfloor];
02831 } else if(xfloor == nx - 1) {
02832
02833 p1 =in[xfloor + yfloor*ny];
02834 p2 =in[ yfloor*ny];
02835 p3 =in[ (yfloor+1)*ny];
02836 p4 =in[xfloor + (yfloor+1)*ny];
02837 } else if(yfloor == ny - 1) {
02838
02839 p1 =in[xfloor + yfloor*ny];
02840 p2 =in[xfloor+1 + yfloor*ny];
02841 p3 =in[xfloor+1 ];
02842 p4 =in[xfloor ];
02843 } else {
02844
02845 p1 =in[xfloor + yfloor*ny];
02846 p2 =in[xfloor+1 + yfloor*ny];
02847 p3 =in[xfloor+1 + (yfloor+1)*ny];
02848 p4 =in[xfloor + (yfloor+1)*ny];
02849 }
02850 (*ret)(ix,iy) = p1 + u * ( p4 - p1) + t * ( p2 - p1 + u *(p3-p2-p4+p1));
02851 }
02852 }
02853 set_array_offsets(saved_offsets);
02854 return ret;
02855 } else {
02856
02857
02858 float delx = translations.at(0);
02859 float dely = translations.at(1);
02860 float delz = translations.at(2);
02861 delx = restrict2(delx, nx);
02862 dely = restrict2(dely, ny);
02863 delz = restrict2(delz, nz);
02864 int xc = nx/2;
02865 int yc = ny/2;
02866 int zc = nz/2;
02867
02868 float shiftxc = xc + delx;
02869 float shiftyc = yc + dely;
02870 float shiftzc = zc + delz;
02871
02872 for (int iz = 0; iz < nz; iz++) {
02873 float z = float(iz) - shiftzc;
02874 float xoldz = z*RAinv[0][2]+xc;
02875 float yoldz = z*RAinv[1][2]+yc;
02876 float zoldz = z*RAinv[2][2]+zc;
02877 for (int iy = 0; iy < ny; iy++) {
02878 float y = float(iy) - shiftyc;
02879 float xoldzy = xoldz + y*RAinv[0][1] ;
02880 float yoldzy = yoldz + y*RAinv[1][1] ;
02881 float zoldzy = zoldz + y*RAinv[2][1] ;
02882 for (int ix = 0; ix < nx; ix++) {
02883 float x = float(ix) - shiftxc;
02884 float xold = xoldzy + x*RAinv[0][0] ;
02885 float yold = yoldzy + x*RAinv[1][0] ;
02886 float zold = zoldzy + x*RAinv[2][0] ;
02887
02888
02889
02890 if ( (xold < 0.0f) || (xold >= (float)(nx)) || (yold < 0.0f) || (yold >= (float)(ny)) || (zold < 0.0f) || (zold >= (float)(nz)) ){
02891 xold = (float)ix;
02892 yold = (float)iy;
02893 zold = (float)iz;
02894 }
02895
02896 int IOX = int(xold);
02897 int IOY = int(yold);
02898 int IOZ = int(zold);
02899
02900 #ifdef _WIN32
02901 int IOXp1 = _cpp_min( nx-1 ,IOX+1);
02902 #else
02903 int IOXp1 = std::min( nx-1 ,IOX+1);
02904 #endif //_WIN32
02905
02906 #ifdef _WIN32
02907 int IOYp1 = _cpp_min( ny-1 ,IOY+1);
02908 #else
02909 int IOYp1 = std::min( ny-1 ,IOY+1);
02910 #endif //_WIN32
02911
02912 #ifdef _WIN32
02913 int IOZp1 = _cpp_min( nz-1 ,IOZ+1);
02914 #else
02915 int IOZp1 = std::min( nz-1 ,IOZ+1);
02916 #endif //_WIN32
02917
02918 float dx = xold-IOX;
02919 float dy = yold-IOY;
02920 float dz = zold-IOZ;
02921
02922 float a1 = in(IOX,IOY,IOZ);
02923 float a2 = in(IOXp1,IOY,IOZ) - in(IOX,IOY,IOZ);
02924 float a3 = in(IOX,IOYp1,IOZ) - in(IOX,IOY,IOZ);
02925 float a4 = in(IOX,IOY,IOZp1) - in(IOX,IOY,IOZ);
02926 float a5 = in(IOX,IOY,IOZ) - in(IOXp1,IOY,IOZ) - in(IOX,IOYp1,IOZ) + in(IOXp1,IOYp1,IOZ);
02927 float a6 = in(IOX,IOY,IOZ) - in(IOXp1,IOY,IOZ) - in(IOX,IOY,IOZp1) + in(IOXp1,IOY,IOZp1);
02928 float a7 = in(IOX,IOY,IOZ) - in(IOX,IOYp1,IOZ) - in(IOX,IOY,IOZp1) + in(IOX,IOYp1,IOZp1);
02929 float a8 = in(IOXp1,IOY,IOZ) + in(IOX,IOYp1,IOZ)+ in(IOX,IOY,IOZp1)
02930 - in(IOX,IOY,IOZ)- in(IOXp1,IOYp1,IOZ) - in(IOXp1,IOY,IOZp1)
02931 - in(IOX,IOYp1,IOZp1) + in(IOXp1,IOYp1,IOZp1);
02932 (*ret)(ix,iy,iz) = a1 + dz*(a4 + a6*dx + (a7 + a8*dx)*dy) + a3*dy + dx*(a2 + a5*dy);
02933 }
02934 }
02935 }
02936
02937 set_array_offsets(saved_offsets);
02938 return ret;
02939
02940 }
02941 }
02942 #undef in
02943
02944
02945
02946
02947
02948
02949
02950
02951
02952
02953
02954
02955
02956
02957
02958
02959
02960
02961
02962
02963
02964
02965
02966
02967
02968
02969
02970
02971
02972
02973
02974
02975
02976
02977
02978
02979
02980
02981
02982
02983
02984
02985
02986
02987
02988
02989
02990
02991
02992
02993
02994
02995
02996
02997
02998
02999
03000
03001
03002
03003
03004
03005
03006
03007
03008
03009
03010
03011
03012
03013
03014
03015
03016
03017
03018
03019
03020 EMData* EMData::rot_scale_conv(float ang, float delx, float dely, Util::KaiserBessel& kb, float scale_input) {
03021 int nxn, nyn, nzn;
03022 if(scale_input == 0.0f) scale_input = 1.0f;
03023
03024 float scale = 0.5f*scale_input;
03025 float sum, w;
03026 if (1 >= ny)
03027 throw ImageDimensionException("Can't rotate 1D image");
03028 if (1 < nz)
03029 throw ImageDimensionException("Volume not currently supported");
03030 nxn=nx/2;nyn=ny/2;nzn=nz/2;
03031
03032 int K = kb.get_window_size();
03033 int kbmin = -K/2;
03034 int kbmax = -kbmin;
03035 int kbc = kbmax+1;
03036 vector<int> saved_offsets = get_array_offsets();
03037 set_array_offsets(0,0,0);
03038 EMData* ret = this->copy_head();
03039 #ifdef _WIN32
03040 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03041 #else
03042 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03043 #endif //_WIN32
03044
03045 delx = restrict2(delx, nx);
03046 dely = restrict2(dely, ny);
03047
03048 int xc = nxn;
03049 int ixs = nxn%2;
03050 int yc = nyn;
03051 int iys = nyn%2;
03052
03053 int xcn = nxn/2;
03054 int ycn = nyn/2;
03055
03056 float shiftxc = xcn + delx;
03057 float shiftyc = ycn + dely;
03058
03059 float ymin = -ny/2.0f;
03060 float xmin = -nx/2.0f;
03061 float ymax = -ymin;
03062 float xmax = -xmin;
03063 if (0 == nx%2) xmax--;
03064 if (0 == ny%2) ymax--;
03065
03066 float *t = (float*)calloc(kbmax-kbmin+1, sizeof(float));
03067
03068
03069 float cang = cos(ang);
03070 float sang = sin(ang);
03071 for (int iy = 0; iy < nyn; iy++) {
03072 float y = float(iy) - shiftyc;
03073 float ycang = y*cang/scale + yc;
03074 float ysang = -y*sang/scale + xc;
03075 for (int ix = 0; ix < nxn; ix++) {
03076 float x = float(ix) - shiftxc;
03077 float xold = x*cang/scale + ysang-ixs;
03078 float yold = x*sang/scale + ycang-iys;
03079
03080 xold = restrict1(xold, nx);
03081 yold = restrict1(yold, ny);
03082
03083 int inxold = int(Util::round(xold)); int inyold = int(Util::round(yold));
03084 sum=0.0f; w=0.0f;
03085 for (int m1 =kbmin; m1 <=kbmax; m1++) t[m1-kbmin] = kb.i0win_tab(xold - inxold-m1);
03086 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 ) {
03087 for (int m2 =kbmin; m2 <=kbmax; m2++) {
03088 float qt = kb.i0win_tab(yold - inyold-m2);
03089 for (int m1 =kbmin; m1 <=kbmax; m1++) {
03090 float q = t[m1-kbmin]*qt;
03091 sum += (*this)((inxold+m1+nx)%nx,(inyold+m2+ny)%ny)*q; w+=q;
03092 }
03093 }
03094 } else {
03095 for (int m2 =kbmin; m2 <=kbmax; m2++) {
03096 float qt = kb.i0win_tab(yold - inyold-m2);
03097 for (int m1 =kbmin; m1 <=kbmax; m1++) {
03098 float q = t[m1-kbmin]*qt;
03099 sum += (*this)(inxold+m1,inyold+m2)*q; w+=q;}
03100 }
03101 }
03102 (*ret)(ix,iy)=sum/w;
03103 }
03104 }
03105 if (t) free(t);
03106 set_array_offsets(saved_offsets);
03107 return ret;
03108 }
03109
03110
03111
03112 EMData* EMData::rot_scale_conv7(float ang, float delx, float dely, Util::KaiserBessel& kb, float scale_input) {
03113 int nxn, nyn, nzn;
03114 float scale = 0.5f*scale_input;
03115 float sum, w;
03116 if (1 >= ny)
03117 throw ImageDimensionException("Can't rotate 1D image");
03118 if (1 < nz)
03119 throw ImageDimensionException("Volume not currently supported");
03120 nxn = nx/2; nyn=ny/2; nzn=nz/2;
03121
03122 int K = kb.get_window_size();
03123 int kbmin = -K/2;
03124 int kbmax = -kbmin;
03125 int kbc = kbmax+1;
03126 vector<int> saved_offsets = get_array_offsets();
03127 set_array_offsets(0,0,0);
03128 EMData* ret = this->copy_head();
03129 #ifdef _WIN32
03130 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03131 #else
03132 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03133 #endif //_WIN32
03134
03135 delx = restrict2(delx, nx);
03136 dely = restrict2(dely, ny);
03137
03138 int xc = nxn;
03139 int ixs = nxn%2;
03140 int yc = nyn;
03141 int iys = nyn%2;
03142
03143 int xcn = nxn/2;
03144 int ycn = nyn/2;
03145
03146 float shiftxc = xcn + delx;
03147 float shiftyc = ycn + dely;
03148
03149 float ymin = -ny/2.0f;
03150 float xmin = -nx/2.0f;
03151 float ymax = -ymin;
03152 float xmax = -xmin;
03153 if (0 == nx%2) xmax--;
03154 if (0 == ny%2) ymax--;
03155
03156 float *t = (float*)calloc(kbmax-kbmin+1, sizeof(float));
03157
03158
03159 float cang = cos(ang);
03160 float sang = sin(ang);
03161 for (int iy = 0; iy < nyn; iy++) {
03162 float y = float(iy) - shiftyc;
03163 float ycang = y*cang/scale + yc;
03164 float ysang = -y*sang/scale + xc;
03165 for (int ix = 0; ix < nxn; ix++) {
03166 float x = float(ix) - shiftxc;
03167 float xold = x*cang/scale + ysang-ixs;
03168 float yold = x*sang/scale + ycang-iys;
03169
03170 xold = restrict1(xold, nx);
03171 yold = restrict1(yold, ny);
03172
03173 int inxold = int(Util::round(xold)); int inyold = int(Util::round(yold));
03174 sum=0.0f; w=0.0f;
03175
03176 float tablex1 = kb.i0win_tab(xold-inxold+3);
03177 float tablex2 = kb.i0win_tab(xold-inxold+2);
03178 float tablex3 = kb.i0win_tab(xold-inxold+1);
03179 float tablex4 = kb.i0win_tab(xold-inxold);
03180 float tablex5 = kb.i0win_tab(xold-inxold-1);
03181 float tablex6 = kb.i0win_tab(xold-inxold-2);
03182 float tablex7 = kb.i0win_tab(xold-inxold-3);
03183
03184 float tabley1 = kb.i0win_tab(yold-inyold+3);
03185 float tabley2 = kb.i0win_tab(yold-inyold+2);
03186 float tabley3 = kb.i0win_tab(yold-inyold+1);
03187 float tabley4 = kb.i0win_tab(yold-inyold);
03188 float tabley5 = kb.i0win_tab(yold-inyold-1);
03189 float tabley6 = kb.i0win_tab(yold-inyold-2);
03190 float tabley7 = kb.i0win_tab(yold-inyold-3);
03191
03192 int x1, x2, x3, x4, x5, x6, x7, y1, y2, y3, y4, y5, y6, y7;
03193
03194 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 ) {
03195 x1 = (inxold-3+nx)%nx;
03196 x2 = (inxold-2+nx)%nx;
03197 x3 = (inxold-1+nx)%nx;
03198 x4 = (inxold +nx)%nx;
03199 x5 = (inxold+1+nx)%nx;
03200 x6 = (inxold+2+nx)%nx;
03201 x7 = (inxold+3+nx)%nx;
03202
03203 y1 = (inyold-3+ny)%ny;
03204 y2 = (inyold-2+ny)%ny;
03205 y3 = (inyold-1+ny)%ny;
03206 y4 = (inyold +ny)%ny;
03207 y5 = (inyold+1+ny)%ny;
03208 y6 = (inyold+2+ny)%ny;
03209 y7 = (inyold+3+ny)%ny;
03210 } else {
03211 x1 = inxold-3;
03212 x2 = inxold-2;
03213 x3 = inxold-1;
03214 x4 = inxold;
03215 x5 = inxold+1;
03216 x6 = inxold+2;
03217 x7 = inxold+3;
03218
03219 y1 = inyold-3;
03220 y2 = inyold-2;
03221 y3 = inyold-1;
03222 y4 = inyold;
03223 y5 = inyold+1;
03224 y6 = inyold+2;
03225 y7 = inyold+3;
03226 }
03227 sum = ( (*this)(x1,y1)*tablex1 + (*this)(x2,y1)*tablex2 + (*this)(x3,y1)*tablex3 +
03228 (*this)(x4,y1)*tablex4 + (*this)(x5,y1)*tablex5 + (*this)(x6,y1)*tablex6 +
03229 (*this)(x7,y1)*tablex7 ) * tabley1 +
03230 ( (*this)(x1,y2)*tablex1 + (*this)(x2,y2)*tablex2 + (*this)(x3,y2)*tablex3 +
03231 (*this)(x4,y2)*tablex4 + (*this)(x5,y2)*tablex5 + (*this)(x6,y2)*tablex6 +
03232 (*this)(x7,y2)*tablex7 ) * tabley2 +
03233 ( (*this)(x1,y3)*tablex1 + (*this)(x2,y3)*tablex2 + (*this)(x3,y3)*tablex3 +
03234 (*this)(x4,y3)*tablex4 + (*this)(x5,y3)*tablex5 + (*this)(x6,y3)*tablex6 +
03235 (*this)(x7,y3)*tablex7 ) * tabley3 +
03236 ( (*this)(x1,y4)*tablex1 + (*this)(x2,y4)*tablex2 + (*this)(x3,y4)*tablex3 +
03237 (*this)(x4,y4)*tablex4 + (*this)(x5,y4)*tablex5 + (*this)(x6,y4)*tablex6 +
03238 (*this)(x7,y4)*tablex7 ) * tabley4 +
03239 ( (*this)(x1,y5)*tablex1 + (*this)(x2,y5)*tablex2 + (*this)(x3,y5)*tablex3 +
03240 (*this)(x4,y5)*tablex4 + (*this)(x5,y5)*tablex5 + (*this)(x6,y5)*tablex6 +
03241 (*this)(x7,y5)*tablex7 ) * tabley5 +
03242 ( (*this)(x1,y6)*tablex1 + (*this)(x2,y6)*tablex2 + (*this)(x3,y6)*tablex3 +
03243 (*this)(x4,y6)*tablex4 + (*this)(x5,y6)*tablex5 + (*this)(x6,y6)*tablex6 +
03244 (*this)(x7,y6)*tablex7 ) * tabley6 +
03245 ( (*this)(x1,y7)*tablex1 + (*this)(x2,y7)*tablex2 + (*this)(x3,y7)*tablex3 +
03246 (*this)(x4,y7)*tablex4 + (*this)(x5,y7)*tablex5 + (*this)(x6,y7)*tablex6 +
03247 (*this)(x7,y7)*tablex7 ) * tabley7;
03248
03249 w = (tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7) *
03250 (tabley1+tabley2+tabley3+tabley4+tabley5+tabley6+tabley7);
03251
03252 (*ret)(ix,iy)=sum/w;
03253 }
03254 }
03255 if (t) free(t);
03256 set_array_offsets(saved_offsets);
03257 return ret;
03258 }
03259
03260 EMData* EMData::downsample(Util::sincBlackman& kb, float scale) {
03261
03262
03263
03264
03265
03266 int nxn, nyn, nzn;
03267 nxn = (int)(nx*scale); nyn = (int)(ny*scale); nzn = (int)(nz*scale);
03268
03269 vector<int> saved_offsets = get_array_offsets();
03270 set_array_offsets(0,0,0);
03271 EMData* ret = this->copy_head();
03272 #ifdef _WIN32
03273 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03274 #else
03275 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03276 #endif //_WIN32
03277 ret->to_zero();
03278
03279
03280 for (int iy =0; iy < nyn; iy++) {
03281 float y = float(iy)/scale;
03282 for (int ix = 0; ix < nxn; ix++) {
03283 float x = float(ix)/scale;
03284 (*ret)(ix,iy) = this->get_pixel_filtered(x, y, 1.0f, kb);
03285 }
03286 }
03287 set_array_offsets(saved_offsets);
03288 return ret;
03289 }
03290
03291
03292 EMData* EMData::rot_scale_conv_new(float ang, float delx, float dely, Util::KaiserBessel& kb, float scale_input) {
03293
03294 if (scale_input == 0.0f) scale_input = 1.0f;
03295 float scale = 0.5f*scale_input;
03296
03297 if (1 >= ny)
03298 throw ImageDimensionException("Can't rotate 1D image");
03299 if (1 < nz)
03300 throw ImageDimensionException("Use rot_scale_conv_new_3D for volumes");
03301 int nxn = nx/2; int nyn = ny/2; int nzn = nz/2;
03302
03303 vector<int> saved_offsets = get_array_offsets();
03304 set_array_offsets(0,0,0);
03305 EMData* ret = this->copy_head();
03306 #ifdef _WIN32
03307 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03308 #else
03309 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03310 #endif //_WIN32
03311
03312 delx = restrict2(delx, nx);
03313 dely = restrict2(dely, ny);
03314
03315 int xc = nxn;
03316 int ixs = nxn%2;
03317 int yc = nyn;
03318 int iys = nyn%2;
03319
03320 int xcn = nxn/2;
03321 int ycn = nyn/2;
03322
03323 float shiftxc = xcn + delx;
03324 float shiftyc = ycn + dely;
03325
03326 float ymin = -ny/2.0f;
03327 float xmin = -nx/2.0f;
03328 float ymax = -ymin;
03329 float xmax = -xmin;
03330 if (0 == nx%2) xmax--;
03331 if (0 == ny%2) ymax--;
03332
03333 float* data = this->get_data();
03334
03335 float cang = cos(ang);
03336 float sang = sin(ang);
03337 for (int iy = 0; iy < nyn; iy++) {
03338 float y = float(iy) - shiftyc;
03339 float ycang = y*cang/scale + yc;
03340 float ysang = -y*sang/scale + xc;
03341 for (int ix = 0; ix < nxn; ix++) {
03342 float x = float(ix) - shiftxc;
03343 float xold = x*cang/scale + ysang-ixs;
03344 float yold = x*sang/scale + ycang-iys;
03345
03346 (*ret)(ix,iy) = Util::get_pixel_conv_new(nx, ny, 1, xold, yold, 1, data, kb);
03347 }
03348 }
03349 set_array_offsets(saved_offsets);
03350 return ret;
03351 }
03352
03353 EMData* EMData::rot_scale_conv_new_3D(float phi, float theta, float psi, float delx, float dely, float delz, Util::KaiserBessel& kb, float scale_input, bool wrap) {
03354
03355 if (scale_input == 0.0f) scale_input = 1.0f;
03356 float scale = 0.5f*scale_input;
03357
03358 if (1 >= ny)
03359 throw ImageDimensionException("Can't rotate 1D image");
03360 int nxn = nx/2; int nyn = ny/2; int nzn = nz/2;
03361
03362 vector<int> saved_offsets = get_array_offsets();
03363 set_array_offsets(0,0,0);
03364 EMData* ret = this->copy_head();
03365 #ifdef _WIN32
03366 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03367 #else
03368 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03369 #endif //_WIN32
03370
03371 if(wrap){
03372 delx = restrict2(delx, nx);
03373 dely = restrict2(dely, ny);
03374 delz = restrict2(delz, nz);
03375 }
03376
03377 int xc = nxn;
03378 int ixs = nxn%2;
03379 int yc = nyn;
03380 int iys = nyn%2;
03381 int zc = nzn;
03382 int izs = nzn%2;
03383
03384 int xcn = nxn/2;
03385 int ycn = nyn/2;
03386 int zcn = nzn/2;
03387
03388 float shiftxc = xcn + delx;
03389 float shiftyc = ycn + dely;
03390 float shiftzc = zcn + delz;
03391
03392 float zmin = -nz/2.0f;
03393 float ymin = -ny/2.0f;
03394 float xmin = -nx/2.0f;
03395 float zmax = -zmin;
03396 float ymax = -ymin;
03397 float xmax = -xmin;
03398 if (0 == nx%2) xmax--;
03399 if (0 == ny%2) ymax--;
03400 if (0 == nz%2) zmax--;
03401
03402 float* data = this->get_data();
03403
03404 float cf = cos(phi); float sf = sin(phi);
03405 float ct = cos(theta); float st = sin(theta);
03406 float cp = cos(psi); float sp = sin(psi);
03407
03408 float a11 = cp*ct*cf-sp*sf; float a12 = cp*ct*sf+sp*cf; float a13 = -cp*st;
03409 float a21 = -sp*ct*cf-cp*sf; float a22 = -sp*ct*sf+cp*cf; float a23 = sp*st;
03410 float a31 = st*cf; float a32 = st*sf; float a33 = ct;
03411 for (int iz = 0; iz < nzn; iz++) {
03412 float z = (float(iz) - shiftzc)/scale;
03413 float zco1 = a31*z+xc;
03414 float zco2 = a32*z+yc;
03415 float zco3 = a33*z+zc;
03416 for (int iy = 0; iy < nyn; iy++) {
03417 float y = (float(iy) - shiftyc)/scale;
03418 float yco1 = zco1+a21*y;
03419 float yco2 = zco2+a22*y;
03420 float yco3 = zco3+a23*y;
03421 for (int ix = 0; ix < nxn; ix++) {
03422 float x = (float(ix) - shiftxc)/scale;
03423 float xold = yco1+a11*x-ixs;
03424 float yold = yco2+a12*x-iys;
03425 float zold = yco3+a13*x-izs;
03426 if(!wrap && (xold<0.0 || xold>nx-1 || yold<0.0 || yold>ny-1 || zold<0.0 || zold>nz-1))
03427 (*ret)(ix,iy,iz) = 0.0;
03428 else
03429 (*ret)(ix,iy,iz) = Util::get_pixel_conv_new(nx, ny, nz, xold, yold, zold, data, kb);
03430 }
03431 }
03432 }
03433 set_array_offsets(saved_offsets);
03434 return ret;
03435 }
03436
03437 EMData* EMData::rot_scale_conv_new_background(float ang, float delx, float dely, Util::KaiserBessel& kb, float scale_input) {
03438
03439 int nxn, nyn, nzn;
03440
03441 if (scale_input == 0.0f) scale_input = 1.0f;
03442 float scale = 0.5f*scale_input;
03443
03444 if (1 >= ny)
03445 throw ImageDimensionException("Can't rotate 1D image");
03446 if (1 < nz)
03447 throw ImageDimensionException("Use rot_scale_conv_new_background_3D for volumes");
03448 nxn = nx/2; nyn = ny/2; nzn = nz/2;
03449
03450 vector<int> saved_offsets = get_array_offsets();
03451 set_array_offsets(0,0,0);
03452 EMData* ret = this->copy_head();
03453 #ifdef _WIN32
03454 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03455 #else
03456 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03457 #endif //_WIN32
03458
03459 delx = restrict2(delx, nx);
03460 dely = restrict2(dely, ny);
03461
03462 int xc = nxn;
03463 int ixs = nxn%2;
03464 int yc = nyn;
03465 int iys = nyn%2;
03466
03467 int xcn = nxn/2;
03468 int ycn = nyn/2;
03469
03470 float shiftxc = xcn + delx;
03471 float shiftyc = ycn + dely;
03472
03473 float ymin = -ny/2.0f;
03474 float xmin = -nx/2.0f;
03475 float ymax = -ymin;
03476 float xmax = -xmin;
03477 if (0 == nx%2) xmax--;
03478 if (0 == ny%2) ymax--;
03479
03480 float* data = this->get_data();
03481
03482
03483 float cang = cos(ang);
03484 float sang = sin(ang);
03485 for (int iy = 0; iy < nyn; iy++) {
03486 float y = float(iy) - shiftyc;
03487 float ycang = y*cang/scale + yc;
03488 float ysang = -y*sang/scale + xc;
03489 for (int ix = 0; ix < nxn; ix++) {
03490 float x = float(ix) - shiftxc;
03491 float xold = x*cang/scale + ysang-ixs;
03492 float yold = x*sang/scale + ycang-iys;
03493
03494 (*ret)(ix,iy) = Util::get_pixel_conv_new_background(nx, ny, 1, xold, yold, 1, data, kb, ix, iy);
03495 }
03496 }
03497 set_array_offsets(saved_offsets);
03498 return ret;
03499 }
03500
03501 EMData* EMData::rot_scale_conv_new_background_3D(float phi, float theta, float psi, float delx, float dely, float delz, Util::KaiserBessel& kb, float scale_input, bool wrap) {
03502
03503 if (scale_input == 0.0f) scale_input = 1.0f;
03504 float scale = 0.5f*scale_input;
03505
03506 if (1 >= ny)
03507 throw ImageDimensionException("Can't rotate 1D image");
03508 int nxn = nx/2; int nyn = ny/2; int nzn = nz/2;
03509
03510 vector<int> saved_offsets = get_array_offsets();
03511 set_array_offsets(0,0,0);
03512 EMData* ret = this->copy_head();
03513 #ifdef _WIN32
03514 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03515 #else
03516 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03517 #endif //_WIN32
03518
03519 if (wrap){
03520 delx = restrict2(delx, nx);
03521 dely = restrict2(dely, ny);
03522 delz = restrict2(delz, nz);
03523 }
03524
03525 int xc = nxn;
03526 int ixs = nxn%2;
03527 int yc = nyn;
03528 int iys = nyn%2;
03529 int zc = nzn;
03530 int izs = nzn%2;
03531
03532 int xcn = nxn/2;
03533 int ycn = nyn/2;
03534 int zcn = nzn/2;
03535
03536 float shiftxc = xcn + delx;
03537 float shiftyc = ycn + dely;
03538 float shiftzc = zcn + delz;
03539
03540 float zmin = -nz/2.0f;
03541 float ymin = -ny/2.0f;
03542 float xmin = -nx/2.0f;
03543 float zmax = -zmin;
03544 float ymax = -ymin;
03545 float xmax = -xmin;
03546 if (0 == nx%2) xmax--;
03547 if (0 == ny%2) ymax--;
03548 if (0 == nz%2) zmax--;
03549
03550 float* data = this->get_data();
03551
03552 float cf = cos(phi); float sf = sin(phi);
03553 float ct = cos(theta); float st = sin(theta);
03554 float cp = cos(psi); float sp = sin(psi);
03555
03556 float a11 = cp*ct*cf-sp*sf; float a12 = cp*ct*sf+sp*cf; float a13 = -cp*st;
03557 float a21 = -sp*ct*cf-cp*sf; float a22 = -sp*ct*sf+cp*cf; float a23 = sp*st;
03558 float a31 = st*cf; float a32 = st*sf; float a33 = ct;
03559 for (int iz = 0; iz < nzn; iz++) {
03560 float z = (float(iz) - shiftzc)/scale;
03561 float zco1 = a31*z+xc;
03562 float zco2 = a32*z+yc;
03563 float zco3 = a33*z+zc;
03564 for (int iy = 0; iy < nyn; iy++) {
03565 float y = (float(iy) - shiftyc)/scale;
03566 float yco1 = zco1+a21*y;
03567 float yco2 = zco2+a22*y;
03568 float yco3 = zco3+a23*y;
03569 for (int ix = 0; ix < nxn; ix++) {
03570 float x = (float(ix) - shiftxc)/scale;
03571 float xold = yco1+a11*x-ixs;
03572 float yold = yco2+a12*x-iys;
03573 float zold = yco3+a13*x-izs;
03574 if(!wrap && (xold<0.0 || xold>nx-1 || yold<0.0 || yold>ny-1 || zold<0.0 || zold>nz-1))
03575 (*ret)(ix,iy,iz) = 0.0;
03576 else
03577 (*ret)(ix,iy,iz) = Util::get_pixel_conv_new_background(nx, ny, nz, xold, yold, zold, data, kb, ix, iy);
03578 }
03579 }
03580 }
03581 set_array_offsets(saved_offsets);
03582 return ret;
03583 }
03584
03585
03586 float EMData::get_pixel_conv(float delx, float dely, float delz, Util::KaiserBessel& kb) {
03587
03588
03589 int K = kb.get_window_size();
03590 int kbmin = -K/2;
03591 int kbmax = -kbmin;
03592 int kbc = kbmax+1;
03593
03594 float pixel =0.0f;
03595 float w=0.0f;
03596
03597 delx = restrict2(delx, nx);
03598 int inxold = int(Util::round(delx));
03599 if(ny<2) {
03600 if(inxold <= kbc || inxold >=nx-kbc-2 ) {
03601
03602 for (int m1 =kbmin; m1 <=kbmax; m1++) {
03603 float q = kb.i0win_tab(delx - inxold-m1);
03604 pixel += (*this)((inxold+m1+nx)%nx)*q; w+=q;
03605 }
03606 } else {
03607 for (int m1 =kbmin; m1 <=kbmax; m1++) {
03608 float q = kb.i0win_tab(delx - inxold-m1);
03609 pixel += (*this)(inxold+m1)*q; w+=q;
03610 }
03611 }
03612
03613 } else if(nz<2) {
03614 dely = restrict2(dely, ny);
03615 int inyold = int(Util::round(dely));
03616 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 ) {
03617
03618 for (int m2 =kbmin; m2 <=kbmax; m2++){ for (int m1 =kbmin; m1 <=kbmax; m1++) {
03619 float q = kb.i0win_tab(delx - inxold-m1)*kb.i0win_tab(dely - inyold-m2);
03620 pixel += (*this)((inxold+m1+nx)%nx,(inyold+m2+ny)%ny)*q; w+=q;}
03621 }
03622 } else {
03623 for (int m2 =kbmin; m2 <=kbmax; m2++){ for (int m1 =kbmin; m1 <=kbmax; m1++) {
03624 float q = kb.i0win_tab(delx - inxold-m1)*kb.i0win_tab(dely - inyold-m2);
03625 pixel += (*this)(inxold+m1,inyold+m2)*q; w+=q;}
03626 }
03627 }
03628 } else {
03629 dely = restrict2(dely, ny);
03630 int inyold = int(Util::round(dely));
03631 delz = restrict2(delz, nz);
03632 int inzold = int(Util::round(delz));
03633
03634 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 || inzold <= kbc || inzold >=nz-kbc-2 ) {
03635
03636 for (int m3 =kbmin; m3 <=kbmax; m3++){ for (int m2 =kbmin; m2 <=kbmax; m2++){ for (int m1 =kbmin; m1 <=kbmax; m1++) {
03637 float q = kb.i0win_tab(delx - inxold-m1)*kb.i0win_tab(dely - inyold-m2)*kb.i0win_tab(delz - inzold-m3);
03638
03639 pixel += (*this)((inxold+m1+nx)%nx,(inyold+m2+ny)%ny,(inzold+m3+nz)%nz)*q ;w+=q;}}
03640 }
03641 } else {
03642 for (int m3 =kbmin; m3 <=kbmax; m3++){ for (int m2 =kbmin; m2 <=kbmax; m2++){ for (int m1 =kbmin; m1 <=kbmax; m1++) {
03643 float q = kb.i0win_tab(delx - inxold-m1)*kb.i0win_tab(dely - inyold-m2)*kb.i0win_tab(delz - inzold-m3);
03644
03645 pixel += (*this)(inxold+m1,inyold+m2,inzold+m3)*q; w+=q;}}
03646 }
03647 }
03648 }
03649 return pixel/w;
03650 }
03651
03652
03653 float EMData::get_pixel_filtered(float delx, float dely, float, Util::sincBlackman& kb) {
03654
03655
03656 int K = kb.get_sB_size();
03657 int kbmin = -K/2;
03658 int kbmax = -kbmin;
03659 int kbc = kbmax+1;
03660
03661 float pixel =0.0f;
03662 float w=0.0f;
03663
03664
03665 int inxold = int(Util::round(delx));
03666
03667
03668
03669
03670
03671
03672
03673
03674
03675
03676
03677
03678
03679
03680
03681
03682 int inyold = int(Util::round(dely));
03683 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 ) {
03684
03685 for (int m2 =kbmin; m2 <=kbmax; m2++){
03686 float t = kb.sBwin_tab(dely - inyold-m2);
03687 for (int m1 =kbmin; m1 <=kbmax; m1++) {
03688 float q = kb.sBwin_tab(delx - inxold-m1)*t;
03689 pixel += (*this)((inxold+m1+nx)%nx,(inyold+m2+ny)%ny)*q;
03690 w += q;
03691 }
03692 }
03693 } else {
03694 for (int m2 =kbmin; m2 <=kbmax; m2++){
03695 float t = kb.sBwin_tab(dely - inyold-m2);
03696 for (int m1 =kbmin; m1 <=kbmax; m1++) {
03697 float q = kb.sBwin_tab(delx - inxold-m1)*t;
03698 pixel += (*this)(inxold+m1,inyold+m2)*q;
03699 w += q;
03700 }
03701 }
03702 }
03703
03704
03705
03706
03707
03708
03709
03710
03711
03712
03713
03714
03715
03716
03717
03718
03719
03720
03721
03722
03723
03724 return pixel/w;
03725 }
03726
03727
03728
03729
03730 float EMData::get_pixel_conv7(float delx, float dely, float delz, Util::KaiserBessel& kb) {
03731
03732
03733 float *image=(this->get_data());
03734 int nx = this->get_xsize();
03735 int ny = this->get_ysize();
03736 int nz = this->get_zsize();
03737
03738 float result;
03739
03740 result = Util::get_pixel_conv_new(nx,ny,nz,delx,dely,delz,image,kb);
03741 return result;
03742 }
03743
03744 float EMData::getconvpt2d_kbi0(float x, float y, Util::KaiserBessel::kbi0_win win, int size) {
03745 const int nxhalf = nx/2;
03746 const int nyhalf = ny/2;
03747 const int bd = size/2;
03748 float* wxarr = new float[size];
03749 float* wyarr = new float[size];
03750 float* wx = wxarr + bd;
03751 float* wy = wyarr + bd;
03752 int ixc = int(x + 0.5f*Util::sgn(x));
03753 int iyc = int(y + 0.5f*Util::sgn(y));
03754 if (abs(ixc) > nxhalf)
03755 throw InvalidValueException(ixc, "getconv: X value out of range");
03756 if (abs(iyc) > nyhalf)
03757 throw InvalidValueException(ixc, "getconv: Y value out of range");
03758 for (int i = -bd; i <= bd; i++) {
03759 int iyp = iyc + i;
03760 wy[i] = win(y - iyp);
03761 int ixp = ixc + i;
03762 wx[i] = win(x - ixp);
03763 }
03764 vector<int> saved_offsets = get_array_offsets();
03765 set_array_offsets(-nxhalf, -nyhalf);
03766 float conv = 0.f, wsum = 0.f;
03767 for (int iy = -bd; iy <= bd; iy++) {
03768 int iyp = iyc + iy;
03769 for (int ix = -bd; ix <= bd; ix++) {
03770 int ixp = ixc + ix;
03771 float wg = wx[ix]*wy[iy];
03772 conv += (*this)(ixp,iyp)*wg;
03773 wsum += wg;
03774 }
03775 }
03776 set_array_offsets(saved_offsets);
03777 delete [] wxarr;
03778 delete [] wyarr;
03779
03780 return conv;
03781 }
03782
03783 std::complex<float> EMData::extractpoint(float nuxnew, float nuynew, Util::KaiserBessel& kb) {
03784 if (2 != get_ndim())
03785 throw ImageDimensionException("extractpoint needs a 2-D image.");
03786 if (!is_complex())
03787 throw ImageFormatException("extractpoint requires a fourier image");
03788 int nxreal = nx - 2;
03789 if (nxreal != ny)
03790 throw ImageDimensionException("extractpoint requires ny == nx");
03791 int nhalf = nxreal/2;
03792 int kbsize = kb.get_window_size();
03793 int kbmin = -kbsize/2;
03794 int kbmax = -kbmin;
03795 bool flip = (nuxnew < 0.f);
03796 if (flip) {
03797 nuxnew *= -1;
03798 nuynew *= -1;
03799 }
03800
03801
03802
03803 int ixn = int(Util::round(nuxnew));
03804 int iyn = int(Util::round(nuynew));
03805
03806 float* wy0 = new float[kbmax - kbmin + 1];
03807 float* wy = wy0 - kbmin;
03808 float* wx0 = new float[kbmax - kbmin + 1];
03809 float* wx = wx0 - kbmin;
03810 for (int i = kbmin; i <= kbmax; i++) {
03811 int iyp = iyn + i;
03812 wy[i] = kb.i0win_tab(nuynew - iyp);
03813 int ixp = ixn + i;
03814 wx[i] = kb.i0win_tab(nuxnew - ixp);
03815 }
03816
03817 int iymin = 0;
03818 for (int iy = kbmin; iy <= -1; iy++) {
03819 if (wy[iy] != 0.f) {
03820 iymin = iy;
03821 break;
03822 }
03823 }
03824 int iymax = 0;
03825 for (int iy = kbmax; iy >= 1; iy--) {
03826 if (wy[iy] != 0.f) {
03827 iymax = iy;
03828 break;
03829 }
03830 }
03831 int ixmin = 0;
03832 for (int ix = kbmin; ix <= -1; ix++) {
03833 if (wx[ix] != 0.f) {
03834 ixmin = ix;
03835 break;
03836 }
03837 }
03838 int ixmax = 0;
03839 for (int ix = kbmax; ix >= 1; ix--) {
03840 if (wx[ix] != 0.f) {
03841 ixmax = ix;
03842 break;
03843 }
03844 }
03845 float wsum = 0.0f;
03846 for (int iy = iymin; iy <= iymax; iy++)
03847 for (int ix = ixmin; ix <= ixmax; ix++)
03848 wsum += wx[ix]*wy[iy];
03849 std::complex<float> result(0.f,0.f);
03850 if ((ixn >= -kbmin) && (ixn <= nhalf-1-kbmax) && (iyn >= -nhalf-kbmin) && (iyn <= nhalf-1-kbmax)) {
03851
03852 for (int iy = iymin; iy <= iymax; iy++) {
03853 int iyp = iyn + iy;
03854 for (int ix = ixmin; ix <= ixmax; ix++) {
03855 int ixp = ixn + ix;
03856 float w = wx[ix]*wy[iy];
03857 std::complex<float> val = cmplx(ixp,iyp);
03858 result += val*w;
03859 }
03860 }
03861 } else {
03862
03863 for (int iy = iymin; iy <= iymax; iy++) {
03864 int iyp = iyn + iy;
03865 for (int ix = ixmin; ix <= ixmax; ix++) {
03866 int ixp = ixn + ix;
03867 bool mirror = false;
03868 int ixt= ixp, iyt= iyp;
03869 if (ixt < 0) {
03870 ixt = -ixt;
03871 iyt = -iyt;
03872 mirror = !mirror;
03873 }
03874 if (ixt > nhalf) {
03875 ixt = nxreal - ixt;
03876 iyt = -iyt;
03877 mirror = !mirror;
03878 }
03879 if (iyt > nhalf-1) iyt -= nxreal;
03880 if (iyt < -nhalf) iyt += nxreal;
03881 float w = wx[ix]*wy[iy];
03882 std::complex<float> val = this->cmplx(ixt,iyt);
03883 if (mirror) result += conj(val)*w;
03884 else result += val*w;
03885 }
03886 }
03887 }
03888 if (flip) result = conj(result)/wsum;
03889 else result /= wsum;
03890 delete [] wx0;
03891 delete [] wy0;
03892 return result;
03893 }
03894
03895 EMData* EMData::extractline(Util::KaiserBessel& kb, float nuxnew, float nuynew)
03896 {
03897 if (!is_complex())
03898 throw ImageFormatException("extractline requires a fourier image");
03899 if (nx%2 != 0)
03900 throw ImageDimensionException("extractline requires nx to be even");
03901 int nxreal = nx - 2;
03902 if (nxreal != ny)
03903 throw ImageDimensionException("extractline requires ny == nx");
03904
03905 EMData* res = new EMData();
03906 res->set_size(nx,1,1);
03907 res->to_zero();
03908 res->set_complex(true);
03909 res->set_fftodd(false);
03910 res->set_fftpad(true);
03911 res->set_ri(true);
03912
03913 int n = nxreal;
03914 int nhalf = n/2;
03915 vector<int> saved_offsets = get_array_offsets();
03916 set_array_offsets(0,-nhalf,-nhalf);
03917
03918
03919 int kbsize = kb.get_window_size();
03920 int kbmin = -kbsize/2;
03921 int kbmax = -kbmin;
03922 float* wy0 = new float[kbmax - kbmin + 1];
03923 float* wy = wy0 - kbmin;
03924 float* wx0 = new float[kbmax - kbmin + 1];
03925 float* wx = wx0 - kbmin;
03926
03927 int count = 0;
03928 float wsum = 0.f;
03929 bool flip = (nuxnew < 0.f);
03930
03931 for (int jx = 0; jx <= nhalf; jx++) {
03932 float xnew = jx*nuxnew, ynew = jx*nuynew;
03933 count++;
03934 std::complex<float> btq(0.f,0.f);
03935 if (flip) {
03936 xnew = -xnew;
03937 ynew = -ynew;
03938 }
03939 int ixn = int(Util::round(xnew));
03940 int iyn = int(Util::round(ynew));
03941
03942 for (int i=kbmin; i <= kbmax; i++) {
03943 int iyp = iyn + i;
03944 wy[i] = kb.i0win_tab(ynew - iyp);
03945 int ixp = ixn + i;
03946 wx[i] = kb.i0win_tab(xnew - ixp);
03947 }
03948
03949
03950 int lnby = 0;
03951 for (int iy = kbmin; iy <= -1; iy++) {
03952 if (wy[iy] != 0.f) {
03953 lnby = iy;
03954 break;
03955 }
03956 }
03957 int lney = 0;
03958 for (int iy = kbmax; iy >= 1; iy--) {
03959 if (wy[iy] != 0.f) {
03960 lney = iy;
03961 break;
03962 }
03963 }
03964 int lnbx = 0;
03965 for (int ix = kbmin; ix <= -1; ix++) {
03966 if (wx[ix] != 0.f) {
03967 lnbx = ix;
03968 break;
03969 }
03970 }
03971 int lnex = 0;
03972 for (int ix = kbmax; ix >= 1; ix--) {
03973 if (wx[ix] != 0.f) {
03974 lnex = ix;
03975 break;
03976 }
03977 }
03978 if (ixn >= -kbmin && ixn <= nhalf-1-kbmax
03979 && iyn >= -nhalf-kbmin && iyn <= nhalf-1-kbmax) {
03980
03981 for (int ly=lnby; ly<=lney; ly++) {
03982 int iyp = iyn + ly;
03983 for (int lx=lnbx; lx<=lnex; lx++) {
03984 int ixp = ixn + lx;
03985 float wg = wx[lx]*wy[ly];
03986 btq += cmplx(ixp,iyp)*wg;
03987 wsum += wg;
03988 }
03989 }
03990 } else {
03991
03992 for (int ly=lnby; ly<=lney; ly++) {
03993 int iyp = iyn + ly;
03994 for (int lx=lnbx; lx<=lnex; lx++) {
03995 int ixp = ixn + lx;
03996 float wg = wx[lx]*wy[ly];
03997 bool mirror = false;
03998 int ixt(ixp), iyt(iyp);
03999 if (ixt > nhalf || ixt < -nhalf) {
04000 ixt = Util::sgn(ixt)*(n - abs(ixt));
04001 iyt = -iyt;
04002 mirror = !mirror;
04003 }
04004 if (iyt >= nhalf || iyt < -nhalf) {
04005 if (ixt != 0) {
04006 ixt = -ixt;
04007 iyt = Util::sgn(iyt)*(n - abs(iyt));
04008 mirror = !mirror;
04009 } else {
04010 iyt -= n*Util::sgn(iyt);
04011 }
04012 }
04013 if (ixt < 0) {
04014 ixt = -ixt;
04015 iyt = -iyt;
04016 mirror = !mirror;
04017 }
04018 if (iyt == nhalf) iyt = -nhalf;
04019 if (mirror) btq += conj(cmplx(ixt,iyt))*wg;
04020 else btq += cmplx(ixt,iyt)*wg;
04021 wsum += wg;
04022 }
04023 }
04024 }
04025 if (flip) res->cmplx(jx) = conj(btq);
04026 else res->cmplx(jx) = btq;
04027 }
04028 for (int jx = 0; jx <= nhalf; jx++) res->cmplx(jx) *= count/wsum;
04029
04030 delete[] wx0; delete[] wy0;
04031 set_array_offsets(saved_offsets);
04032 res->set_array_offsets(0,0,0);
04033 return res;
04034 }
04035
04036
04038 inline void swapx(float* a, float* b, float* temp, size_t nbytes) {
04039 memcpy(temp, a, nbytes);
04040 memcpy(a, b, nbytes);
04041 memcpy(b, temp, nbytes);
04042 }
04043
04044 void EMData::fft_shuffle() {
04045 if (!is_complex())
04046 throw ImageFormatException("fft_shuffle requires a fourier image");
04047 vector<int> offsets = get_array_offsets();
04048 set_array_offsets();
04049 EMData& self = *this;
04050 int nyhalf = ny/2;
04051 int nzhalf = nz/2;
04052 int nbytes = nx*sizeof(float);
04053 float* temp = new float[nx];
04054 for (int iz=0; iz < nz; iz++)
04055 for (int iy=0; iy < nyhalf; iy++)
04056 swapx(&self(0,iy,iz),&self(0,iy+nyhalf,iz),temp,nbytes);
04057 if (nz > 1) {
04058 for (int iy=0; iy < ny; iy++)
04059 for (int iz=0; iz < nzhalf; iz++)
04060 swapx(&self(0,iy,iz),&self(0,iy,iz+nzhalf),temp,nbytes);
04061 }
04062 set_shuffled(!is_shuffled());
04063 set_array_offsets(offsets);
04064 update();
04065 delete[] temp;
04066 }
04067
04068 void EMData::pad_corner(float *pad_image) {
04069 size_t nbytes = nx*sizeof(float);
04070 for (int iy=0; iy<ny; iy++)
04071 memcpy(&(*this)(0,iy), pad_image+3+(iy+3)*nx, nbytes);
04072 }
04073
04074 void EMData::shuffle_pad_corner(float *pad_image) {
04075 int nyhalf = ny/2;
04076 size_t nbytes = nx*sizeof(float);
04077 for (int iy = 0; iy < nyhalf; iy++)
04078 memcpy(&(*this)(0,iy), pad_image+6+(iy+nyhalf+3)*nx, nbytes);
04079 for (int iy = nyhalf; iy < ny; iy++)
04080 memcpy(&(*this)(0,iy), pad_image+6+(iy-nyhalf+3)*nx, nbytes);
04081 }
04082
04083 #define QUADPI 3.141592653589793238462643383279502884197
04084 #define DGR_TO_RAD QUADPI/180
04085
04086
04087
04088
04089
04090
04091
04092
04093
04094
04095
04096
04097
04098
04099
04100
04101
04102
04103
04104
04105
04106
04107
04108
04109
04110
04111
04112
04113
04114
04115
04116
04117
04118
04119
04120
04121
04122
04123
04124
04125
04126
04127
04128
04129
04130
04131
04132
04133
04134
04135
04136
04137
04138
04139 EMData* EMData::fouriergridrot2d(float ang, float scale, Util::KaiserBessel& kb) {
04140 if (2 != get_ndim())
04141 throw ImageDimensionException("fouriergridrot2d needs a 2-D image.");
04142 if (!is_complex())
04143 throw ImageFormatException("fouriergridrot2d requires a fourier image");
04144 int nxreal = nx - 2 + int(is_fftodd());
04145 if (nxreal != ny)
04146 throw ImageDimensionException("fouriergridrot2d requires ny == nx(real)");
04147 if (0 != nxreal%2)
04148 throw ImageDimensionException("fouriergridrot2d needs an even image.");
04149 if (scale == 0.0f) scale = 1.0f;
04150 int nxhalf = nxreal/2;
04151 int nyhalf = ny/2;
04152 float cir = (float)((nxhalf-1)*(nxhalf-1));
04153
04154 if (!is_shuffled()) fft_shuffle();
04155
04156 EMData* result = copy_head();
04157 set_array_offsets(0,-nyhalf);
04158 result->set_array_offsets(0,-nyhalf);
04159
04160
04161
04162 ang = ang*(float)DGR_TO_RAD;
04163 float cang = cos(ang);
04164 float sang = sin(ang);
04165 for (int iy = -nyhalf; iy < nyhalf; iy++) {
04166 float ycang = iy*cang;
04167 float ysang = iy*sang;
04168 for (int ix = 0; ix <= nxhalf; ix++) {
04169 float nuxold = (ix*cang - ysang)*scale;
04170 float nuyold = (ix*sang + ycang)*scale;
04171 if (nuxold*nuxold+nuyold*nuyold<cir) result->cmplx(ix,iy) = Util::extractpoint2(nx, ny, nuxold, nuyold, this, kb);
04172
04173 }
04174 }
04175 result->set_array_offsets();
04176 result->fft_shuffle();
04177 result->update();
04178 set_array_offsets();
04179 fft_shuffle();
04180 return result;
04181 }
04182
04183 EMData* EMData::fouriergridrot_shift2d(float ang, float sx, float sy, Util::KaiserBessel& kb) {
04184 if (2 != get_ndim())
04185 throw ImageDimensionException("fouriergridrot_shift2d needs a 2-D image.");
04186 if (!is_complex())
04187 throw ImageFormatException("fouriergridrot_shift2d requires a fourier image");
04188 int nxreal = nx - 2 + int(is_fftodd());
04189 if (nxreal != ny)
04190 throw ImageDimensionException("fouriergridrot_shift2d requires ny == nx(real)");
04191 if (0 != nxreal%2)
04192 throw ImageDimensionException("fouriergridrot_shift2d needs an even image.");
04193 int nxhalf = nxreal/2;
04194 int nyhalf = ny/2;
04195
04196 if (!is_shuffled()) fft_shuffle();
04197
04198 EMData* result = copy_head();
04199 set_array_offsets(0, -nyhalf);
04200 result->set_array_offsets(0, -nyhalf);
04201
04202 ang = ang*(float)DGR_TO_RAD;
04203 float cang = cos(ang);
04204 float sang = sin(ang);
04205 float temp = -2.0f*M_PI/nxreal;
04206 for (int iy = -nyhalf; iy < nyhalf; iy++) {
04207 float ycang = iy*cang;
04208 float ysang = iy*sang;
04209 for (int ix = 0; ix <= nxhalf; ix++) {
04210 float nuxold = ix*cang - ysang;
04211 float nuyold = ix*sang + ycang;
04212 result->cmplx(ix,iy) = Util::extractpoint2(nx, ny, nuxold, nuyold, this, kb);
04213
04214 float phase_ang = temp*(sx*ix+sy*iy);
04215 result->cmplx(ix,iy) *= complex<float>(cos(phase_ang), sin(phase_ang));
04216 }
04217 }
04218 result->set_array_offsets();
04219 result->fft_shuffle();
04220 result->update();
04221 set_array_offsets();
04222 fft_shuffle();
04223 return result;
04224 }
04225
04226 #undef QUADPI
04227 #undef DGR_TO_RAD
04228
04229 void EMData::divkbsinh(const Util::KaiserBessel& kb) {
04230
04231 if (is_complex())
04232 throw ImageFormatException("divkbsinh requires a real image.");
04233 vector<int> saved_offsets = get_array_offsets();
04234 set_array_offsets(0,0,0);
04235
04236
04237
04238
04239 for (int iz=0; iz < nz; iz++) {
04240 float wz = kb.sinhwin(static_cast<float>(iz-nz/2));
04241 for (int iy=0; iy < ny; iy++) {
04242 float wy = kb.sinhwin(static_cast<float>(iy-ny/2));
04243 for (int ix=0; ix < nx; ix++) {
04244 float wx = kb.sinhwin(static_cast<float>(ix-nx/2));
04245 float w = wx*wy*wz;
04246 (*this)(ix,iy,iz) /= w;
04247 }
04248 }
04249 }
04250 set_array_offsets(saved_offsets);
04251 }
04252
04253 void EMData::divkbsinh_rect(const Util::KaiserBessel& kbx, const Util::KaiserBessel& kby, const Util::KaiserBessel& kbz) {
04254
04255 if (is_complex())
04256 throw ImageFormatException("divkbsinh requires a real image.");
04257 vector<int> saved_offsets = get_array_offsets();
04258 set_array_offsets(0,0,0);
04259
04260
04261
04262
04263 for (int iz=0; iz < nz; iz++) {
04264 float wz = kbz.sinhwin(static_cast<float>(iz-nz/2));
04265 for (int iy=0; iy < ny; iy++) {
04266 float wy = kby.sinhwin(static_cast<float>(iy-ny/2));
04267 for (int ix=0; ix < nx; ix++) {
04268 float wx = kbx.sinhwin(static_cast<float>(ix-nx/2));
04269 float w = wx*wy*wz;
04270 (*this)(ix,iy,iz) /= w;
04271 }
04272 }
04273 }
04274
04275 set_array_offsets(saved_offsets);
04276 }
04277
04278
04279
04280
04281
04282
04283
04284
04285
04286
04287
04288
04289
04290
04291
04292
04293
04294
04295
04296
04297
04298
04299
04300
04301
04302
04303
04304 EMData* EMData::extract_plane(const Transform& tf, Util::KaiserBessel& kb) {
04305 if (!is_complex())
04306 throw ImageFormatException("extractplane requires a complex image");
04307 if (nx%2 != 0)
04308 throw ImageDimensionException("extractplane requires nx to be even");
04309 int nxreal = nx - 2;
04310 if (nxreal != ny || nxreal != nz)
04311 throw ImageDimensionException("extractplane requires ny == nx == nz");
04312
04313 EMData* res = new EMData();
04314 res->set_size(nx,ny,1);
04315 res->to_zero();
04316 res->set_complex(true);
04317 res->set_fftodd(false);
04318 res->set_fftpad(true);
04319 res->set_ri(true);
04320
04321 int n = nxreal;
04322 int nhalf = n/2;
04323 vector<int> saved_offsets = get_array_offsets();
04324 set_array_offsets(0,-nhalf,-nhalf);
04325 res->set_array_offsets(0,-nhalf,0);
04326
04327 int kbsize = kb.get_window_size();
04328 int kbmin = -kbsize/2;
04329 int kbmax = -kbmin;
04330 float* wy0 = new float[kbmax - kbmin + 1];
04331 float* wy = wy0 - kbmin;
04332 float* wx0 = new float[kbmax - kbmin + 1];
04333 float* wx = wx0 - kbmin;
04334 float* wz0 = new float[kbmax - kbmin + 1];
04335 float* wz = wz0 - kbmin;
04336 float rim = nhalf*float(nhalf);
04337 int count = 0;
04338 float wsum = 0.f;
04339 Transform tftrans = tf;
04340 tftrans.invert();
04341 for (int jy = -nhalf; jy < nhalf; jy++)
04342 {
04343 for (int jx = 0; jx <= nhalf; jx++)
04344 {
04345 Vec3f nucur((float)jx, (float)jy, 0.f);
04346 Vec3f nunew = tftrans*nucur;
04347 float xnew = nunew[0], ynew = nunew[1], znew = nunew[2];
04348 if (xnew*xnew+ynew*ynew+znew*znew <= rim)
04349 {
04350 count++;
04351 std::complex<float> btq(0.f,0.f);
04352 bool flip = false;
04353 if (xnew < 0.f) {
04354 flip = true;
04355 xnew = -xnew;
04356 ynew = -ynew;
04357 znew = -znew;
04358 }
04359 int ixn = int(Util::round(xnew));
04360 int iyn = int(Util::round(ynew));
04361 int izn = int(Util::round(znew));
04362
04363 for (int i=kbmin; i <= kbmax; i++) {
04364 int izp = izn + i;
04365 wz[i] = kb.i0win_tab(znew - izp);
04366 int iyp = iyn + i;
04367 wy[i] = kb.i0win_tab(ynew - iyp);
04368 int ixp = ixn + i;
04369 wx[i] = kb.i0win_tab(xnew - ixp);
04370
04371 }
04372
04373 int lnbz = 0;
04374 for (int iz = kbmin; iz <= -1; iz++) {
04375 if (wz[iz] != 0.f) {
04376 lnbz = iz;
04377 break;
04378 }
04379 }
04380 int lnez = 0;
04381 for (int iz = kbmax; iz >= 1; iz--) {
04382 if (wz[iz] != 0.f) {
04383 lnez = iz;
04384 break;
04385 }
04386 }
04387 int lnby = 0;
04388 for (int iy = kbmin; iy <= -1; iy++) {
04389 if (wy[iy] != 0.f) {
04390 lnby = iy;
04391 break;
04392 }
04393 }
04394 int lney = 0;
04395 for (int iy = kbmax; iy >= 1; iy--) {
04396 if (wy[iy] != 0.f) {
04397 lney = iy;
04398 break;
04399 }
04400 }
04401 int lnbx = 0;
04402 for (int ix = kbmin; ix <= -1; ix++) {
04403 if (wx[ix] != 0.f) {
04404 lnbx = ix;
04405 break;
04406 }
04407 }
04408 int lnex = 0;
04409 for (int ix = kbmax; ix >= 1; ix--) {
04410 if (wx[ix] != 0.f) {
04411 lnex = ix;
04412 break;
04413 }
04414 }
04415 if (ixn >= -kbmin && ixn <= nhalf-1-kbmax
04416 && iyn >= -nhalf-kbmin && iyn <= nhalf-1-kbmax
04417 && izn >= -nhalf-kbmin && izn <= nhalf-1-kbmax) {
04418
04419 for (int lz = lnbz; lz <= lnez; lz++) {
04420 int izp = izn + lz;
04421 for (int ly=lnby; ly<=lney; ly++) {
04422 int iyp = iyn + ly;
04423 float ty = wz[lz]*wy[ly];
04424 for (int lx=lnbx; lx<=lnex; lx++) {
04425 int ixp = ixn + lx;
04426 float wg = wx[lx]*ty;
04427 btq += cmplx(ixp,iyp,izp)*wg;
04428 wsum += wg;
04429 }
04430 }
04431 }
04432 } else {
04433
04434 for (int lz = lnbz; lz <= lnez; lz++) {
04435 int izp = izn + lz;
04436 for (int ly=lnby; ly<=lney; ly++) {
04437 int iyp = iyn + ly;
04438 float ty = wz[lz]*wy[ly];
04439 for (int lx=lnbx; lx<=lnex; lx++) {
04440 int ixp = ixn + lx;
04441 float wg = wx[lx]*ty;
04442 bool mirror = false;
04443 int ixt(ixp), iyt(iyp), izt(izp);
04444 if (ixt > nhalf || ixt < -nhalf) {
04445 ixt = Util::sgn(ixt)
04446 *(n - abs(ixt));
04447 iyt = -iyt;
04448 izt = -izt;
04449 mirror = !mirror;
04450 }
04451 if (iyt >= nhalf || iyt < -nhalf) {
04452 if (ixt != 0) {
04453 ixt = -ixt;
04454 iyt = Util::sgn(iyt)
04455 *(n - abs(iyt));
04456 izt = -izt;
04457 mirror = !mirror;
04458 } else {
04459 iyt -= n*Util::sgn(iyt);
04460 }
04461 }
04462 if (izt >= nhalf || izt < -nhalf) {
04463 if (ixt != 0) {
04464 ixt = -ixt;
04465 iyt = -iyt;
04466 izt = Util::sgn(izt)
04467 *(n - abs(izt));
04468 mirror = !mirror;
04469 } else {
04470 izt -= Util::sgn(izt)*n;
04471 }
04472 }
04473 if (ixt < 0) {
04474 ixt = -ixt;
04475 iyt = -iyt;
04476 izt = -izt;
04477 mirror = !mirror;
04478 }
04479 if (iyt == nhalf) iyt = -nhalf;
04480 if (izt == nhalf) izt = -nhalf;
04481 if (mirror) btq += conj(cmplx(ixt,iyt,izt))*wg;
04482 else btq += cmplx(ixt,iyt,izt)*wg;
04483 wsum += wg;
04484 }
04485 }
04486 }
04487 }
04488 if (flip) res->cmplx(jx,jy) = conj(btq);
04489 else res->cmplx(jx,jy) = btq;
04490 }
04491 }
04492 }
04493 for (int jy = -nhalf; jy < nhalf; jy++)
04494 for (int jx = 0; jx <= nhalf; jx++)
04495 res->cmplx(jx,jy) *= count/wsum;
04496 delete[] wx0; delete[] wy0; delete[] wz0;
04497 set_array_offsets(saved_offsets);
04498 res->set_array_offsets(0,0,0);
04499 res->set_shuffled(true);
04500 return res;
04501 }
04502
04503
04505
04506
04507
04508 EMData* EMData::extract_plane_rect(const Transform& tf, Util::KaiserBessel& kbx,Util::KaiserBessel& kby, Util::KaiserBessel& kbz) {
04509
04510
04511 if (!is_complex())
04512 throw ImageFormatException("extractplane requires a complex image");
04513 if (nx%2 != 0)
04514 throw ImageDimensionException("extractplane requires nx to be even");
04515
04516 int nxfromz = nz+2;
04517 int nxcircal = nxfromz - 2;
04518 EMData* res = new EMData();
04519 res->set_size(nxfromz,nz,1);
04520 res->to_zero();
04521 res->set_complex(true);
04522 res->set_fftodd(false);
04523 res->set_fftpad(true);
04524 res->set_ri(true);
04525
04526 int n = nxcircal;
04527 int nhalf = n/2;
04528 int nxhalf = (nx-2)/2;
04529 int nyhalf = ny/2;
04530 int nzhalf = nz/2;
04531
04532 vector<int> saved_offsets = get_array_offsets();
04533 set_array_offsets(0, -nyhalf, -nzhalf);
04534 res->set_array_offsets(0,-nhalf,0);
04535
04536 int kbxsize = kbx.get_window_size();
04537 int kbxmin = -kbxsize/2;
04538 int kbxmax = -kbxmin;
04539
04540 int kbysize = kby.get_window_size();
04541 int kbymin = -kbysize/2;
04542 int kbymax = -kbymin;
04543
04544 int kbzsize = kbz.get_window_size();
04545 int kbzmin = -kbzsize/2;
04546 int kbzmax = -kbzmin;
04547
04548
04549 float* wy0 = new float[kbymax - kbymin + 1];
04550 float* wy = wy0 - kbymin;
04551 float* wx0 = new float[kbxmax - kbxmin + 1];
04552 float* wx = wx0 - kbxmin;
04553 float* wz0 = new float[kbzmax - kbzmin + 1];
04554 float* wz = wz0 - kbzmin;
04555 float rim = nhalf*float(nhalf);
04556 int count = 0;
04557 float wsum = 0.f;
04558 Transform tftrans = tf;
04559 tftrans.invert();
04560 float xratio=float(nx-2)/float(nz);
04561 float yratio=float(ny)/float(nz);
04562
04563 for (int jy = -nhalf; jy < nhalf; jy++)
04564 {
04565 for (int jx = 0; jx <= nhalf; jx++)
04566 {
04567 Vec3f nucur((float)jx, (float)jy, 0.f);
04568 Vec3f nunew = tftrans*nucur;
04569 float xnew = nunew[0]*xratio, ynew = nunew[1]*yratio, znew = nunew[2];
04570
04571 if (nunew[0]*nunew[0]+nunew[1]*nunew[1]+nunew[2]*nunew[2] <= rim)
04572 {
04573 count++;
04574 std::complex<float> btq(0.f,0.f);
04575 bool flip = false;
04576 if (xnew < 0.f) {
04577 flip = true;
04578 xnew = -xnew;
04579 ynew = -ynew;
04580 znew = -znew;
04581 }
04582 int ixn = int(Util::round(xnew));
04583 int iyn = int(Util::round(ynew));
04584 int izn = int(Util::round(znew));
04585
04586 for (int i=kbzmin; i <= kbzmax; i++) {
04587 int izp = izn + i;
04588 wz[i] = kbz.i0win_tab(znew - izp);
04589 }
04590 for (int i=kbymin; i <= kbymax; i++) {
04591 int iyp = iyn + i;
04592 wy[i] = kby.i0win_tab(ynew - iyp);
04593 }
04594 for (int i=kbxmin; i <= kbxmax; i++) {
04595 int ixp = ixn + i;
04596 wx[i] = kbx.i0win_tab(xnew - ixp);
04597 }
04598
04599
04600
04601
04602 int lnbz = 0;
04603 for (int iz = kbzmin; iz <= -1; iz++) {
04604 if (wz[iz] != 0.f) {
04605 lnbz = iz;
04606 break;
04607 }
04608 }
04609 int lnez = 0;
04610 for (int iz = kbzmax; iz >= 1; iz--) {
04611 if (wz[iz] != 0.f) {
04612 lnez = iz;
04613 break;
04614 }
04615 }
04616 int lnby = 0;
04617 for (int iy = kbymin; iy <= -1; iy++) {
04618 if (wy[iy] != 0.f) {
04619 lnby = iy;
04620 break;
04621 }
04622 }
04623 int lney = 0;
04624 for (int iy = kbymax; iy >= 1; iy--) {
04625 if (wy[iy] != 0.f) {
04626 lney = iy;
04627 break;
04628 }
04629 }
04630 int lnbx = 0;
04631 for (int ix = kbxmin; ix <= -1; ix++) {
04632 if (wx[ix] != 0.f) {
04633 lnbx = ix;
04634 break;
04635 }
04636 }
04637 int lnex = 0;
04638 for (int ix = kbxmax; ix >= 1; ix--) {
04639 if (wx[ix] != 0.f) {
04640 lnex = ix;
04641 break;
04642 }
04643 }
04644 if (ixn >= -kbxmin && ixn <= nxhalf-1-kbxmax
04645 && iyn >= -nyhalf-kbymin && iyn <= nyhalf-1-kbymax
04646 && izn >= -nzhalf-kbzmin && izn <= nzhalf-1-kbzmax) {
04647
04648 for (int lz = lnbz; lz <= lnez; lz++) {
04649 int izp = izn + lz;
04650 for (int ly=lnby; ly<=lney; ly++) {
04651 int iyp = iyn + ly;
04652 float ty = wz[lz]*wy[ly];
04653 for (int lx=lnbx; lx<=lnex; lx++) {
04654 int ixp = ixn + lx;
04655 float wg = wx[lx]*ty;
04656 btq += cmplx(ixp,iyp,izp)*wg;
04657 wsum += wg;
04658 }
04659 }
04660 }
04661 }
04662 else {
04663
04664 for (int lz = lnbz; lz <= lnez; lz++) {
04665 int izp = izn + lz;
04666 for (int ly=lnby; ly<=lney; ly++) {
04667 int iyp = iyn + ly;
04668 float ty = wz[lz]*wy[ly];
04669 for (int lx=lnbx; lx<=lnex; lx++) {
04670 int ixp = ixn + lx;
04671 float wg = wx[lx]*ty;
04672 bool mirror = false;
04673 int ixt(ixp), iyt(iyp), izt(izp);
04674 if (ixt > nxhalf || ixt < -nxhalf) {
04675 ixt = Util::sgn(ixt)
04676 *(nx-2-abs(ixt));
04677 iyt = -iyt;
04678 izt = -izt;
04679 mirror = !mirror;
04680 }
04681 if (iyt >= nyhalf || iyt < -nyhalf) {
04682 if (ixt != 0) {
04683 ixt = -ixt;
04684 iyt = Util::sgn(iyt)
04685 *(ny - abs(iyt));
04686 izt = -izt;
04687 mirror = !mirror;
04688 } else {
04689 iyt -= ny*Util::sgn(iyt);
04690 }
04691 }
04692 if (izt >= nzhalf || izt < -nzhalf) {
04693 if (ixt != 0) {
04694 ixt = -ixt;
04695 iyt = -iyt;
04696 izt = Util::sgn(izt)
04697 *(nz - abs(izt));
04698 mirror = !mirror;
04699 } else {
04700 izt -= Util::sgn(izt)*nz;
04701 }
04702 }
04703 if (ixt < 0) {
04704 ixt = -ixt;
04705 iyt = -iyt;
04706 izt = -izt;
04707 mirror = !mirror;
04708 }
04709 if (iyt == nyhalf) iyt = -nyhalf;
04710 if (izt == nzhalf) izt = -nzhalf;
04711 if (mirror) btq += conj(cmplx(ixt,iyt,izt))*wg;
04712 else btq += cmplx(ixt,iyt,izt)*wg;
04713 wsum += wg;
04714 }
04715 }
04716 }
04717 }
04718 if (flip) res->cmplx(jx,jy) = conj(btq);
04719 else res->cmplx(jx,jy) = btq;
04720 }
04721 }
04722 }
04723 for (int jy = -nhalf; jy < nhalf; jy++)
04724 for (int jx = 0; jx <= nhalf; jx++)
04725 res->cmplx(jx,jy) *= count/wsum;
04726 delete[] wx0; delete[] wy0; delete[] wz0;
04727 set_array_offsets(saved_offsets);
04728 res->set_array_offsets(0,0,0);
04729 res->set_shuffled(true);
04730 return res;
04731 }
04732
04733
04734
04735 EMData* EMData::extract_plane_rect_fast(const Transform& tf, Util::KaiserBessel& kbx,Util::KaiserBessel& kby, Util::KaiserBessel& kbz) {
04736
04737
04738
04739 if (!is_complex())
04740 throw ImageFormatException("extractplane requires a complex image");
04741 if (nx%2 != 0)
04742 throw ImageDimensionException("extractplane requires nx to be even");
04743
04744 int nxfromz=nz+2;
04745 int nxcircal = nxfromz - 2;
04746
04747
04748 float xratio=float(nx-2)/float(nz);
04749 float yratio=float(ny)/float(nz);
04750 Vec3f axis_newx,axis_newy;
04751 axis_newx[0] = xratio*0.5f*nz*tf[0][0];
04752 axis_newx[1] = yratio*0.5f*nz*tf[0][1];
04753 axis_newx[2] = 0.5f*nz*tf[0][2];
04754
04755
04756 float ellipse_length_x=std::sqrt(axis_newx[0]*axis_newx[0]+axis_newx[1]*axis_newx[1]+axis_newx[2]*axis_newx[2]);
04757
04758 int ellipse_length_x_int=int(ellipse_length_x);
04759 float ellipse_step_x=0.5f*nz/float(ellipse_length_x_int);
04760 float xscale=ellipse_step_x;
04761
04762 axis_newy[0] = xratio*0.5f*nz*tf[1][0];
04763 axis_newy[1] = yratio*0.5f*nz*tf[1][1];
04764 axis_newy[2] = 0.5f*nz*tf[1][2];
04765
04766
04767 float ellipse_length_y=std::sqrt(axis_newy[0]*axis_newy[0]+axis_newy[1]*axis_newy[1]+axis_newy[2]*axis_newy[2]);
04768 int ellipse_length_y_int=int(ellipse_length_y);
04769 float ellipse_step_y=0.5f*nz/float(ellipse_length_y_int);
04770 float yscale=ellipse_step_y;
04771
04772 int nx_e=ellipse_length_x_int*2;
04773 int ny_e=ellipse_length_y_int*2;
04774 int nx_ec=nx_e+2;
04775
04776 EMData* res = new EMData();
04777 res->set_size(nx_ec,ny_e,1);
04778 res->to_zero();
04779 res->set_complex(true);
04780 res->set_fftodd(false);
04781 res->set_fftpad(true);
04782 res->set_ri(true);
04783
04784
04785
04786 int n = nxcircal;
04787 int nhalf = n/2;
04788 int nhalfx_e = nx_e/2;
04789 int nhalfy_e = ny_e/2;
04790 int nxhalf=(nx-2)/2;
04791 int nyhalf=ny/2;
04792 int nzhalf=nz/2;
04793
04794 vector<int> saved_offsets = get_array_offsets();
04795 set_array_offsets(0,-nyhalf,-nzhalf);
04796 res->set_array_offsets(0,-nhalfy_e,0);
04797
04798 int kbxsize = kbx.get_window_size();
04799 int kbxmin = -kbxsize/2;
04800 int kbxmax = -kbxmin;
04801
04802 int kbysize = kby.get_window_size();
04803 int kbymin = -kbysize/2;
04804 int kbymax = -kbymin;
04805
04806 int kbzsize = kbz.get_window_size();
04807 int kbzmin = -kbzsize/2;
04808 int kbzmax = -kbzmin;
04809
04810
04811 float* wy0 = new float[kbymax - kbymin + 1];
04812 float* wy = wy0 - kbymin;
04813 float* wx0 = new float[kbxmax - kbxmin + 1];
04814 float* wx = wx0 - kbxmin;
04815 float* wz0 = new float[kbzmax - kbzmin + 1];
04816 float* wz = wz0 - kbzmin;
04817 float rim = nhalf*float(nhalf);
04818 int count = 0;
04819 float wsum = 0.f;
04820 Transform tftrans = tf;
04821 tftrans.invert();
04822
04823
04824 for (int jy = -nhalfy_e; jy < nhalfy_e; jy++)
04825 {
04826 for (int jx = 0; jx <= nhalfx_e; jx++)
04827 {
04828 Vec3f nucur((float)jx, (float)jy, 0.f);
04829 nucur[0]=nucur[0]*xscale;nucur[1]=nucur[1]*yscale;;
04830 Vec3f nunew = tftrans*nucur;
04831 float xnew = nunew[0]*xratio, ynew = nunew[1]*yratio, znew = nunew[2];
04832
04833 if (nunew[0]*nunew[0]+nunew[1]*nunew[1]+nunew[2]*nunew[2] <= rim)
04834 {
04835 count++;
04836 std::complex<float> btq(0.f,0.f);
04837 bool flip = false;
04838 if (xnew < 0.f) {
04839 flip = true;
04840 xnew = -xnew;
04841 ynew = -ynew;
04842 znew = -znew;
04843 }
04844 int ixn = int(Util::round(xnew));
04845 int iyn = int(Util::round(ynew));
04846 int izn = int(Util::round(znew));
04847
04848 for (int i=kbzmin; i <= kbzmax; i++) {
04849 int izp = izn + i;
04850 wz[i] = kbz.i0win_tab(znew - izp);
04851 }
04852 for (int i=kbymin; i <= kbymax; i++) {
04853 int iyp = iyn + i;
04854 wy[i] = kby.i0win_tab(ynew - iyp);
04855 }
04856 for (int i=kbxmin; i <= kbxmax; i++) {
04857 int ixp = ixn + i;
04858 wx[i] = kbx.i0win_tab(xnew - ixp);
04859 }
04860
04861
04862
04863
04864 int lnbz = 0;
04865 for (int iz = kbzmin; iz <= -1; iz++) {
04866 if (wz[iz] != 0.f) {
04867 lnbz = iz;
04868 break;
04869 }
04870 }
04871 int lnez = 0;
04872 for (int iz = kbzmax; iz >= 1; iz--) {
04873 if (wz[iz] != 0.f) {
04874 lnez = iz;
04875 break;
04876 }
04877 }
04878 int lnby = 0;
04879 for (int iy = kbymin; iy <= -1; iy++) {
04880 if (wy[iy] != 0.f) {
04881 lnby = iy;
04882 break;
04883 }
04884 }
04885 int lney = 0;
04886 for (int iy = kbymax; iy >= 1; iy--) {
04887 if (wy[iy] != 0.f) {
04888 lney = iy;
04889 break;
04890 }
04891 }
04892 int lnbx = 0;
04893 for (int ix = kbxmin; ix <= -1; ix++) {
04894 if (wx[ix] != 0.f) {
04895 lnbx = ix;
04896 break;
04897 }
04898 }
04899 int lnex = 0;
04900 for (int ix = kbxmax; ix >= 1; ix--) {
04901 if (wx[ix] != 0.f) {
04902 lnex = ix;
04903 break;
04904 }
04905 }
04906 if (ixn >= -kbxmin && ixn <= nxhalf-1-kbxmax
04907 && iyn >= -nyhalf-kbymin && iyn <= nyhalf-1-kbymax
04908 && izn >= -nzhalf-kbzmin && izn <= nzhalf-1-kbzmax) {
04909
04910 for (int lz = lnbz; lz <= lnez; lz++) {
04911 int izp = izn + lz;
04912 for (int ly=lnby; ly<=lney; ly++) {
04913 int iyp = iyn + ly;
04914 float ty = wz[lz]*wy[ly];
04915 for (int lx=lnbx; lx<=lnex; lx++) {
04916 int ixp = ixn + lx;
04917 float wg = wx[lx]*ty;
04918 btq += cmplx(ixp,iyp,izp)*wg;
04919 wsum += wg;
04920 }
04921 }
04922 }
04923 }
04924 else {
04925
04926 for (int lz = lnbz; lz <= lnez; lz++) {
04927 int izp = izn + lz;
04928 for (int ly=lnby; ly<=lney; ly++) {
04929 int iyp = iyn + ly;
04930 float ty = wz[lz]*wy[ly];
04931 for (int lx=lnbx; lx<=lnex; lx++) {
04932 int ixp = ixn + lx;
04933 float wg = wx[lx]*ty;
04934 bool mirror = false;
04935 int ixt(ixp), iyt(iyp), izt(izp);
04936 if (ixt > nxhalf || ixt < -nxhalf) {
04937 ixt = Util::sgn(ixt)
04938 *(nx-2-abs(ixt));
04939 iyt = -iyt;
04940 izt = -izt;
04941 mirror = !mirror;
04942 }
04943 if (iyt >= nyhalf || iyt < -nyhalf) {
04944 if (ixt != 0) {
04945 ixt = -ixt;
04946 iyt = Util::sgn(iyt)
04947 *(ny - abs(iyt));
04948 izt = -izt;
04949 mirror = !mirror;
04950 } else {
04951 iyt -= ny*Util::sgn(iyt);
04952 }
04953 }
04954 if (izt >= nzhalf || izt < -nzhalf) {
04955 if (ixt != 0) {
04956 ixt = -ixt;
04957 iyt = -iyt;
04958 izt = Util::sgn(izt)
04959 *(nz - abs(izt));
04960 mirror = !mirror;
04961 } else {
04962 izt -= Util::sgn(izt)*nz;
04963 }
04964 }
04965 if (ixt < 0) {
04966 ixt = -ixt;
04967 iyt = -iyt;
04968 izt = -izt;
04969 mirror = !mirror;
04970 }
04971 if (iyt == nyhalf) iyt = -nyhalf;
04972 if (izt == nzhalf) izt = -nzhalf;
04973 if (mirror) btq += conj(cmplx(ixt,iyt,izt))*wg;
04974 else btq += cmplx(ixt,iyt,izt)*wg;
04975 wsum += wg;
04976 }
04977 }
04978 }
04979 }
04980 if (flip) res->cmplx(jx,jy) = conj(btq);
04981 else res->cmplx(jx,jy) = btq;
04982 }
04983 }
04984 }
04985 for (int jy = -nhalfy_e; jy < nhalfy_e; jy++)
04986 for (int jx = 0; jx <= nhalfx_e; jx++)
04987 res->cmplx(jx,jy) *= count/wsum;
04988 delete[] wx0; delete[] wy0; delete[] wz0;
04989 set_array_offsets(saved_offsets);
04990 res->set_array_offsets(0,0,0);
04991 res->set_shuffled(true);
04992 return res;
04993 }
04994
04995
04996
04997
04998 bool EMData::peakcmp(const Pixel& p1, const Pixel& p2) {
04999 return (p1.value > p2.value);
05000 }
05001
05002 ostream& operator<< (ostream& os, const Pixel& peak) {
05003 os << peak.x << peak.y << peak.z << peak.value;
05004 return os;
05005 }
05006
05007
05008
05009
05010
05011
05012
05013
05014
05015
05016
05017
05018
05019
05020
05021
05022
05023
05024
05025
05026
05027
05028
05029
05030
05031
05032
05033
05034
05035
05036
05037
05038
05039
05040
05041
05042
05043
05044
05045
05046
05047
05048
05049
05050
05051
05052
05053
05054
05055
05056
05057
05058
05059
05060
05061
05062
05063
05064
05065
05066
05067
05068
05069
05070
05071 vector<float> EMData::peak_search(int ml, float invert) {
05072
05073 EMData& buf = *this;
05074 vector<Pixel> peaks;
05075 int img_dim;
05076 int i, j, k;
05077 int i__1, i__2;
05078 int j__1, j__2;
05079
05080 bool peak_check;
05081 img_dim=buf.get_ndim();
05082 vector<int> ix, jy, kz;
05083 vector<float>res;
05084 int nx = buf.get_xsize();
05085 int ny = buf.get_ysize();
05086 int nz = buf.get_zsize();
05087 if(invert <= 0.0f) invert=-1.0f;
05088 else invert=1.0f ;
05089 int count = 0;
05090 switch (img_dim) {
05091 case(1):
05092 for(i=0;i<=nx-1;++i) {
05093 i__1 = (i-1+nx)%nx;
05094 i__2 = (i+1)%nx;
05095
05096
05097
05098 float qbf = buf(i)*invert;
05099 peak_check = qbf > buf(i__1)*invert && qbf > buf(i__2)*invert;
05100 if(peak_check) {
05101 if(count < ml) {
05102 count++;
05103 peaks.push_back( Pixel(i, 0, 0, qbf) );
05104 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
05105 } else {
05106 if( qbf > (peaks.back()).value ) {
05107
05108 peaks.pop_back();
05109 peaks.push_back( Pixel(i, 0, 0, qbf) );
05110 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
05111 }
05112 }
05113 }
05114 }
05115 break;
05116 case(2):
05117
05118
05119
05120
05121
05122
05123
05124
05125 for(j=1;j<=ny-2;++j) {
05126 j__1 = j-1;
05127 j__2 = j+1;
05128 for(i=1;i<=nx-2;++i) {
05129 i__1 = i-1;
05130 i__2 = i+1;
05131 float qbf = buf(i,j)*invert;
05132 peak_check = (qbf > buf(i,j__1)*invert) && (qbf > buf(i,j__2)*invert);
05133 if(peak_check) {
05134 peak_check = (qbf > buf(i__1,j)*invert) && (qbf > buf(i__2,j)*invert);
05135 if(peak_check) {
05136 peak_check = (qbf > buf(i__1,j__1)*invert) && (qbf > buf(i__1,j__2)*invert);
05137 if(peak_check) {
05138 peak_check = (qbf > buf(i__2,j__1)*invert) && (qbf > buf(i__2,j__2)*invert);
05139 if(peak_check) {
05140 if(count < ml) {
05141 count++;
05142 peaks.push_back( Pixel(i, j, 0, qbf) );
05143 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
05144 } else {
05145 if( qbf > (peaks.back()).value ) {
05146
05147 peaks.pop_back();
05148 peaks.push_back( Pixel(i, j, 0, qbf) );
05149 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
05150 }
05151 }
05152 }
05153 }
05154 }
05155 }
05156 }
05157 }
05158 break;
05159 case(3):
05160
05161
05162
05163
05164
05165
05166
05167
05168
05169
05170
05171
05172
05173
05174
05175
05176
05177
05178
05179
05180 for(k=1; k<=nz-2; ++k) {
05181 kz.clear();
05182 kz.push_back(k-1);
05183 kz.push_back(k);
05184 kz.push_back(k+1);
05185 for(j=1; j<=ny-2; ++j) {
05186 jy.clear();
05187 jy.push_back(j-1);
05188 jy.push_back(j);
05189 jy.push_back(j+1);
05190 for(i=1; i<=nx-2; ++i) {
05191 ix.clear();
05192 ix.push_back(i-1);
05193 ix.push_back(i);
05194 ix.push_back(i+1);
05195 float qbf = buf(i,j,k)*invert;
05196 peak_check = qbf > buf(ix[0],jy[0],kz[0])*invert;
05197 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[0])*invert;
05198 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[0])*invert;
05199 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[0])*invert;
05200 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[0])*invert;
05201 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[0])*invert;
05202 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[0])*invert;
05203 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[0])*invert;
05204 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[0])*invert;
05205 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[1])*invert;
05206 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[1])*invert;
05207 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[1])*invert;
05208 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[1])*invert;
05209
05210 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[1])*invert;
05211 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[1])*invert;
05212 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[1])*invert;
05213 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[1])*invert;
05214 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[2])*invert;
05215 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[2])*invert;
05216 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[2])*invert;
05217 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[2])*invert;
05218 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[2])*invert;
05219 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[2])*invert;
05220 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[2])*invert;
05221 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[2])*invert;
05222 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[2])*invert;
05223 if(peak_check) {
05224 if(count < ml) {
05225 count++;
05226 peaks.push_back( Pixel(i, j, k, qbf) );
05227 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
05228 } else {
05229 if( qbf > (peaks.back()).value ) {
05230
05231
05232 peaks.pop_back();
05233 peaks.push_back( Pixel(i, j, k, qbf) );
05234 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
05235 }
05236 }
05237 }
05238 }}}}}}}}}}}}}}}}}}}}}}}}}
05239 }
05240 }
05241 }
05242
05243
05244 for(k=1; k<=nz-2; ++k) {
05245 kz.clear();
05246 kz.push_back(k-1);
05247 kz.push_back(k);
05248 kz.push_back(k+1);
05249 for(j=1; j<=ny-2; ++j) {
05250 jy.clear();
05251 jy.push_back(j-1);
05252 jy.push_back(j);
05253 jy.push_back(j+1);
05254 for(i=0; i<=0; ++i) {
05255 ix.clear();
05256 ix.push_back(nx-1);
05257 ix.push_back(i);
05258 ix.push_back(i+1);
05259 float qbf = buf(i,j,k)*invert;
05260 peak_check = qbf > buf(ix[0],jy[0],kz[0])*invert;
05261 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[0])*invert;
05262 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[0])*invert;
05263 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[0])*invert;
05264 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[0])*invert;
05265 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[0])*invert;
05266 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[0])*invert;
05267 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[0])*invert;
05268 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[0])*invert;
05269 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[1])*invert;
05270 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[1])*invert;
05271 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[1])*invert;
05272 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[1])*invert;
05273
05274 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[1])*invert;
05275 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[1])*invert;
05276 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[1])*invert;
05277 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[1])*invert;
05278 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[2])*invert;
05279 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[2])*invert;
05280 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[2])*invert;
05281 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[2])*invert;
05282 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[2])*invert;
05283 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[2])*invert;
05284 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[2])*invert;
05285 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[2])*invert;
05286 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[2])*invert;
05287 if(peak_check) {
05288 if(count < ml) {
05289 count++;
05290 peaks.push_back( Pixel(i, j, k, qbf) );
05291 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
05292 } else {
05293 if( qbf > (peaks.back()).value ) {
05294
05295
05296 peaks.pop_back();
05297 peaks.push_back( Pixel(i, j, k, qbf) );
05298 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
05299 }
05300 }
05301 }
05302 }}}}}}}}}}}}}}}}}}}}}}}}}
05303 }
05304 for(i=nx-1; i<=nx-1; ++i) {
05305 ix.clear();
05306 ix.push_back(i-1);
05307 ix.push_back(i);
05308 ix.push_back(0);
05309 float qbf = buf(i,j,k)*invert;
05310 peak_check = qbf > buf(ix[0],jy[0],kz[0])*invert;
05311 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[0])*invert;
05312 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[0])*invert;
05313 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[0])*invert;
05314 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[0])*invert;
05315 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[0])*invert;
05316 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[0])*invert;
05317 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[0])*invert;
05318 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[0])*invert;
05319 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[1])*invert;
05320 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[1])*invert;
05321 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[1])*invert;
05322 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[1])*invert;
05323
05324 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[1])*invert;
05325 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[1])*invert;
05326 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[1])*invert;
05327 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[1])*invert;
05328 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[2])*invert;
05329 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[2])*invert;
05330 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[2])*invert;
05331 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[2])*invert;
05332 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[2])*invert;
05333 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[2])*invert;
05334 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[2])*invert;
05335 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[2])*invert;
05336 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[2])*invert;
05337 if(peak_check) {
05338 if(count < ml) {
05339 count++;
05340 peaks.push_back( Pixel(i, j, k, qbf) );
05341 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
05342 } else {
05343 if( qbf > (peaks.back()).value ) {
05344
05345
05346 peaks.pop_back();
05347 peaks.push_back( Pixel(i, j, k, qbf) );
05348 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
05349 }
05350 }
05351 }
05352 }}}}}}}}}}}}}}}}}}}}}}}}}
05353 }
05354 }
05355 }
05356 break;
05357
05358
05359
05360
05361
05362
05363
05364
05365
05366
05367
05368
05369
05370
05371
05372
05373
05374
05375
05376
05377
05378
05379
05380
05381
05382
05383
05384
05385
05386
05387
05388
05389
05390
05391
05392
05393
05394
05395
05396
05397
05398
05399
05400
05401
05402
05403
05404
05405
05406
05407
05408
05409
05410
05411
05412
05413
05414
05415
05416
05417
05418
05419
05420
05421
05422
05423
05424
05425
05426
05427
05428
05429
05430
05431
05432
05433
05434
05435
05436
05437
05438
05439
05440
05441
05442
05443
05444
05445
05446
05447
05448
05449
05450
05451
05452
05453
05454
05455
05456
05457
05458
05459
05460
05461
05462
05463
05464
05465
05466
05467
05468
05469
05470
05471
05472
05473
05474
05475
05476
05477
05478
05479
05480
05481
05482
05483
05484
05485
05486
05487
05488
05489
05490
05491
05492
05493
05494
05495
05496
05497
05498
05499
05500
05501
05502
05503
05504
05505
05506
05507
05508
05509
05510
05511
05512
05513
05514
05515
05516
05517
05518
05519
05520
05521
05522
05523
05524
05525
05526
05527
05528
05529
05530
05531
05532
05533
05534
05535
05536
05537
05538
05539
05540
05541
05542
05543
05544
05545
05546
05547
05548
05549
05550
05551
05552
05553
05554
05555
05556
05557
05558
05559
05560
05561
05562
05563
05564
05565
05566
05567
05568
05569
05570
05571
05572
05573
05574
05575
05576
05577
05578
05579
05580
05581
05582
05583
05584
05585
05586
05587
05588
05589
05590
05591
05592
05593
05594
05595
05596
05597
05598
05599
05600
05601
05602
05603
05604
05605
05606
05607
05608
05609
05610
05611
05612
05613
05614
05615
05616
05617
05618
05619
05620
05621
05622
05623
05624
05625
05626
05627
05628
05629
05630
05631
05632
05633
05634
05635
05636
05637
05638
05639
05640
05641
05642
05643
05644
05645
05646
05647
05648
05649
05650
05651
05652
05653
05654
05655
05656
05657
05658
05659
05660
05661
05662
05663
05664
05665
05666
05667
05668
05669
05670
05671
05672
05673
05674
05675
05676
05677
05678
05679
05680
05681
05682
05683
05684
05685
05686
05687
05688
05689
05690
05691
05692
05693
05694
05695
05696
05697
05698
05699
05700
05701
05702
05703
05704
05705
05706
05707
05708
05709
05710
05711
05712
05713
05714
05715
05716
05717
05718
05719
05720
05721
05722
05723
05724
05725
05726
05727
05728
05729
05730
05731
05732
05733
05734
05735
05736
05737
05738
05739
05740
05741
05742
05743
05744
05745
05746
05747
05748
05749
05750
05751
05752
05753
05754
05755
05756
05757
05758
05759
05760
05761
05762
05763
05764
05765
05766
05767
05768
05769
05770
05771
05772
05773
05774
05775
05776
05777
05778
05779
05780
05781
05782
05783
05784
05785
05786
05787
05788
05789
05790
05791
05792
05793
05794
05795
05796 }
05797
05798 if (peaks.begin() != peaks.end()) {
05799
05800 sort(peaks.begin(), peaks.end(), peakcmp);
05801
05802 int count = 0;
05803
05804 float xval = (*peaks.begin()).value;
05805
05806 for (vector<Pixel>::iterator it = peaks.begin(); it != peaks.end(); it++) {
05807
05808 count++;
05809
05810 if(count <= ml) {
05811
05812 res.push_back((*it).value);
05813 res.push_back(static_cast<float>((*it).x));
05814
05815 if(img_dim > 1) {
05816 res.push_back(static_cast<float>((*it).y));
05817 if(nz > 1) res.push_back(static_cast<float>((*it).z));
05818 }
05819
05820 if(xval != 0.0) res.push_back((*it).value/xval);
05821 else res.push_back((*it).value);
05822 res.push_back((*it).x-float(int(nx/2)));
05823 if(img_dim >1) {
05824 res.push_back((*it).y-float(int(ny/2)));
05825 if(nz>1) res.push_back((*it).z-float(nz/2));
05826 }
05827 }
05828 }
05829 res.insert(res.begin(),1,img_dim);
05830 } else {
05831
05832 res.push_back(buf(0,0,0));
05833 res.insert(res.begin(),1,0.0);
05834 }
05835
05836
05837 return res;
05838 }
05839
05840 #define rdata(i,j,k) rdata[(i-1)+((j-1)+(k-1)*ny)*(size_t)nx]
05841 #define X(i) X[i-1]
05842 #define Y(j) Y[j-1]
05843 #define Z(k) Z[k-1]
05844 vector<float> EMData::phase_cog()
05845 {
05846 vector<float> ph_cntog;
05847 int i=1,j=1,k=1;
05848 float C=0.f,S=0.f,P=0.f,F1=0.f,SNX;
05849 if (get_ndim()==1) {
05850 P = 8*atan(1.0f)/nx;
05851 for (i=1;i<=nx;i++) {
05852 C += cos(P * (i-1)) * rdata(i,j,k);
05853 S += sin(P * (i-1)) * rdata(i,j,k);
05854 }
05855 F1 = atan2(S,C);
05856 if (F1 < 0.0) F1 += 8*atan(1.0f);
05857 SNX = F1/P +1.0f;
05858 SNX = SNX - ((nx/2)+1);
05859 ph_cntog.push_back(SNX);
05860 #ifdef _WIN32
05861 ph_cntog.push_back((float)Util::round(SNX));
05862 #else
05863 ph_cntog.push_back(round(SNX));
05864 #endif //_WIN32
05865 } else if (get_ndim()==2) {
05866 #ifdef _WIN32
05867 float SNY;
05868 float T=0.0f;
05869 vector<float> X;
05870 X.resize(nx);
05871 #else
05872 float SNY,X[nx],T=0.f;
05873 #endif //_WIN32
05874 for ( i=1;i<=nx;i++) X(i)=0.0;
05875 P = 8*atan(1.0f)/ny;
05876 for(j=1;j<=ny;j++) {
05877 T=0.f;
05878 for(i=1;i<=nx;i++) {
05879 T += rdata(i,j,k);
05880 X(i)+=rdata(i,j,k);
05881 }
05882 C += cos(P*(j-1))*T;
05883 S += sin(P*(j-1))*T;
05884 }
05885 F1=atan2(S,C);
05886 if(F1<0.0) F1 += 8*atan(1.0f);
05887 SNY = F1/P +1.0f;
05888 C=0.f; S=0.f;
05889 P = 8*atan(1.0f)/nx;
05890 for(i=1;i<=nx;i++) {
05891 C += cos(P*(i-1))*X(i);
05892 S += sin(P*(i-1))*X(i);
05893 }
05894 F1=atan2(S,C);
05895 if(F1<0.0) F1 += 8*atan(1.0f);
05896 SNX = F1/P +1.0f;
05897 SNX = SNX - ((nx/2)+1);
05898 SNY = SNY - ((ny/2)+1);
05899 ph_cntog.push_back(SNX); ph_cntog.push_back(SNY);
05900 #ifdef _WIN32
05901 ph_cntog.push_back((float)Util::round(SNX)); ph_cntog.push_back((float)Util::round(SNY));
05902 #else
05903 ph_cntog.push_back(round(SNX)); ph_cntog.push_back(round(SNY));
05904 #endif //_WIN32
05905 } else {
05906 #ifdef _WIN32
05907 float val=0.f,sum1=0.f, SNY,SNZ;
05908 vector<float> X;
05909 X.resize(nx);
05910 vector<float> Y;
05911 Y.resize(ny);
05912 vector<float> Z;
05913 Z.resize(nz);
05914 #else
05915 float val=0.f, sum1=0.f, X[nx], Y[ny], Z[nz], SNY, SNZ;
05916 #endif //_WIN32
05917 for (i=1;i<=nx;i++) X(i)=0.0;
05918 for (j=1;j<=ny;j++) Y(j)=0.0;
05919 for (k=1;k<=nz;k++) Z(k)=0.0;
05920 for(k=1;k<=nz;k++) {
05921 for(j=1;j<=ny;j++) {
05922 sum1=0.f;
05923 for(i=1;i<=nx;i++) {
05924 val = rdata(i,j,k);
05925 sum1 += val;
05926 X(i) += val;
05927 }
05928 Y(j) += sum1;
05929 Z(k) += sum1;
05930 }
05931 }
05932 P = 8*atan(1.0f)/nx;
05933 for (i=1;i<=nx;i++) {
05934 C += cos(P*(i-1))*X(i);
05935 S += sin(P*(i-1))*X(i);
05936 }
05937 F1=atan2(S,C);
05938 if(F1<0.0) F1 += 8*atan(1.0f);
05939 SNX = F1/P +1.0f;
05940 C=0.f; S=0.f;
05941 P = 8*atan(1.0f)/ny;
05942 for(j=1;j<=ny;j++) {
05943 C += cos(P*(j-1))*Y(j);
05944 S += sin(P*(j-1))*Y(j);
05945 }
05946 F1=atan2(S,C);
05947 if(F1<0.0) F1 += 8*atan(1.0f);
05948 SNY = F1/P +1.0f;
05949 C=0.f; S=0.f;
05950 P = 8*atan(1.0f)/nz;
05951 for(k=1;k<=nz;k++) {
05952 C += cos(P*(k-1))*Z(k);
05953 S += sin(P*(k-1))*Z(k);
05954 }
05955 F1=atan2(S,C);
05956 if(F1<0.0) F1 += 8*atan(1.0f);
05957 SNZ = F1/P +1.0f;
05958 SNX = SNX - ((nx/2)+1);
05959 SNY = SNY - ((ny/2)+1);
05960 SNZ = SNZ - ((nz/2)+1);
05961 ph_cntog.push_back(SNX); ph_cntog.push_back(SNY); ph_cntog.push_back(SNZ);
05962 #ifdef _WIN32
05963 ph_cntog.push_back((float)Util::round(SNX)); ph_cntog.push_back((float)Util::round(SNY)); ph_cntog.push_back((float)Util::round(SNZ));
05964 #else
05965 ph_cntog.push_back(round(SNX)); ph_cntog.push_back(round(SNY));ph_cntog.push_back(round(SNZ));
05966 #endif
05967 }
05968 return ph_cntog;
05969 }
05970 #undef rdata
05971 #undef X
05972 #undef Y
05973 #undef Z
05974
05975 #define avagadro (6.023*(double)pow(10.0,23.0))
05976 #define density_protein (1.36)
05977 #define R (0.61803399f)
05978 #define C (1.f-R)
05979 float EMData::find_3d_threshold(float mass, float pixel_size)
05980 {
05981
05982 if(get_ndim()!=3)
05983 throw ImageDimensionException("The image should be 3D");
05984
05985
05986
05987 float density_1_mole, vol_1_mole, vol_angstrom;
05988 int vol_voxels;
05989 density_1_mole = static_cast<float>( (mass*1000.0f)/avagadro );
05990 vol_1_mole = static_cast<float>( density_1_mole/density_protein );
05991 vol_angstrom = static_cast<float>( vol_1_mole*(double)pow((double)pow(10.0,8),3) );
05992 vol_voxels = static_cast<int> (vol_angstrom/(double)pow(pixel_size,3));
05993
05994
05995
05996 float thr1 = get_attr("maximum");
05997 float thr3 = get_attr("minimum");
05998 float thr2 = (thr1-thr3)/2 + thr3;
05999 size_t size = (size_t)nx*ny*nz;
06000 float x0 = thr1,x3 = thr3,x1,x2,THR=0;
06001
06002 #ifdef _WIN32
06003 int ILE = _cpp_min(nx*ny*nx,_cpp_max(1,vol_voxels));
06004 #else
06005 int ILE = std::min(nx*ny*nx,std::max(1,vol_voxels));
06006 #endif //_WIN32
06007
06008 if (abs(thr3-thr2)>abs(thr2-thr1)) {
06009 x1=thr2;
06010 x2=thr2+C*(thr3-thr2);
06011 } else {
06012 x2=thr2;
06013 x1=thr2-C*(thr2-thr1);
06014 }
06015
06016 int cnt1=0,cnt2=0;
06017 for (size_t i=0;i<size;++i) {
06018 if(rdata[i]>=x1) cnt1++;
06019 if(rdata[i]>=x2) cnt2++;
06020 }
06021 float LF1 = static_cast<float>( cnt1 - ILE );
06022 float F1 = LF1*LF1;
06023 float LF2 = static_cast<float>( cnt2 - ILE );
06024 float F2 = LF2*LF2;
06025
06026 while ((LF1 != 0 || LF2 != 0) && (fabs(LF1-LF2) >= 1.f) && (abs(x1-x2) > (double)pow(10.0,-5) && abs(x1-x3) > (double)pow(10.0,-5) && abs(x2-x3) > (double)pow(10.0,-5)))
06027 {
06028 if(F2 < F1) {
06029 x0=x1;
06030 x1=x2;
06031 x2 = R*x1 + C*x3;
06032 F1=F2;
06033 int cnt=0;
06034 for(size_t i=0;i<size;++i)
06035 if(rdata[i]>=x2)
06036 cnt++;
06037 LF2 = static_cast<float>( cnt - ILE );
06038 F2 = LF2*LF2;
06039 } else {
06040 x3=x2;
06041 x2=x1;
06042 x1=R*x2 + C*x0;
06043 F2=F1;
06044 int cnt=0;
06045 for(size_t i=0;i<size;++i)
06046 if(rdata[i]>=x1)
06047 cnt++;
06048 LF1 = static_cast<float>( cnt - ILE );
06049 F1 = LF1*LF1;
06050 }
06051 }
06052
06053 if(F1 < F2) {
06054 ILE = static_cast<int> (LF1 + ILE);
06055 THR = x1;
06056 } else {
06057 ILE = static_cast<int> (LF2 + ILE);
06058 THR = x2;
06059 }
06060 return THR;
06061
06062 }
06063 #undef avagadro
06064 #undef density_protein
06065 #undef R
06066 #undef C
06067
06068
06069
06070
06071
06072
06073
06074 vector<float> EMData::peak_ccf(float hf_p)
06075 {
06076
06077
06078
06079 EMData & buf = *this;
06080 vector<Pixel> peaks;
06081 int half=int(hf_p);
06082 float hf_p2 = hf_p*hf_p;
06083 int i,j;
06084 int i__1,i__2;
06085 int j__1,j__2;
06086 vector<float>res;
06087 int nx = buf.get_xsize()-half;
06088 int ny = buf.get_ysize()-half;
06089
06090 for(i=half; i<=nx; ++i) {
06091
06092 i__1 = i-1;
06093 i__2 = i+1;
06094 for (j=half;j<=ny;++j) {
06095 j__1 = j-1;
06096 j__2 = j+1;
06097
06098 if((buf(i,j)>0.0f)&&buf(i,j)>buf(i,j__1)) {
06099 if(buf(i,j)>buf(i,j__2)) {
06100 if(buf(i,j)>buf(i__1,j)) {
06101 if(buf(i,j)>buf(i__2,j)) {
06102 if(buf(i,j)>buf(i__1,j__1)) {
06103 if((buf(i,j))> buf(i__1,j__2)) {
06104 if(buf(i,j)>buf(i__2,j__1)) {
06105 if(buf(i,j)> buf(i__2,j__2)) {
06106
06107
06108
06109 if (peaks.size()==0) {
06110
06111 peaks.push_back(Pixel(i,j,0,buf(i,j)));
06112
06113 } else {
06114
06115
06116 bool overlap = false;
06117
06118
06119
06120
06121
06122 std::stack <vector<Pixel>::iterator> delete_stack;
06123
06124
06125 for (vector<Pixel>::iterator it=peaks.begin();it!=peaks.end();++it) {
06126
06127
06128
06129
06130 float radius=((*it).x-float(i))*((*it).x-float(i))+((*it).y-float(j))*((*it).y-float(j));
06131 if (radius <= hf_p2 ) {
06132
06133 if( buf(i,j) > (*it).value) {
06134
06135
06136
06137
06138
06139
06140 delete_stack.push(it);
06141
06142
06143
06144
06145 } else {
06146 overlap = true;
06147
06148
06149 break;
06150 }
06151 }
06152 }
06153
06154
06155
06156 if (false == overlap) {
06157 vector<Pixel>::iterator delete_iterator;
06158 while (!delete_stack.empty()) {
06159
06160
06161 delete_iterator = delete_stack.top();
06162 peaks.erase(delete_iterator);
06163 delete_stack.pop();
06164 }
06165
06166
06167
06168 if (! (peaks.size() < 2000 )) {
06169
06170
06171
06172
06173 sort(peaks.begin(), peaks.end(), peakcmp);
06174
06175
06176 peaks.pop_back();
06177 }
06178
06179
06180 peaks.push_back(Pixel(i,j,0,buf(i,j)));
06181
06182
06183 } else {
06184
06185
06186 while (!delete_stack.empty()) delete_stack.pop();
06187 }
06188 }
06189 }
06190 }}}}}}}
06191 }
06192 }
06193
06194
06195 if(peaks.size()>0) {
06196
06197 sort(peaks.begin(),peaks.end(), peakcmp);
06198
06199 for (vector<Pixel>::iterator it = peaks.begin(); it != peaks.end(); it++) {
06200
06201 res.push_back((*it).value);
06202 res.push_back(static_cast<float>((*it).x));
06203 res.push_back(static_cast<float>((*it).y));
06204 }
06205 } else {
06206
06207 res.push_back(buf(0,0,0));
06208 res.insert(res.begin(),1,0.0);
06209 }
06210 return res;
06211 }
06212
06213 EMData* EMData::get_pow(float n_pow)
06214 {
06215 EMData* buf_new = this->copy_head();
06216 float *in = this->get_data();
06217 float *out = buf_new->get_data();
06218 for(size_t i=0; i<(size_t)nx*ny*nz; ++i) out[i] = pow(in[i],n_pow);
06219 return buf_new;
06220 }
06221
06222 EMData* EMData::conjg()
06223 {
06224 if(this->is_complex()) {
06225 EMData* buf_new = this->copy_head();
06226 float *in = this->get_data();
06227 float *out = buf_new->get_data();
06228 for(size_t i=0; i<(size_t)nx*ny*nz; i+=2) {out[i] = in[i]; out[i+1] = -in[i+1];}
06229 return buf_new;
06230 } else throw ImageFormatException("image has to be complex");
06231 }
06232
06233 EMData* EMData::delete_disconnected_regions(int ix, int iy, int iz) {
06234 if (3 != get_ndim())
06235 throw ImageDimensionException("delete_disconnected_regions needs a 3-D image.");
06236 if (is_complex())
06237 throw ImageFormatException("delete_disconnected_regions requires a real image");
06238 if ((*this)(ix+nx/2,iy+ny/2,iz+nz/2) == 0)
06239 throw ImageDimensionException("delete_disconnected_regions starting point is zero.");
06240
06241 EMData* result = this->copy_head();
06242 result->to_zero();
06243 (*result)(ix+nx/2,iy+ny/2,iz+nz/2) = (*this)(ix+nx/2,iy+ny/2,iz+nz/2);
06244 bool kpt = true;
06245
06246 while(kpt) {
06247 kpt = false;
06248 for (int cz = 1; cz < nz-1; cz++) {
06249 for (int cy = 1; cy < ny-1; cy++) {
06250 for (int cx = 1; cx < nx-1; cx++) {
06251 if((*result)(cx,cy,cz) == 1) {
06252 for (int lz = -1; lz <= 1; lz++) {
06253 for (int ly = -1; ly <= 1; ly++) {
06254 for (int lx = -1; lx <= 1; lx++) {
06255 if(((*this)(cx+lx,cy+ly,cz+lz) == 1) && ((*result)(cx+lx,cy+ly,cz+lz) == 0)) {
06256 (*result)(cx+lx,cy+ly,cz+lz) = 1;
06257 kpt = true;
06258 }
06259 }
06260 }
06261 }
06262 }
06263 }
06264 }
06265 }
06266 }
06267 result->update();
06268 return result;
06269 }
06270
06271 #define QUADPI 3.141592653589793238462643383279502884197
06272 #define DGR_TO_RAD QUADPI/180
06273
06274 EMData* EMData::helicise(float pixel_size, float dp, float dphi, float section_use, float radius, float minrad) {
06275 if (3 != get_ndim())
06276 throw ImageDimensionException("helicise needs a 3-D image.");
06277 if (is_complex())
06278 throw ImageFormatException("helicise requires a real image");
06279
06280 EMData* result = this->copy_head();
06281 result->to_zero();
06282 int nyc = ny/2;
06283 int nxc = nx/2;
06284 int nb = int(nz*(1.0f - section_use)/2.);
06285 int ne = nz - nb -1;
06286 int numst = int((ne - nb)/dp*pixel_size + 0.5);
06287
06288 int nst = int(nz*pixel_size/dp+0.5);
06289 float r2, ir;
06290 if(radius < 0.0f) r2 = (float)((nxc-1)*(nxc-1));
06291 else r2 = radius*radius;
06292 if(minrad < 0.0f) ir = 0.0f;
06293 else ir = minrad*minrad;
06294 for (int k = 0; k<nz; k++) {
06295 for (int j = 0; j<ny; j++) {
06296 int jy = j - nyc;
06297 int jj = jy*jy;
06298 for (int i = 0; i<nx; i++) {
06299 int ix = i - nxc;
06300 float d2 = (float)(ix*ix + jj);
06301 if(d2 <= r2 && d2>=ir) {
06302 int nq = 1;
06303 for ( int ist = -nst; ist <= nst; ist++) {
06304 float zold = (k*pixel_size + ist*dp)/pixel_size;
06305 int IOZ = int(zold);
06306 if(IOZ >= nb && IOZ <= ne) {
06307
06308 float cphi = ist*dphi*(float)DGR_TO_RAD;
06309 float ca = cos(cphi);
06310 float sa = sin(cphi);
06311 float xold = ix*ca - jy*sa + nxc;
06312 float yold = ix*sa + jy*ca + nyc;
06313 nq++;
06314
06315
06316
06317 int IOX = int(xold);
06318 int IOY = int(yold);
06319
06320
06321 #ifdef _WIN32
06322 int IOXp1 = _cpp_min( nx-1 ,IOX+1);
06323 #else
06324 int IOXp1 = std::min( nx-1 ,IOX+1);
06325 #endif //_WIN32
06326
06327 #ifdef _WIN32
06328 int IOYp1 = _cpp_min( ny-1 ,IOY+1);
06329 #else
06330 int IOYp1 = std::min( ny-1 ,IOY+1);
06331 #endif //_WIN32
06332
06333 #ifdef _WIN32
06334 int IOZp1 = _cpp_min( nz-1 ,IOZ+1);
06335 #else
06336 int IOZp1 = std::min( nz-1 ,IOZ+1);
06337 #endif //_WIN32
06338
06339 float dx = xold-IOX;
06340 float dy = yold-IOY;
06341 float dz = zold-IOZ;
06342
06343 float a1 = (*this)(IOX,IOY,IOZ);
06344 float a2 = (*this)(IOXp1,IOY,IOZ) - (*this)(IOX,IOY,IOZ);
06345 float a3 = (*this)(IOX,IOYp1,IOZ) - (*this)(IOX,IOY,IOZ);
06346 float a4 = (*this)(IOX,IOY,IOZp1) - (*this)(IOX,IOY,IOZ);
06347 float a5 = (*this)(IOX,IOY,IOZ) - (*this)(IOXp1,IOY,IOZ) - (*this)(IOX,IOYp1,IOZ) + (*this)(IOXp1,IOYp1,IOZ);
06348 float a6 = (*this)(IOX,IOY,IOZ) - (*this)(IOXp1,IOY,IOZ) - (*this)(IOX,IOY,IOZp1) + (*this)(IOXp1,IOY,IOZp1);
06349 float a7 = (*this)(IOX,IOY,IOZ) - (*this)(IOX,IOYp1,IOZ) - (*this)(IOX,IOY,IOZp1) + (*this)(IOX,IOYp1,IOZp1);
06350 float a8 = (*this)(IOXp1,IOY,IOZ) + (*this)(IOX,IOYp1,IOZ)+ (*this)(IOX,IOY,IOZp1)
06351 - (*this)(IOX,IOY,IOZ)- (*this)(IOXp1,IOYp1,IOZ) - (*this)(IOXp1,IOY,IOZp1)
06352 - (*this)(IOX,IOYp1,IOZp1) + (*this)(IOXp1,IOYp1,IOZp1);
06353 (*result)(i,j,k) += a1 + dz*(a4 + a6*dx + (a7 + a8*dx)*dy) + a3*dy + dx*(a2 + a5*dy);
06354
06355
06356 if(nq == numst) break;
06357 }
06358 }
06359 if(nq != numst)
06360 throw InvalidValueException(nq, "incorrect number of repeats encoutered.");
06361 }
06362 }
06363 }
06364 }
06365 for (int k = 0; k<nz; k++) for (int j = 0; j<ny; j++) for (int i = 0; i<nx; i++) (*result)(i,j,k) /= numst ;
06366
06367 result->update();
06368 return result;
06369 }
06370
06371
06372 EMData* EMData::helicise_rect(float pixel_size, float dp, float dphi, float section_use, float radius, float minrad) {
06373 if (3 != get_ndim())
06374 throw ImageDimensionException("helicise needs a 3-D image.");
06375 if (is_complex())
06376 throw ImageFormatException("helicise requires a real image");
06377
06378 EMData* result = this->copy_head();
06379 result->to_zero();
06380 int nyc = ny/2;
06381 int nxc = nx/2;
06382 int nb = int(nz*(1.0f - section_use)/2.);
06383 int ne = nz - nb -1;
06384 int numst = int((ne - nb)/dp*pixel_size + 0.5);
06385
06386 int nst = int(nz*pixel_size/dp+0.5);
06387 float r2, ir;
06388 if(radius < 0.0f) r2 = (float)((nxc-1)*(nxc-1));
06389 else r2 = radius*radius;
06390 if(minrad < 0.0f) ir = 0.0f;
06391 else ir = minrad*minrad;
06392 for (int k = 0; k<nz; k++) {
06393 for (int j = 0; j<ny; j++) {
06394 int jy = j - nyc;
06395 int jj = jy*jy;
06396 for (int i = 0; i<nx; i++) {
06397 int ix = i - nxc;
06398 float d2 = (float)(ix*ix + jj);
06399 if(d2 <= r2 && d2>=ir) {
06400 int nq = 1;
06401 for ( int ist = -nst; ist <= nst; ist++) {
06402 float zold = (k*pixel_size + ist*dp)/pixel_size;
06403 int IOZ = int(zold);
06404 if(IOZ >= nb && IOZ <= ne) {
06405
06406 float cphi = ist*dphi*(float)DGR_TO_RAD;
06407 float ca = cos(cphi);
06408 float sa = sin(cphi);
06409 float xold = ix*ca - jy*sa + nxc;
06410 float yold = ix*sa + jy*ca + nyc;
06411 nq++;
06412
06413
06414
06415 int IOX = int(xold);
06416 int IOY = int(yold);
06417
06418
06419 #ifdef _WIN32
06420 int IOXp1 = _cpp_min( nx-1 ,IOX+1);
06421 #else
06422 int IOXp1 = std::min( nx-1 ,IOX+1);
06423 #endif //_WIN32
06424
06425 #ifdef _WIN32
06426 int IOYp1 = _cpp_min( ny-1 ,IOY+1);
06427 #else
06428 int IOYp1 = std::min( ny-1 ,IOY+1);
06429 #endif //_WIN32
06430
06431 #ifdef _WIN32
06432 int IOZp1 = _cpp_min( nz-1 ,IOZ+1);
06433 #else
06434 int IOZp1 = std::min( nz-1 ,IOZ+1);
06435 #endif //_WIN32
06436
06437 float dx = xold-IOX;
06438 float dy = yold-IOY;
06439 float dz = zold-IOZ;
06440
06441 float a1 = (*this)(IOX,IOY,IOZ);
06442 float a2 = (*this)(IOXp1,IOY,IOZ) - (*this)(IOX,IOY,IOZ);
06443 float a3 = (*this)(IOX,IOYp1,IOZ) - (*this)(IOX,IOY,IOZ);
06444 float a4 = (*this)(IOX,IOY,IOZp1) - (*this)(IOX,IOY,IOZ);
06445 float a5 = (*this)(IOX,IOY,IOZ) - (*this)(IOXp1,IOY,IOZ) - (*this)(IOX,IOYp1,IOZ) + (*this)(IOXp1,IOYp1,IOZ);
06446 float a6 = (*this)(IOX,IOY,IOZ) - (*this)(IOXp1,IOY,IOZ) - (*this)(IOX,IOY,IOZp1) + (*this)(IOXp1,IOY,IOZp1);
06447 float a7 = (*this)(IOX,IOY,IOZ) - (*this)(IOX,IOYp1,IOZ) - (*this)(IOX,IOY,IOZp1) + (*this)(IOX,IOYp1,IOZp1);
06448 float a8 = (*this)(IOXp1,IOY,IOZ) + (*this)(IOX,IOYp1,IOZ)+ (*this)(IOX,IOY,IOZp1)
06449 - (*this)(IOX,IOY,IOZ)- (*this)(IOXp1,IOYp1,IOZ) - (*this)(IOXp1,IOY,IOZp1)
06450 - (*this)(IOX,IOYp1,IOZp1) + (*this)(IOXp1,IOYp1,IOZp1);
06451 (*result)(i,j,k) += a1 + dz*(a4 + a6*dx + (a7 + a8*dx)*dy) + a3*dy + dx*(a2 + a5*dy);
06452
06453
06454 if(nq == numst) break;
06455 }
06456 }
06457 if(nq != numst)
06458 throw InvalidValueException(nq, "incorrect number of repeats encoutered.");
06459 }
06460 }
06461 }
06462 }
06463 for (int k = 0; k<nz; k++) for (int j = 0; j<ny; j++) for (int i = 0; i<nx; i++) (*result)(i,j,k) /= numst ;
06464
06465 result->update();
06466 return result;
06467 }
06468
06469
06470
06471
06472
06473
06474
06475
06476
06477
06478 void EMData::depad() {
06479 if (is_complex())
06480 throw ImageFormatException("Depadding of complex images not supported");
06481 vector<int> saved_offsets = get_array_offsets();
06482 set_array_offsets(0,0,0);
06483 int npad = attr_dict["npad"];
06484 if (0 == npad) npad = 1;
06485 int offset = is_fftodd() ? 1 : 2;
06486 int nxold = (nx - offset)/npad;
06487 #ifdef _WIN32
06488 int nyold = _cpp_max(ny/npad, 1);
06489 int nzold = _cpp_max(nz/npad, 1);
06490 #else
06491 int nyold = std::max<int>(ny/npad, 1);
06492 int nzold = std::max<int>(nz/npad, 1);
06493 #endif //_WIN32
06494 int xstart = 0, ystart = 0, zstart = 0;
06495 if( npad > 1) {
06496 xstart = (nx - offset - nxold)/2 + nxold%2;
06497 if(ny > 1) {
06498 ystart = (ny - nyold)/2 + nyold%2;
06499 if(nz > 1) {
06500 zstart = (nz - nzold)/2 + nzold%2;
06501 }
06502 }
06503 }
06504 int bytes = nxold*sizeof(float);
06505 float* dest = get_data();
06506 for (int iz=0; iz < nzold; iz++) {
06507 for (int iy = 0; iy < nyold; iy++) {
06508 memmove(dest, &(*this)(xstart,iy+ystart,iz+zstart), bytes);
06509 dest += nxold;
06510 }
06511 }
06512 set_size(nxold, nyold, nzold);
06513 set_attr("npad", 1);
06514 set_fftpad(false);
06515 set_fftodd(false);
06516 set_complex(false);
06517 if(ny==1 && nz==1) set_complex_x(false);
06518 set_array_offsets(saved_offsets);
06519 update();
06520 EXITFUNC;
06521 }
06522
06523
06524
06525
06526
06527
06528
06529
06530 void EMData::depad_corner() {
06531 if(is_complex())
06532 throw ImageFormatException("Depadding of complex images not allowed");
06533 vector<int> saved_offsets = get_array_offsets();
06534 set_array_offsets(0,0,0);
06535 int npad = attr_dict["npad"];
06536 if(0 == npad) npad = 1;
06537 int offset = is_fftodd() ? 1 : 2;
06538 int nxold = (nx - offset)/npad;
06539 #ifdef _WIN32
06540 int nyold = _cpp_max(ny/npad, 1);
06541 int nzold = _cpp_max(nz/npad, 1);
06542 #else
06543 int nyold = std::max<int>(ny/npad, 1);
06544 int nzold = std::max<int>(nz/npad, 1);
06545 #endif //_WIN32
06546 size_t bytes = nxold*sizeof(float);
06547 float* dest = get_data();
06548 for (int iz=0; iz < nzold; iz++) {
06549 for (int iy = 0; iy < nyold; iy++) {
06550 memmove(dest, &(*this)(0,iy,iz), bytes);
06551 dest += nxold;
06552 }
06553 }
06554 set_size(nxold, nyold, nzold);
06555 set_attr("npad", 1);
06556 set_fftpad(false);
06557 set_fftodd(false);
06558 set_complex(false);
06559 if(ny==1 && nz==1) set_complex_x(false);
06560 set_array_offsets(saved_offsets);
06561 update();
06562 EXITFUNC;
06563 }
06564
06565
06566
06567 float circumference( EMData* emdata, int npixel )
06568 {
06569 int nx = emdata->get_xsize();
06570 int ny = emdata->get_ysize();
06571 int nz = emdata->get_zsize();
06572
06573 float* data = emdata->get_data();
06574 if( ny==1 && nz==1 ) {
06575
06576 float sumf=0.0;
06577 int sumn=0;
06578 for( int i=0; i < npixel; ++i ) {
06579 sumf += data[i];
06580 sumf += data[nx-1-i];
06581 sumn += 2;
06582 }
06583 return sumf/sumn;
06584 }
06585
06586 if( nz==1 ) {
06587 float sumf=0.0;
06588 int sumn=0;
06589 int id=0;
06590 for( int iy=0; iy < ny; ++iy ) {
06591 for( int ix=0; ix < nx; ++ix ) {
06592 if( iy<npixel || iy>ny-1-npixel || ix<npixel || ix>nx-1-npixel ) {
06593 sumf += data[id];
06594 sumn += 1;
06595 }
06596 id++;
06597 }
06598 }
06599
06600 Assert( id==nx*ny );
06601 Assert( sumn == nx*ny - (nx-2*npixel)*(ny-2*npixel) );
06602 return sumf/sumn;
06603 }
06604
06605
06606
06607 float sumf = 0.0;
06608 size_t sumn = 0;
06609 size_t id = 0;
06610 for( int iz=0; iz < nz; ++iz) {
06611 for( int iy=0; iy < ny; ++iy) {
06612 for( int ix=0; ix < nx; ++ix ) {
06613 if( iz<npixel||iz>nz-1-npixel||iy<npixel||iy>ny-1-npixel||ix<npixel||ix>nx-1-npixel) {
06614 sumf += data[id];
06615 sumn += 1;
06616 }
06617 id++;
06618 }
06619 }
06620 }
06621
06622
06623 Assert( id==(size_t)nx*ny*nz);
06624 Assert( sumn==(size_t)nx*ny*nz-(size_t)(nx-2*npixel)*(ny-2*npixel)*(nz-2*npixel) );
06625 return sumf/sumn;
06626 }
06627
06628
06629
06630
06631
06632
06633
06634
06635 EMData* EMData::norm_pad(bool donorm, int npad, int valtype) {
06636 if (this->is_complex())
06637 throw ImageFormatException("Padding of complex images not supported");
06638 int nx = this->get_xsize();
06639 int ny = this->get_ysize();
06640 int nz = this->get_zsize();
06641 float mean = 0., stddev = 1.;
06642 if(donorm) {
06643 mean = this->get_attr("mean");
06644 stddev = this->get_attr("sigma");
06645 }
06646
06647 if (npad < 1) npad = 1;
06648 int nxpad = npad*nx;
06649 int nypad = npad*ny;
06650 int nzpad = npad*nz;
06651 if (1 == ny) {
06652
06653
06654 nypad = ny;
06655 nzpad = nz;
06656 } else if (nz == 1) {
06657
06658 nzpad = nz;
06659 }
06660 size_t bytes;
06661 size_t offset;
06662
06663 offset = 2 - nxpad%2;
06664 bytes = nx*sizeof(float);
06665 EMData* fpimage = copy_head();
06666 fpimage->set_size(nxpad+offset, nypad, nzpad);
06667 int xstart = 0, ystart = 0, zstart = 0;
06668 if( npad > 1) {
06669 if( valtype==0 ) {
06670 fpimage->to_zero();
06671 } else {
06672 float val = circumference(this, 1);
06673 float* data = fpimage->get_data();
06674 int nxyz = (nxpad+offset)*nypad*nzpad;
06675 for( int i=0; i < nxyz; ++i ) data[i] = val;
06676 }
06677
06678 xstart = (nxpad - nx)/2 + nx%2;
06679 if(ny > 1) {
06680 ystart = (nypad - ny)/2 + ny%2;
06681 if(nz > 1) {
06682 zstart = (nzpad - nz)/2 + nz%2;
06683 }
06684 }
06685 }
06686
06687
06688 vector<int> saved_offsets = this->get_array_offsets();
06689 this->set_array_offsets( 0, 0, 0 );
06690 for (int iz = 0; iz < nz; iz++) {
06691 for (int iy = 0; iy < ny; iy++) {
06692 memcpy(&(*fpimage)(xstart,iy+ystart,iz+zstart), &(*this)(0,iy,iz), bytes);
06693 }
06694 }
06695 this->set_array_offsets( saved_offsets );
06696
06697
06698
06699
06700 if (donorm) {
06701 for (int iz = zstart; iz < nz+zstart; iz++)
06702 for (int iy = ystart; iy < ny+ystart; iy++)
06703 for (int ix = xstart; ix < nx+xstart; ix++)
06704 (*fpimage)(ix,iy,iz) = ((*fpimage)(ix,iy,iz)-mean)/stddev;
06705 }
06706
06707 fpimage->set_fftpad(true);
06708 fpimage->set_attr("npad", npad);
06709 if (offset == 1) fpimage->set_fftodd(true);
06710 else fpimage->set_fftodd(false);
06711 return fpimage;
06712 }
06713
06714 void EMData::center_origin()
06715 {
06716 ENTERFUNC;
06717 if (is_complex()) {
06718 LOGERR("Real image expected. Input image is complex.");
06719 throw ImageFormatException("Real image expected. Input image is complex.");
06720 }
06721 for (int iz = 0; iz < nz; iz++) {
06722 for (int iy = 0; iy < ny; iy++) {
06723 for (int ix = 0; ix < nx; ix++) {
06724
06725 (*this)(ix,iy,iz) *= -2*((ix+iy+iz)%2) + 1;
06726 }
06727 }
06728 }
06729 update();
06730 EXITFUNC;
06731 }
06732
06733 void EMData::center_origin_yz()
06734 {
06735 ENTERFUNC;
06736 if (is_complex()) {
06737 LOGERR("Real image expected. Input image is complex.");
06738 throw ImageFormatException("Real image expected. Input image is complex.");
06739 }
06740 for (int iz = 0; iz < nz; iz++) {
06741 for (int iy = (iz+1)%2; iy < ny; iy+=2) {
06742 for (int ix = 0; ix < nx; ix++) {
06743 (*this)(ix,iy,iz) *= -1;
06744 }
06745 }
06746 }
06747 update();
06748 EXITFUNC;
06749 }
06750
06751 void EMData::center_origin_fft()
06752 {
06753 ENTERFUNC;
06754 if (!is_complex()) {
06755 LOGERR("complex image expected. Input image is real image.");
06756 throw ImageFormatException("complex image expected. Input image is real image.");
06757 }
06758
06759 if (!is_ri()) {
06760 LOGWARN("Only RI should be used. ");
06761 }
06762 vector<int> saved_offsets = get_array_offsets();
06763
06764
06765
06766 set_array_offsets(0,1,1);
06767 int nxc = nx/2;
06768
06769 if (is_fftodd()) {
06770 for (int iz = 1; iz <= nz; iz++) {
06771 for (int iy = 1; iy <= ny; iy++) {
06772 for (int ix = 0; ix < nxc; ix++) {
06773 cmplx(ix,iy,iz) *= float(-2*((ix+iy+iz)%2) + 1);
06774 float temp = float(iz-1+iy-1+ix)/float(ny)*M_PI;
06775 complex<float> temp2 = complex<float>(cos(temp), -sin(temp));
06776 cmplx(ix,iy,iz) *= temp2;
06777 }
06778 }
06779 }
06780 } else {
06781 for (int iz = 1; iz <= nz; iz++) {
06782 for (int iy = 1; iy <= ny; iy++) {
06783 for (int ix = 0; ix < nxc; ix++) {
06784
06785 cmplx(ix,iy,iz) *= float(-2*((ix+iy+iz)%2) + 1);
06786 }
06787 }
06788 }
06789 }
06790 set_array_offsets(saved_offsets);
06791 update();
06792 EXITFUNC;
06793 }
06794
06795
06796 #define fint(i,j,k) fint[(i-1) + ((j-1) + (k-1)*ny)*(size_t)lsd]
06797 #define fout(i,j,k) fout[(i-1) + ((j-1) + (k-1)*nyn)*(size_t)lsdn]
06798 EMData *EMData::FourInterpol(int nxn, int nyni, int nzni, bool RetReal) {
06799
06800 int nyn, nzn, lsd, lsdn, inx, iny, inz;
06801 int i, j, k;
06802 if (is_complex())
06803 throw ImageFormatException("Input image has to be real");
06804
06805 if(ny > 1) {
06806 nyn = nyni;
06807 if(nz > 1) {
06808 nzn = nzni;
06809 } else {
06810 nzn = 1;
06811 }
06812 } else {
06813 nyn = 1; nzn = 1;
06814 }
06815 if(nxn<nx || nyn<ny || nzn<nz) throw ImageDimensionException("Cannot reduce the image size");
06816 lsd = nx + 2 - nx%2;
06817 lsdn = nxn + 2 - nxn%2;
06818
06819 EMData *temp_ft = do_fft();
06820 EMData *ret = this->copy();
06821 ret->set_size(lsdn, nyn, nzn);
06822 ret->to_zero();
06823 float *fout = ret->get_data();
06824 float *fint = temp_ft->get_data();
06825
06826
06827 float sq2 = 1.0f/std::sqrt(2.0f);
06828 float anorm = (float) nxn* (float) nyn* (float) nzn/(float) nx/ (float) ny/ (float) nz;
06829 for (i = 0; i < lsd*ny*nz; i++) fint[i] *= anorm;
06830 inx = nxn-nx; iny = nyn - ny; inz = nzn - nz;
06831 for (k=1; k<=nz/2+1; k++) for (j=1; j<=ny/2+1; j++) for (i=1; i<=lsd; i++) fout(i,j,k)=fint(i,j,k);
06832 if(nyn>1) {
06833
06834 for (k=1; k<=nz/2+1; k++) for (j=ny/2+2+iny; j<=nyn; j++) for (i=1; i<=lsd; i++) fout(i,j,k)=fint(i,j-iny,k);
06835 if(nzn>1) {
06836 for (k=nz/2+2+inz; k<=nzn; k++) {
06837 for (j=1; j<=ny/2+1; j++) {
06838 for (i=1; i<=lsd; i++) {
06839 fout(i,j,k)=fint(i,j,k-inz);
06840 }
06841 }
06842 for (j=ny/2+2+iny; j<=nyn; j++) {
06843 for (i=1; i<=lsd; i++) {
06844 fout(i,j,k)=fint(i,j-iny,k-inz);
06845 }
06846 }
06847 }
06848 }
06849 }
06850
06851
06852
06853 if(nx%2 == 0 && inx !=0) {
06854 for (k=1; k<=nzn; k++) {
06855 for (j=1; j<=nyn; j++) {
06856 fout(nx+1,j,k) *= sq2;
06857 fout(nx+2,j,k) *= sq2;
06858 }
06859 }
06860 if(nyn>1) {
06861 for (k=1; k<=nzn; k++) {
06862 for (i=1; i<=lsd; i++) {
06863 fout(i,ny/2+1+iny,k) = sq2*fout(i,ny/2+1,k);
06864 fout(i,ny/2+1,k) *= sq2;
06865 }
06866 }
06867 if(nzn>1) {
06868 for (j=1; j<=nyn; j++) {
06869 for (i=1; i<=lsd; i++) {
06870 fout(i,j,nz/2+1+inz) = sq2*fout(i,j,nz/2+1);
06871 fout(i,j,nz/2+1) *= sq2;
06872 }
06873 }
06874 }
06875 }
06876 }
06877 ret->set_complex(true);
06878
06879
06880
06881
06882
06883
06884
06885
06886
06887
06888
06889
06890
06891
06892
06893
06894
06895
06896
06897
06898
06899
06900
06901
06902
06903
06904
06905
06906
06907
06908
06909
06910
06911 ret->set_ri(1);
06912 ret->set_fftpad(true);
06913 ret->set_attr("npad", 1);
06914 if (nxn%2 == 1) {ret->set_fftodd(true);} else {ret->set_fftodd(false);}
06915 if(RetReal) {
06916 ret->do_ift_inplace();
06917 ret->depad();
06918 }
06919 ret->update();
06920
06921
06922
06923
06924
06925
06926
06927 delete temp_ft;
06928 temp_ft = 0;
06929 return ret;
06930 }
06931
06932 EMData *EMData::FourTruncate(int nxn, int nyni, int nzni, bool RetReal) {
06933
06934 int nyn, nzn, lsd, lsdn, inx, iny, inz;
06935 int i, j, k;
06936 float *fint;
06937 EMData *temp_ft = NULL;
06938
06939
06940
06941 if(ny > 1) {
06942 nyn = nyni;
06943 if(nz > 1) {
06944 nzn = nzni;
06945 } else {
06946 nzn = 1;
06947 }
06948 } else {
06949 nyn = 1; nzn = 1;
06950 }
06951 if (is_complex()) {
06952 nx = nx - 2 + nx%2;
06953 fint = get_data();
06954 } else {
06955
06956 temp_ft = do_fft();
06957 fint = temp_ft->get_data();
06958 }
06959 if(nxn>nx || nyn>ny || nzn>nz) throw ImageDimensionException("Cannot increase the image size");
06960 lsd = nx + 2 - nx%2;
06961 lsdn = nxn + 2 - nxn%2;
06962 EMData *ret = this->copy_head();
06963 ret->set_size(lsdn, nyn, nzn);
06964 float *fout = ret->get_data();
06965
06966
06967
06968 float anorm = (float) nxn* (float) nyn* (float) nzn/(float) nx/ (float) ny/ (float) nz;
06969 for (i = 0; i < lsd*ny*nz; i++) fint[i] *= anorm;
06970 inx = nx - nxn; iny = ny - nyn; inz = nz - nzn;
06971 for (k=1; k<=nzn/2+1; k++) for (j=1; j<=nyn/2+1; j++) for (i=1; i<=lsdn; i++) fout(i,j,k)=fint(i,j,k);
06972 if(nyn>1) {
06973 for (k=1; k<=nzn/2+1; k++) for (j=nyn/2+2; j<=nyn; j++) for (i=1; i<=lsdn; i++) fout(i,j,k)=fint(i,j+iny,k);
06974 if(nzn>1) {
06975 for (k=nzn/2+2; k<=nzn; k++) {
06976 for (j=1; j<=nyn/2+1; j++) {
06977 for (i=1; i<=lsdn; i++) {
06978 fout(i,j,k)=fint(i,j,k+inz);
06979 }
06980 }
06981 for (j=nyn/2+2; j<=nyn; j++) {
06982 for (i=1; i<=lsdn; i++) {
06983 fout(i,j,k)=fint(i,j+iny,k+inz);
06984 }
06985 }
06986 }
06987 }
06988 }
06989
06990
06991
06992
06993
06994
06995
06996
06997
06998
06999
07000
07001
07002
07003
07004
07005
07006
07007
07008
07009
07010
07011
07012
07013
07014
07015
07016
07017 ret->set_complex(true);
07018 ret->set_ri(1);
07019 ret->set_fftpad(true);
07020 ret->set_attr("npad", 1);
07021 if (nxn%2 == 1) {ret->set_fftodd(true);} else {ret->set_fftodd(false);}
07022 if(RetReal) {
07023 ret->do_ift_inplace();
07024 ret->depad();
07025 }
07026 ret->update();
07027
07028
07029
07030
07031
07032
07033
07034 if (!is_complex()) {
07035 delete temp_ft;
07036 temp_ft = 0;
07037 }
07038 return ret;
07039 }
07040
07041
07042
07043
07044
07045
07046
07047
07048
07049
07050
07051
07052
07053
07054
07055
07056
07057
07058
07059
07060
07061
07062
07063
07064
07065
07066
07067
07068
07069
07070
07071
07072
07073
07074
07075
07076
07077
07078
07079
07080
07081
07082
07083
07084
07085
07086
07087
07088
07089
07090
07091
07092
07093
07094
07095
07096
07097
07098
07099
07100
07101
07102
07103
07104
07105
07106
07107
07108
07109
07110
07111
07112
07113
07114
07115
07116
07117
07118
07119
07120
07121
07122
07123
07124
07125
07126
07127
07128
07129
07130
07131
07132
07133
07134
07135 EMData *EMData::Four_ds(int nxn, int nyni, int nzni, bool RetReal) {
07136
07137 int nyn, nzn, lsd, lsdn, inx, iny, inz;
07138 int i, j;
07139
07140 if(ny > 1) {
07141 nyn = nyni;
07142 if(nz > 1) {
07143 nzn = nzni;
07144 } else {
07145 nzn = 1;
07146 }
07147 } else {
07148 nyn = 1; nzn = 1;
07149 }
07150 lsd = nx-2 + 2 - nx%2;
07151 lsdn = nxn + 2 - nxn%2;
07152
07153 EMData *temp_ft = this->copy();
07154 EMData *ret = this->copy();
07155 ret->set_size(lsdn, nyn, nzn);
07156 ret->to_zero();
07157 float *fout = ret->get_data();
07158 float *fint = temp_ft->get_data();
07159
07160
07161
07162 float anorm = (float) nxn* (float) nyn* (float) nzn/(float) nx/ (float) ny/ (float) nz;
07163 for (i = 0; i < lsd*ny*nz; i++) fint[i] *= anorm;
07164 inx = nxn-(nx-2); iny = nyn - ny; inz = nzn - nz;
07165 for (j=1; j<=nyn; j++)
07166 for (i=1; i<=lsdn; i++)
07167 fout(i,j,1)=fint((i-1)/2*4+2-i%2,j*2-1,1);
07168 ret->set_complex(true);
07169 ret->set_ri(1);
07170
07171
07172 if (nxn%2 == 1) {ret->set_fftodd(true);} else {ret->set_fftodd(false);}
07173 if(RetReal) {
07174 ret->do_ift_inplace();
07175 ret->depad();
07176 }
07177 ret->update();
07178
07179 delete temp_ft;
07180 temp_ft = 0;
07181 return ret;
07182 }
07183
07184 EMData *EMData::Four_shuf_ds_cen_us(int nxn, int nyni, int, bool RetReal) {
07185
07186 int nyn, nzn, lsd, lsdn, inx, iny, inz;
07187 int i, j;
07188
07189 nyn = nyni;
07190 nzn = 1;
07191 lsd = nx;
07192 lsdn = nxn + 2 - nxn%2;
07193
07194 EMData *temp_ft = this->copy();
07195 EMData *ret = this->copy();
07196 ret->set_size(lsdn, nyn, nzn);
07197 ret->to_zero();
07198 float *fout = ret->get_data();
07199 float *fint = temp_ft->get_data();
07200
07201
07202 float sq2 = 1.0f/std::sqrt(2.0f);
07203
07204 for (size_t i = 0; i < (size_t)lsd*ny*nz; i++) fint[i] *= 4;
07205
07206 inx = nxn-(nx-2); iny = nyn - ny; inz = nzn - nz;
07207 for (j=1; j<=ny/4; j++)
07208 for (i=1; i<=(nx-2)/2+2; i++) {
07209 int g = (i-1)/2+1;
07210 if ((g+j)%2 == 0) {
07211 fout(i,j,1)=fint(g*4-2-i%2,j*2-1+ny/2,1);
07212 } else {
07213 fout(i,j,1)=-fint(g*4-2-i%2,j*2-1+ny/2,1);
07214 }
07215 }
07216
07217 for (j=ny/4+1; j<=ny/4+1; j++)
07218 for (i=1; i<=(nx-2)/2+2; i++) {
07219 int g = (i-1)/2+1;
07220 if ((g+j)%2 == 0) {
07221 fout(i,j,1)=fint(g*4-2-i%2,j*2-1-ny/2,1);
07222 } else {
07223 fout(i,j,1)=-fint(g*4-2-i%2,j*2-1-ny/2,1);
07224 }
07225 }
07226
07227 for (j=ny/4+2; j<=ny/2; j++)
07228 for (i=1; i<=(nx-2)/2+2; i++) {
07229 int g = (i-1)/2+1;
07230 if ((g+j)%2 == 0) {
07231 fout(i,j+ny/2,1)=fint(g*4-2-i%2,j*2-1-ny/2,1);
07232 } else {
07233 fout(i,j+ny/2,1)=-fint(g*4-2-i%2,j*2-1-ny/2,1);
07234 }
07235 }
07236
07237 if (nx%2 == 0) {
07238 for (j=1; j<=nyn; j++) {
07239 fout((nx-2)/2+1,j,1) *= sq2;
07240 fout((nx-2)/2+2,j,1) *= sq2;
07241 }
07242 for (i=1; i<=lsd/2+1; i++) {
07243 fout(i,ny/4+1+ny/2,1) = sq2*fout(i,ny/4+1,1);
07244 fout(i,ny/4+1,1) *= sq2;
07245 }
07246 }
07247
07248 ret->set_complex(true);
07249 ret->set_ri(1);
07250
07251 if (nxn%2 == 1) {ret->set_fftodd(true);} else {ret->set_fftodd(false);}
07252 if(RetReal) {
07253 ret->do_ift_inplace();
07254 ret->depad();
07255 }
07256 ret->update();
07257
07258 delete temp_ft;
07259 temp_ft = 0;
07260 return ret;
07261 }
07262
07263 #undef fint
07264 #undef fout
07265
07266 #define fint(jx,jy,jz) fint[jx + (jy + jz*ny)*(size_t)nox]
07267 EMData *EMData::filter_by_image(EMData* image, bool RetReal) {
07268
07269
07270 bool complex_input = this->is_complex();
07271 nx = this->get_xsize();
07272 ny = this->get_ysize();
07273 nz = this->get_zsize();
07274 int nox;
07275 if (complex_input) nox = (nx - 2 + this->is_fftodd()); else nox = nx;
07276
07277 int lsd2 = (nox + 2 - nox%2) / 2;
07278
07279 EMData* fp = NULL;
07280 if(complex_input) {
07281
07282 fp = this->copy();
07283 } else {
07284 fp = this->norm_pad( false, 1);
07285 fp->do_fft_inplace();
07286 }
07287 fp->set_array_offsets(1,1,1);
07288 int nx2 = nox/2;
07289 int ny2 = ny/2;
07290 int nz2 = nz/2;
07291 float *fint = image->get_data();
07292 for ( int iz = 1; iz <= nz; iz++) {
07293 int jz=nz2-iz+1; if(jz<0) jz += nz;
07294 for ( int iy = 1; iy <= ny; iy++) {
07295 int jy=ny2-iy+1; if(jy<0) jy += ny;
07296 for ( int ix = 1; ix <= lsd2; ix++) {
07297 int jx = nx2-ix+1;
07298 fp->cmplx(ix,iy,iz) *= fint(jx,jy,jz);
07299 }
07300 }
07301 }
07302
07303 fp->set_ri(1);
07304 fp->set_fftpad(true);
07305 fp->set_attr("npad", 1);
07306 if (nx%2 == 1) fp->set_fftodd(true);
07307 else fp->set_fftodd(false);
07308 if(RetReal) {
07309 fp->do_ift_inplace();
07310 fp->depad();
07311 }
07312 fp->set_array_offsets(0,0,0);
07313 fp->update();
07314
07315 return fp;
07316 }
07317 #undef fint
07318 #define fint(jx,jy,jz) fint[jx + (jy + jz*ny)*(size_t)nx]
07319 #define fout(jx,jy,jz) fout[jx + (jy + jz*ny)*(size_t)nx]
07320 EMData *EMData::replace_amplitudes(EMData* image, bool RetReal) {
07321
07322
07323 bool complex_input = this->is_complex();
07324 nx = this->get_xsize();
07325 ny = this->get_ysize();
07326 nz = this->get_zsize();
07327 int nox;
07328 if (complex_input) nox = (nx - 2 + this->is_fftodd()); else nox = nx;
07329
07330 EMData* fp = NULL;
07331 if(complex_input) {
07332
07333 fp = this->copy();
07334 } else {
07335 fp = this->norm_pad( false, 1);
07336 fp->do_fft_inplace();
07337 }
07338 float *fout = fp->get_data();
07339 float *fint = image->get_data();
07340 for ( int iz = 0; iz < nz; iz++) {
07341 for ( int iy = 0; iy < ny; iy++) {
07342 for ( int ix = 0; ix < nx; ix+=2) {
07343 float qt = fint(ix,iy,iz)*fint(ix,iy,iz)+fint(ix+1,iy,iz)*fint(ix+1,iy,iz);
07344 float rt = fout(ix,iy,iz)*fout(ix,iy,iz)+fout(ix+1,iy,iz)*fout(ix+1,iy,iz);
07345 if(rt > 1.0e-20) {
07346 fout(ix,iy,iz) *= (qt/rt);
07347 fout(ix+1,iy,iz) *= (qt/rt);
07348 } else {
07349 qt = std::sqrt(qt/2.0f);
07350 fout(ix,iy,iz) = qt;
07351 fout(ix+1,iy,iz) = qt;
07352 }
07353 }
07354 }
07355 }
07356
07357 fp->set_ri(1);
07358 fp->set_fftpad(true);
07359 fp->set_attr("npad", 1);
07360 if (nx%2 == 1) fp->set_fftodd(true);
07361 else fp->set_fftodd(false);
07362 if(RetReal) {
07363 fp->do_ift_inplace();
07364 fp->depad();
07365 }
07366 fp->set_array_offsets(0,0,0);
07367 fp->update();
07368
07369 return fp;
07370 }
07371 #undef fint
07372 #undef fout
07373
07374
07375 #undef QUADPI
07376 #undef DGR_TO_RAD