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 EMData* EMData::symvol(string symString) {
00984 ENTERFUNC;
00985 int nsym = Transform::get_nsym(symString);
00986 Transform sym;
00987
00988 EMData *svol = new EMData;
00989 svol->set_size(nx, ny, nz);
00990 svol->to_zero();
00991
00992 for (int isym = 0; isym < nsym; isym++) {
00993 Transform rm = sym.get_sym(symString, isym);
00994 EMData* symcopy = this -> rot_scale_trans(rm);
00995 *svol += (*symcopy);
00996 delete symcopy;
00997 }
00998 *svol /= ((float) nsym);
00999 svol->update();
01000 EXITFUNC;
01001 return svol;
01002 }
01003
01004 #define proj(ix,iy,iz) proj[ix-1+(iy-1+(iz-1)*ny)*(size_t)nx]
01005 #define pnewimg(ix,iy,iz) pnewimg[ix-1+(iy-1+(iz-1)*ny)*(size_t)nx]
01006 EMData* EMData::average_circ_sub() const
01007 {
01008
01009
01010 ENTERFUNC;
01011 const EMData* const image = this;
01012 EMData* newimg = copy_head();
01013 float *proj = image->get_data();
01014 float *pnewimg = newimg->get_data();
01015
01016 float r2 = static_cast<float>( (nx/2)*(nx/2) );
01017 float qs=0.0f;
01018 int m=0;
01019 int ncz = nz/2 + 1;
01020 int ncy = ny/2 + 1;
01021 int ncx = nx/2 + 1;
01022 for (int iz = 1; iz <= nz; iz++) {
01023 float yy = static_cast<float>( (iz-ncz)*(iz-ncz) );
01024 for (int iy = 1; iy <=ny; iy++) { float xx = yy + (iy-ncy)*(iy-ncy);
01025 for (int ix = 1; ix <= nx; ix++) {
01026 if ( xx+float((ix-ncx)*(ix-ncx)) > r2 ) {
01027 qs += proj(ix,iy,iz);
01028 m++;
01029 }
01030 }
01031 }
01032 }
01033
01034
01035 if( m > 0 ) qs /= m;
01036
01037 for (int iz = 1; iz <= nz; iz++)
01038 for (int iy = 1; iy <= ny; iy++)
01039 for (int ix = 1; ix <= nx; ix++)
01040 pnewimg(ix,iy,iz) = proj(ix,iy,iz) - qs;
01041 newimg->update();
01042 return newimg;
01043 EXITFUNC;
01044 }
01045
01046
01047
01048
01049
01050 void EMData::onelinenn(int j, int n, int n2, EMData* wptr, EMData* bi, const Transform& tf)
01051 {
01052
01053 int jp = (j >= 0) ? j+1 : n+j+1;
01054
01055
01056 for (int i = 0; i <= n2; i++) {
01057 if (((i*i+j*j) < n*n/4) && !((0 == i) && (j < 0))) {
01058
01059 float xnew = i*tf[0][0] + j*tf[1][0];
01060 float ynew = i*tf[0][1] + j*tf[1][1];
01061 float znew = i*tf[0][2] + j*tf[1][2];
01062 std::complex<float> btq;
01063 if (xnew < 0.) {
01064 xnew = -xnew;
01065 ynew = -ynew;
01066 znew = -znew;
01067 btq = conj(bi->cmplx(i,jp));
01068 } else {
01069 btq = bi->cmplx(i,jp);
01070 }
01071 int ixn = int(xnew + 0.5 + n) - n;
01072 int iyn = int(ynew + 0.5 + n) - n;
01073 int izn = int(znew + 0.5 + n) - n;
01074
01075 int iza, iya;
01076 if (izn >= 0) iza = izn + 1;
01077 else iza = n + izn + 1;
01078
01079 if (iyn >= 0) iya = iyn + 1;
01080 else iya = n + iyn + 1;
01081
01082 cmplx(ixn,iya,iza) += btq;
01083
01084 (*wptr)(ixn,iya,iza)++;
01085
01086
01087
01088
01089
01090
01091
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 void EMData::onelinenn_mult(int j, int n, int n2, EMData* wptr, EMData* bi, const Transform& tf, int mult)
01117 {
01118
01119 int jp = (j >= 0) ? j+1 : n+j+1;
01120
01121
01122 for (int i = 0; i <= n2; i++) {
01123 if (((i*i+j*j) < n*n/4) && !((0 == i) && (j < 0))) {
01124
01125 float xnew = i*tf[0][0] + j*tf[1][0];
01126 float ynew = i*tf[0][1] + j*tf[1][1];
01127 float znew = i*tf[0][2] + j*tf[1][2];
01128 std::complex<float> btq;
01129 if (xnew < 0.) {
01130 xnew = -xnew;
01131 ynew = -ynew;
01132 znew = -znew;
01133 btq = conj(bi->cmplx(i,jp));
01134 } else {
01135 btq = bi->cmplx(i,jp);
01136 }
01137 int ixn = int(xnew + 0.5 + n) - n;
01138 int iyn = int(ynew + 0.5 + n) - n;
01139 int izn = int(znew + 0.5 + n) - n;
01140
01141
01142 int iza, iya;
01143 if (izn >= 0) iza = izn + 1;
01144 else iza = n + izn + 1;
01145
01146 if (iyn >= 0) iya = iyn + 1;
01147 else iya = n + iyn + 1;
01148
01149 cmplx(ixn,iya,iza) += btq*float(mult);
01150 (*wptr)(ixn,iya,iza)+=float(mult);
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175 }
01176 }
01177 }
01178
01179 void EMData::nn(EMData* wptr, EMData* myfft, const Transform& tf, int mult)
01180 {
01181 ENTERFUNC;
01182 int nxc = attr_dict["nxc"];
01183
01184 vector<int> saved_offsets = get_array_offsets();
01185 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01186 set_array_offsets(0,1,1);
01187 myfft->set_array_offsets(0,1);
01188
01189
01190
01191
01192 if( mult == 1 ) {
01193 for (int iy = -ny/2 + 1; iy <= ny/2; iy++) onelinenn(iy, ny, nxc, wptr, myfft, tf);
01194 } else {
01195 for (int iy = -ny/2 + 1; iy <= ny/2; iy++) onelinenn_mult(iy, ny, nxc, wptr, myfft, tf, mult);
01196 }
01197
01198 set_array_offsets(saved_offsets);
01199 myfft->set_array_offsets(myfft_saved_offsets);
01200 EXITFUNC;
01201 }
01202
01203
01204 void EMData::insert_rect_slice(EMData* w, EMData* myfft, const Transform& trans, int sizeofprojection, float xratio, float yratio, int npad, int mult)
01205 {
01206 ENTERFUNC;
01207 vector<int> saved_offsets = get_array_offsets();
01208 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01209 set_array_offsets(0,1,1);
01210 myfft->set_array_offsets(0,1);
01211
01212
01213
01214 Vec2f coordinate_2d_square;
01215 Vec3f coordinate_3dnew;
01216 Vec3f axis_newx;
01217 Vec3f axis_newy;
01218 Vec3f tempv;
01219
01220
01221
01222 axis_newx[0] = xratio*0.5f*(sizeofprojection*npad)*trans[0][0];
01223 axis_newx[1] = yratio*0.5f*(sizeofprojection*npad)*trans[0][1];
01224 axis_newx[2] = 0.5f*(sizeofprojection*npad)*trans[0][2];
01225
01226 float ellipse_length_x = std::sqrt(axis_newx[0]*axis_newx[0]+axis_newx[1]*axis_newx[1]+axis_newx[2]*axis_newx[2]);
01227
01228 int ellipse_length_x_int = int(ellipse_length_x);
01229 float ellipse_step_x = 0.5f*(sizeofprojection*npad)/float(ellipse_length_x_int);
01230 float xscale = ellipse_step_x;
01231
01232 axis_newy[0] = xratio*0.5f*(sizeofprojection*npad)*trans[1][0];
01233 axis_newy[1] = yratio*0.5f*(sizeofprojection*npad)*trans[1][1];
01234 axis_newy[2] = 0.5f*(sizeofprojection*npad)*trans[1][2];
01235
01236
01237
01238 float ellipse_length_y = std::sqrt(axis_newy[0]*axis_newy[0]+axis_newy[1]*axis_newy[1]+axis_newy[2]*axis_newy[2]);
01239 int ellipse_length_y_int = int(ellipse_length_y);
01240 float ellipse_step_y = 0.5f*(sizeofprojection*npad)/float(ellipse_length_y_int);
01241 float yscale = ellipse_step_y;
01242
01243 std::complex<float> c1;
01244 nz = get_zsize();
01245
01246 float r2=0.25f*sizeofprojection*npad*sizeofprojection*npad;
01247 float r2_at_point;
01248
01249 for(int i=0;i<ellipse_length_x_int;i++) {
01250 for(int j=-1*ellipse_length_y_int+1; j<=ellipse_length_y_int; j++) {
01251
01252 r2_at_point=i*xscale*i*xscale+j*yscale*j*yscale;
01253 if(r2_at_point<=r2 ) {
01254
01255
01256 coordinate_2d_square[0] = xscale*float(i);
01257 coordinate_2d_square[1] = yscale*float(j);
01258 float xnew = coordinate_2d_square[0]*trans[0][0] + coordinate_2d_square[1]*trans[1][0];
01259 float ynew = coordinate_2d_square[0]*trans[0][1] + coordinate_2d_square[1]*trans[1][1];
01260 float znew = coordinate_2d_square[0]*trans[0][2] + coordinate_2d_square[1]*trans[1][2];
01261 coordinate_3dnew[0] =xnew*xratio;
01262 coordinate_3dnew[1] = ynew*yratio;
01263 coordinate_3dnew[2] = znew;
01264
01265
01266 float xp = coordinate_2d_square[0];
01267 float yp = ( coordinate_2d_square[1] >= 0) ? coordinate_2d_square[1]+1 : nz+coordinate_2d_square[1]+1;
01268 std::complex<float> lin_interpolated(0,0);
01269 int xlow=int(xp),xhigh=int(xp)+1;
01270 int ylow=int(yp),yhigh=int(yp)+1;
01271 float tx=xp-xlow,ty=yp-ylow;
01272
01273
01274 if(j == -1) {
01275
01276 if(ylow<yp)
01277 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)*(1-ty) + myfft->cmplx(xlow,yhigh)*(1-tx)*ty
01278 + myfft->cmplx(xhigh,ylow)*tx*(1-ty) + myfft->cmplx(xhigh,yhigh)*tx*ty;
01279 else
01280 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)
01281 + myfft->cmplx(xhigh,ylow)*tx;
01282
01283 }
01284 else {
01285 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)*(1-ty) + myfft->cmplx(xlow,yhigh)*(1-tx)*ty
01286 + myfft->cmplx(xhigh,ylow)*tx*(1-ty)+ myfft->cmplx(xhigh,yhigh)*tx*ty;
01287
01288 }
01289
01290 c1 = lin_interpolated;
01291
01292
01293
01294 std::complex<float> btq;
01295 if ( coordinate_3dnew[0] < 0.) {
01296 coordinate_3dnew[0] = -coordinate_3dnew[0];
01297 coordinate_3dnew[1] = -coordinate_3dnew[1];
01298 coordinate_3dnew[2] = -coordinate_3dnew[2];
01299 btq = conj(c1);
01300 } else {
01301 btq = c1;
01302 }
01303 int ixn = int(coordinate_3dnew[0] + 0.5 + nx) - nx;
01304 int iyn = int(coordinate_3dnew[1] + 0.5 + ny) - ny;
01305 int izn = int(coordinate_3dnew[2] + 0.5 + nz) - nz;
01306
01307 int iza, iya;
01308 if (izn >= 0) iza = izn + 1;
01309 else iza = nz + izn + 1;
01310
01311 if (iyn >= 0) iya = iyn + 1;
01312 else iya = ny + iyn + 1;
01313
01314 cmplx(ixn,iya,iza) += btq*float(mult);
01315 (*w)(ixn,iya,iza) += mult;
01316
01317 }
01318 }
01319
01320 }
01321
01322
01323
01324
01325 set_array_offsets(saved_offsets);
01326 myfft->set_array_offsets(myfft_saved_offsets);
01327 EXITFUNC;
01328
01329 }
01330
01331
01332
01333 void EMData::nn_SSNR(EMData* wptr, EMData* wptr2, EMData* myfft, const Transform& tf, int)
01334 {
01335 ENTERFUNC;
01336 int nxc = attr_dict["nxc"];
01337
01338 vector<int> saved_offsets = get_array_offsets();
01339 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01340
01341 set_array_offsets(0,1,1);
01342 myfft->set_array_offsets(0,1);
01343
01344 int iymin = is_fftodd() ? -ny/2 : -ny/2 + 1 ;
01345 int iymax = ny/2;
01346 int izmin = is_fftodd() ? -nz/2 : -nz/2 + 1 ;
01347 int izmax = nz/2;
01348
01349 for (int iy = iymin; iy <= iymax; iy++) {
01350 int jp = iy >= 0 ? iy+1 : ny+iy+1;
01351 for (int ix = 0; ix <= nxc; ix++) {
01352 if (( 4*(ix*ix+iy*iy) < ny*ny ) && !( ix == 0 && iy < 0 ) ) {
01353 float xnew = ix*tf[0][0] + iy*tf[1][0];
01354 float ynew = ix*tf[0][1] + iy*tf[1][1];
01355 float znew = ix*tf[0][2] + iy*tf[1][2];
01356 std::complex<float> btq;
01357 if (xnew < 0.0) {
01358 xnew = -xnew;
01359 ynew = -ynew;
01360 znew = -znew;
01361 btq = conj(myfft->cmplx(ix,jp));
01362 } else {
01363 btq = myfft->cmplx(ix,jp);
01364 }
01365 int ixn = int(xnew + 0.5 + nx) - nx;
01366 int iyn = int(ynew + 0.5 + ny) - ny;
01367 int izn = int(znew + 0.5 + nz) - nz;
01368 if ((ixn <= nxc) && (iyn >= iymin) && (iyn <= iymax) && (izn >= izmin) && (izn <= izmax)) {
01369 if (ixn >= 0) {
01370 int iza, iya;
01371 if (izn >= 0) iza = izn + 1;
01372 else iza = nz + izn + 1;
01373
01374 if (iyn >= 0) iya = iyn + 1;
01375 else iya = ny + iyn + 1;
01376
01377 cmplx(ixn,iya,iza) += btq;
01378 (*wptr)(ixn,iya,iza)++;
01379 (*wptr2)(ixn,iya,iza) += norm(btq);
01380 } else {
01381 int izt, iyt;
01382 if (izn > 0) izt = nz - izn + 1;
01383 else izt = -izn + 1;
01384
01385 if (iyn > 0) iyt = ny - iyn + 1;
01386 else iyt = -iyn + 1;
01387
01388 cmplx(-ixn,iyt,izt) += conj(btq);
01389 (*wptr)(-ixn,iyt,izt)++;
01390 (*wptr2)(-ixn,iyt,izt) += norm(btq);
01391 }
01392 }
01393 }
01394 }
01395 }
01396 set_array_offsets(saved_offsets);
01397 myfft->set_array_offsets(myfft_saved_offsets);
01398 EXITFUNC;
01399 }
01400
01401
01402
01403 void EMData::symplane0(EMData* wptr) {
01404 ENTERFUNC;
01405 int nxc = attr_dict["nxc"];
01406 int n = nxc*2;
01407
01408 vector<int> saved_offsets = get_array_offsets();
01409 set_array_offsets(0,1,1);
01410 for (int iza = 2; iza <= nxc; iza++) {
01411 for (int iya = 2; iya <= nxc; iya++) {
01412 cmplx(0,iya,iza) += conj(cmplx(0,n-iya+2,n-iza+2));
01413 (*wptr)(0,iya,iza) += (*wptr)(0,n-iya+2,n-iza+2);
01414 cmplx(0,n-iya+2,n-iza+2) = conj(cmplx(0,iya,iza));
01415 (*wptr)(0,n-iya+2,n-iza+2) = (*wptr)(0,iya,iza);
01416 cmplx(0,n-iya+2,iza) += conj(cmplx(0,iya,n-iza+2));
01417 (*wptr)(0,n-iya+2,iza) += (*wptr)(0,iya,n-iza+2);
01418 cmplx(0,iya,n-iza+2) = conj(cmplx(0,n-iya+2,iza));
01419 (*wptr)(0,iya,n-iza+2) = (*wptr)(0,n-iya+2,iza);
01420 }
01421 }
01422 for (int iya = 2; iya <= nxc; iya++) {
01423 cmplx(0,iya,1) += conj(cmplx(0,n-iya+2,1));
01424 (*wptr)(0,iya,1) += (*wptr)(0,n-iya+2,1);
01425 cmplx(0,n-iya+2,1) = conj(cmplx(0,iya,1));
01426 (*wptr)(0,n-iya+2,1) = (*wptr)(0,iya,1);
01427 }
01428 for (int iza = 2; iza <= nxc; iza++) {
01429 cmplx(0,1,iza) += conj(cmplx(0,1,n-iza+2));
01430 (*wptr)(0,1,iza) += (*wptr)(0,1,n-iza+2);
01431 cmplx(0,1,n-iza+2) = conj(cmplx(0,1,iza));
01432 (*wptr)(0,1,n-iza+2) = (*wptr)(0,1,iza);
01433 }
01434 EXITFUNC;
01435 }
01436
01437 void EMData::symplane1(EMData* wptr, EMData* wptr2) {
01438 ENTERFUNC;
01439 int nxc = attr_dict["nxc"];
01440 int n = nxc*2;
01441 vector<int> saved_offsets = get_array_offsets();
01442 set_array_offsets(0,1,1);
01443 for (int iza = 2; iza <= nxc; iza++) {
01444 for (int iya = 2; iya <= nxc; iya++) {
01445 cmplx(0,iya,iza) += conj(cmplx(0,n-iya+2,n-iza+2));
01446 (*wptr)(0,iya,iza) += (*wptr)(0,n-iya+2,n-iza+2);
01447 (*wptr2)(0,iya,iza) += (*wptr2)(0,n-iya+2,n-iza+2);
01448 cmplx(0,n-iya+2,n-iza+2) = conj(cmplx(0,iya,iza));
01449 (*wptr)(0,n-iya+2,n-iza+2) = (*wptr)(0,iya,iza);
01450 (*wptr2)(0,n-iya+2,n-iza+2) = (*wptr2)(0,iya,iza);
01451 cmplx(0,n-iya+2,iza) += conj(cmplx(0,iya,n-iza+2));
01452 (*wptr)(0,n-iya+2,iza) += (*wptr)(0,iya,n-iza+2);
01453 (*wptr2)(0,n-iya+2,iza) += (*wptr2)(0,iya,n-iza+2);
01454 cmplx(0,iya,n-iza+2) = conj(cmplx(0,n-iya+2,iza));
01455 (*wptr)(0,iya,n-iza+2) = (*wptr)(0,n-iya+2,iza);
01456 (*wptr2)(0,iya,n-iza+2) = (*wptr2)(0,n-iya+2,iza);
01457 }
01458 }
01459 for (int iya = 2; iya <= nxc; iya++) {
01460 cmplx(0,iya,1) += conj(cmplx(0,n-iya+2,1));
01461 (*wptr)(0,iya,1) += (*wptr)(0,n-iya+2,1);
01462 (*wptr2)(0,iya,1) += (*wptr2)(0,n-iya+2,1);
01463 cmplx(0,n-iya+2,1) = conj(cmplx(0,iya,1));
01464 (*wptr)(0,n-iya+2,1) = (*wptr)(0,iya,1);
01465 (*wptr2)(0,n-iya+2,1) = (*wptr2)(0,iya,1);
01466 }
01467 for (int iza = 2; iza <= nxc; iza++) {
01468 cmplx(0,1,iza) += conj(cmplx(0,1,n-iza+2));
01469 (*wptr)(0,1,iza) += (*wptr)(0,1,n-iza+2);
01470 (*wptr2)(0,1,iza) += (*wptr2)(0,1,n-iza+2);
01471 cmplx(0,1,n-iza+2) = conj(cmplx(0,1,iza));
01472 (*wptr)(0,1,n-iza+2) = (*wptr)(0,1,iza);
01473 (*wptr2)(0,1,n-iza+2) = (*wptr2)(0,1,iza);
01474 }
01475 EXITFUNC;
01476 }
01477
01478 void EMData::symplane2(EMData* wptr, EMData* wptr2, EMData* wptr3) {
01479 ENTERFUNC;
01480 int nxc = attr_dict["nxc"];
01481 int n = nxc*2;
01482 vector<int> saved_offsets = get_array_offsets();
01483 set_array_offsets(0,1,1);
01484 for (int iza = 2; iza <= nxc; iza++) {
01485 for (int iya = 2; iya <= nxc; iya++) {
01486 cmplx(0,iya,iza) += conj(cmplx(0,n-iya+2,n-iza+2));
01487 (*wptr)(0,iya,iza) += (*wptr)(0,n-iya+2,n-iza+2);
01488 (*wptr2)(0,iya,iza) += (*wptr2)(0,n-iya+2,n-iza+2);
01489 (*wptr3)(0,iya,iza) += (*wptr3)(0,n-iya+2,n-iza+2);
01490
01491 cmplx(0,n-iya+2,n-iza+2) = conj(cmplx(0,iya,iza));
01492 (*wptr)(0,n-iya+2,n-iza+2) = (*wptr)(0,iya,iza);
01493 (*wptr2)(0,n-iya+2,n-iza+2) = (*wptr2)(0,iya,iza);
01494 (*wptr3)(0,n-iya+2,n-iza+2) = (*wptr3)(0,iya,iza);
01495
01496 cmplx(0,n-iya+2,iza) += conj(cmplx(0,iya,n-iza+2));
01497 (*wptr)(0,n-iya+2,iza) += (*wptr)(0,iya,n-iza+2);
01498 (*wptr2)(0,n-iya+2,iza) += (*wptr2)(0,iya,n-iza+2);
01499 (*wptr3)(0,n-iya+2,iza) += (*wptr3)(0,iya,n-iza+2);
01500
01501 cmplx(0,iya,n-iza+2) = conj(cmplx(0,n-iya+2,iza));
01502 (*wptr)(0,iya,n-iza+2) = (*wptr)(0,n-iya+2,iza);
01503 (*wptr2)(0,iya,n-iza+2) = (*wptr2)(0,n-iya+2,iza);
01504 (*wptr3)(0,iya,n-iza+2) = (*wptr3)(0,n-iya+2,iza);
01505 }
01506 }
01507 for (int iya = 2; iya <= nxc; iya++) {
01508 cmplx(0,iya,1) += conj(cmplx(0,n-iya+2,1));
01509 (*wptr)(0,iya,1) += (*wptr)(0,n-iya+2,1);
01510 (*wptr2)(0,iya,1) += (*wptr2)(0,n-iya+2,1);
01511 (*wptr3)(0,iya,1) += (*wptr3)(0,n-iya+2,1);
01512
01513 cmplx(0,n-iya+2,1) = conj(cmplx(0,iya,1));
01514 (*wptr)(0,n-iya+2,1) = (*wptr)(0,iya,1);
01515 (*wptr2)(0,n-iya+2,1) = (*wptr2)(0,iya,1);
01516 (*wptr3)(0,n-iya+2,1) = (*wptr3)(0,iya,1);
01517 }
01518 for (int iza = 2; iza <= nxc; iza++) {
01519 cmplx(0,1,iza) += conj(cmplx(0,1,n-iza+2));
01520 (*wptr)(0,1,iza) += (*wptr)(0,1,n-iza+2);
01521 (*wptr2)(0,1,iza) += (*wptr2)(0,1,n-iza+2);
01522 (*wptr3)(0,1,iza) += (*wptr3)(0,1,n-iza+2);
01523
01524 cmplx(0,1,n-iza+2) = conj(cmplx(0,1,iza));
01525 (*wptr)(0,1,n-iza+2) = (*wptr)(0,1,iza);
01526 (*wptr2)(0,1,n-iza+2) = (*wptr2)(0,1,iza);
01527 (*wptr3)(0,1,n-iza+2) = (*wptr3)(0,1,iza);
01528 }
01529 EXITFUNC;
01530 }
01531
01532
01533 class ctf_store
01534 {
01535 public:
01536
01537 static void init( int winsize, const Ctf* ctf )
01538 {
01539 Dict params = ctf->to_dict();
01540
01541 m_winsize = winsize;
01542
01543 m_voltage = params["voltage"];
01544 m_pixel = params["apix"];
01545 m_cs = params["cs"];
01546 m_ampcont = params["ampcont"];
01547 m_bfactor = params["bfactor"];
01548 m_defocus = params["defocus"];
01549 m_winsize2= m_winsize*m_winsize;
01550 m_vecsize = m_winsize2/4;
01551 }
01552
01553 static float get_ctf( int r2 )
01554 {
01555 float ak = std::sqrt( r2/float(m_winsize2) )/m_pixel;
01556 return Util::tf( m_defocus, ak, m_voltage, m_cs, m_ampcont, m_bfactor, 1 );
01557 }
01558
01559 private:
01560
01561 static int m_winsize, m_winsize2, m_vecsize;
01562 static float m_cs;
01563 static float m_voltage;
01564 static float m_pixel;
01565 static float m_ampcont;
01566 static float m_bfactor;
01567 static float m_defocus;
01568 };
01569
01570
01571 int ctf_store::m_winsize, ctf_store::m_winsize2, ctf_store::m_vecsize;
01572
01573 float ctf_store::m_cs, ctf_store::m_voltage, ctf_store::m_pixel;
01574 float ctf_store::m_ampcont, ctf_store::m_bfactor;
01575 float ctf_store::m_defocus;
01576
01577
01578 class ctf_store_new
01579 {
01580 public:
01581
01582 static void init( int winsize, const Ctf* ctf )
01583 {
01584 Dict params = ctf->to_dict();
01585
01586 m_winsize = winsize;
01587
01588 m_voltage = params["voltage"];
01589 m_pixel = params["apix"];
01590 m_cs = params["cs"];
01591 m_ampcont = params["ampcont"];
01592 m_bfactor = params["bfactor"];
01593 m_defocus = params["defocus"];
01594 m_winsize2= m_winsize*m_winsize;
01595 m_vecsize = m_winsize2/4;
01596 }
01597
01598 static float get_ctf( float r2 )
01599 {
01600 float ak = std::sqrt( r2/float(m_winsize2) )/m_pixel;
01601 return Util::tf( m_defocus, ak, m_voltage, m_cs, m_ampcont, m_bfactor, 1 );
01602 }
01603
01604 private:
01605
01606 static int m_winsize, m_winsize2, m_vecsize;
01607 static float m_cs;
01608 static float m_voltage;
01609 static float m_pixel;
01610 static float m_ampcont;
01611 static float m_bfactor;
01612 static float m_defocus;
01613 };
01614
01615
01616 int ctf_store_new::m_winsize, ctf_store_new::m_winsize2, ctf_store_new::m_vecsize;
01617
01618 float ctf_store_new::m_cs, ctf_store_new::m_voltage, ctf_store_new::m_pixel;
01619 float ctf_store_new::m_ampcont, ctf_store_new::m_bfactor;
01620 float ctf_store_new::m_defocus;
01621
01622
01623
01624
01625 void EMData::onelinenn_ctf(int j, int n, int n2,
01626 EMData* w, EMData* bi, const Transform& tf, int mult) {
01627
01628 int remove = bi->get_attr_default( "remove", 0 );
01629
01630 int jp = (j >= 0) ? j+1 : n+j+1;
01631
01632 for (int i = 0; i <= n2; i++) {
01633 int r2 = i*i+j*j;
01634 if ( (r2<n*n/4) && !((0==i) && (j<0)) ) {
01635 float ctf = ctf_store::get_ctf( r2 );
01636 float xnew = i*tf[0][0] + j*tf[1][0];
01637 float ynew = i*tf[0][1] + j*tf[1][1];
01638 float znew = i*tf[0][2] + j*tf[1][2];
01639 std::complex<float> btq;
01640 if (xnew < 0.) {
01641 xnew = -xnew;
01642 ynew = -ynew;
01643 znew = -znew;
01644 btq = conj(bi->cmplx(i,jp));
01645 } else btq = bi->cmplx(i,jp);
01646 int ixn = int(xnew + 0.5 + n) - n;
01647 int iyn = int(ynew + 0.5 + n) - n;
01648 int izn = int(znew + 0.5 + n) - n;
01649
01650 int iza, iya;
01651 if (izn >= 0) iza = izn + 1;
01652 else iza = n + izn + 1;
01653
01654 if (iyn >= 0) iya = iyn + 1;
01655 else iya = n + iyn + 1;
01656
01657 if(remove > 0 ) {
01658 cmplx(ixn,iya,iza) -= btq*ctf*float(mult);
01659 (*w)(ixn,iya,iza) -= ctf*ctf*mult;
01660 } else {
01661 cmplx(ixn,iya,iza) += btq*ctf*float(mult);
01662 (*w)(ixn,iya,iza) += ctf*ctf*mult;
01663 }
01664
01665
01666
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 void EMData::onelinenn_ctf_applied(int j, int n, int n2,
01709 EMData* w, EMData* bi, const Transform& tf, int mult) {
01710
01711 int remove = bi->get_attr_default( "remove", 0 );
01712
01713 int jp = (j >= 0) ? j+1 : n+j+1;
01714
01715 for (int i = 0; i <= n2; i++) {
01716 int r2 = i*i + j*j;
01717 if ( (r2< n*n/4) && !((0==i) && (j< 0)) ) {
01718 float ctf = ctf_store::get_ctf(r2);
01719
01720
01721 float xnew = i*tf[0][0] + j*tf[1][0];
01722 float ynew = i*tf[0][1] + j*tf[1][1];
01723 float znew = i*tf[0][2] + j*tf[1][2];
01724 std::complex<float> btq;
01725 if (xnew < 0.) {
01726 xnew = -xnew;
01727 ynew = -ynew;
01728 znew = -znew;
01729 btq = conj(bi->cmplx(i,jp));
01730 } else btq = bi->cmplx(i,jp);
01731 int ixn = int(xnew + 0.5 + n) - n;
01732 int iyn = int(ynew + 0.5 + n) - n;
01733 int izn = int(znew + 0.5 + n) - n;
01734
01735 int iza, iya;
01736 if (izn >= 0) iza = izn + 1;
01737 else iza = n + izn + 1;
01738
01739 if (iyn >= 0) iya = iyn + 1;
01740 else iya = n + iyn + 1;
01741
01742 if( remove > 0 ) {
01743 cmplx(ixn,iya,iza) -= btq*float(mult);
01744 (*w)(ixn,iya,iza) -= mult*ctf*ctf;
01745 } else {
01746 cmplx(ixn,iya,iza) += btq*float(mult);
01747 (*w)(ixn,iya,iza) += mult*ctf*ctf;
01748 }
01749
01750
01751
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 void
01792 EMData::nn_ctf(EMData* w, EMData* myfft, const Transform& tf, int mult) {
01793 ENTERFUNC;
01794 int nxc = attr_dict["nxc"];
01795
01796 vector<int> saved_offsets = get_array_offsets();
01797 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01798 set_array_offsets(0,1,1);
01799 myfft->set_array_offsets(0,1);
01800
01801 Ctf* ctf = myfft->get_attr("ctf");
01802 ctf_store::init( ny, ctf );
01803 if(ctf) {delete ctf; ctf=0;}
01804
01805
01806 for (int iy = -ny/2 + 1; iy <= ny/2; iy++) onelinenn_ctf(iy, ny, nxc, w, myfft, tf, mult);
01807 set_array_offsets(saved_offsets);
01808 myfft->set_array_offsets(myfft_saved_offsets);
01809 EXITFUNC;
01810 }
01811
01812 void
01813 EMData::nn_ctf_applied(EMData* w, EMData* myfft, const Transform& tf, int mult) {
01814 ENTERFUNC;
01815 int nxc = attr_dict["nxc"];
01816
01817 vector<int> saved_offsets = get_array_offsets();
01818 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01819 set_array_offsets(0,1,1);
01820 myfft->set_array_offsets(0,1);
01821
01822 Ctf* ctf = myfft->get_attr( "ctf" );
01823 ctf_store::init( ny, ctf );
01824 if(ctf) {delete ctf; ctf=0;}
01825
01826
01827
01828 for (int iy = -ny/2 + 1; iy <= ny/2; iy++) onelinenn_ctf_applied(iy, ny, nxc, w, myfft, tf, mult);
01829 set_array_offsets(saved_offsets);
01830 myfft->set_array_offsets(myfft_saved_offsets);
01831 EXITFUNC;
01832 }
01833
01834
01835 void EMData::insert_rect_slice_ctf(EMData* w, EMData* myfft, const Transform& trans, int sizeofprojection, float xratio, float yratio, int npad, int mult)
01836 {
01837 ENTERFUNC;
01838 vector<int> saved_offsets = get_array_offsets();
01839 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01840 set_array_offsets(0,1,1);
01841 myfft->set_array_offsets(0,1);
01842
01843
01844
01845 Vec2f coordinate_2d_square;
01846 Vec3f coordinate_3dnew;
01847 Vec3f axis_newx;
01848 Vec3f axis_newy;
01849 Vec3f tempv;
01850
01851
01852
01853 axis_newx[0] = xratio*0.5f*(sizeofprojection*npad)*trans[0][0];
01854 axis_newx[1] = yratio*0.5f*(sizeofprojection*npad)*trans[0][1];
01855 axis_newx[2] = 0.5f*(sizeofprojection*npad)*trans[0][2];
01856
01857 float ellipse_length_x = std::sqrt(axis_newx[0]*axis_newx[0]+axis_newx[1]*axis_newx[1]+axis_newx[2]*axis_newx[2]);
01858
01859 int ellipse_length_x_int = int(ellipse_length_x);
01860 float ellipse_step_x = 0.5f*(sizeofprojection*npad)/float(ellipse_length_x_int);
01861 float xscale = ellipse_step_x;
01862
01863 axis_newy[0] = xratio*0.5f*(sizeofprojection*npad)*trans[1][0];
01864 axis_newy[1] = yratio*0.5f*(sizeofprojection*npad)*trans[1][1];
01865 axis_newy[2] = 0.5f*(sizeofprojection*npad)*trans[1][2];
01866
01867
01868
01869 float ellipse_length_y = std::sqrt(axis_newy[0]*axis_newy[0]+axis_newy[1]*axis_newy[1]+axis_newy[2]*axis_newy[2]);
01870 int ellipse_length_y_int = int(ellipse_length_y);
01871 float ellipse_step_y = 0.5f*(sizeofprojection*npad)/float(ellipse_length_y_int);
01872 float yscale = ellipse_step_y;
01873
01874 std::complex<float> c1;
01875 nz = get_zsize();
01876 Ctf* ctf = myfft->get_attr( "ctf" );
01877 ctf_store_new::init( nz, ctf );
01878 if(ctf) {delete ctf; ctf=0;}
01879 int remove = myfft->get_attr_default( "remove", 0 );
01880
01881 float r2=0.25f*sizeofprojection*npad*sizeofprojection*npad;
01882 float r2_at_point;
01883
01884 for(int i=0;i<ellipse_length_x_int;i++) {
01885 for(int j=-1*ellipse_length_y_int+1; j<=ellipse_length_y_int; j++) {
01886
01887 r2_at_point=i*xscale*i*xscale+j*yscale*j*yscale;
01888 if(r2_at_point<=r2 && ! ((0==i) && (j<0))) {
01889
01890 float ctf_value = ctf_store_new::get_ctf( r2_at_point );
01891 coordinate_2d_square[0] = xscale*float(i);
01892 coordinate_2d_square[1] = yscale*float(j);
01893 float xnew = coordinate_2d_square[0]*trans[0][0] + coordinate_2d_square[1]*trans[1][0];
01894 float ynew = coordinate_2d_square[0]*trans[0][1] + coordinate_2d_square[1]*trans[1][1];
01895 float znew = coordinate_2d_square[0]*trans[0][2] + coordinate_2d_square[1]*trans[1][2];
01896 coordinate_3dnew[0] =xnew*xratio;
01897 coordinate_3dnew[1] = ynew*yratio;
01898 coordinate_3dnew[2] = znew;
01899
01900
01901 float xp = coordinate_2d_square[0];
01902 float yp = ( coordinate_2d_square[1] >= 0) ? coordinate_2d_square[1]+1 : nz+coordinate_2d_square[1]+1;
01903 std::complex<float> lin_interpolated(0,0);
01904 int xlow=int(xp),xhigh=int(xp)+1;
01905 int ylow=int(yp),yhigh=int(yp)+1;
01906 float tx=xp-xlow,ty=yp-ylow;
01907
01908
01909 if(j == -1) {
01910
01911 if(ylow<yp)
01912 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)*(1-ty) + myfft->cmplx(xlow,yhigh)*(1-tx)*ty
01913 + myfft->cmplx(xhigh,ylow)*tx*(1-ty) + myfft->cmplx(xhigh,yhigh)*tx*ty;
01914 else
01915 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)
01916 + myfft->cmplx(xhigh,ylow)*tx;
01917
01918 }
01919 else {
01920 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)*(1-ty) + myfft->cmplx(xlow,yhigh)*(1-tx)*ty
01921 + myfft->cmplx(xhigh,ylow)*tx*(1-ty)+ myfft->cmplx(xhigh,yhigh)*tx*ty;
01922
01923 }
01924
01925 c1 = lin_interpolated;
01926
01927
01928
01929 std::complex<float> btq;
01930 if ( coordinate_3dnew[0] < 0.) {
01931 coordinate_3dnew[0] = -coordinate_3dnew[0];
01932 coordinate_3dnew[1] = -coordinate_3dnew[1];
01933 coordinate_3dnew[2] = -coordinate_3dnew[2];
01934 btq = conj(c1);
01935 } else {
01936 btq = c1;
01937 }
01938 int ixn = int(coordinate_3dnew[0] + 0.5 + nx) - nx;
01939 int iyn = int(coordinate_3dnew[1] + 0.5 + ny) - ny;
01940 int izn = int(coordinate_3dnew[2] + 0.5 + nz) - nz;
01941
01942 int iza, iya;
01943 if (izn >= 0) iza = izn + 1;
01944 else iza = nz + izn + 1;
01945
01946 if (iyn >= 0) iya = iyn + 1;
01947 else iya = ny + iyn + 1;
01948
01949 if(remove > 0 ) {
01950 cmplx(ixn,iya,iza) -= btq*ctf_value*float(mult);
01951 (*w)(ixn,iya,iza) -= ctf_value*ctf_value*mult;
01952 } else {
01953 cmplx(ixn,iya,iza) += btq*ctf_value*float(mult);
01954 (*w)(ixn,iya,iza) += ctf_value*ctf_value*mult;
01955 }
01956
01957 }
01958 }
01959
01960 }
01961
01962
01963
01964
01965 set_array_offsets(saved_offsets);
01966 myfft->set_array_offsets(myfft_saved_offsets);
01967 EXITFUNC;
01968
01969 }
01970
01971
01972 void EMData::insert_rect_slice_ctf_applied(EMData* w, EMData* myfft,const Transform& trans,int sizeofprojection,float xratio,float yratio,int npad,int mult)
01973 {
01974 ENTERFUNC;
01975 vector<int> saved_offsets = get_array_offsets();
01976 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01977 set_array_offsets(0,1,1);
01978 myfft->set_array_offsets(0,1);
01979
01980
01981
01982 Vec2f coordinate_2d_square;
01983 Vec3f coordinate_3dnew;
01984 Vec3f axis_newx;
01985 Vec3f axis_newy;
01986 Vec3f tempv;
01987
01988
01989
01990 axis_newx[0] = xratio*0.5f*(sizeofprojection*npad)*trans[0][0];
01991 axis_newx[1] = yratio*0.5f*(sizeofprojection*npad)*trans[0][1];
01992 axis_newx[2] = 0.5f*(sizeofprojection*npad)*trans[0][2];
01993
01994 float ellipse_length_x = std::sqrt(axis_newx[0]*axis_newx[0]+axis_newx[1]*axis_newx[1]+axis_newx[2]*axis_newx[2]);
01995
01996 int ellipse_length_x_int = int(ellipse_length_x);
01997 float ellipse_step_x = 0.5f*(sizeofprojection*npad)/float(ellipse_length_x_int);
01998 float xscale = ellipse_step_x;
01999
02000 axis_newy[0] = xratio*0.5f*(sizeofprojection*npad)*trans[1][0];
02001 axis_newy[1] = yratio*0.5f*(sizeofprojection*npad)*trans[1][1];
02002 axis_newy[2] = 0.5f*(sizeofprojection*npad)*trans[1][2];
02003
02004
02005
02006 float ellipse_length_y = std::sqrt(axis_newy[0]*axis_newy[0]+axis_newy[1]*axis_newy[1]+axis_newy[2]*axis_newy[2]);
02007 int ellipse_length_y_int = int(ellipse_length_y);
02008 float ellipse_step_y = 0.5f*(sizeofprojection*npad)/float(ellipse_length_y_int);
02009 float yscale = ellipse_step_y;
02010
02011 std::complex<float> c1;
02012 nz = get_zsize();
02013 Ctf* ctf = myfft->get_attr( "ctf" );
02014 ctf_store_new::init( nz, ctf );
02015 if(ctf) {delete ctf; ctf=0;}
02016 int remove = myfft->get_attr_default( "remove", 0 );
02017
02018 float r2=0.25f*sizeofprojection*npad*sizeofprojection*npad;
02019 float r2_at_point;
02020
02021 for(int i=0;i<ellipse_length_x_int;i++) {
02022 for(int j=-1*ellipse_length_y_int+1; j<=ellipse_length_y_int; j++) {
02023
02024 r2_at_point=i*xscale*i*xscale+j*yscale*j*yscale;
02025 if(r2_at_point<=r2 && ! ((0==i) && (j<0))) {
02026
02027 float ctf_value = ctf_store_new::get_ctf( r2_at_point );
02028 coordinate_2d_square[0] = xscale*float(i);
02029 coordinate_2d_square[1] = yscale*float(j);
02030 float xnew = coordinate_2d_square[0]*trans[0][0] + coordinate_2d_square[1]*trans[1][0];
02031 float ynew = coordinate_2d_square[0]*trans[0][1] + coordinate_2d_square[1]*trans[1][1];
02032 float znew = coordinate_2d_square[0]*trans[0][2] + coordinate_2d_square[1]*trans[1][2];
02033 coordinate_3dnew[0] =xnew*xratio;
02034 coordinate_3dnew[1] = ynew*yratio;
02035 coordinate_3dnew[2] = znew;
02036
02037
02038 float xp = coordinate_2d_square[0];
02039 float yp = ( coordinate_2d_square[1] >= 0) ? coordinate_2d_square[1]+1 : nz+coordinate_2d_square[1]+1;
02040 std::complex<float> lin_interpolated(0,0);
02041 int xlow=int(xp),xhigh=int(xp)+1;
02042 int ylow=int(yp),yhigh=int(yp)+1;
02043 float tx=xp-xlow,ty=yp-ylow;
02044
02045
02046 if(j == -1) {
02047
02048 if(ylow<yp)
02049 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)*(1-ty) + myfft->cmplx(xlow,yhigh)*(1-tx)*ty
02050 + myfft->cmplx(xhigh,ylow)*tx*(1-ty) + myfft->cmplx(xhigh,yhigh)*tx*ty;
02051 else
02052 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)
02053 + myfft->cmplx(xhigh,ylow)*tx;
02054
02055 }
02056 else {
02057 lin_interpolated=myfft->cmplx(xlow,ylow)*(1-tx)*(1-ty) + myfft->cmplx(xlow,yhigh)*(1-tx)*ty
02058 + myfft->cmplx(xhigh,ylow)*tx*(1-ty)+ myfft->cmplx(xhigh,yhigh)*tx*ty;
02059
02060 }
02061
02062 c1 = lin_interpolated;
02063
02064
02065
02066 std::complex<float> btq;
02067 if ( coordinate_3dnew[0] < 0.) {
02068 coordinate_3dnew[0] = -coordinate_3dnew[0];
02069 coordinate_3dnew[1] = -coordinate_3dnew[1];
02070 coordinate_3dnew[2] = -coordinate_3dnew[2];
02071 btq = conj(c1);
02072 } else {
02073 btq = c1;
02074 }
02075 int ixn = int(coordinate_3dnew[0] + 0.5 + nx) - nx;
02076 int iyn = int(coordinate_3dnew[1] + 0.5 + ny) - ny;
02077 int izn = int(coordinate_3dnew[2] + 0.5 + nz) - nz;
02078
02079 int iza, iya;
02080 if (izn >= 0) iza = izn + 1;
02081 else iza = nz + izn + 1;
02082
02083 if (iyn >= 0) iya = iyn + 1;
02084 else iya = ny + iyn + 1;
02085
02086 if(remove > 0 ) {
02087 cmplx(ixn,iya,iza) -= btq*float(mult);
02088 (*w)(ixn,iya,iza) -= ctf_value*ctf_value*mult;
02089 } else {
02090 cmplx(ixn,iya,iza) += btq*float(mult);
02091 (*w)(ixn,iya,iza) += ctf_value*ctf_value*mult;
02092 }
02093
02094 }
02095 }
02096
02097 }
02098
02099
02100
02101
02102 set_array_offsets(saved_offsets);
02103 myfft->set_array_offsets(myfft_saved_offsets);
02104 EXITFUNC;
02105
02106 }
02107
02108
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 void EMData::nn_SSNR_ctf(EMData* wptr, EMData* wptr2, EMData* wptr3, EMData* myfft, const Transform& tf, int)
02178 {
02179
02180
02181
02182
02183
02184
02185
02186
02187
02188 ENTERFUNC;
02189 int nxc = attr_dict["nxc"];
02190 vector<int> saved_offsets = get_array_offsets();
02191 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
02192 set_array_offsets(0,1,1);
02193 myfft->set_array_offsets(0,1);
02194
02195 Ctf* ctf = myfft->get_attr("ctf");
02196 ctf_store::init( ny, ctf );
02197 int iymin = is_fftodd() ? -ny/2 : -ny/2 + 1;
02198 int iymax = ny/2;
02199 int izmin = is_fftodd() ? -nz/2 : -nz/2 + 1;
02200 int izmax = nz/2;
02201
02202 for (int iy = iymin; iy <= iymax; iy++) {
02203 int jp = iy >= 0 ? iy+1 : ny+iy+1;
02204 for (int ix = 0; ix <= nxc; ix++) {
02205 int r2 = ix*ix+iy*iy;
02206 if (( 4*r2 < ny*ny ) && !( ix == 0 && iy < 0 ) ) {
02207 float ctf = ctf_store::get_ctf( r2 )*10.f;
02208 float xnew = ix*tf[0][0] + iy*tf[1][0];
02209 float ynew = ix*tf[0][1] + iy*tf[1][1];
02210 float znew = ix*tf[0][2] + iy*tf[1][2];
02211 std::complex<float> btq;
02212 if (xnew < 0.0) {
02213 xnew = -xnew;
02214 ynew = -ynew;
02215 znew = -znew;
02216 btq = conj(myfft->cmplx(ix,jp));
02217 } else {
02218 btq = myfft->cmplx(ix,jp);
02219 }
02220 int ixn = int(xnew + 0.5 + nx) - nx;
02221 int iyn = int(ynew + 0.5 + ny) - ny;
02222 int izn = int(znew + 0.5 + nz) - nz;
02223 if ((ixn <= nxc) && (iyn >= iymin) && (iyn <= iymax) && (izn >= izmin) && (izn <= izmax)) {
02224 if (ixn >= 0) {
02225 int iza, iya;
02226 if (izn >= 0) iza = izn + 1;
02227 else iza = nz + izn + 1;
02228
02229 if (iyn >= 0) iya = iyn + 1;
02230 else iya = ny + iyn + 1;
02231
02232 cmplx(ixn,iya,iza) += btq*ctf;
02233 (*wptr)(ixn,iya,iza) += ctf*ctf;
02234 (*wptr2)(ixn,iya,iza) += std::norm(btq);
02235 (*wptr3)(ixn,iya,iza) += 1;
02236 } else {
02237 int izt, iyt;
02238 if (izn > 0) izt = nz - izn + 1;
02239 else izt = -izn + 1;
02240
02241 if (iyn > 0) iyt = ny - iyn + 1;
02242 else iyt = -iyn + 1;
02243
02244 cmplx(-ixn,iyt,izt) += std::conj(btq)*ctf;
02245 (*wptr) (-ixn,iyt,izt) += ctf*ctf;
02246 (*wptr2)(-ixn,iyt,izt) += std::norm(btq);
02247 (*wptr3)(-ixn,iyt,izt) += 1;
02248 }
02249 }
02250 }
02251 }
02252 }
02253 set_array_offsets(saved_offsets);
02254 myfft->set_array_offsets(myfft_saved_offsets);
02255 if(ctf) {delete ctf; ctf=0;}
02256 EXITFUNC;
02257 }
02258
02259
02260
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 void EMData::symplane0_ctf(EMData* w) {
02358 ENTERFUNC;
02359 int nxc = attr_dict["nxc"];
02360 int n = nxc*2;
02361
02362 vector<int> saved_offsets = get_array_offsets();
02363 set_array_offsets(0,1,1);
02364 for (int iza = 2; iza <= nxc; iza++) {
02365 for (int iya = 2; iya <= nxc; iya++) {
02366 cmplx(0,iya,iza) += conj(cmplx(0,n-iya+2,n-iza+2));
02367 (*w)(0,iya,iza) += (*w)(0,n-iya+2,n-iza+2);
02368 cmplx(0,n-iya+2,n-iza+2) = conj(cmplx(0,iya,iza));
02369 (*w)(0,n-iya+2,n-iza+2) = (*w)(0,iya,iza);
02370 cmplx(0,n-iya+2,iza) += conj(cmplx(0,iya,n-iza+2));
02371 (*w)(0,n-iya+2,iza) += (*w)(0,iya,n-iza+2);
02372 cmplx(0,iya,n-iza+2) = conj(cmplx(0,n-iya+2,iza));
02373 (*w)(0,iya,n-iza+2) = (*w)(0,n-iya+2,iza);
02374 }
02375 }
02376 for (int iya = 2; iya <= nxc; iya++) {
02377 cmplx(0,iya,1) += conj(cmplx(0,n-iya+2,1));
02378 (*w)(0,iya,1) += (*w)(0,n-iya+2,1);
02379 cmplx(0,n-iya+2,1) = conj(cmplx(0,iya,1));
02380 (*w)(0,n-iya+2,1) = (*w)(0,iya,1);
02381 }
02382 for (int iza = 2; iza <= nxc; iza++) {
02383 cmplx(0,1,iza) += conj(cmplx(0,1,n-iza+2));
02384 (*w)(0,1,iza) += (*w)(0,1,n-iza+2);
02385 cmplx(0,1,n-iza+2) = conj(cmplx(0,1,iza));
02386 (*w)(0,1,n-iza+2) = (*w)(0,1,iza);
02387 }
02388 EXITFUNC;
02389 }
02390
02391 void EMData::symplane0_rect(EMData* w) {
02392 ENTERFUNC;
02393 nx=get_xsize();
02394 ny=get_ysize();
02395 nz=get_zsize();
02396 int nzc=nz/2;
02397 int nyc=ny/2;
02398
02399
02400
02401 vector<int> saved_offsets = get_array_offsets();
02402 set_array_offsets(0,1,1);
02403 for (int iza = 2; iza <= nzc; iza++) {
02404 for (int iya = 2; iya <= nyc; iya++) {
02405 cmplx(0,iya,iza) += conj(cmplx(0,ny-iya+2,nz-iza+2));
02406 (*w)(0,iya,iza) += (*w)(0,ny-iya+2,nz-iza+2);
02407 cmplx(0,ny-iya+2,nz-iza+2) = conj(cmplx(0,iya,iza));
02408 (*w)(0,ny-iya+2,nz-iza+2) = (*w)(0,iya,iza);
02409 cmplx(0,ny-iya+2,iza) += conj(cmplx(0,iya,nz-iza+2));
02410 (*w)(0,ny-iya+2,iza) += (*w)(0,iya,nz-iza+2);
02411 cmplx(0,iya,nz-iza+2) = conj(cmplx(0,ny-iya+2,iza));
02412 (*w)(0,iya,nz-iza+2) = (*w)(0,ny-iya+2,iza);
02413 }
02414 }
02415 for (int iya = 2; iya <= nyc; iya++) {
02416 cmplx(0,iya,1) += conj(cmplx(0,ny-iya+2,1));
02417 (*w)(0,iya,1) += (*w)(0,ny-iya+2,1);
02418 cmplx(0,ny-iya+2,1) = conj(cmplx(0,iya,1));
02419 (*w)(0,ny-iya+2,1) = (*w)(0,iya,1);
02420 }
02421 for (int iza = 2; iza <= nzc; iza++) {
02422 cmplx(0,1,iza) += conj(cmplx(0,1,nz-iza+2));
02423 (*w)(0,1,iza) += (*w)(0,1,nz-iza+2);
02424 cmplx(0,1,nz-iza+2) = conj(cmplx(0,1,iza));
02425 (*w)(0,1,nz-iza+2) = (*w)(0,1,iza);
02426 }
02427 EXITFUNC;
02428 }
02429
02430 EMData* EMData::rot_scale_trans2D(float angDeg, float delx, float dely, float scale) {
02431 float ang=angDeg*M_PI/180.0f;
02432 if (1 >= ny)
02433 throw ImageDimensionException("Can't rotate 1D image");
02434 if (nz<2) {
02435 vector<int> saved_offsets = get_array_offsets();
02436 set_array_offsets(0,0,0);
02437 if (0.0f == scale) scale = 1.0f;
02438 EMData* ret = copy_head();
02439 delx = restrict2(delx, nx);
02440 dely = restrict2(dely, ny);
02441
02442 int xc = nx/2;
02443 int yc = ny/2;
02444
02445 float shiftxc = xc + delx;
02446 float shiftyc = yc + dely;
02447
02448 float cang = cos(ang);
02449 float sang = sin(ang);
02450 for (int iy = 0; iy < ny; iy++) {
02451 float y = float(iy) - shiftyc;
02452 float ycang = y*cang/scale + yc;
02453 float ysang = -y*sang/scale + xc;
02454 for (int ix = 0; ix < nx; ix++) {
02455 float x = float(ix) - shiftxc;
02456 float xold = x*cang/scale + ysang ;
02457 float yold = x*sang/scale + ycang ;
02458
02459 (*ret)(ix,iy) = Util::quadri(xold+1.0f, yold+1.0f, nx, ny, get_data());
02460
02461 }
02462 }
02463 set_array_offsets(saved_offsets);
02464 return ret;
02465 } else {
02466 throw ImageDimensionException("Volume not currently supported");
02467 }
02468 }
02469
02470 EMData* EMData::rot_scale_trans2D_background(float angDeg, float delx, float dely, float scale) {
02471 float ang=angDeg*M_PI/180.0f;
02472 if (1 >= ny)
02473 throw ImageDimensionException("Can't rotate 1D image");
02474 if (nz<2) {
02475 vector<int> saved_offsets = get_array_offsets();
02476 set_array_offsets(0,0,0);
02477 if (0.0f == scale) scale = 1.0f;
02478 EMData* ret = copy_head();
02479 delx = restrict2(delx, nx);
02480 dely = restrict2(dely, ny);
02481
02482 int xc = nx/2;
02483 int yc = ny/2;
02484
02485 float shiftxc = xc + delx;
02486 float shiftyc = yc + dely;
02487
02488 float cang = cos(ang);
02489 float sang = sin(ang);
02490 for (int iy = 0; iy < ny; iy++) {
02491 float y = float(iy) - shiftyc;
02492 float ycang = y*cang/scale + yc;
02493 float ysang = -y*sang/scale + xc;
02494 for (int ix = 0; ix < nx; ix++) {
02495 float x = float(ix) - shiftxc;
02496 float xold = x*cang/scale + ysang ;
02497 float yold = x*sang/scale + ycang ;
02498
02499 (*ret)(ix,iy) = Util::quadri_background(xold+1.0f, yold+1.0f, nx, ny, get_data(),ix+1,iy+1);
02500
02501 }
02502 }
02503 set_array_offsets(saved_offsets);
02504 return ret;
02505 } else {
02506 throw ImageDimensionException("Volume not currently supported");
02507 }
02508 }
02509
02510 #define in(i,j,k) in[i+(j+(k*ny))*(size_t)nx]
02511 EMData*
02512 EMData::rot_scale_trans(const Transform &RA) {
02513
02514 EMData* ret = copy_head();
02515 float *in = this->get_data();
02516 vector<int> saved_offsets = get_array_offsets();
02517 set_array_offsets(0,0,0);
02518 Vec3f translations = RA.get_trans();
02519 Transform RAinv = RA.inverse();
02520
02521 if (1 >= ny) throw ImageDimensionException("Can't rotate 1D image");
02522 if (nz < 2) {
02523 float p1, p2, p3, p4;
02524 float delx = translations.at(0);
02525 float dely = translations.at(1);
02526 delx = restrict2(delx, nx);
02527 dely = restrict2(dely, ny);
02528 int xc = nx/2;
02529 int yc = ny/2;
02530
02531 float shiftxc = xc + delx;
02532 float shiftyc = yc + dely;
02533 for (int iy = 0; iy < ny; iy++) {
02534 float y = float(iy) - shiftyc;
02535 float ysang = y*RAinv[0][1]+xc;
02536 float ycang = y*RAinv[1][1]+yc;
02537 for (int ix = 0; ix < nx; ix++) {
02538 float x = float(ix) - shiftxc;
02539 float xold = x*RAinv[0][0] + ysang;
02540 float yold = x*RAinv[1][0] + ycang;
02541
02542 xold = restrict1(xold, nx);
02543 yold = restrict1(yold, ny);
02544
02545 int xfloor = int(xold);
02546 int yfloor = int(yold);
02547 float t = xold-xfloor;
02548 float u = yold-yfloor;
02549 if(xfloor == nx -1 && yfloor == ny -1) {
02550
02551 p1 =in[xfloor + yfloor*ny];
02552 p2 =in[ yfloor*ny];
02553 p3 =in[0];
02554 p4 =in[xfloor];
02555 } else if(xfloor == nx - 1) {
02556
02557 p1 =in[xfloor + yfloor*ny];
02558 p2 =in[ yfloor*ny];
02559 p3 =in[ (yfloor+1)*ny];
02560 p4 =in[xfloor + (yfloor+1)*ny];
02561 } else if(yfloor == ny - 1) {
02562
02563 p1 =in[xfloor + yfloor*ny];
02564 p2 =in[xfloor+1 + yfloor*ny];
02565 p3 =in[xfloor+1 ];
02566 p4 =in[xfloor ];
02567 } else {
02568 p1 =in[xfloor + yfloor*ny];
02569 p2 =in[xfloor+1 + yfloor*ny];
02570 p3 =in[xfloor+1 + (yfloor+1)*ny];
02571 p4 =in[xfloor + (yfloor+1)*ny];
02572 }
02573 (*ret)(ix,iy) = p1 + u * ( p4 - p1) + t * ( p2 - p1 + u *(p3-p2-p4+p1));
02574 }
02575 }
02576 set_array_offsets(saved_offsets);
02577 return ret;
02578 } else {
02579
02580
02581 float delx = translations.at(0);
02582 float dely = translations.at(1);
02583 float delz = translations.at(2);
02584 delx = restrict2(delx, nx);
02585 dely = restrict2(dely, ny);
02586 delz = restrict2(delz, nz);
02587 int xc = nx/2;
02588 int yc = ny/2;
02589 int zc = nz/2;
02590
02591 float shiftxc = xc + delx;
02592 float shiftyc = yc + dely;
02593 float shiftzc = zc + delz;
02594
02595 for (int iz = 0; iz < nz; iz++) {
02596 float z = float(iz) - shiftzc;
02597 float xoldz = z*RAinv[0][2]+xc;
02598 float yoldz = z*RAinv[1][2]+yc;
02599 float zoldz = z*RAinv[2][2]+zc;
02600 for (int iy = 0; iy < ny; iy++) {
02601 float y = float(iy) - shiftyc;
02602 float xoldzy = xoldz + y*RAinv[0][1] ;
02603 float yoldzy = yoldz + y*RAinv[1][1] ;
02604 float zoldzy = zoldz + y*RAinv[2][1] ;
02605 for (int ix = 0; ix < nx; ix++) {
02606 float x = float(ix) - shiftxc;
02607 float xold = xoldzy + x*RAinv[0][0] ;
02608 float yold = yoldzy + x*RAinv[1][0] ;
02609 float zold = zoldzy + x*RAinv[2][0] ;
02610
02611 xold = restrict1(xold, nx);
02612 yold = restrict1(yold, ny);
02613 zold = restrict1(zold, nz);
02614
02615
02616 int IOX = int(xold);
02617 int IOY = int(yold);
02618 int IOZ = int(zold);
02619
02620 #ifdef _WIN32
02621 int IOXp1 = _cpp_min( nx-1 ,IOX+1);
02622 #else
02623 int IOXp1 = std::min( nx-1 ,IOX+1);
02624 #endif //_WIN32
02625
02626 #ifdef _WIN32
02627 int IOYp1 = _cpp_min( ny-1 ,IOY+1);
02628 #else
02629 int IOYp1 = std::min( ny-1 ,IOY+1);
02630 #endif //_WIN32
02631
02632 #ifdef _WIN32
02633 int IOZp1 = _cpp_min( nz-1 ,IOZ+1);
02634 #else
02635 int IOZp1 = std::min( nz-1 ,IOZ+1);
02636 #endif //_WIN32
02637
02638 float dx = xold-IOX;
02639 float dy = yold-IOY;
02640 float dz = zold-IOZ;
02641
02642 float a1 = in(IOX,IOY,IOZ);
02643 float a2 = in(IOXp1,IOY,IOZ) - in(IOX,IOY,IOZ);
02644 float a3 = in(IOX,IOYp1,IOZ) - in(IOX,IOY,IOZ);
02645 float a4 = in(IOX,IOY,IOZp1) - in(IOX,IOY,IOZ);
02646 float a5 = in(IOX,IOY,IOZ) - in(IOXp1,IOY,IOZ) - in(IOX,IOYp1,IOZ) + in(IOXp1,IOYp1,IOZ);
02647 float a6 = in(IOX,IOY,IOZ) - in(IOXp1,IOY,IOZ) - in(IOX,IOY,IOZp1) + in(IOXp1,IOY,IOZp1);
02648 float a7 = in(IOX,IOY,IOZ) - in(IOX,IOYp1,IOZ) - in(IOX,IOY,IOZp1) + in(IOX,IOYp1,IOZp1);
02649 float a8 = in(IOXp1,IOY,IOZ) + in(IOX,IOYp1,IOZ)+ in(IOX,IOY,IOZp1)
02650 - in(IOX,IOY,IOZ)- in(IOXp1,IOYp1,IOZ) - in(IOXp1,IOY,IOZp1)
02651 - in(IOX,IOYp1,IOZp1) + in(IOXp1,IOYp1,IOZp1);
02652 (*ret)(ix,iy,iz) = a1 + dz*(a4 + a6*dx + (a7 + a8*dx)*dy) + a3*dy + dx*(a2 + a5*dy);
02653 }
02654 }
02655 }
02656
02657 set_array_offsets(saved_offsets);
02658 return ret;
02659
02660
02661
02662
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 #undef in
02777
02778
02779 #define in(i,j,k) in[i+(j+(k*ny))*(size_t)nx]
02780 EMData*
02781 EMData::rot_scale_trans_background(const Transform &RA) {
02782 EMData* ret = copy_head();
02783 float *in = this->get_data();
02784 vector<int> saved_offsets = get_array_offsets();
02785 set_array_offsets(0,0,0);
02786 Vec3f translations = RA.get_trans();
02787 Transform RAinv = RA.inverse();
02788
02789 if (1 >= ny) throw ImageDimensionException("Can't rotate 1D image");
02790 if (nz < 2) {
02791 float p1, p2, p3, p4;
02792 float delx = translations.at(0);
02793 float dely = translations.at(1);
02794 delx = restrict2(delx, nx);
02795 dely = restrict2(dely, ny);
02796 int xc = nx/2;
02797 int yc = ny/2;
02798
02799 float shiftxc = xc + delx;
02800 float shiftyc = yc + dely;
02801 for (int iy = 0; iy < ny; iy++) {
02802 float y = float(iy) - shiftyc;
02803 float ysang = y*RAinv[0][1]+xc;
02804 float ycang = y*RAinv[1][1]+yc;
02805 for (int ix = 0; ix < nx; ix++) {
02806 float x = float(ix) - shiftxc;
02807 float xold = x*RAinv[0][0] + ysang;
02808 float yold = x*RAinv[1][0] + ycang;
02809
02810
02811
02812 if ( (xold < 0.0f) || (xold >= (float)(nx)) || (yold < 0.0f) || (yold >= (float)(ny)) ){
02813 xold = (float)ix;
02814 yold = (float)iy;
02815 }
02816
02817 int xfloor = int(xold);
02818 int yfloor = int(yold);
02819 float t = xold-xfloor;
02820 float u = yold-yfloor;
02821 if(xfloor == nx -1 && yfloor == ny -1) {
02822
02823 p1 =in[xfloor + yfloor*ny];
02824 p2 =in[ yfloor*ny];
02825 p3 =in[0];
02826 p4 =in[xfloor];
02827 } else if(xfloor == nx - 1) {
02828
02829 p1 =in[xfloor + yfloor*ny];
02830 p2 =in[ yfloor*ny];
02831 p3 =in[ (yfloor+1)*ny];
02832 p4 =in[xfloor + (yfloor+1)*ny];
02833 } else if(yfloor == ny - 1) {
02834
02835 p1 =in[xfloor + yfloor*ny];
02836 p2 =in[xfloor+1 + yfloor*ny];
02837 p3 =in[xfloor+1 ];
02838 p4 =in[xfloor ];
02839 } else {
02840
02841 p1 =in[xfloor + yfloor*ny];
02842 p2 =in[xfloor+1 + yfloor*ny];
02843 p3 =in[xfloor+1 + (yfloor+1)*ny];
02844 p4 =in[xfloor + (yfloor+1)*ny];
02845 }
02846 (*ret)(ix,iy) = p1 + u * ( p4 - p1) + t * ( p2 - p1 + u *(p3-p2-p4+p1));
02847 }
02848 }
02849 set_array_offsets(saved_offsets);
02850 return ret;
02851 } else {
02852
02853
02854 float delx = translations.at(0);
02855 float dely = translations.at(1);
02856 float delz = translations.at(2);
02857 delx = restrict2(delx, nx);
02858 dely = restrict2(dely, ny);
02859 delz = restrict2(delz, nz);
02860 int xc = nx/2;
02861 int yc = ny/2;
02862 int zc = nz/2;
02863
02864 float shiftxc = xc + delx;
02865 float shiftyc = yc + dely;
02866 float shiftzc = zc + delz;
02867
02868 for (int iz = 0; iz < nz; iz++) {
02869 float z = float(iz) - shiftzc;
02870 float xoldz = z*RAinv[0][2]+xc;
02871 float yoldz = z*RAinv[1][2]+yc;
02872 float zoldz = z*RAinv[2][2]+zc;
02873 for (int iy = 0; iy < ny; iy++) {
02874 float y = float(iy) - shiftyc;
02875 float xoldzy = xoldz + y*RAinv[0][1] ;
02876 float yoldzy = yoldz + y*RAinv[1][1] ;
02877 float zoldzy = zoldz + y*RAinv[2][1] ;
02878 for (int ix = 0; ix < nx; ix++) {
02879 float x = float(ix) - shiftxc;
02880 float xold = xoldzy + x*RAinv[0][0] ;
02881 float yold = yoldzy + x*RAinv[1][0] ;
02882 float zold = zoldzy + x*RAinv[2][0] ;
02883
02884
02885
02886 if ( (xold < 0.0f) || (xold >= (float)(nx)) || (yold < 0.0f) || (yold >= (float)(ny)) || (zold < 0.0f) || (zold >= (float)(nz)) ){
02887 xold = (float)ix;
02888 yold = (float)iy;
02889 zold = (float)iz;
02890 }
02891
02892 int IOX = int(xold);
02893 int IOY = int(yold);
02894 int IOZ = int(zold);
02895
02896 #ifdef _WIN32
02897 int IOXp1 = _cpp_min( nx-1 ,IOX+1);
02898 #else
02899 int IOXp1 = std::min( nx-1 ,IOX+1);
02900 #endif //_WIN32
02901
02902 #ifdef _WIN32
02903 int IOYp1 = _cpp_min( ny-1 ,IOY+1);
02904 #else
02905 int IOYp1 = std::min( ny-1 ,IOY+1);
02906 #endif //_WIN32
02907
02908 #ifdef _WIN32
02909 int IOZp1 = _cpp_min( nz-1 ,IOZ+1);
02910 #else
02911 int IOZp1 = std::min( nz-1 ,IOZ+1);
02912 #endif //_WIN32
02913
02914 float dx = xold-IOX;
02915 float dy = yold-IOY;
02916 float dz = zold-IOZ;
02917
02918 float a1 = in(IOX,IOY,IOZ);
02919 float a2 = in(IOXp1,IOY,IOZ) - in(IOX,IOY,IOZ);
02920 float a3 = in(IOX,IOYp1,IOZ) - in(IOX,IOY,IOZ);
02921 float a4 = in(IOX,IOY,IOZp1) - in(IOX,IOY,IOZ);
02922 float a5 = in(IOX,IOY,IOZ) - in(IOXp1,IOY,IOZ) - in(IOX,IOYp1,IOZ) + in(IOXp1,IOYp1,IOZ);
02923 float a6 = in(IOX,IOY,IOZ) - in(IOXp1,IOY,IOZ) - in(IOX,IOY,IOZp1) + in(IOXp1,IOY,IOZp1);
02924 float a7 = in(IOX,IOY,IOZ) - in(IOX,IOYp1,IOZ) - in(IOX,IOY,IOZp1) + in(IOX,IOYp1,IOZp1);
02925 float a8 = in(IOXp1,IOY,IOZ) + in(IOX,IOYp1,IOZ)+ in(IOX,IOY,IOZp1)
02926 - in(IOX,IOY,IOZ)- in(IOXp1,IOYp1,IOZ) - in(IOXp1,IOY,IOZp1)
02927 - in(IOX,IOYp1,IOZp1) + in(IOXp1,IOYp1,IOZp1);
02928 (*ret)(ix,iy,iz) = a1 + dz*(a4 + a6*dx + (a7 + a8*dx)*dy) + a3*dy + dx*(a2 + a5*dy);
02929 }
02930 }
02931 }
02932
02933 set_array_offsets(saved_offsets);
02934 return ret;
02935
02936 }
02937 }
02938 #undef in
02939
02940
02941
02942
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 EMData* EMData::rot_scale_conv(float ang, float delx, float dely, Util::KaiserBessel& kb, float scale_input) {
03017 int nxn, nyn, nzn;
03018 if(scale_input == 0.0f) scale_input = 1.0f;
03019
03020 float scale = 0.5f*scale_input;
03021 float sum, w;
03022 if (1 >= ny)
03023 throw ImageDimensionException("Can't rotate 1D image");
03024 if (1 < nz)
03025 throw ImageDimensionException("Volume not currently supported");
03026 nxn=nx/2;nyn=ny/2;nzn=nz/2;
03027
03028 int K = kb.get_window_size();
03029 int kbmin = -K/2;
03030 int kbmax = -kbmin;
03031 int kbc = kbmax+1;
03032 vector<int> saved_offsets = get_array_offsets();
03033 set_array_offsets(0,0,0);
03034 EMData* ret = this->copy_head();
03035 #ifdef _WIN32
03036 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03037 #else
03038 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03039 #endif //_WIN32
03040
03041 delx = restrict2(delx, nx);
03042 dely = restrict2(dely, ny);
03043
03044 int xc = nxn;
03045 int ixs = nxn%2;
03046 int yc = nyn;
03047 int iys = nyn%2;
03048
03049 int xcn = nxn/2;
03050 int ycn = nyn/2;
03051
03052 float shiftxc = xcn + delx;
03053 float shiftyc = ycn + dely;
03054
03055 float ymin = -ny/2.0f;
03056 float xmin = -nx/2.0f;
03057 float ymax = -ymin;
03058 float xmax = -xmin;
03059 if (0 == nx%2) xmax--;
03060 if (0 == ny%2) ymax--;
03061
03062 float *t = (float*)calloc(kbmax-kbmin+1, sizeof(float));
03063
03064
03065 float cang = cos(ang);
03066 float sang = sin(ang);
03067 for (int iy = 0; iy < nyn; iy++) {
03068 float y = float(iy) - shiftyc;
03069 float ycang = y*cang/scale + yc;
03070 float ysang = -y*sang/scale + xc;
03071 for (int ix = 0; ix < nxn; ix++) {
03072 float x = float(ix) - shiftxc;
03073 float xold = x*cang/scale + ysang-ixs;
03074 float yold = x*sang/scale + ycang-iys;
03075
03076 xold = restrict1(xold, nx);
03077 yold = restrict1(yold, ny);
03078
03079 int inxold = int(Util::round(xold)); int inyold = int(Util::round(yold));
03080 sum=0.0f; w=0.0f;
03081 for (int m1 =kbmin; m1 <=kbmax; m1++) t[m1-kbmin] = kb.i0win_tab(xold - inxold-m1);
03082 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 ) {
03083 for (int m2 =kbmin; m2 <=kbmax; m2++) {
03084 float qt = kb.i0win_tab(yold - inyold-m2);
03085 for (int m1 =kbmin; m1 <=kbmax; m1++) {
03086 float q = t[m1-kbmin]*qt;
03087 sum += (*this)((inxold+m1+nx)%nx,(inyold+m2+ny)%ny)*q; w+=q;
03088 }
03089 }
03090 } else {
03091 for (int m2 =kbmin; m2 <=kbmax; m2++) {
03092 float qt = kb.i0win_tab(yold - inyold-m2);
03093 for (int m1 =kbmin; m1 <=kbmax; m1++) {
03094 float q = t[m1-kbmin]*qt;
03095 sum += (*this)(inxold+m1,inyold+m2)*q; w+=q;}
03096 }
03097 }
03098 (*ret)(ix,iy)=sum/w;
03099 }
03100 }
03101 if (t) free(t);
03102 set_array_offsets(saved_offsets);
03103 return ret;
03104 }
03105
03106
03107
03108 EMData* EMData::rot_scale_conv7(float ang, float delx, float dely, Util::KaiserBessel& kb, float scale_input) {
03109 int nxn, nyn, nzn;
03110 float scale = 0.5f*scale_input;
03111 float sum, w;
03112 if (1 >= ny)
03113 throw ImageDimensionException("Can't rotate 1D image");
03114 if (1 < nz)
03115 throw ImageDimensionException("Volume not currently supported");
03116 nxn = nx/2; nyn=ny/2; nzn=nz/2;
03117
03118 int K = kb.get_window_size();
03119 int kbmin = -K/2;
03120 int kbmax = -kbmin;
03121 int kbc = kbmax+1;
03122 vector<int> saved_offsets = get_array_offsets();
03123 set_array_offsets(0,0,0);
03124 EMData* ret = this->copy_head();
03125 #ifdef _WIN32
03126 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03127 #else
03128 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03129 #endif //_WIN32
03130
03131 delx = restrict2(delx, nx);
03132 dely = restrict2(dely, ny);
03133
03134 int xc = nxn;
03135 int ixs = nxn%2;
03136 int yc = nyn;
03137 int iys = nyn%2;
03138
03139 int xcn = nxn/2;
03140 int ycn = nyn/2;
03141
03142 float shiftxc = xcn + delx;
03143 float shiftyc = ycn + dely;
03144
03145 float ymin = -ny/2.0f;
03146 float xmin = -nx/2.0f;
03147 float ymax = -ymin;
03148 float xmax = -xmin;
03149 if (0 == nx%2) xmax--;
03150 if (0 == ny%2) ymax--;
03151
03152 float *t = (float*)calloc(kbmax-kbmin+1, sizeof(float));
03153
03154
03155 float cang = cos(ang);
03156 float sang = sin(ang);
03157 for (int iy = 0; iy < nyn; iy++) {
03158 float y = float(iy) - shiftyc;
03159 float ycang = y*cang/scale + yc;
03160 float ysang = -y*sang/scale + xc;
03161 for (int ix = 0; ix < nxn; ix++) {
03162 float x = float(ix) - shiftxc;
03163 float xold = x*cang/scale + ysang-ixs;
03164 float yold = x*sang/scale + ycang-iys;
03165
03166 xold = restrict1(xold, nx);
03167 yold = restrict1(yold, ny);
03168
03169 int inxold = int(Util::round(xold)); int inyold = int(Util::round(yold));
03170 sum=0.0f; w=0.0f;
03171
03172 float tablex1 = kb.i0win_tab(xold-inxold+3);
03173 float tablex2 = kb.i0win_tab(xold-inxold+2);
03174 float tablex3 = kb.i0win_tab(xold-inxold+1);
03175 float tablex4 = kb.i0win_tab(xold-inxold);
03176 float tablex5 = kb.i0win_tab(xold-inxold-1);
03177 float tablex6 = kb.i0win_tab(xold-inxold-2);
03178 float tablex7 = kb.i0win_tab(xold-inxold-3);
03179
03180 float tabley1 = kb.i0win_tab(yold-inyold+3);
03181 float tabley2 = kb.i0win_tab(yold-inyold+2);
03182 float tabley3 = kb.i0win_tab(yold-inyold+1);
03183 float tabley4 = kb.i0win_tab(yold-inyold);
03184 float tabley5 = kb.i0win_tab(yold-inyold-1);
03185 float tabley6 = kb.i0win_tab(yold-inyold-2);
03186 float tabley7 = kb.i0win_tab(yold-inyold-3);
03187
03188 int x1, x2, x3, x4, x5, x6, x7, y1, y2, y3, y4, y5, y6, y7;
03189
03190 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 ) {
03191 x1 = (inxold-3+nx)%nx;
03192 x2 = (inxold-2+nx)%nx;
03193 x3 = (inxold-1+nx)%nx;
03194 x4 = (inxold +nx)%nx;
03195 x5 = (inxold+1+nx)%nx;
03196 x6 = (inxold+2+nx)%nx;
03197 x7 = (inxold+3+nx)%nx;
03198
03199 y1 = (inyold-3+ny)%ny;
03200 y2 = (inyold-2+ny)%ny;
03201 y3 = (inyold-1+ny)%ny;
03202 y4 = (inyold +ny)%ny;
03203 y5 = (inyold+1+ny)%ny;
03204 y6 = (inyold+2+ny)%ny;
03205 y7 = (inyold+3+ny)%ny;
03206 } else {
03207 x1 = inxold-3;
03208 x2 = inxold-2;
03209 x3 = inxold-1;
03210 x4 = inxold;
03211 x5 = inxold+1;
03212 x6 = inxold+2;
03213 x7 = inxold+3;
03214
03215 y1 = inyold-3;
03216 y2 = inyold-2;
03217 y3 = inyold-1;
03218 y4 = inyold;
03219 y5 = inyold+1;
03220 y6 = inyold+2;
03221 y7 = inyold+3;
03222 }
03223 sum = ( (*this)(x1,y1)*tablex1 + (*this)(x2,y1)*tablex2 + (*this)(x3,y1)*tablex3 +
03224 (*this)(x4,y1)*tablex4 + (*this)(x5,y1)*tablex5 + (*this)(x6,y1)*tablex6 +
03225 (*this)(x7,y1)*tablex7 ) * tabley1 +
03226 ( (*this)(x1,y2)*tablex1 + (*this)(x2,y2)*tablex2 + (*this)(x3,y2)*tablex3 +
03227 (*this)(x4,y2)*tablex4 + (*this)(x5,y2)*tablex5 + (*this)(x6,y2)*tablex6 +
03228 (*this)(x7,y2)*tablex7 ) * tabley2 +
03229 ( (*this)(x1,y3)*tablex1 + (*this)(x2,y3)*tablex2 + (*this)(x3,y3)*tablex3 +
03230 (*this)(x4,y3)*tablex4 + (*this)(x5,y3)*tablex5 + (*this)(x6,y3)*tablex6 +
03231 (*this)(x7,y3)*tablex7 ) * tabley3 +
03232 ( (*this)(x1,y4)*tablex1 + (*this)(x2,y4)*tablex2 + (*this)(x3,y4)*tablex3 +
03233 (*this)(x4,y4)*tablex4 + (*this)(x5,y4)*tablex5 + (*this)(x6,y4)*tablex6 +
03234 (*this)(x7,y4)*tablex7 ) * tabley4 +
03235 ( (*this)(x1,y5)*tablex1 + (*this)(x2,y5)*tablex2 + (*this)(x3,y5)*tablex3 +
03236 (*this)(x4,y5)*tablex4 + (*this)(x5,y5)*tablex5 + (*this)(x6,y5)*tablex6 +
03237 (*this)(x7,y5)*tablex7 ) * tabley5 +
03238 ( (*this)(x1,y6)*tablex1 + (*this)(x2,y6)*tablex2 + (*this)(x3,y6)*tablex3 +
03239 (*this)(x4,y6)*tablex4 + (*this)(x5,y6)*tablex5 + (*this)(x6,y6)*tablex6 +
03240 (*this)(x7,y6)*tablex7 ) * tabley6 +
03241 ( (*this)(x1,y7)*tablex1 + (*this)(x2,y7)*tablex2 + (*this)(x3,y7)*tablex3 +
03242 (*this)(x4,y7)*tablex4 + (*this)(x5,y7)*tablex5 + (*this)(x6,y7)*tablex6 +
03243 (*this)(x7,y7)*tablex7 ) * tabley7;
03244
03245 w = (tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7) *
03246 (tabley1+tabley2+tabley3+tabley4+tabley5+tabley6+tabley7);
03247
03248 (*ret)(ix,iy)=sum/w;
03249 }
03250 }
03251 if (t) free(t);
03252 set_array_offsets(saved_offsets);
03253 return ret;
03254 }
03255
03256 EMData* EMData::downsample(Util::sincBlackman& kb, float scale) {
03257
03258
03259
03260
03261
03262 int nxn, nyn, nzn;
03263 nxn = (int)(nx*scale); nyn = (int)(ny*scale); nzn = (int)(nz*scale);
03264
03265 vector<int> saved_offsets = get_array_offsets();
03266 set_array_offsets(0,0,0);
03267 EMData* ret = this->copy_head();
03268 #ifdef _WIN32
03269 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03270 #else
03271 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03272 #endif //_WIN32
03273 ret->to_zero();
03274
03275
03276 for (int iy =0; iy < nyn; iy++) {
03277 float y = float(iy)/scale;
03278 for (int ix = 0; ix < nxn; ix++) {
03279 float x = float(ix)/scale;
03280 (*ret)(ix,iy) = this->get_pixel_filtered(x, y, 1.0f, kb);
03281 }
03282 }
03283 set_array_offsets(saved_offsets);
03284 return ret;
03285 }
03286
03287
03288 EMData* EMData::rot_scale_conv_new(float ang, float delx, float dely, Util::KaiserBessel& kb, float scale_input) {
03289
03290 if (scale_input == 0.0f) scale_input = 1.0f;
03291 float scale = 0.5f*scale_input;
03292
03293 if (1 >= ny)
03294 throw ImageDimensionException("Can't rotate 1D image");
03295 if (1 < nz)
03296 throw ImageDimensionException("Use rot_scale_conv_new_3D for volumes");
03297 int nxn = nx/2; int nyn = ny/2; int nzn = nz/2;
03298
03299 vector<int> saved_offsets = get_array_offsets();
03300 set_array_offsets(0,0,0);
03301 EMData* ret = this->copy_head();
03302 #ifdef _WIN32
03303 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03304 #else
03305 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03306 #endif //_WIN32
03307
03308 delx = restrict2(delx, nx);
03309 dely = restrict2(dely, ny);
03310
03311 int xc = nxn;
03312 int ixs = nxn%2;
03313 int yc = nyn;
03314 int iys = nyn%2;
03315
03316 int xcn = nxn/2;
03317 int ycn = nyn/2;
03318
03319 float shiftxc = xcn + delx;
03320 float shiftyc = ycn + dely;
03321
03322 float ymin = -ny/2.0f;
03323 float xmin = -nx/2.0f;
03324 float ymax = -ymin;
03325 float xmax = -xmin;
03326 if (0 == nx%2) xmax--;
03327 if (0 == ny%2) ymax--;
03328
03329 float* data = this->get_data();
03330
03331 float cang = cos(ang);
03332 float sang = sin(ang);
03333 for (int iy = 0; iy < nyn; iy++) {
03334 float y = float(iy) - shiftyc;
03335 float ycang = y*cang/scale + yc;
03336 float ysang = -y*sang/scale + xc;
03337 for (int ix = 0; ix < nxn; ix++) {
03338 float x = float(ix) - shiftxc;
03339 float xold = x*cang/scale + ysang-ixs;
03340 float yold = x*sang/scale + ycang-iys;
03341
03342 (*ret)(ix,iy) = Util::get_pixel_conv_new(nx, ny, 1, xold, yold, 1, data, kb);
03343 }
03344 }
03345 set_array_offsets(saved_offsets);
03346 return ret;
03347 }
03348
03349 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) {
03350
03351 if (scale_input == 0.0f) scale_input = 1.0f;
03352 float scale = 0.5f*scale_input;
03353
03354 if (1 >= ny)
03355 throw ImageDimensionException("Can't rotate 1D image");
03356 int nxn = nx/2; int nyn = ny/2; int nzn = nz/2;
03357
03358 vector<int> saved_offsets = get_array_offsets();
03359 set_array_offsets(0,0,0);
03360 EMData* ret = this->copy_head();
03361 #ifdef _WIN32
03362 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03363 #else
03364 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03365 #endif //_WIN32
03366
03367 if(wrap){
03368 delx = restrict2(delx, nx);
03369 dely = restrict2(dely, ny);
03370 delz = restrict2(delz, nz);
03371 }
03372
03373 int xc = nxn;
03374 int ixs = nxn%2;
03375 int yc = nyn;
03376 int iys = nyn%2;
03377 int zc = nzn;
03378 int izs = nzn%2;
03379
03380 int xcn = nxn/2;
03381 int ycn = nyn/2;
03382 int zcn = nzn/2;
03383
03384 float shiftxc = xcn + delx;
03385 float shiftyc = ycn + dely;
03386 float shiftzc = zcn + delz;
03387
03388 float zmin = -nz/2.0f;
03389 float ymin = -ny/2.0f;
03390 float xmin = -nx/2.0f;
03391 float zmax = -zmin;
03392 float ymax = -ymin;
03393 float xmax = -xmin;
03394 if (0 == nx%2) xmax--;
03395 if (0 == ny%2) ymax--;
03396 if (0 == nz%2) zmax--;
03397
03398 float* data = this->get_data();
03399
03400 float cf = cos(phi); float sf = sin(phi);
03401 float ct = cos(theta); float st = sin(theta);
03402 float cp = cos(psi); float sp = sin(psi);
03403
03404 float a11 = cp*ct*cf-sp*sf; float a12 = cp*ct*sf+sp*cf; float a13 = -cp*st;
03405 float a21 = -sp*ct*cf-cp*sf; float a22 = -sp*ct*sf+cp*cf; float a23 = sp*st;
03406 float a31 = st*cf; float a32 = st*sf; float a33 = ct;
03407 for (int iz = 0; iz < nzn; iz++) {
03408 float z = (float(iz) - shiftzc)/scale;
03409 float zco1 = a31*z+xc;
03410 float zco2 = a32*z+yc;
03411 float zco3 = a33*z+zc;
03412 for (int iy = 0; iy < nyn; iy++) {
03413 float y = (float(iy) - shiftyc)/scale;
03414 float yco1 = zco1+a21*y;
03415 float yco2 = zco2+a22*y;
03416 float yco3 = zco3+a23*y;
03417 for (int ix = 0; ix < nxn; ix++) {
03418 float x = (float(ix) - shiftxc)/scale;
03419 float xold = yco1+a11*x-ixs;
03420 float yold = yco2+a12*x-iys;
03421 float zold = yco3+a13*x-izs;
03422 if(!wrap && (xold<0.0 || xold>nx-1 || yold<0.0 || yold>ny-1 || zold<0.0 || zold>nz-1))
03423 (*ret)(ix,iy,iz) = 0.0;
03424 else
03425 (*ret)(ix,iy,iz) = Util::get_pixel_conv_new(nx, ny, nz, xold, yold, zold, data, kb);
03426 }
03427 }
03428 }
03429 set_array_offsets(saved_offsets);
03430 return ret;
03431 }
03432
03433 EMData* EMData::rot_scale_conv_new_background(float ang, float delx, float dely, Util::KaiserBessel& kb, float scale_input) {
03434
03435 int nxn, nyn, nzn;
03436
03437 if (scale_input == 0.0f) scale_input = 1.0f;
03438 float scale = 0.5f*scale_input;
03439
03440 if (1 >= ny)
03441 throw ImageDimensionException("Can't rotate 1D image");
03442 if (1 < nz)
03443 throw ImageDimensionException("Use rot_scale_conv_new_background_3D for volumes");
03444 nxn = nx/2; nyn = ny/2; nzn = nz/2;
03445
03446 vector<int> saved_offsets = get_array_offsets();
03447 set_array_offsets(0,0,0);
03448 EMData* ret = this->copy_head();
03449 #ifdef _WIN32
03450 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03451 #else
03452 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03453 #endif //_WIN32
03454
03455 delx = restrict2(delx, nx);
03456 dely = restrict2(dely, ny);
03457
03458 int xc = nxn;
03459 int ixs = nxn%2;
03460 int yc = nyn;
03461 int iys = nyn%2;
03462
03463 int xcn = nxn/2;
03464 int ycn = nyn/2;
03465
03466 float shiftxc = xcn + delx;
03467 float shiftyc = ycn + dely;
03468
03469 float ymin = -ny/2.0f;
03470 float xmin = -nx/2.0f;
03471 float ymax = -ymin;
03472 float xmax = -xmin;
03473 if (0 == nx%2) xmax--;
03474 if (0 == ny%2) ymax--;
03475
03476 float* data = this->get_data();
03477
03478
03479 float cang = cos(ang);
03480 float sang = sin(ang);
03481 for (int iy = 0; iy < nyn; iy++) {
03482 float y = float(iy) - shiftyc;
03483 float ycang = y*cang/scale + yc;
03484 float ysang = -y*sang/scale + xc;
03485 for (int ix = 0; ix < nxn; ix++) {
03486 float x = float(ix) - shiftxc;
03487 float xold = x*cang/scale + ysang-ixs;
03488 float yold = x*sang/scale + ycang-iys;
03489
03490 (*ret)(ix,iy) = Util::get_pixel_conv_new_background(nx, ny, 1, xold, yold, 1, data, kb, ix, iy);
03491 }
03492 }
03493 set_array_offsets(saved_offsets);
03494 return ret;
03495 }
03496
03497 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) {
03498
03499 if (scale_input == 0.0f) scale_input = 1.0f;
03500 float scale = 0.5f*scale_input;
03501
03502 if (1 >= ny)
03503 throw ImageDimensionException("Can't rotate 1D image");
03504 int nxn = nx/2; int nyn = ny/2; int nzn = nz/2;
03505
03506 vector<int> saved_offsets = get_array_offsets();
03507 set_array_offsets(0,0,0);
03508 EMData* ret = this->copy_head();
03509 #ifdef _WIN32
03510 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
03511 #else
03512 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
03513 #endif //_WIN32
03514
03515 if (wrap){
03516 delx = restrict2(delx, nx);
03517 dely = restrict2(dely, ny);
03518 delz = restrict2(delz, nz);
03519 }
03520
03521 int xc = nxn;
03522 int ixs = nxn%2;
03523 int yc = nyn;
03524 int iys = nyn%2;
03525 int zc = nzn;
03526 int izs = nzn%2;
03527
03528 int xcn = nxn/2;
03529 int ycn = nyn/2;
03530 int zcn = nzn/2;
03531
03532 float shiftxc = xcn + delx;
03533 float shiftyc = ycn + dely;
03534 float shiftzc = zcn + delz;
03535
03536 float zmin = -nz/2.0f;
03537 float ymin = -ny/2.0f;
03538 float xmin = -nx/2.0f;
03539 float zmax = -zmin;
03540 float ymax = -ymin;
03541 float xmax = -xmin;
03542 if (0 == nx%2) xmax--;
03543 if (0 == ny%2) ymax--;
03544 if (0 == nz%2) zmax--;
03545
03546 float* data = this->get_data();
03547
03548 float cf = cos(phi); float sf = sin(phi);
03549 float ct = cos(theta); float st = sin(theta);
03550 float cp = cos(psi); float sp = sin(psi);
03551
03552 float a11 = cp*ct*cf-sp*sf; float a12 = cp*ct*sf+sp*cf; float a13 = -cp*st;
03553 float a21 = -sp*ct*cf-cp*sf; float a22 = -sp*ct*sf+cp*cf; float a23 = sp*st;
03554 float a31 = st*cf; float a32 = st*sf; float a33 = ct;
03555 for (int iz = 0; iz < nzn; iz++) {
03556 float z = (float(iz) - shiftzc)/scale;
03557 float zco1 = a31*z+xc;
03558 float zco2 = a32*z+yc;
03559 float zco3 = a33*z+zc;
03560 for (int iy = 0; iy < nyn; iy++) {
03561 float y = (float(iy) - shiftyc)/scale;
03562 float yco1 = zco1+a21*y;
03563 float yco2 = zco2+a22*y;
03564 float yco3 = zco3+a23*y;
03565 for (int ix = 0; ix < nxn; ix++) {
03566 float x = (float(ix) - shiftxc)/scale;
03567 float xold = yco1+a11*x-ixs;
03568 float yold = yco2+a12*x-iys;
03569 float zold = yco3+a13*x-izs;
03570 if(!wrap && (xold<0.0 || xold>nx-1 || yold<0.0 || yold>ny-1 || zold<0.0 || zold>nz-1))
03571 (*ret)(ix,iy,iz) = 0.0;
03572 else
03573 (*ret)(ix,iy,iz) = Util::get_pixel_conv_new_background(nx, ny, nz, xold, yold, zold, data, kb, ix, iy);
03574 }
03575 }
03576 }
03577 set_array_offsets(saved_offsets);
03578 return ret;
03579 }
03580
03581
03582 float EMData::get_pixel_conv(float delx, float dely, float delz, Util::KaiserBessel& kb) {
03583
03584
03585 int K = kb.get_window_size();
03586 int kbmin = -K/2;
03587 int kbmax = -kbmin;
03588 int kbc = kbmax+1;
03589
03590 float pixel =0.0f;
03591 float w=0.0f;
03592
03593 delx = restrict2(delx, nx);
03594 int inxold = int(Util::round(delx));
03595 if(ny<2) {
03596 if(inxold <= kbc || inxold >=nx-kbc-2 ) {
03597
03598 for (int m1 =kbmin; m1 <=kbmax; m1++) {
03599 float q = kb.i0win_tab(delx - inxold-m1);
03600 pixel += (*this)((inxold+m1+nx)%nx)*q; w+=q;
03601 }
03602 } else {
03603 for (int m1 =kbmin; m1 <=kbmax; m1++) {
03604 float q = kb.i0win_tab(delx - inxold-m1);
03605 pixel += (*this)(inxold+m1)*q; w+=q;
03606 }
03607 }
03608
03609 } else if(nz<2) {
03610 dely = restrict2(dely, ny);
03611 int inyold = int(Util::round(dely));
03612 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 ) {
03613
03614 for (int m2 =kbmin; m2 <=kbmax; m2++){ for (int m1 =kbmin; m1 <=kbmax; m1++) {
03615 float q = kb.i0win_tab(delx - inxold-m1)*kb.i0win_tab(dely - inyold-m2);
03616 pixel += (*this)((inxold+m1+nx)%nx,(inyold+m2+ny)%ny)*q; w+=q;}
03617 }
03618 } else {
03619 for (int m2 =kbmin; m2 <=kbmax; m2++){ for (int m1 =kbmin; m1 <=kbmax; m1++) {
03620 float q = kb.i0win_tab(delx - inxold-m1)*kb.i0win_tab(dely - inyold-m2);
03621 pixel += (*this)(inxold+m1,inyold+m2)*q; w+=q;}
03622 }
03623 }
03624 } else {
03625 dely = restrict2(dely, ny);
03626 int inyold = int(Util::round(dely));
03627 delz = restrict2(delz, nz);
03628 int inzold = int(Util::round(delz));
03629
03630 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 || inzold <= kbc || inzold >=nz-kbc-2 ) {
03631
03632 for (int m3 =kbmin; m3 <=kbmax; m3++){ for (int m2 =kbmin; m2 <=kbmax; m2++){ for (int m1 =kbmin; m1 <=kbmax; m1++) {
03633 float q = kb.i0win_tab(delx - inxold-m1)*kb.i0win_tab(dely - inyold-m2)*kb.i0win_tab(delz - inzold-m3);
03634
03635 pixel += (*this)((inxold+m1+nx)%nx,(inyold+m2+ny)%ny,(inzold+m3+nz)%nz)*q ;w+=q;}}
03636 }
03637 } else {
03638 for (int m3 =kbmin; m3 <=kbmax; m3++){ for (int m2 =kbmin; m2 <=kbmax; m2++){ for (int m1 =kbmin; m1 <=kbmax; m1++) {
03639 float q = kb.i0win_tab(delx - inxold-m1)*kb.i0win_tab(dely - inyold-m2)*kb.i0win_tab(delz - inzold-m3);
03640
03641 pixel += (*this)(inxold+m1,inyold+m2,inzold+m3)*q; w+=q;}}
03642 }
03643 }
03644 }
03645 return pixel/w;
03646 }
03647
03648
03649 float EMData::get_pixel_filtered(float delx, float dely, float, Util::sincBlackman& kb) {
03650
03651
03652 int K = kb.get_sB_size();
03653 int kbmin = -K/2;
03654 int kbmax = -kbmin;
03655 int kbc = kbmax+1;
03656
03657 float pixel =0.0f;
03658 float w=0.0f;
03659
03660
03661 int inxold = int(Util::round(delx));
03662
03663
03664
03665
03666
03667
03668
03669
03670
03671
03672
03673
03674
03675
03676
03677
03678 int inyold = int(Util::round(dely));
03679 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 ) {
03680
03681 for (int m2 =kbmin; m2 <=kbmax; m2++){
03682 float t = kb.sBwin_tab(dely - inyold-m2);
03683 for (int m1 =kbmin; m1 <=kbmax; m1++) {
03684 float q = kb.sBwin_tab(delx - inxold-m1)*t;
03685 pixel += (*this)((inxold+m1+nx)%nx,(inyold+m2+ny)%ny)*q;
03686 w += q;
03687 }
03688 }
03689 } else {
03690 for (int m2 =kbmin; m2 <=kbmax; m2++){
03691 float t = kb.sBwin_tab(dely - inyold-m2);
03692 for (int m1 =kbmin; m1 <=kbmax; m1++) {
03693 float q = kb.sBwin_tab(delx - inxold-m1)*t;
03694 pixel += (*this)(inxold+m1,inyold+m2)*q;
03695 w += q;
03696 }
03697 }
03698 }
03699
03700
03701
03702
03703
03704
03705
03706
03707
03708
03709
03710
03711
03712
03713
03714
03715
03716
03717
03718
03719
03720 return pixel/w;
03721 }
03722
03723
03724
03725
03726 float EMData::get_pixel_conv7(float delx, float dely, float delz, Util::KaiserBessel& kb) {
03727
03728
03729 float *image=(this->get_data());
03730 int nx = this->get_xsize();
03731 int ny = this->get_ysize();
03732 int nz = this->get_zsize();
03733
03734 float result;
03735
03736 result = Util::get_pixel_conv_new(nx,ny,nz,delx,dely,delz,image,kb);
03737 return result;
03738 }
03739
03740 float EMData::getconvpt2d_kbi0(float x, float y, Util::KaiserBessel::kbi0_win win, int size) {
03741 const int nxhalf = nx/2;
03742 const int nyhalf = ny/2;
03743 const int bd = size/2;
03744 float* wxarr = new float[size];
03745 float* wyarr = new float[size];
03746 float* wx = wxarr + bd;
03747 float* wy = wyarr + bd;
03748 int ixc = int(x + 0.5f*Util::sgn(x));
03749 int iyc = int(y + 0.5f*Util::sgn(y));
03750 if (abs(ixc) > nxhalf)
03751 throw InvalidValueException(ixc, "getconv: X value out of range");
03752 if (abs(iyc) > nyhalf)
03753 throw InvalidValueException(ixc, "getconv: Y value out of range");
03754 for (int i = -bd; i <= bd; i++) {
03755 int iyp = iyc + i;
03756 wy[i] = win(y - iyp);
03757 int ixp = ixc + i;
03758 wx[i] = win(x - ixp);
03759 }
03760 vector<int> saved_offsets = get_array_offsets();
03761 set_array_offsets(-nxhalf, -nyhalf);
03762 float conv = 0.f, wsum = 0.f;
03763 for (int iy = -bd; iy <= bd; iy++) {
03764 int iyp = iyc + iy;
03765 for (int ix = -bd; ix <= bd; ix++) {
03766 int ixp = ixc + ix;
03767 float wg = wx[ix]*wy[iy];
03768 conv += (*this)(ixp,iyp)*wg;
03769 wsum += wg;
03770 }
03771 }
03772 set_array_offsets(saved_offsets);
03773 delete [] wxarr;
03774 delete [] wyarr;
03775
03776 return conv;
03777 }
03778
03779 std::complex<float> EMData::extractpoint(float nuxnew, float nuynew, Util::KaiserBessel& kb) {
03780 if (2 != get_ndim())
03781 throw ImageDimensionException("extractpoint needs a 2-D image.");
03782 if (!is_complex())
03783 throw ImageFormatException("extractpoint requires a fourier image");
03784 int nxreal = nx - 2;
03785 if (nxreal != ny)
03786 throw ImageDimensionException("extractpoint requires ny == nx");
03787 int nhalf = nxreal/2;
03788 int kbsize = kb.get_window_size();
03789 int kbmin = -kbsize/2;
03790 int kbmax = -kbmin;
03791 bool flip = (nuxnew < 0.f);
03792 if (flip) {
03793 nuxnew *= -1;
03794 nuynew *= -1;
03795 }
03796
03797
03798
03799 int ixn = int(Util::round(nuxnew));
03800 int iyn = int(Util::round(nuynew));
03801
03802 float* wy0 = new float[kbmax - kbmin + 1];
03803 float* wy = wy0 - kbmin;
03804 float* wx0 = new float[kbmax - kbmin + 1];
03805 float* wx = wx0 - kbmin;
03806 for (int i = kbmin; i <= kbmax; i++) {
03807 int iyp = iyn + i;
03808 wy[i] = kb.i0win_tab(nuynew - iyp);
03809 int ixp = ixn + i;
03810 wx[i] = kb.i0win_tab(nuxnew - ixp);
03811 }
03812
03813 int iymin = 0;
03814 for (int iy = kbmin; iy <= -1; iy++) {
03815 if (wy[iy] != 0.f) {
03816 iymin = iy;
03817 break;
03818 }
03819 }
03820 int iymax = 0;
03821 for (int iy = kbmax; iy >= 1; iy--) {
03822 if (wy[iy] != 0.f) {
03823 iymax = iy;
03824 break;
03825 }
03826 }
03827 int ixmin = 0;
03828 for (int ix = kbmin; ix <= -1; ix++) {
03829 if (wx[ix] != 0.f) {
03830 ixmin = ix;
03831 break;
03832 }
03833 }
03834 int ixmax = 0;
03835 for (int ix = kbmax; ix >= 1; ix--) {
03836 if (wx[ix] != 0.f) {
03837 ixmax = ix;
03838 break;
03839 }
03840 }
03841 float wsum = 0.0f;
03842 for (int iy = iymin; iy <= iymax; iy++)
03843 for (int ix = ixmin; ix <= ixmax; ix++)
03844 wsum += wx[ix]*wy[iy];
03845 std::complex<float> result(0.f,0.f);
03846 if ((ixn >= -kbmin) && (ixn <= nhalf-1-kbmax) && (iyn >= -nhalf-kbmin) && (iyn <= nhalf-1-kbmax)) {
03847
03848 for (int iy = iymin; iy <= iymax; iy++) {
03849 int iyp = iyn + iy;
03850 for (int ix = ixmin; ix <= ixmax; ix++) {
03851 int ixp = ixn + ix;
03852 float w = wx[ix]*wy[iy];
03853 std::complex<float> val = cmplx(ixp,iyp);
03854 result += val*w;
03855 }
03856 }
03857 } else {
03858
03859 for (int iy = iymin; iy <= iymax; iy++) {
03860 int iyp = iyn + iy;
03861 for (int ix = ixmin; ix <= ixmax; ix++) {
03862 int ixp = ixn + ix;
03863 bool mirror = false;
03864 int ixt= ixp, iyt= iyp;
03865 if (ixt < 0) {
03866 ixt = -ixt;
03867 iyt = -iyt;
03868 mirror = !mirror;
03869 }
03870 if (ixt > nhalf) {
03871 ixt = nxreal - ixt;
03872 iyt = -iyt;
03873 mirror = !mirror;
03874 }
03875 if (iyt > nhalf-1) iyt -= nxreal;
03876 if (iyt < -nhalf) iyt += nxreal;
03877 float w = wx[ix]*wy[iy];
03878 std::complex<float> val = this->cmplx(ixt,iyt);
03879 if (mirror) result += conj(val)*w;
03880 else result += val*w;
03881 }
03882 }
03883 }
03884 if (flip) result = conj(result)/wsum;
03885 else result /= wsum;
03886 delete [] wx0;
03887 delete [] wy0;
03888 return result;
03889 }
03890
03891 EMData* EMData::extractline(Util::KaiserBessel& kb, float nuxnew, float nuynew)
03892 {
03893 if (!is_complex())
03894 throw ImageFormatException("extractline requires a fourier image");
03895 if (nx%2 != 0)
03896 throw ImageDimensionException("extractline requires nx to be even");
03897 int nxreal = nx - 2;
03898 if (nxreal != ny)
03899 throw ImageDimensionException("extractline requires ny == nx");
03900
03901 EMData* res = new EMData();
03902 res->set_size(nx,1,1);
03903 res->to_zero();
03904 res->set_complex(true);
03905 res->set_fftodd(false);
03906 res->set_fftpad(true);
03907 res->set_ri(true);
03908
03909 int n = nxreal;
03910 int nhalf = n/2;
03911 vector<int> saved_offsets = get_array_offsets();
03912 set_array_offsets(0,-nhalf,-nhalf);
03913
03914
03915 int kbsize = kb.get_window_size();
03916 int kbmin = -kbsize/2;
03917 int kbmax = -kbmin;
03918 float* wy0 = new float[kbmax - kbmin + 1];
03919 float* wy = wy0 - kbmin;
03920 float* wx0 = new float[kbmax - kbmin + 1];
03921 float* wx = wx0 - kbmin;
03922
03923 int count = 0;
03924 float wsum = 0.f;
03925 bool flip = (nuxnew < 0.f);
03926
03927 for (int jx = 0; jx <= nhalf; jx++) {
03928 float xnew = jx*nuxnew, ynew = jx*nuynew;
03929 count++;
03930 std::complex<float> btq(0.f,0.f);
03931 if (flip) {
03932 xnew = -xnew;
03933 ynew = -ynew;
03934 }
03935 int ixn = int(Util::round(xnew));
03936 int iyn = int(Util::round(ynew));
03937
03938 for (int i=kbmin; i <= kbmax; i++) {
03939 int iyp = iyn + i;
03940 wy[i] = kb.i0win_tab(ynew - iyp);
03941 int ixp = ixn + i;
03942 wx[i] = kb.i0win_tab(xnew - ixp);
03943 }
03944
03945
03946 int lnby = 0;
03947 for (int iy = kbmin; iy <= -1; iy++) {
03948 if (wy[iy] != 0.f) {
03949 lnby = iy;
03950 break;
03951 }
03952 }
03953 int lney = 0;
03954 for (int iy = kbmax; iy >= 1; iy--) {
03955 if (wy[iy] != 0.f) {
03956 lney = iy;
03957 break;
03958 }
03959 }
03960 int lnbx = 0;
03961 for (int ix = kbmin; ix <= -1; ix++) {
03962 if (wx[ix] != 0.f) {
03963 lnbx = ix;
03964 break;
03965 }
03966 }
03967 int lnex = 0;
03968 for (int ix = kbmax; ix >= 1; ix--) {
03969 if (wx[ix] != 0.f) {
03970 lnex = ix;
03971 break;
03972 }
03973 }
03974 if (ixn >= -kbmin && ixn <= nhalf-1-kbmax
03975 && iyn >= -nhalf-kbmin && iyn <= nhalf-1-kbmax) {
03976
03977 for (int ly=lnby; ly<=lney; ly++) {
03978 int iyp = iyn + ly;
03979 for (int lx=lnbx; lx<=lnex; lx++) {
03980 int ixp = ixn + lx;
03981 float wg = wx[lx]*wy[ly];
03982 btq += cmplx(ixp,iyp)*wg;
03983 wsum += wg;
03984 }
03985 }
03986 } else {
03987
03988 for (int ly=lnby; ly<=lney; ly++) {
03989 int iyp = iyn + ly;
03990 for (int lx=lnbx; lx<=lnex; lx++) {
03991 int ixp = ixn + lx;
03992 float wg = wx[lx]*wy[ly];
03993 bool mirror = false;
03994 int ixt(ixp), iyt(iyp);
03995 if (ixt > nhalf || ixt < -nhalf) {
03996 ixt = Util::sgn(ixt)*(n - abs(ixt));
03997 iyt = -iyt;
03998 mirror = !mirror;
03999 }
04000 if (iyt >= nhalf || iyt < -nhalf) {
04001 if (ixt != 0) {
04002 ixt = -ixt;
04003 iyt = Util::sgn(iyt)*(n - abs(iyt));
04004 mirror = !mirror;
04005 } else {
04006 iyt -= n*Util::sgn(iyt);
04007 }
04008 }
04009 if (ixt < 0) {
04010 ixt = -ixt;
04011 iyt = -iyt;
04012 mirror = !mirror;
04013 }
04014 if (iyt == nhalf) iyt = -nhalf;
04015 if (mirror) btq += conj(cmplx(ixt,iyt))*wg;
04016 else btq += cmplx(ixt,iyt)*wg;
04017 wsum += wg;
04018 }
04019 }
04020 }
04021 if (flip) res->cmplx(jx) = conj(btq);
04022 else res->cmplx(jx) = btq;
04023 }
04024 for (int jx = 0; jx <= nhalf; jx++) res->cmplx(jx) *= count/wsum;
04025
04026 delete[] wx0; delete[] wy0;
04027 set_array_offsets(saved_offsets);
04028 res->set_array_offsets(0,0,0);
04029 return res;
04030 }
04031
04032
04034 inline void swapx(float* a, float* b, float* temp, size_t nbytes) {
04035 memcpy(temp, a, nbytes);
04036 memcpy(a, b, nbytes);
04037 memcpy(b, temp, nbytes);
04038 }
04039
04040 void EMData::fft_shuffle() {
04041 if (!is_complex())
04042 throw ImageFormatException("fft_shuffle requires a fourier image");
04043 vector<int> offsets = get_array_offsets();
04044 set_array_offsets();
04045 EMData& self = *this;
04046 int nyhalf = ny/2;
04047 int nzhalf = nz/2;
04048 int nbytes = nx*sizeof(float);
04049 float* temp = new float[nx];
04050 for (int iz=0; iz < nz; iz++)
04051 for (int iy=0; iy < nyhalf; iy++)
04052 swapx(&self(0,iy,iz),&self(0,iy+nyhalf,iz),temp,nbytes);
04053 if (nz > 1) {
04054 for (int iy=0; iy < ny; iy++)
04055 for (int iz=0; iz < nzhalf; iz++)
04056 swapx(&self(0,iy,iz),&self(0,iy,iz+nzhalf),temp,nbytes);
04057 }
04058 set_shuffled(!is_shuffled());
04059 set_array_offsets(offsets);
04060 update();
04061 delete[] temp;
04062 }
04063
04064 void EMData::pad_corner(float *pad_image) {
04065 size_t nbytes = nx*sizeof(float);
04066 for (int iy=0; iy<ny; iy++)
04067 memcpy(&(*this)(0,iy), pad_image+3+(iy+3)*nx, nbytes);
04068 }
04069
04070 void EMData::shuffle_pad_corner(float *pad_image) {
04071 int nyhalf = ny/2;
04072 size_t nbytes = nx*sizeof(float);
04073 for (int iy = 0; iy < nyhalf; iy++)
04074 memcpy(&(*this)(0,iy), pad_image+6+(iy+nyhalf+3)*nx, nbytes);
04075 for (int iy = nyhalf; iy < ny; iy++)
04076 memcpy(&(*this)(0,iy), pad_image+6+(iy-nyhalf+3)*nx, nbytes);
04077 }
04078
04079 #define QUADPI 3.141592653589793238462643383279502884197
04080 #define DGR_TO_RAD QUADPI/180
04081
04082
04083
04084
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 EMData* EMData::fouriergridrot2d(float ang, float scale, Util::KaiserBessel& kb) {
04136 if (2 != get_ndim())
04137 throw ImageDimensionException("fouriergridrot2d needs a 2-D image.");
04138 if (!is_complex())
04139 throw ImageFormatException("fouriergridrot2d requires a fourier image");
04140 int nxreal = nx - 2 + int(is_fftodd());
04141 if (nxreal != ny)
04142 throw ImageDimensionException("fouriergridrot2d requires ny == nx(real)");
04143 if (0 != nxreal%2)
04144 throw ImageDimensionException("fouriergridrot2d needs an even image.");
04145 if (scale == 0.0f) scale = 1.0f;
04146 int nxhalf = nxreal/2;
04147 int nyhalf = ny/2;
04148 float cir = (float)((nxhalf-1)*(nxhalf-1));
04149
04150 if (!is_shuffled()) fft_shuffle();
04151
04152 EMData* result = copy_head();
04153 set_array_offsets(0,-nyhalf);
04154 result->set_array_offsets(0,-nyhalf);
04155
04156
04157
04158 ang = ang*(float)DGR_TO_RAD;
04159 float cang = cos(ang);
04160 float sang = sin(ang);
04161 for (int iy = -nyhalf; iy < nyhalf; iy++) {
04162 float ycang = iy*cang;
04163 float ysang = iy*sang;
04164 for (int ix = 0; ix <= nxhalf; ix++) {
04165 float nuxold = (ix*cang - ysang)*scale;
04166 float nuyold = (ix*sang + ycang)*scale;
04167 if (nuxold*nuxold+nuyold*nuyold<cir) result->cmplx(ix,iy) = Util::extractpoint2(nx, ny, nuxold, nuyold, this, kb);
04168
04169 }
04170 }
04171 result->set_array_offsets();
04172 result->fft_shuffle();
04173 result->update();
04174 set_array_offsets();
04175 fft_shuffle();
04176 return result;
04177 }
04178
04179 EMData* EMData::fouriergridrot_shift2d(float ang, float sx, float sy, Util::KaiserBessel& kb) {
04180 if (2 != get_ndim())
04181 throw ImageDimensionException("fouriergridrot_shift2d needs a 2-D image.");
04182 if (!is_complex())
04183 throw ImageFormatException("fouriergridrot_shift2d requires a fourier image");
04184 int nxreal = nx - 2 + int(is_fftodd());
04185 if (nxreal != ny)
04186 throw ImageDimensionException("fouriergridrot_shift2d requires ny == nx(real)");
04187 if (0 != nxreal%2)
04188 throw ImageDimensionException("fouriergridrot_shift2d needs an even image.");
04189 int nxhalf = nxreal/2;
04190 int nyhalf = ny/2;
04191
04192 if (!is_shuffled()) fft_shuffle();
04193
04194 EMData* result = copy_head();
04195 set_array_offsets(0, -nyhalf);
04196 result->set_array_offsets(0, -nyhalf);
04197
04198 ang = ang*(float)DGR_TO_RAD;
04199 float cang = cos(ang);
04200 float sang = sin(ang);
04201 float temp = -2.0f*M_PI/nxreal;
04202 for (int iy = -nyhalf; iy < nyhalf; iy++) {
04203 float ycang = iy*cang;
04204 float ysang = iy*sang;
04205 for (int ix = 0; ix <= nxhalf; ix++) {
04206 float nuxold = ix*cang - ysang;
04207 float nuyold = ix*sang + ycang;
04208 result->cmplx(ix,iy) = Util::extractpoint2(nx, ny, nuxold, nuyold, this, kb);
04209
04210 float phase_ang = temp*(sx*ix+sy*iy);
04211 result->cmplx(ix,iy) *= complex<float>(cos(phase_ang), sin(phase_ang));
04212 }
04213 }
04214 result->set_array_offsets();
04215 result->fft_shuffle();
04216 result->update();
04217 set_array_offsets();
04218 fft_shuffle();
04219 return result;
04220 }
04221
04222 #undef QUADPI
04223 #undef DGR_TO_RAD
04224
04225 void EMData::divkbsinh(const Util::KaiserBessel& kb) {
04226
04227 if (is_complex())
04228 throw ImageFormatException("divkbsinh requires a real image.");
04229 vector<int> saved_offsets = get_array_offsets();
04230 set_array_offsets(0,0,0);
04231
04232
04233
04234
04235 for (int iz=0; iz < nz; iz++) {
04236 float wz = kb.sinhwin(static_cast<float>(iz-nz/2));
04237 for (int iy=0; iy < ny; iy++) {
04238 float wy = kb.sinhwin(static_cast<float>(iy-ny/2));
04239 for (int ix=0; ix < nx; ix++) {
04240 float wx = kb.sinhwin(static_cast<float>(ix-nx/2));
04241 float w = wx*wy*wz;
04242 (*this)(ix,iy,iz) /= w;
04243 }
04244 }
04245 }
04246 set_array_offsets(saved_offsets);
04247 }
04248
04249 void EMData::divkbsinh_rect(const Util::KaiserBessel& kbx, const Util::KaiserBessel& kby, const Util::KaiserBessel& kbz) {
04250
04251 if (is_complex())
04252 throw ImageFormatException("divkbsinh requires a real image.");
04253 vector<int> saved_offsets = get_array_offsets();
04254 set_array_offsets(0,0,0);
04255
04256
04257
04258
04259 for (int iz=0; iz < nz; iz++) {
04260 float wz = kbz.sinhwin(static_cast<float>(iz-nz/2));
04261 for (int iy=0; iy < ny; iy++) {
04262 float wy = kby.sinhwin(static_cast<float>(iy-ny/2));
04263 for (int ix=0; ix < nx; ix++) {
04264 float wx = kbx.sinhwin(static_cast<float>(ix-nx/2));
04265 float w = wx*wy*wz;
04266 (*this)(ix,iy,iz) /= w;
04267 }
04268 }
04269 }
04270
04271 set_array_offsets(saved_offsets);
04272 }
04273
04274
04275
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 EMData* EMData::extract_plane(const Transform& tf, Util::KaiserBessel& kb) {
04301 if (!is_complex())
04302 throw ImageFormatException("extractplane requires a complex image");
04303 if (nx%2 != 0)
04304 throw ImageDimensionException("extractplane requires nx to be even");
04305 int nxreal = nx - 2;
04306 if (nxreal != ny || nxreal != nz)
04307 throw ImageDimensionException("extractplane requires ny == nx == nz");
04308
04309 EMData* res = new EMData();
04310 res->set_size(nx,ny,1);
04311 res->to_zero();
04312 res->set_complex(true);
04313 res->set_fftodd(false);
04314 res->set_fftpad(true);
04315 res->set_ri(true);
04316
04317 int n = nxreal;
04318 int nhalf = n/2;
04319 vector<int> saved_offsets = get_array_offsets();
04320 set_array_offsets(0,-nhalf,-nhalf);
04321 res->set_array_offsets(0,-nhalf,0);
04322
04323 int kbsize = kb.get_window_size();
04324 int kbmin = -kbsize/2;
04325 int kbmax = -kbmin;
04326 float* wy0 = new float[kbmax - kbmin + 1];
04327 float* wy = wy0 - kbmin;
04328 float* wx0 = new float[kbmax - kbmin + 1];
04329 float* wx = wx0 - kbmin;
04330 float* wz0 = new float[kbmax - kbmin + 1];
04331 float* wz = wz0 - kbmin;
04332 float rim = nhalf*float(nhalf);
04333 int count = 0;
04334 float wsum = 0.f;
04335 Transform tftrans = tf;
04336 tftrans.invert();
04337 for (int jy = -nhalf; jy < nhalf; jy++)
04338 {
04339 for (int jx = 0; jx <= nhalf; jx++)
04340 {
04341 Vec3f nucur((float)jx, (float)jy, 0.f);
04342 Vec3f nunew = tftrans*nucur;
04343 float xnew = nunew[0], ynew = nunew[1], znew = nunew[2];
04344 if (xnew*xnew+ynew*ynew+znew*znew <= rim)
04345 {
04346 count++;
04347 std::complex<float> btq(0.f,0.f);
04348 bool flip = false;
04349 if (xnew < 0.f) {
04350 flip = true;
04351 xnew = -xnew;
04352 ynew = -ynew;
04353 znew = -znew;
04354 }
04355 int ixn = int(Util::round(xnew));
04356 int iyn = int(Util::round(ynew));
04357 int izn = int(Util::round(znew));
04358
04359 for (int i=kbmin; i <= kbmax; i++) {
04360 int izp = izn + i;
04361 wz[i] = kb.i0win_tab(znew - izp);
04362 int iyp = iyn + i;
04363 wy[i] = kb.i0win_tab(ynew - iyp);
04364 int ixp = ixn + i;
04365 wx[i] = kb.i0win_tab(xnew - ixp);
04366
04367 }
04368
04369 int lnbz = 0;
04370 for (int iz = kbmin; iz <= -1; iz++) {
04371 if (wz[iz] != 0.f) {
04372 lnbz = iz;
04373 break;
04374 }
04375 }
04376 int lnez = 0;
04377 for (int iz = kbmax; iz >= 1; iz--) {
04378 if (wz[iz] != 0.f) {
04379 lnez = iz;
04380 break;
04381 }
04382 }
04383 int lnby = 0;
04384 for (int iy = kbmin; iy <= -1; iy++) {
04385 if (wy[iy] != 0.f) {
04386 lnby = iy;
04387 break;
04388 }
04389 }
04390 int lney = 0;
04391 for (int iy = kbmax; iy >= 1; iy--) {
04392 if (wy[iy] != 0.f) {
04393 lney = iy;
04394 break;
04395 }
04396 }
04397 int lnbx = 0;
04398 for (int ix = kbmin; ix <= -1; ix++) {
04399 if (wx[ix] != 0.f) {
04400 lnbx = ix;
04401 break;
04402 }
04403 }
04404 int lnex = 0;
04405 for (int ix = kbmax; ix >= 1; ix--) {
04406 if (wx[ix] != 0.f) {
04407 lnex = ix;
04408 break;
04409 }
04410 }
04411 if (ixn >= -kbmin && ixn <= nhalf-1-kbmax
04412 && iyn >= -nhalf-kbmin && iyn <= nhalf-1-kbmax
04413 && izn >= -nhalf-kbmin && izn <= nhalf-1-kbmax) {
04414
04415 for (int lz = lnbz; lz <= lnez; lz++) {
04416 int izp = izn + lz;
04417 for (int ly=lnby; ly<=lney; ly++) {
04418 int iyp = iyn + ly;
04419 float ty = wz[lz]*wy[ly];
04420 for (int lx=lnbx; lx<=lnex; lx++) {
04421 int ixp = ixn + lx;
04422 float wg = wx[lx]*ty;
04423 btq += cmplx(ixp,iyp,izp)*wg;
04424 wsum += wg;
04425 }
04426 }
04427 }
04428 } else {
04429
04430 for (int lz = lnbz; lz <= lnez; lz++) {
04431 int izp = izn + lz;
04432 for (int ly=lnby; ly<=lney; ly++) {
04433 int iyp = iyn + ly;
04434 float ty = wz[lz]*wy[ly];
04435 for (int lx=lnbx; lx<=lnex; lx++) {
04436 int ixp = ixn + lx;
04437 float wg = wx[lx]*ty;
04438 bool mirror = false;
04439 int ixt(ixp), iyt(iyp), izt(izp);
04440 if (ixt > nhalf || ixt < -nhalf) {
04441 ixt = Util::sgn(ixt)
04442 *(n - abs(ixt));
04443 iyt = -iyt;
04444 izt = -izt;
04445 mirror = !mirror;
04446 }
04447 if (iyt >= nhalf || iyt < -nhalf) {
04448 if (ixt != 0) {
04449 ixt = -ixt;
04450 iyt = Util::sgn(iyt)
04451 *(n - abs(iyt));
04452 izt = -izt;
04453 mirror = !mirror;
04454 } else {
04455 iyt -= n*Util::sgn(iyt);
04456 }
04457 }
04458 if (izt >= nhalf || izt < -nhalf) {
04459 if (ixt != 0) {
04460 ixt = -ixt;
04461 iyt = -iyt;
04462 izt = Util::sgn(izt)
04463 *(n - abs(izt));
04464 mirror = !mirror;
04465 } else {
04466 izt -= Util::sgn(izt)*n;
04467 }
04468 }
04469 if (ixt < 0) {
04470 ixt = -ixt;
04471 iyt = -iyt;
04472 izt = -izt;
04473 mirror = !mirror;
04474 }
04475 if (iyt == nhalf) iyt = -nhalf;
04476 if (izt == nhalf) izt = -nhalf;
04477 if (mirror) btq += conj(cmplx(ixt,iyt,izt))*wg;
04478 else btq += cmplx(ixt,iyt,izt)*wg;
04479 wsum += wg;
04480 }
04481 }
04482 }
04483 }
04484 if (flip) res->cmplx(jx,jy) = conj(btq);
04485 else res->cmplx(jx,jy) = btq;
04486 }
04487 }
04488 }
04489 for (int jy = -nhalf; jy < nhalf; jy++)
04490 for (int jx = 0; jx <= nhalf; jx++)
04491 res->cmplx(jx,jy) *= count/wsum;
04492 delete[] wx0; delete[] wy0; delete[] wz0;
04493 set_array_offsets(saved_offsets);
04494 res->set_array_offsets(0,0,0);
04495 res->set_shuffled(true);
04496 return res;
04497 }
04498
04499
04501
04502
04503
04504 EMData* EMData::extract_plane_rect(const Transform& tf, Util::KaiserBessel& kbx,Util::KaiserBessel& kby, Util::KaiserBessel& kbz) {
04505
04506
04507 if (!is_complex())
04508 throw ImageFormatException("extractplane requires a complex image");
04509 if (nx%2 != 0)
04510 throw ImageDimensionException("extractplane requires nx to be even");
04511
04512 int nxfromz = nz+2;
04513 int nxcircal = nxfromz - 2;
04514 EMData* res = new EMData();
04515 res->set_size(nxfromz,nz,1);
04516 res->to_zero();
04517 res->set_complex(true);
04518 res->set_fftodd(false);
04519 res->set_fftpad(true);
04520 res->set_ri(true);
04521
04522 int n = nxcircal;
04523 int nhalf = n/2;
04524 int nxhalf = (nx-2)/2;
04525 int nyhalf = ny/2;
04526 int nzhalf = nz/2;
04527
04528 vector<int> saved_offsets = get_array_offsets();
04529 set_array_offsets(0, -nyhalf, -nzhalf);
04530 res->set_array_offsets(0,-nhalf,0);
04531
04532 int kbxsize = kbx.get_window_size();
04533 int kbxmin = -kbxsize/2;
04534 int kbxmax = -kbxmin;
04535
04536 int kbysize = kby.get_window_size();
04537 int kbymin = -kbysize/2;
04538 int kbymax = -kbymin;
04539
04540 int kbzsize = kbz.get_window_size();
04541 int kbzmin = -kbzsize/2;
04542 int kbzmax = -kbzmin;
04543
04544
04545 float* wy0 = new float[kbymax - kbymin + 1];
04546 float* wy = wy0 - kbymin;
04547 float* wx0 = new float[kbxmax - kbxmin + 1];
04548 float* wx = wx0 - kbxmin;
04549 float* wz0 = new float[kbzmax - kbzmin + 1];
04550 float* wz = wz0 - kbzmin;
04551 float rim = nhalf*float(nhalf);
04552 int count = 0;
04553 float wsum = 0.f;
04554 Transform tftrans = tf;
04555 tftrans.invert();
04556 float xratio=float(nx-2)/float(nz);
04557 float yratio=float(ny)/float(nz);
04558
04559 for (int jy = -nhalf; jy < nhalf; jy++)
04560 {
04561 for (int jx = 0; jx <= nhalf; jx++)
04562 {
04563 Vec3f nucur((float)jx, (float)jy, 0.f);
04564 Vec3f nunew = tftrans*nucur;
04565 float xnew = nunew[0]*xratio, ynew = nunew[1]*yratio, znew = nunew[2];
04566
04567 if (nunew[0]*nunew[0]+nunew[1]*nunew[1]+nunew[2]*nunew[2] <= rim)
04568 {
04569 count++;
04570 std::complex<float> btq(0.f,0.f);
04571 bool flip = false;
04572 if (xnew < 0.f) {
04573 flip = true;
04574 xnew = -xnew;
04575 ynew = -ynew;
04576 znew = -znew;
04577 }
04578 int ixn = int(Util::round(xnew));
04579 int iyn = int(Util::round(ynew));
04580 int izn = int(Util::round(znew));
04581
04582 for (int i=kbzmin; i <= kbzmax; i++) {
04583 int izp = izn + i;
04584 wz[i] = kbz.i0win_tab(znew - izp);
04585 }
04586 for (int i=kbymin; i <= kbymax; i++) {
04587 int iyp = iyn + i;
04588 wy[i] = kby.i0win_tab(ynew - iyp);
04589 }
04590 for (int i=kbxmin; i <= kbxmax; i++) {
04591 int ixp = ixn + i;
04592 wx[i] = kbx.i0win_tab(xnew - ixp);
04593 }
04594
04595
04596
04597
04598 int lnbz = 0;
04599 for (int iz = kbzmin; iz <= -1; iz++) {
04600 if (wz[iz] != 0.f) {
04601 lnbz = iz;
04602 break;
04603 }
04604 }
04605 int lnez = 0;
04606 for (int iz = kbzmax; iz >= 1; iz--) {
04607 if (wz[iz] != 0.f) {
04608 lnez = iz;
04609 break;
04610 }
04611 }
04612 int lnby = 0;
04613 for (int iy = kbymin; iy <= -1; iy++) {
04614 if (wy[iy] != 0.f) {
04615 lnby = iy;
04616 break;
04617 }
04618 }
04619 int lney = 0;
04620 for (int iy = kbymax; iy >= 1; iy--) {
04621 if (wy[iy] != 0.f) {
04622 lney = iy;
04623 break;
04624 }
04625 }
04626 int lnbx = 0;
04627 for (int ix = kbxmin; ix <= -1; ix++) {
04628 if (wx[ix] != 0.f) {
04629 lnbx = ix;
04630 break;
04631 }
04632 }
04633 int lnex = 0;
04634 for (int ix = kbxmax; ix >= 1; ix--) {
04635 if (wx[ix] != 0.f) {
04636 lnex = ix;
04637 break;
04638 }
04639 }
04640 if (ixn >= -kbxmin && ixn <= nxhalf-1-kbxmax
04641 && iyn >= -nyhalf-kbymin && iyn <= nyhalf-1-kbymax
04642 && izn >= -nzhalf-kbzmin && izn <= nzhalf-1-kbzmax) {
04643
04644 for (int lz = lnbz; lz <= lnez; lz++) {
04645 int izp = izn + lz;
04646 for (int ly=lnby; ly<=lney; ly++) {
04647 int iyp = iyn + ly;
04648 float ty = wz[lz]*wy[ly];
04649 for (int lx=lnbx; lx<=lnex; lx++) {
04650 int ixp = ixn + lx;
04651 float wg = wx[lx]*ty;
04652 btq += cmplx(ixp,iyp,izp)*wg;
04653 wsum += wg;
04654 }
04655 }
04656 }
04657 }
04658 else {
04659
04660 for (int lz = lnbz; lz <= lnez; lz++) {
04661 int izp = izn + lz;
04662 for (int ly=lnby; ly<=lney; ly++) {
04663 int iyp = iyn + ly;
04664 float ty = wz[lz]*wy[ly];
04665 for (int lx=lnbx; lx<=lnex; lx++) {
04666 int ixp = ixn + lx;
04667 float wg = wx[lx]*ty;
04668 bool mirror = false;
04669 int ixt(ixp), iyt(iyp), izt(izp);
04670 if (ixt > nxhalf || ixt < -nxhalf) {
04671 ixt = Util::sgn(ixt)
04672 *(nx-2-abs(ixt));
04673 iyt = -iyt;
04674 izt = -izt;
04675 mirror = !mirror;
04676 }
04677 if (iyt >= nyhalf || iyt < -nyhalf) {
04678 if (ixt != 0) {
04679 ixt = -ixt;
04680 iyt = Util::sgn(iyt)
04681 *(ny - abs(iyt));
04682 izt = -izt;
04683 mirror = !mirror;
04684 } else {
04685 iyt -= ny*Util::sgn(iyt);
04686 }
04687 }
04688 if (izt >= nzhalf || izt < -nzhalf) {
04689 if (ixt != 0) {
04690 ixt = -ixt;
04691 iyt = -iyt;
04692 izt = Util::sgn(izt)
04693 *(nz - abs(izt));
04694 mirror = !mirror;
04695 } else {
04696 izt -= Util::sgn(izt)*nz;
04697 }
04698 }
04699 if (ixt < 0) {
04700 ixt = -ixt;
04701 iyt = -iyt;
04702 izt = -izt;
04703 mirror = !mirror;
04704 }
04705 if (iyt == nyhalf) iyt = -nyhalf;
04706 if (izt == nzhalf) izt = -nzhalf;
04707 if (mirror) btq += conj(cmplx(ixt,iyt,izt))*wg;
04708 else btq += cmplx(ixt,iyt,izt)*wg;
04709 wsum += wg;
04710 }
04711 }
04712 }
04713 }
04714 if (flip) res->cmplx(jx,jy) = conj(btq);
04715 else res->cmplx(jx,jy) = btq;
04716 }
04717 }
04718 }
04719 for (int jy = -nhalf; jy < nhalf; jy++)
04720 for (int jx = 0; jx <= nhalf; jx++)
04721 res->cmplx(jx,jy) *= count/wsum;
04722 delete[] wx0; delete[] wy0; delete[] wz0;
04723 set_array_offsets(saved_offsets);
04724 res->set_array_offsets(0,0,0);
04725 res->set_shuffled(true);
04726 return res;
04727 }
04728
04729
04730
04731 EMData* EMData::extract_plane_rect_fast(const Transform& tf, Util::KaiserBessel& kbx,Util::KaiserBessel& kby, Util::KaiserBessel& kbz) {
04732
04733
04734
04735 if (!is_complex())
04736 throw ImageFormatException("extractplane requires a complex image");
04737 if (nx%2 != 0)
04738 throw ImageDimensionException("extractplane requires nx to be even");
04739
04740 int nxfromz=nz+2;
04741 int nxcircal = nxfromz - 2;
04742
04743
04744 float xratio=float(nx-2)/float(nz);
04745 float yratio=float(ny)/float(nz);
04746 Vec3f axis_newx,axis_newy;
04747 axis_newx[0] = xratio*0.5f*nz*tf[0][0];
04748 axis_newx[1] = yratio*0.5f*nz*tf[0][1];
04749 axis_newx[2] = 0.5f*nz*tf[0][2];
04750
04751
04752 float ellipse_length_x=std::sqrt(axis_newx[0]*axis_newx[0]+axis_newx[1]*axis_newx[1]+axis_newx[2]*axis_newx[2]);
04753
04754 int ellipse_length_x_int=int(ellipse_length_x);
04755 float ellipse_step_x=0.5f*nz/float(ellipse_length_x_int);
04756 float xscale=ellipse_step_x;
04757
04758 axis_newy[0] = xratio*0.5f*nz*tf[1][0];
04759 axis_newy[1] = yratio*0.5f*nz*tf[1][1];
04760 axis_newy[2] = 0.5f*nz*tf[1][2];
04761
04762
04763 float ellipse_length_y=std::sqrt(axis_newy[0]*axis_newy[0]+axis_newy[1]*axis_newy[1]+axis_newy[2]*axis_newy[2]);
04764 int ellipse_length_y_int=int(ellipse_length_y);
04765 float ellipse_step_y=0.5f*nz/float(ellipse_length_y_int);
04766 float yscale=ellipse_step_y;
04767
04768 int nx_e=ellipse_length_x_int*2;
04769 int ny_e=ellipse_length_y_int*2;
04770 int nx_ec=nx_e+2;
04771
04772 EMData* res = new EMData();
04773 res->set_size(nx_ec,ny_e,1);
04774 res->to_zero();
04775 res->set_complex(true);
04776 res->set_fftodd(false);
04777 res->set_fftpad(true);
04778 res->set_ri(true);
04779
04780
04781
04782 int n = nxcircal;
04783 int nhalf = n/2;
04784 int nhalfx_e = nx_e/2;
04785 int nhalfy_e = ny_e/2;
04786 int nxhalf=(nx-2)/2;
04787 int nyhalf=ny/2;
04788 int nzhalf=nz/2;
04789
04790 vector<int> saved_offsets = get_array_offsets();
04791 set_array_offsets(0,-nyhalf,-nzhalf);
04792 res->set_array_offsets(0,-nhalfy_e,0);
04793
04794 int kbxsize = kbx.get_window_size();
04795 int kbxmin = -kbxsize/2;
04796 int kbxmax = -kbxmin;
04797
04798 int kbysize = kby.get_window_size();
04799 int kbymin = -kbysize/2;
04800 int kbymax = -kbymin;
04801
04802 int kbzsize = kbz.get_window_size();
04803 int kbzmin = -kbzsize/2;
04804 int kbzmax = -kbzmin;
04805
04806
04807 float* wy0 = new float[kbymax - kbymin + 1];
04808 float* wy = wy0 - kbymin;
04809 float* wx0 = new float[kbxmax - kbxmin + 1];
04810 float* wx = wx0 - kbxmin;
04811 float* wz0 = new float[kbzmax - kbzmin + 1];
04812 float* wz = wz0 - kbzmin;
04813 float rim = nhalf*float(nhalf);
04814 int count = 0;
04815 float wsum = 0.f;
04816 Transform tftrans = tf;
04817 tftrans.invert();
04818
04819
04820 for (int jy = -nhalfy_e; jy < nhalfy_e; jy++)
04821 {
04822 for (int jx = 0; jx <= nhalfx_e; jx++)
04823 {
04824 Vec3f nucur((float)jx, (float)jy, 0.f);
04825 nucur[0]=nucur[0]*xscale;nucur[1]=nucur[1]*yscale;;
04826 Vec3f nunew = tftrans*nucur;
04827 float xnew = nunew[0]*xratio, ynew = nunew[1]*yratio, znew = nunew[2];
04828
04829 if (nunew[0]*nunew[0]+nunew[1]*nunew[1]+nunew[2]*nunew[2] <= rim)
04830 {
04831 count++;
04832 std::complex<float> btq(0.f,0.f);
04833 bool flip = false;
04834 if (xnew < 0.f) {
04835 flip = true;
04836 xnew = -xnew;
04837 ynew = -ynew;
04838 znew = -znew;
04839 }
04840 int ixn = int(Util::round(xnew));
04841 int iyn = int(Util::round(ynew));
04842 int izn = int(Util::round(znew));
04843
04844 for (int i=kbzmin; i <= kbzmax; i++) {
04845 int izp = izn + i;
04846 wz[i] = kbz.i0win_tab(znew - izp);
04847 }
04848 for (int i=kbymin; i <= kbymax; i++) {
04849 int iyp = iyn + i;
04850 wy[i] = kby.i0win_tab(ynew - iyp);
04851 }
04852 for (int i=kbxmin; i <= kbxmax; i++) {
04853 int ixp = ixn + i;
04854 wx[i] = kbx.i0win_tab(xnew - ixp);
04855 }
04856
04857
04858
04859
04860 int lnbz = 0;
04861 for (int iz = kbzmin; iz <= -1; iz++) {
04862 if (wz[iz] != 0.f) {
04863 lnbz = iz;
04864 break;
04865 }
04866 }
04867 int lnez = 0;
04868 for (int iz = kbzmax; iz >= 1; iz--) {
04869 if (wz[iz] != 0.f) {
04870 lnez = iz;
04871 break;
04872 }
04873 }
04874 int lnby = 0;
04875 for (int iy = kbymin; iy <= -1; iy++) {
04876 if (wy[iy] != 0.f) {
04877 lnby = iy;
04878 break;
04879 }
04880 }
04881 int lney = 0;
04882 for (int iy = kbymax; iy >= 1; iy--) {
04883 if (wy[iy] != 0.f) {
04884 lney = iy;
04885 break;
04886 }
04887 }
04888 int lnbx = 0;
04889 for (int ix = kbxmin; ix <= -1; ix++) {
04890 if (wx[ix] != 0.f) {
04891 lnbx = ix;
04892 break;
04893 }
04894 }
04895 int lnex = 0;
04896 for (int ix = kbxmax; ix >= 1; ix--) {
04897 if (wx[ix] != 0.f) {
04898 lnex = ix;
04899 break;
04900 }
04901 }
04902 if (ixn >= -kbxmin && ixn <= nxhalf-1-kbxmax
04903 && iyn >= -nyhalf-kbymin && iyn <= nyhalf-1-kbymax
04904 && izn >= -nzhalf-kbzmin && izn <= nzhalf-1-kbzmax) {
04905
04906 for (int lz = lnbz; lz <= lnez; lz++) {
04907 int izp = izn + lz;
04908 for (int ly=lnby; ly<=lney; ly++) {
04909 int iyp = iyn + ly;
04910 float ty = wz[lz]*wy[ly];
04911 for (int lx=lnbx; lx<=lnex; lx++) {
04912 int ixp = ixn + lx;
04913 float wg = wx[lx]*ty;
04914 btq += cmplx(ixp,iyp,izp)*wg;
04915 wsum += wg;
04916 }
04917 }
04918 }
04919 }
04920 else {
04921
04922 for (int lz = lnbz; lz <= lnez; lz++) {
04923 int izp = izn + lz;
04924 for (int ly=lnby; ly<=lney; ly++) {
04925 int iyp = iyn + ly;
04926 float ty = wz[lz]*wy[ly];
04927 for (int lx=lnbx; lx<=lnex; lx++) {
04928 int ixp = ixn + lx;
04929 float wg = wx[lx]*ty;
04930 bool mirror = false;
04931 int ixt(ixp), iyt(iyp), izt(izp);
04932 if (ixt > nxhalf || ixt < -nxhalf) {
04933 ixt = Util::sgn(ixt)
04934 *(nx-2-abs(ixt));
04935 iyt = -iyt;
04936 izt = -izt;
04937 mirror = !mirror;
04938 }
04939 if (iyt >= nyhalf || iyt < -nyhalf) {
04940 if (ixt != 0) {
04941 ixt = -ixt;
04942 iyt = Util::sgn(iyt)
04943 *(ny - abs(iyt));
04944 izt = -izt;
04945 mirror = !mirror;
04946 } else {
04947 iyt -= ny*Util::sgn(iyt);
04948 }
04949 }
04950 if (izt >= nzhalf || izt < -nzhalf) {
04951 if (ixt != 0) {
04952 ixt = -ixt;
04953 iyt = -iyt;
04954 izt = Util::sgn(izt)
04955 *(nz - abs(izt));
04956 mirror = !mirror;
04957 } else {
04958 izt -= Util::sgn(izt)*nz;
04959 }
04960 }
04961 if (ixt < 0) {
04962 ixt = -ixt;
04963 iyt = -iyt;
04964 izt = -izt;
04965 mirror = !mirror;
04966 }
04967 if (iyt == nyhalf) iyt = -nyhalf;
04968 if (izt == nzhalf) izt = -nzhalf;
04969 if (mirror) btq += conj(cmplx(ixt,iyt,izt))*wg;
04970 else btq += cmplx(ixt,iyt,izt)*wg;
04971 wsum += wg;
04972 }
04973 }
04974 }
04975 }
04976 if (flip) res->cmplx(jx,jy) = conj(btq);
04977 else res->cmplx(jx,jy) = btq;
04978 }
04979 }
04980 }
04981 for (int jy = -nhalfy_e; jy < nhalfy_e; jy++)
04982 for (int jx = 0; jx <= nhalfx_e; jx++)
04983 res->cmplx(jx,jy) *= count/wsum;
04984 delete[] wx0; delete[] wy0; delete[] wz0;
04985 set_array_offsets(saved_offsets);
04986 res->set_array_offsets(0,0,0);
04987 res->set_shuffled(true);
04988 return res;
04989 }
04990
04991
04992
04993
04994 bool EMData::peakcmp(const Pixel& p1, const Pixel& p2) {
04995 return (p1.value > p2.value);
04996 }
04997
04998 ostream& operator<< (ostream& os, const Pixel& peak) {
04999 os << peak.x << peak.y << peak.z << peak.value;
05000 return os;
05001 }
05002
05003
05004
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 vector<float> EMData::peak_search(int ml, float invert) {
05068
05069 EMData& buf = *this;
05070 vector<Pixel> peaks;
05071 int img_dim;
05072 int i, j, k;
05073 int i__1, i__2;
05074 int j__1, j__2;
05075
05076 bool peak_check;
05077 img_dim=buf.get_ndim();
05078 vector<int> ix, jy, kz;
05079 vector<float>res;
05080 int nx = buf.get_xsize();
05081 int ny = buf.get_ysize();
05082 int nz = buf.get_zsize();
05083 if(invert <= 0.0f) invert=-1.0f;
05084 else invert=1.0f ;
05085 int count = 0;
05086 switch (img_dim) {
05087 case(1):
05088 for(i=0;i<=nx-1;++i) {
05089 i__1 = (i-1+nx)%nx;
05090 i__2 = (i+1)%nx;
05091
05092
05093
05094 float qbf = buf(i)*invert;
05095 peak_check = qbf > buf(i__1)*invert && qbf > buf(i__2)*invert;
05096 if(peak_check) {
05097 if(count < ml) {
05098 count++;
05099 peaks.push_back( Pixel(i, 0, 0, qbf) );
05100 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
05101 } else {
05102 if( qbf > (peaks.back()).value ) {
05103
05104 peaks.pop_back();
05105 peaks.push_back( Pixel(i, 0, 0, qbf) );
05106 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
05107 }
05108 }
05109 }
05110 }
05111 break;
05112 case(2):
05113
05114
05115
05116
05117
05118
05119
05120
05121 for(j=1;j<=ny-2;++j) {
05122 j__1 = j-1;
05123 j__2 = j+1;
05124 for(i=1;i<=nx-2;++i) {
05125 i__1 = i-1;
05126 i__2 = i+1;
05127 float qbf = buf(i,j)*invert;
05128 peak_check = (qbf > buf(i,j__1)*invert) && (qbf > buf(i,j__2)*invert);
05129 if(peak_check) {
05130 peak_check = (qbf > buf(i__1,j)*invert) && (qbf > buf(i__2,j)*invert);
05131 if(peak_check) {
05132 peak_check = (qbf > buf(i__1,j__1)*invert) && (qbf > buf(i__1,j__2)*invert);
05133 if(peak_check) {
05134 peak_check = (qbf > buf(i__2,j__1)*invert) && (qbf > buf(i__2,j__2)*invert);
05135 if(peak_check) {
05136 if(count < ml) {
05137 count++;
05138 peaks.push_back( Pixel(i, j, 0, qbf) );
05139 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
05140 } else {
05141 if( qbf > (peaks.back()).value ) {
05142
05143 peaks.pop_back();
05144 peaks.push_back( Pixel(i, j, 0, qbf) );
05145 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
05146 }
05147 }
05148 }
05149 }
05150 }
05151 }
05152 }
05153 }
05154 break;
05155 case(3):
05156
05157
05158
05159
05160
05161
05162
05163
05164
05165
05166
05167
05168
05169
05170
05171
05172
05173
05174
05175
05176 for(k=1; k<=nz-2; ++k) {
05177 kz.clear();
05178 kz.push_back(k-1);
05179 kz.push_back(k);
05180 kz.push_back(k+1);
05181 for(j=1; j<=ny-2; ++j) {
05182 jy.clear();
05183 jy.push_back(j-1);
05184 jy.push_back(j);
05185 jy.push_back(j+1);
05186 for(i=1; i<=nx-2; ++i) {
05187 ix.clear();
05188 ix.push_back(i-1);
05189 ix.push_back(i);
05190 ix.push_back(i+1);
05191 float qbf = buf(i,j,k)*invert;
05192 peak_check = qbf > buf(ix[0],jy[0],kz[0])*invert;
05193 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[0])*invert;
05194 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[0])*invert;
05195 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[0])*invert;
05196 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[0])*invert;
05197 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[0])*invert;
05198 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[0])*invert;
05199 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[0])*invert;
05200 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[0])*invert;
05201 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[1])*invert;
05202 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[1])*invert;
05203 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[1])*invert;
05204 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[1])*invert;
05205
05206 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[1])*invert;
05207 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[1])*invert;
05208 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[1])*invert;
05209 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[1])*invert;
05210 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[2])*invert;
05211 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[2])*invert;
05212 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[2])*invert;
05213 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[2])*invert;
05214 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[2])*invert;
05215 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[2])*invert;
05216 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[2])*invert;
05217 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[2])*invert;
05218 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[2])*invert;
05219 if(peak_check) {
05220 if(count < ml) {
05221 count++;
05222 peaks.push_back( Pixel(i, j, k, qbf) );
05223 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
05224 } else {
05225 if( qbf > (peaks.back()).value ) {
05226
05227
05228 peaks.pop_back();
05229 peaks.push_back( Pixel(i, j, k, qbf) );
05230 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
05231 }
05232 }
05233 }
05234 }}}}}}}}}}}}}}}}}}}}}}}}}
05235 }
05236 }
05237 }
05238
05239
05240 for(k=1; k<=nz-2; ++k) {
05241 kz.clear();
05242 kz.push_back(k-1);
05243 kz.push_back(k);
05244 kz.push_back(k+1);
05245 for(j=1; j<=ny-2; ++j) {
05246 jy.clear();
05247 jy.push_back(j-1);
05248 jy.push_back(j);
05249 jy.push_back(j+1);
05250 for(i=0; i<=0; ++i) {
05251 ix.clear();
05252 ix.push_back(nx-1);
05253 ix.push_back(i);
05254 ix.push_back(i+1);
05255 float qbf = buf(i,j,k)*invert;
05256 peak_check = qbf > buf(ix[0],jy[0],kz[0])*invert;
05257 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[0])*invert;
05258 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[0])*invert;
05259 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[0])*invert;
05260 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[0])*invert;
05261 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[0])*invert;
05262 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[0])*invert;
05263 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[0])*invert;
05264 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[0])*invert;
05265 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[1])*invert;
05266 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[1])*invert;
05267 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[1])*invert;
05268 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[1])*invert;
05269
05270 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[1])*invert;
05271 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[1])*invert;
05272 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[1])*invert;
05273 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[1])*invert;
05274 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[2])*invert;
05275 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[2])*invert;
05276 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[2])*invert;
05277 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[2])*invert;
05278 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[2])*invert;
05279 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[2])*invert;
05280 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[2])*invert;
05281 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[2])*invert;
05282 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[2])*invert;
05283 if(peak_check) {
05284 if(count < ml) {
05285 count++;
05286 peaks.push_back( Pixel(i, j, k, qbf) );
05287 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
05288 } else {
05289 if( qbf > (peaks.back()).value ) {
05290
05291
05292 peaks.pop_back();
05293 peaks.push_back( Pixel(i, j, k, qbf) );
05294 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
05295 }
05296 }
05297 }
05298 }}}}}}}}}}}}}}}}}}}}}}}}}
05299 }
05300 for(i=nx-1; i<=nx-1; ++i) {
05301 ix.clear();
05302 ix.push_back(i-1);
05303 ix.push_back(i);
05304 ix.push_back(0);
05305 float qbf = buf(i,j,k)*invert;
05306 peak_check = qbf > buf(ix[0],jy[0],kz[0])*invert;
05307 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[0])*invert;
05308 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[0])*invert;
05309 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[0])*invert;
05310 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[0])*invert;
05311 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[0])*invert;
05312 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[0])*invert;
05313 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[0])*invert;
05314 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[0])*invert;
05315 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[1])*invert;
05316 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[1])*invert;
05317 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[1])*invert;
05318 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[1])*invert;
05319
05320 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[1])*invert;
05321 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[1])*invert;
05322 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[1])*invert;
05323 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[1])*invert;
05324 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[2])*invert;
05325 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[2])*invert;
05326 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[2])*invert;
05327 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[2])*invert;
05328 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[2])*invert;
05329 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[2])*invert;
05330 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[2])*invert;
05331 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[2])*invert;
05332 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[2])*invert;
05333 if(peak_check) {
05334 if(count < ml) {
05335 count++;
05336 peaks.push_back( Pixel(i, j, k, qbf) );
05337 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
05338 } else {
05339 if( qbf > (peaks.back()).value ) {
05340
05341
05342 peaks.pop_back();
05343 peaks.push_back( Pixel(i, j, k, qbf) );
05344 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
05345 }
05346 }
05347 }
05348 }}}}}}}}}}}}}}}}}}}}}}}}}
05349 }
05350 }
05351 }
05352 break;
05353
05354
05355
05356
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 if (peaks.begin() != peaks.end()) {
05795
05796 sort(peaks.begin(), peaks.end(), peakcmp);
05797
05798 int count = 0;
05799
05800 float xval = (*peaks.begin()).value;
05801
05802 for (vector<Pixel>::iterator it = peaks.begin(); it != peaks.end(); it++) {
05803
05804 count++;
05805
05806 if(count <= ml) {
05807
05808 res.push_back((*it).value);
05809 res.push_back(static_cast<float>((*it).x));
05810
05811 if(img_dim > 1) {
05812 res.push_back(static_cast<float>((*it).y));
05813 if(nz > 1) res.push_back(static_cast<float>((*it).z));
05814 }
05815
05816 if(xval != 0.0) res.push_back((*it).value/xval);
05817 else res.push_back((*it).value);
05818 res.push_back((*it).x-float(int(nx/2)));
05819 if(img_dim >1) {
05820 res.push_back((*it).y-float(int(ny/2)));
05821 if(nz>1) res.push_back((*it).z-float(nz/2));
05822 }
05823 }
05824 }
05825 res.insert(res.begin(),1,img_dim);
05826 } else {
05827
05828 res.push_back(buf(0,0,0));
05829 res.insert(res.begin(),1,0.0);
05830 }
05831
05832
05833 return res;
05834 }
05835
05836 #define rdata(i,j,k) rdata[(i-1)+((j-1)+(k-1)*ny)*(size_t)nx]
05837 #define X(i) X[i-1]
05838 #define Y(j) Y[j-1]
05839 #define Z(k) Z[k-1]
05840 vector<float> EMData::phase_cog()
05841 {
05842 vector<float> ph_cntog;
05843 int i=1,j=1,k=1;
05844 float C=0.f,S=0.f,P=0.f,F1=0.f,SNX;
05845 if (get_ndim()==1) {
05846 P = 8*atan(1.0f)/nx;
05847 for (i=1;i<=nx;i++) {
05848 C += cos(P * (i-1)) * rdata(i,j,k);
05849 S += sin(P * (i-1)) * rdata(i,j,k);
05850 }
05851 F1 = atan2(S,C);
05852 if (F1 < 0.0) F1 += 8*atan(1.0f);
05853 SNX = F1/P +1.0f;
05854 SNX = SNX - ((nx/2)+1);
05855 ph_cntog.push_back(SNX);
05856 #ifdef _WIN32
05857 ph_cntog.push_back((float)Util::round(SNX));
05858 #else
05859 ph_cntog.push_back(round(SNX));
05860 #endif //_WIN32
05861 } else if (get_ndim()==2) {
05862 #ifdef _WIN32
05863 float SNY;
05864 float T=0.0f;
05865 vector<float> X;
05866 X.resize(nx);
05867 #else
05868 float SNY,X[nx],T=0.f;
05869 #endif //_WIN32
05870 for ( i=1;i<=nx;i++) X(i)=0.0;
05871 P = 8*atan(1.0f)/ny;
05872 for(j=1;j<=ny;j++) {
05873 T=0.f;
05874 for(i=1;i<=nx;i++) {
05875 T += rdata(i,j,k);
05876 X(i)+=rdata(i,j,k);
05877 }
05878 C += cos(P*(j-1))*T;
05879 S += sin(P*(j-1))*T;
05880 }
05881 F1=atan2(S,C);
05882 if(F1<0.0) F1 += 8*atan(1.0f);
05883 SNY = F1/P +1.0f;
05884 C=0.f; S=0.f;
05885 P = 8*atan(1.0f)/nx;
05886 for(i=1;i<=nx;i++) {
05887 C += cos(P*(i-1))*X(i);
05888 S += sin(P*(i-1))*X(i);
05889 }
05890 F1=atan2(S,C);
05891 if(F1<0.0) F1 += 8*atan(1.0f);
05892 SNX = F1/P +1.0f;
05893 SNX = SNX - ((nx/2)+1);
05894 SNY = SNY - ((ny/2)+1);
05895 ph_cntog.push_back(SNX); ph_cntog.push_back(SNY);
05896 #ifdef _WIN32
05897 ph_cntog.push_back((float)Util::round(SNX)); ph_cntog.push_back((float)Util::round(SNY));
05898 #else
05899 ph_cntog.push_back(round(SNX)); ph_cntog.push_back(round(SNY));
05900 #endif //_WIN32
05901 } else {
05902 #ifdef _WIN32
05903 float val=0.f,sum1=0.f, SNY,SNZ;
05904 vector<float> X;
05905 X.resize(nx);
05906 vector<float> Y;
05907 Y.resize(ny);
05908 vector<float> Z;
05909 Z.resize(nz);
05910 #else
05911 float val=0.f, sum1=0.f, X[nx], Y[ny], Z[nz], SNY, SNZ;
05912 #endif //_WIN32
05913 for (i=1;i<=nx;i++) X(i)=0.0;
05914 for (j=1;j<=ny;j++) Y(j)=0.0;
05915 for (k=1;k<=nz;k++) Z(k)=0.0;
05916 for(k=1;k<=nz;k++) {
05917 for(j=1;j<=ny;j++) {
05918 sum1=0.f;
05919 for(i=1;i<=nx;i++) {
05920 val = rdata(i,j,k);
05921 sum1 += val;
05922 X(i) += val;
05923 }
05924 Y(j) += sum1;
05925 Z(k) += sum1;
05926 }
05927 }
05928 P = 8*atan(1.0f)/nx;
05929 for (i=1;i<=nx;i++) {
05930 C += cos(P*(i-1))*X(i);
05931 S += sin(P*(i-1))*X(i);
05932 }
05933 F1=atan2(S,C);
05934 if(F1<0.0) F1 += 8*atan(1.0f);
05935 SNX = F1/P +1.0f;
05936 C=0.f; S=0.f;
05937 P = 8*atan(1.0f)/ny;
05938 for(j=1;j<=ny;j++) {
05939 C += cos(P*(j-1))*Y(j);
05940 S += sin(P*(j-1))*Y(j);
05941 }
05942 F1=atan2(S,C);
05943 if(F1<0.0) F1 += 8*atan(1.0f);
05944 SNY = F1/P +1.0f;
05945 C=0.f; S=0.f;
05946 P = 8*atan(1.0f)/nz;
05947 for(k=1;k<=nz;k++) {
05948 C += cos(P*(k-1))*Z(k);
05949 S += sin(P*(k-1))*Z(k);
05950 }
05951 F1=atan2(S,C);
05952 if(F1<0.0) F1 += 8*atan(1.0f);
05953 SNZ = F1/P +1.0f;
05954 SNX = SNX - ((nx/2)+1);
05955 SNY = SNY - ((ny/2)+1);
05956 SNZ = SNZ - ((nz/2)+1);
05957 ph_cntog.push_back(SNX); ph_cntog.push_back(SNY); ph_cntog.push_back(SNZ);
05958 #ifdef _WIN32
05959 ph_cntog.push_back((float)Util::round(SNX)); ph_cntog.push_back((float)Util::round(SNY)); ph_cntog.push_back((float)Util::round(SNZ));
05960 #else
05961 ph_cntog.push_back(round(SNX)); ph_cntog.push_back(round(SNY));ph_cntog.push_back(round(SNZ));
05962 #endif
05963 }
05964 return ph_cntog;
05965 }
05966 #undef rdata
05967 #undef X
05968 #undef Y
05969 #undef Z
05970
05971 #define avagadro (6.023*(double)pow(10.0,23.0))
05972 #define density_protein (1.36)
05973 #define R (0.61803399f)
05974 #define C (1.f-R)
05975 float EMData::find_3d_threshold(float mass, float pixel_size)
05976 {
05977
05978 if(get_ndim()!=3)
05979 throw ImageDimensionException("The image should be 3D");
05980
05981
05982
05983 float density_1_mole, vol_1_mole, vol_angstrom;
05984 int vol_voxels;
05985 density_1_mole = static_cast<float>( (mass*1000.0f)/avagadro );
05986 vol_1_mole = static_cast<float>( density_1_mole/density_protein );
05987 vol_angstrom = static_cast<float>( vol_1_mole*(double)pow((double)pow(10.0,8),3) );
05988 vol_voxels = static_cast<int> (vol_angstrom/(double)pow(pixel_size,3));
05989
05990
05991
05992 float thr1 = get_attr("maximum");
05993 float thr3 = get_attr("minimum");
05994 float thr2 = (thr1-thr3)/2 + thr3;
05995 size_t size = (size_t)nx*ny*nz;
05996 float x0 = thr1,x3 = thr3,x1,x2,THR=0;
05997
05998 #ifdef _WIN32
05999 int ILE = _cpp_min(nx*ny*nx,_cpp_max(1,vol_voxels));
06000 #else
06001 int ILE = std::min(nx*ny*nx,std::max(1,vol_voxels));
06002 #endif //_WIN32
06003
06004 if (abs(thr3-thr2)>abs(thr2-thr1)) {
06005 x1=thr2;
06006 x2=thr2+C*(thr3-thr2);
06007 } else {
06008 x2=thr2;
06009 x1=thr2-C*(thr2-thr1);
06010 }
06011
06012 int cnt1=0,cnt2=0;
06013 for (size_t i=0;i<size;++i) {
06014 if(rdata[i]>=x1) cnt1++;
06015 if(rdata[i]>=x2) cnt2++;
06016 }
06017 float LF1 = static_cast<float>( cnt1 - ILE );
06018 float F1 = LF1*LF1;
06019 float LF2 = static_cast<float>( cnt2 - ILE );
06020 float F2 = LF2*LF2;
06021
06022 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)))
06023 {
06024 if(F2 < F1) {
06025 x0=x1;
06026 x1=x2;
06027 x2 = R*x1 + C*x3;
06028 F1=F2;
06029 int cnt=0;
06030 for(size_t i=0;i<size;++i)
06031 if(rdata[i]>=x2)
06032 cnt++;
06033 LF2 = static_cast<float>( cnt - ILE );
06034 F2 = LF2*LF2;
06035 } else {
06036 x3=x2;
06037 x2=x1;
06038 x1=R*x2 + C*x0;
06039 F2=F1;
06040 int cnt=0;
06041 for(size_t i=0;i<size;++i)
06042 if(rdata[i]>=x1)
06043 cnt++;
06044 LF1 = static_cast<float>( cnt - ILE );
06045 F1 = LF1*LF1;
06046 }
06047 }
06048
06049 if(F1 < F2) {
06050 ILE = static_cast<int> (LF1 + ILE);
06051 THR = x1;
06052 } else {
06053 ILE = static_cast<int> (LF2 + ILE);
06054 THR = x2;
06055 }
06056 return THR;
06057
06058 }
06059 #undef avagadro
06060 #undef density_protein
06061 #undef R
06062 #undef C
06063
06064
06065
06066
06067
06068
06069
06070 vector<float> EMData::peak_ccf(float hf_p)
06071 {
06072
06073
06074
06075 EMData & buf = *this;
06076 vector<Pixel> peaks;
06077 int half=int(hf_p);
06078 float hf_p2 = hf_p*hf_p;
06079 int i,j;
06080 int i__1,i__2;
06081 int j__1,j__2;
06082 vector<float>res;
06083 int nx = buf.get_xsize()-half;
06084 int ny = buf.get_ysize()-half;
06085
06086 for(i=half; i<=nx; ++i) {
06087
06088 i__1 = i-1;
06089 i__2 = i+1;
06090 for (j=half;j<=ny;++j) {
06091 j__1 = j-1;
06092 j__2 = j+1;
06093
06094 if((buf(i,j)>0.0f)&&buf(i,j)>buf(i,j__1)) {
06095 if(buf(i,j)>buf(i,j__2)) {
06096 if(buf(i,j)>buf(i__1,j)) {
06097 if(buf(i,j)>buf(i__2,j)) {
06098 if(buf(i,j)>buf(i__1,j__1)) {
06099 if((buf(i,j))> buf(i__1,j__2)) {
06100 if(buf(i,j)>buf(i__2,j__1)) {
06101 if(buf(i,j)> buf(i__2,j__2)) {
06102
06103
06104
06105 if (peaks.size()==0) {
06106
06107 peaks.push_back(Pixel(i,j,0,buf(i,j)));
06108
06109 } else {
06110
06111
06112 bool overlap = false;
06113
06114
06115
06116
06117
06118 std::stack <vector<Pixel>::iterator> delete_stack;
06119
06120
06121 for (vector<Pixel>::iterator it=peaks.begin();it!=peaks.end();++it) {
06122
06123
06124
06125
06126 float radius=((*it).x-float(i))*((*it).x-float(i))+((*it).y-float(j))*((*it).y-float(j));
06127 if (radius <= hf_p2 ) {
06128
06129 if( buf(i,j) > (*it).value) {
06130
06131
06132
06133
06134
06135
06136 delete_stack.push(it);
06137
06138
06139
06140
06141 } else {
06142 overlap = true;
06143
06144
06145 break;
06146 }
06147 }
06148 }
06149
06150
06151
06152 if (false == overlap) {
06153 vector<Pixel>::iterator delete_iterator;
06154 while (!delete_stack.empty()) {
06155
06156
06157 delete_iterator = delete_stack.top();
06158 peaks.erase(delete_iterator);
06159 delete_stack.pop();
06160 }
06161
06162
06163
06164 if (! (peaks.size() < 2000 )) {
06165
06166
06167
06168
06169 sort(peaks.begin(), peaks.end(), peakcmp);
06170
06171
06172 peaks.pop_back();
06173 }
06174
06175
06176 peaks.push_back(Pixel(i,j,0,buf(i,j)));
06177
06178
06179 } else {
06180
06181
06182 while (!delete_stack.empty()) delete_stack.pop();
06183 }
06184 }
06185 }
06186 }}}}}}}
06187 }
06188 }
06189
06190
06191 if(peaks.size()>0) {
06192
06193 sort(peaks.begin(),peaks.end(), peakcmp);
06194
06195 for (vector<Pixel>::iterator it = peaks.begin(); it != peaks.end(); it++) {
06196
06197 res.push_back((*it).value);
06198 res.push_back(static_cast<float>((*it).x));
06199 res.push_back(static_cast<float>((*it).y));
06200 }
06201 } else {
06202
06203 res.push_back(buf(0,0,0));
06204 res.insert(res.begin(),1,0.0);
06205 }
06206 return res;
06207 }
06208
06209 EMData* EMData::get_pow(float n_pow)
06210 {
06211 EMData* buf_new = this->copy_head();
06212 float *in = this->get_data();
06213 float *out = buf_new->get_data();
06214 for(size_t i=0; i<(size_t)nx*ny*nz; ++i) out[i] = pow(in[i],n_pow);
06215 return buf_new;
06216 }
06217
06218 EMData* EMData::conjg()
06219 {
06220 if(this->is_complex()) {
06221 EMData* buf_new = this->copy_head();
06222 float *in = this->get_data();
06223 float *out = buf_new->get_data();
06224 for(size_t i=0; i<(size_t)nx*ny*nz; i+=2) {out[i] = in[i]; out[i+1] = -in[i+1];}
06225 return buf_new;
06226 } else throw ImageFormatException("image has to be complex");
06227 }
06228
06229 EMData* EMData::delete_disconnected_regions(int ix, int iy, int iz) {
06230 if (3 != get_ndim())
06231 throw ImageDimensionException("delete_disconnected_regions needs a 3-D image.");
06232 if (is_complex())
06233 throw ImageFormatException("delete_disconnected_regions requires a real image");
06234 if ((*this)(ix+nx/2,iy+ny/2,iz+nz/2) == 0)
06235 throw ImageDimensionException("delete_disconnected_regions starting point is zero.");
06236
06237 EMData* result = this->copy_head();
06238 result->to_zero();
06239 (*result)(ix+nx/2,iy+ny/2,iz+nz/2) = (*this)(ix+nx/2,iy+ny/2,iz+nz/2);
06240 bool kpt = true;
06241
06242 while(kpt) {
06243 kpt = false;
06244 for (int cz = 1; cz < nz-1; cz++) {
06245 for (int cy = 1; cy < ny-1; cy++) {
06246 for (int cx = 1; cx < nx-1; cx++) {
06247 if((*result)(cx,cy,cz) == 1) {
06248 for (int lz = -1; lz <= 1; lz++) {
06249 for (int ly = -1; ly <= 1; ly++) {
06250 for (int lx = -1; lx <= 1; lx++) {
06251 if(((*this)(cx+lx,cy+ly,cz+lz) == 1) && ((*result)(cx+lx,cy+ly,cz+lz) == 0)) {
06252 (*result)(cx+lx,cy+ly,cz+lz) = 1;
06253 kpt = true;
06254 }
06255 }
06256 }
06257 }
06258 }
06259 }
06260 }
06261 }
06262 }
06263 result->update();
06264 return result;
06265 }
06266
06267 #define QUADPI 3.141592653589793238462643383279502884197
06268 #define DGR_TO_RAD QUADPI/180
06269
06270 EMData* EMData::helicise(float pixel_size, float dp, float dphi, float section_use, float radius, float minrad) {
06271 if (3 != get_ndim())
06272 throw ImageDimensionException("helicise needs a 3-D image.");
06273 if (is_complex())
06274 throw ImageFormatException("helicise requires a real image");
06275 EMData* result = this->copy_head();
06276 result->to_zero();
06277 int nyc = ny/2;
06278 int nxc = nx/2;
06279 int vl = nz-1;
06280 if ( section_use < dp/int(vl*pixel_size) )
06281 section_use = (dp)/int(vl*pixel_size);
06282
06283 float nb = vl*(1.0f - section_use)/2.0f;
06284
06285 float ne = nb+vl*section_use;
06286 int numst = int( (ne-nb)*pixel_size/dp );
06287
06288
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 int nst1 = int ( (nb-k)*pixel_size/dp) -1;
06296 int nst2 = int ( (ne-k)*pixel_size/dp) +1;
06297 for (int j = 0; j<ny; j++) {
06298 int jy = j - nyc;
06299 int jj = jy*jy;
06300 for (int i = 0; i<nx; i++) {
06301 int ix = i - nxc;
06302 float d2 = (float)(ix*ix + jj);
06303 if(d2 <= r2 && d2>=ir) {
06304 int nq = 0;
06305 for ( int ist = nst1; ist < nst2; ist++) {
06306 float zold = (k*pixel_size + ist*dp)/pixel_size;
06307
06308 if(zold >= nb && zold <= ne) {
06309
06310 float cphi = ist*dphi*(float)DGR_TO_RAD;
06311 float ca = cos(cphi);
06312 float sa = sin(cphi);
06313 float xold = ix*ca - jy*sa + nxc;
06314 float yold = ix*sa + jy*ca + nyc;
06315 nq++;
06316
06317 int IOZ = int(zold);
06318
06319 int IOX = int(xold);
06320 int IOY = int(yold);
06321
06322
06323 #ifdef _WIN32
06324 int IOXp1 = _cpp_min( nx-1 ,IOX+1);
06325 #else
06326 int IOXp1 = std::min( nx-1 ,IOX+1);
06327 #endif //_WIN32
06328
06329 #ifdef _WIN32
06330 int IOYp1 = _cpp_min( ny-1 ,IOY+1);
06331 #else
06332 int IOYp1 = std::min( ny-1 ,IOY+1);
06333 #endif //_WIN32
06334
06335 #ifdef _WIN32
06336 int IOZp1 = _cpp_min( nz-1 ,IOZ+1);
06337 #else
06338 int IOZp1 = std::min( nz-1 ,IOZ+1);
06339 #endif //_WIN32
06340
06341 float dx = xold-IOX;
06342 float dy = yold-IOY;
06343 float dz = zold-IOZ;
06344
06345 float a1 = (*this)(IOX,IOY,IOZ);
06346 float a2 = (*this)(IOXp1,IOY,IOZ) - (*this)(IOX,IOY,IOZ);
06347 float a3 = (*this)(IOX,IOYp1,IOZ) - (*this)(IOX,IOY,IOZ);
06348 float a4 = (*this)(IOX,IOY,IOZp1) - (*this)(IOX,IOY,IOZ);
06349 float a5 = (*this)(IOX,IOY,IOZ) - (*this)(IOXp1,IOY,IOZ) - (*this)(IOX,IOYp1,IOZ) + (*this)(IOXp1,IOYp1,IOZ);
06350 float a6 = (*this)(IOX,IOY,IOZ) - (*this)(IOXp1,IOY,IOZ) - (*this)(IOX,IOY,IOZp1) + (*this)(IOXp1,IOY,IOZp1);
06351 float a7 = (*this)(IOX,IOY,IOZ) - (*this)(IOX,IOYp1,IOZ) - (*this)(IOX,IOY,IOZp1) + (*this)(IOX,IOYp1,IOZp1);
06352 float a8 = (*this)(IOXp1,IOY,IOZ) + (*this)(IOX,IOYp1,IOZ)+ (*this)(IOX,IOY,IOZp1)
06353 - (*this)(IOX,IOY,IOZ)- (*this)(IOXp1,IOYp1,IOZ) - (*this)(IOXp1,IOY,IOZp1)
06354 - (*this)(IOX,IOYp1,IOZp1) + (*this)(IOXp1,IOYp1,IOZp1);
06355
06356
06357
06358 (*result)(i,j,k) += a1 + dz*(a4 + a6*dx + (a7 + a8*dx)*dy) + a3*dy + dx*(a2 + a5*dy);
06359 if(nq == numst) break;
06360 }
06361 }
06362 if(nq != numst)
06363 throw InvalidValueException(nq, "Helicise: incorrect number of repeats encoutered.");
06364 }
06365 }
06366 }
06367 }
06368 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 ;
06369
06370 result->update();
06371 return result;
06372 }
06373
06374
06375
06376 EMData* EMData::helicise_grid(float pixel_size, float dp, float dphi, Util::KaiserBessel& kb, float section_use, float radius, float minrad) {
06377 std::cout<<"111111"<<std::endl;
06378 if (3 != get_ndim())
06379 throw ImageDimensionException("helicise needs a 3-D image.");
06380 if (is_complex())
06381 throw ImageFormatException("helicise requires a real image");
06382
06383
06384 float scale = 0.5f;
06385
06386
06387 int nxn = nx/2; int nyn = ny/2; int nzn = nz/2;
06388
06389 vector<int> saved_offsets = get_array_offsets();
06390 set_array_offsets(0,0,0);
06391 EMData* ret = this->copy_head();
06392 #ifdef _WIN32
06393 ret->set_size(nxn, _cpp_max(nyn,1), _cpp_max(nzn,1));
06394 #else
06395 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
06396 #endif //_WIN32
06397 ret->to_zero();
06398
06399
06400 int xc = nxn;
06401 int ixs = nxn%2;
06402 int yc = nyn;
06403 int iys = nyn%2;
06404 int zc = nzn;
06405 int izs = nzn%2;
06406
06407 int xcn = nxn/2;
06408 int ycn = nyn/2;
06409 int zcn = nzn/2;
06410
06411 float shiftxc = xcn;
06412 float shiftyc = ycn;
06413 float shiftzc = zcn;
06414
06415 float zmin = -nz/2.0f;
06416 float ymin = -ny/2.0f;
06417 float xmin = -nx/2.0f;
06418 float zmax = -zmin;
06419 float ymax = -ymin;
06420 float xmax = -xmin;
06421 if (0 == nx%2) xmax--;
06422 if (0 == ny%2) ymax--;
06423 if (0 == nz%2) zmax--;
06424
06425 float* data = this->get_data();
06426
06427
06428
06429
06430
06431
06432
06433
06434
06435
06436 int nyc = nyn/2;
06437 int nxc = nxn/2;
06438 int nb = int(nzn*(1.0f - section_use)/2.);
06439 int ne = nzn - nb -1;
06440 int numst = int(nzn*section_use*pixel_size/dp);
06441
06442 int nst = int(nzn*pixel_size/dp);
06443 float r2, ir;
06444 if(radius < 0.0f) r2 = (float)((nxc-1)*(nxc-1));
06445 else r2 = radius*radius;
06446 if(minrad < 0.0f) ir = 0.0f;
06447 else ir = minrad*minrad;
06448
06449 for (int k = 0; k<nzn; k++) {
06450 for (int j = 0; j<nyn; j++) {
06451 int jy = j - nyc;
06452 int jj = jy*jy;
06453 for (int i = 0; i<nxn; i++) {
06454 int ix = i - nxc;
06455 float d2 = (float)(ix*ix + jj);
06456 if(d2 <= r2 && d2>=ir) {
06457 int nq = 0;
06458 for ( int ist = -nst; ist <= nst; ist++) {
06459 float zold = (k*pixel_size + ist*dp)/pixel_size;
06460 int IOZ = int(zold);
06461 if(IOZ >= nb && IOZ <= ne) {
06462
06463 float cphi = ist*dphi*(float)DGR_TO_RAD;
06464 float ca = cos(cphi);
06465 float sa = sin(cphi);
06466
06467 float xold = ix*ca - jy*sa + nxc;
06468 float yold = ix*sa + jy*ca + nyc;
06469
06470 float xold_big = (xold-shiftxc)/scale - ixs + xc;
06471 float yold_big = (yold-shiftyc)/scale - iys + yc;
06472 float zold_big = (zold-shiftzc)/scale - izs + zc;
06473
06474
06475
06476
06477
06478
06479
06480
06481
06482
06483
06484
06485
06486
06487
06488
06489
06490
06491
06492
06493 nq++;
06494
06495
06496 (*ret)(i,j,k) += Util::get_pixel_conv_new(nx, ny, nz, xold_big, yold_big, zold_big, data, kb);
06497
06498
06499 if(nq == numst) break;
06500 }
06501 }
06502 if(nq != numst)
06503 throw InvalidValueException(nq, "Helicise: incorrect number of repeats encoutered.");
06504 }
06505 }
06506 }
06507 }
06508
06509 for (int k = 0; k<nzn; k++) for (int j = 0; j<nyn; j++) for (int i = 0; i<nxn; i++) (*ret)(i,j,k) /= numst ;
06510 set_array_offsets(saved_offsets);
06511 ret->update();
06512 return ret;
06513 }
06514
06515
06516
06517
06518
06519
06520
06521
06522
06523 void EMData::depad() {
06524 if (is_complex())
06525 throw ImageFormatException("Depadding of complex images not supported");
06526 vector<int> saved_offsets = get_array_offsets();
06527 set_array_offsets(0,0,0);
06528 int npad = attr_dict["npad"];
06529 if (0 == npad) npad = 1;
06530 int offset = is_fftodd() ? 1 : 2;
06531 int nxold = (nx - offset)/npad;
06532 #ifdef _WIN32
06533 int nyold = _cpp_max(ny/npad, 1);
06534 int nzold = _cpp_max(nz/npad, 1);
06535 #else
06536 int nyold = std::max<int>(ny/npad, 1);
06537 int nzold = std::max<int>(nz/npad, 1);
06538 #endif //_WIN32
06539 int xstart = 0, ystart = 0, zstart = 0;
06540 if( npad > 1) {
06541 xstart = (nx - offset - nxold)/2 + nxold%2;
06542 if(ny > 1) {
06543 ystart = (ny - nyold)/2 + nyold%2;
06544 if(nz > 1) {
06545 zstart = (nz - nzold)/2 + nzold%2;
06546 }
06547 }
06548 }
06549 int bytes = nxold*sizeof(float);
06550 float* dest = get_data();
06551 for (int iz=0; iz < nzold; iz++) {
06552 for (int iy = 0; iy < nyold; iy++) {
06553 memmove(dest, &(*this)(xstart,iy+ystart,iz+zstart), bytes);
06554 dest += nxold;
06555 }
06556 }
06557 set_size(nxold, nyold, nzold);
06558 set_attr("npad", 1);
06559 set_fftpad(false);
06560 set_fftodd(false);
06561 set_complex(false);
06562 if(ny==1 && nz==1) set_complex_x(false);
06563 set_array_offsets(saved_offsets);
06564 update();
06565 EXITFUNC;
06566 }
06567
06568
06569
06570
06571
06572
06573
06574
06575 void EMData::depad_corner() {
06576 if(is_complex())
06577 throw ImageFormatException("Depadding of complex images not allowed");
06578 vector<int> saved_offsets = get_array_offsets();
06579 set_array_offsets(0,0,0);
06580 int npad = attr_dict["npad"];
06581 if(0 == npad) npad = 1;
06582 int offset = is_fftodd() ? 1 : 2;
06583 int nxold = (nx - offset)/npad;
06584 #ifdef _WIN32
06585 int nyold = _cpp_max(ny/npad, 1);
06586 int nzold = _cpp_max(nz/npad, 1);
06587 #else
06588 int nyold = std::max<int>(ny/npad, 1);
06589 int nzold = std::max<int>(nz/npad, 1);
06590 #endif //_WIN32
06591 size_t bytes = nxold*sizeof(float);
06592 float* dest = get_data();
06593 for (int iz=0; iz < nzold; iz++) {
06594 for (int iy = 0; iy < nyold; iy++) {
06595 memmove(dest, &(*this)(0,iy,iz), bytes);
06596 dest += nxold;
06597 }
06598 }
06599 set_size(nxold, nyold, nzold);
06600 set_attr("npad", 1);
06601 set_fftpad(false);
06602 set_fftodd(false);
06603 set_complex(false);
06604 if(ny==1 && nz==1) set_complex_x(false);
06605 set_array_offsets(saved_offsets);
06606 update();
06607 EXITFUNC;
06608 }
06609
06610
06611
06612 float circumference( EMData* emdata, int npixel )
06613 {
06614 int nx = emdata->get_xsize();
06615 int ny = emdata->get_ysize();
06616 int nz = emdata->get_zsize();
06617
06618 float* data = emdata->get_data();
06619 if( ny==1 && nz==1 ) {
06620
06621 float sumf=0.0;
06622 int sumn=0;
06623 for( int i=0; i < npixel; ++i ) {
06624 sumf += data[i];
06625 sumf += data[nx-1-i];
06626 sumn += 2;
06627 }
06628 return sumf/sumn;
06629 }
06630
06631 if( nz==1 ) {
06632 float sumf=0.0;
06633 int sumn=0;
06634 int id=0;
06635 for( int iy=0; iy < ny; ++iy ) {
06636 for( int ix=0; ix < nx; ++ix ) {
06637 if( iy<npixel || iy>ny-1-npixel || ix<npixel || ix>nx-1-npixel ) {
06638 sumf += data[id];
06639 sumn += 1;
06640 }
06641 id++;
06642 }
06643 }
06644
06645 Assert( id==nx*ny );
06646 Assert( sumn == nx*ny - (nx-2*npixel)*(ny-2*npixel) );
06647 return sumf/sumn;
06648 }
06649
06650
06651
06652 float sumf = 0.0;
06653 size_t sumn = 0;
06654 size_t id = 0;
06655 for( int iz=0; iz < nz; ++iz) {
06656 for( int iy=0; iy < ny; ++iy) {
06657 for( int ix=0; ix < nx; ++ix ) {
06658 if( iz<npixel||iz>nz-1-npixel||iy<npixel||iy>ny-1-npixel||ix<npixel||ix>nx-1-npixel) {
06659 sumf += data[id];
06660 sumn += 1;
06661 }
06662 id++;
06663 }
06664 }
06665 }
06666
06667
06668 Assert( id==(size_t)nx*ny*nz);
06669 Assert( sumn==(size_t)nx*ny*nz-(size_t)(nx-2*npixel)*(ny-2*npixel)*(nz-2*npixel) );
06670 return sumf/sumn;
06671 }
06672
06673
06674
06675
06676
06677
06678
06679
06680 EMData* EMData::norm_pad(bool donorm, int npad, int valtype) {
06681 if (this->is_complex())
06682 throw ImageFormatException("Padding of complex images not supported");
06683 int nx = this->get_xsize();
06684 int ny = this->get_ysize();
06685 int nz = this->get_zsize();
06686 float mean = 0., stddev = 1.;
06687 if(donorm) {
06688 mean = this->get_attr("mean");
06689 stddev = this->get_attr("sigma");
06690 }
06691
06692 if (npad < 1) npad = 1;
06693 int nxpad = npad*nx;
06694 int nypad = npad*ny;
06695 int nzpad = npad*nz;
06696 if (1 == ny) {
06697
06698
06699 nypad = ny;
06700 nzpad = nz;
06701 } else if (nz == 1) {
06702
06703 nzpad = nz;
06704 }
06705 size_t bytes;
06706 size_t offset;
06707
06708 offset = 2 - nxpad%2;
06709 bytes = nx*sizeof(float);
06710 EMData* fpimage = copy_head();
06711 fpimage->set_size(nxpad+offset, nypad, nzpad);
06712 int xstart = 0, ystart = 0, zstart = 0;
06713 if( npad > 1) {
06714 if( valtype==0 ) {
06715 fpimage->to_zero();
06716 } else {
06717 float val = circumference(this, 1);
06718 float* data = fpimage->get_data();
06719 int nxyz = (nxpad+offset)*nypad*nzpad;
06720 for( int i=0; i < nxyz; ++i ) data[i] = val;
06721 }
06722
06723 xstart = (nxpad - nx)/2 + nx%2;
06724 if(ny > 1) {
06725 ystart = (nypad - ny)/2 + ny%2;
06726 if(nz > 1) {
06727 zstart = (nzpad - nz)/2 + nz%2;
06728 }
06729 }
06730 }
06731
06732
06733 vector<int> saved_offsets = this->get_array_offsets();
06734 this->set_array_offsets( 0, 0, 0 );
06735 for (int iz = 0; iz < nz; iz++) {
06736 for (int iy = 0; iy < ny; iy++) {
06737 memcpy(&(*fpimage)(xstart,iy+ystart,iz+zstart), &(*this)(0,iy,iz), bytes);
06738 }
06739 }
06740 this->set_array_offsets( saved_offsets );
06741
06742
06743
06744
06745 if (donorm) {
06746 for (int iz = zstart; iz < nz+zstart; iz++)
06747 for (int iy = ystart; iy < ny+ystart; iy++)
06748 for (int ix = xstart; ix < nx+xstart; ix++)
06749 (*fpimage)(ix,iy,iz) = ((*fpimage)(ix,iy,iz)-mean)/stddev;
06750 }
06751
06752 fpimage->set_fftpad(true);
06753 fpimage->set_attr("npad", npad);
06754 if (offset == 1) fpimage->set_fftodd(true);
06755 else fpimage->set_fftodd(false);
06756 return fpimage;
06757 }
06758
06759 void EMData::center_origin()
06760 {
06761 ENTERFUNC;
06762 if (is_complex()) {
06763 LOGERR("Real image expected. Input image is complex.");
06764 throw ImageFormatException("Real image expected. Input image is complex.");
06765 }
06766 for (int iz = 0; iz < nz; iz++) {
06767 for (int iy = 0; iy < ny; iy++) {
06768 for (int ix = 0; ix < nx; ix++) {
06769
06770 (*this)(ix,iy,iz) *= -2*((ix+iy+iz)%2) + 1;
06771 }
06772 }
06773 }
06774 update();
06775 EXITFUNC;
06776 }
06777
06778 void EMData::center_origin_yz()
06779 {
06780 ENTERFUNC;
06781 if (is_complex()) {
06782 LOGERR("Real image expected. Input image is complex.");
06783 throw ImageFormatException("Real image expected. Input image is complex.");
06784 }
06785 for (int iz = 0; iz < nz; iz++) {
06786 for (int iy = (iz+1)%2; iy < ny; iy+=2) {
06787 for (int ix = 0; ix < nx; ix++) {
06788 (*this)(ix,iy,iz) *= -1;
06789 }
06790 }
06791 }
06792 update();
06793 EXITFUNC;
06794 }
06795
06796 void EMData::center_origin_fft()
06797 {
06798 ENTERFUNC;
06799 if (!is_complex()) {
06800 LOGERR("complex image expected. Input image is real image.");
06801 throw ImageFormatException("complex image expected. Input image is real image.");
06802 }
06803
06804 if (!is_ri()) {
06805 LOGWARN("Only RI should be used. ");
06806 }
06807 vector<int> saved_offsets = get_array_offsets();
06808
06809
06810
06811 set_array_offsets(0,1,1);
06812 int nxc = nx/2;
06813
06814 if (is_fftodd()) {
06815 for (int iz = 1; iz <= nz; iz++) {
06816 for (int iy = 1; iy <= ny; iy++) {
06817 for (int ix = 0; ix < nxc; ix++) {
06818 cmplx(ix,iy,iz) *= float(-2*((ix+iy+iz)%2) + 1);
06819 float temp = float(iz-1+iy-1+ix)/float(ny)*M_PI;
06820 complex<float> temp2 = complex<float>(cos(temp), -sin(temp));
06821 cmplx(ix,iy,iz) *= temp2;
06822 }
06823 }
06824 }
06825 } else {
06826 for (int iz = 1; iz <= nz; iz++) {
06827 for (int iy = 1; iy <= ny; iy++) {
06828 for (int ix = 0; ix < nxc; ix++) {
06829
06830 cmplx(ix,iy,iz) *= float(-2*((ix+iy+iz)%2) + 1);
06831 }
06832 }
06833 }
06834 }
06835 set_array_offsets(saved_offsets);
06836 update();
06837 EXITFUNC;
06838 }
06839
06840
06841 #define fint(i,j,k) fint[(i-1) + ((j-1) + (k-1)*ny)*(size_t)lsd]
06842 #define fout(i,j,k) fout[(i-1) + ((j-1) + (k-1)*nyn)*(size_t)lsdn]
06843 EMData *EMData::FourInterpol(int nxn, int nyni, int nzni, bool RetReal) {
06844
06845 int nyn, nzn, lsd, lsdn, inx, iny, inz;
06846 int i, j, k;
06847 if (is_complex())
06848 throw ImageFormatException("Input image has to be real");
06849
06850 if(ny > 1) {
06851 nyn = nyni;
06852 if(nz > 1) {
06853 nzn = nzni;
06854 } else {
06855 nzn = 1;
06856 }
06857 } else {
06858 nyn = 1; nzn = 1;
06859 }
06860 if(nxn<nx || nyn<ny || nzn<nz) throw ImageDimensionException("Cannot reduce the image size");
06861 lsd = nx + 2 - nx%2;
06862 lsdn = nxn + 2 - nxn%2;
06863
06864 EMData *temp_ft = do_fft();
06865 EMData *ret = this->copy();
06866 ret->set_size(lsdn, nyn, nzn);
06867 ret->to_zero();
06868 float *fout = ret->get_data();
06869 float *fint = temp_ft->get_data();
06870
06871
06872 float sq2 = 1.0f/std::sqrt(2.0f);
06873 float anorm = (float) nxn* (float) nyn* (float) nzn/(float) nx/ (float) ny/ (float) nz;
06874 for (i = 0; i < lsd*ny*nz; i++) fint[i] *= anorm;
06875 inx = nxn-nx; iny = nyn - ny; inz = nzn - nz;
06876 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);
06877 if(nyn>1) {
06878
06879 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);
06880 if(nzn>1) {
06881 for (k=nz/2+2+inz; k<=nzn; k++) {
06882 for (j=1; j<=ny/2+1; j++) {
06883 for (i=1; i<=lsd; i++) {
06884 fout(i,j,k)=fint(i,j,k-inz);
06885 }
06886 }
06887 for (j=ny/2+2+iny; j<=nyn; j++) {
06888 for (i=1; i<=lsd; i++) {
06889 fout(i,j,k)=fint(i,j-iny,k-inz);
06890 }
06891 }
06892 }
06893 }
06894 }
06895
06896
06897
06898 if(nx%2 == 0 && inx !=0) {
06899 for (k=1; k<=nzn; k++) {
06900 for (j=1; j<=nyn; j++) {
06901 fout(nx+1,j,k) *= sq2;
06902 fout(nx+2,j,k) *= sq2;
06903 }
06904 }
06905 if(nyn>1) {
06906 for (k=1; k<=nzn; k++) {
06907 for (i=1; i<=lsd; i++) {
06908 fout(i,ny/2+1+iny,k) = sq2*fout(i,ny/2+1,k);
06909 fout(i,ny/2+1,k) *= sq2;
06910 }
06911 }
06912 if(nzn>1) {
06913 for (j=1; j<=nyn; j++) {
06914 for (i=1; i<=lsd; i++) {
06915 fout(i,j,nz/2+1+inz) = sq2*fout(i,j,nz/2+1);
06916 fout(i,j,nz/2+1) *= sq2;
06917 }
06918 }
06919 }
06920 }
06921 }
06922 ret->set_complex(true);
06923
06924
06925
06926
06927
06928
06929
06930
06931
06932
06933
06934
06935
06936
06937
06938
06939
06940
06941
06942
06943
06944
06945
06946
06947
06948
06949
06950
06951
06952
06953
06954
06955
06956 ret->set_ri(1);
06957 ret->set_fftpad(true);
06958 ret->set_attr("npad", 1);
06959 if (nxn%2 == 1) {ret->set_fftodd(true);} else {ret->set_fftodd(false);}
06960 if(RetReal) {
06961 ret->do_ift_inplace();
06962 ret->depad();
06963 }
06964 ret->update();
06965
06966
06967
06968
06969
06970
06971
06972 delete temp_ft;
06973 temp_ft = 0;
06974 return ret;
06975 }
06976
06977 EMData *EMData::FourTruncate(int nxn, int nyni, int nzni, bool RetReal) {
06978
06979 int nyn, nzn, lsd, lsdn, inx, iny, inz;
06980 int i, j, k;
06981 float *fint;
06982 EMData *temp_ft = NULL;
06983
06984
06985
06986 if(ny > 1) {
06987 nyn = nyni;
06988 if(nz > 1) {
06989 nzn = nzni;
06990 } else {
06991 nzn = 1;
06992 }
06993 } else {
06994 nyn = 1; nzn = 1;
06995 }
06996 if (is_complex()) {
06997 nx = nx - 2 + nx%2;
06998 fint = get_data();
06999 } else {
07000
07001 temp_ft = do_fft();
07002 fint = temp_ft->get_data();
07003 }
07004 if(nxn>nx || nyn>ny || nzn>nz) throw ImageDimensionException("Cannot increase the image size");
07005 lsd = nx + 2 - nx%2;
07006 lsdn = nxn + 2 - nxn%2;
07007 EMData *ret = this->copy_head();
07008 ret->set_size(lsdn, nyn, nzn);
07009 float *fout = ret->get_data();
07010
07011
07012
07013 float anorm = (float) nxn* (float) nyn* (float) nzn/(float) nx/ (float) ny/ (float) nz;
07014 for (i = 0; i < lsd*ny*nz; i++) fint[i] *= anorm;
07015 inx = nx - nxn; iny = ny - nyn; inz = nz - nzn;
07016 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);
07017 if(nyn>1) {
07018 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);
07019 if(nzn>1) {
07020 for (k=nzn/2+2; k<=nzn; k++) {
07021 for (j=1; j<=nyn/2+1; j++) {
07022 for (i=1; i<=lsdn; i++) {
07023 fout(i,j,k)=fint(i,j,k+inz);
07024 }
07025 }
07026 for (j=nyn/2+2; j<=nyn; j++) {
07027 for (i=1; i<=lsdn; i++) {
07028 fout(i,j,k)=fint(i,j+iny,k+inz);
07029 }
07030 }
07031 }
07032 }
07033 }
07034
07035
07036
07037
07038
07039
07040
07041
07042
07043
07044
07045
07046
07047
07048
07049
07050
07051
07052
07053
07054
07055
07056
07057
07058
07059
07060
07061
07062 ret->set_complex(true);
07063 ret->set_ri(1);
07064 ret->set_fftpad(true);
07065 ret->set_attr("npad", 1);
07066 if (nxn%2 == 1) {ret->set_fftodd(true);} else {ret->set_fftodd(false);}
07067 if(RetReal) {
07068 ret->do_ift_inplace();
07069 ret->depad();
07070 }
07071 ret->update();
07072
07073
07074
07075
07076
07077
07078
07079 if (!is_complex()) {
07080 delete temp_ft;
07081 temp_ft = 0;
07082 }
07083 return ret;
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
07136
07137
07138
07139
07140
07141
07142
07143
07144
07145
07146
07147
07148
07149
07150
07151
07152
07153
07154
07155
07156
07157
07158
07159
07160
07161
07162
07163
07164
07165
07166
07167
07168
07169
07170
07171
07172
07173
07174
07175
07176
07177
07178
07179
07180 EMData *EMData::Four_ds(int nxn, int nyni, int nzni, bool RetReal) {
07181
07182 int nyn, nzn, lsd, lsdn, inx, iny, inz;
07183 int i, j;
07184
07185 if(ny > 1) {
07186 nyn = nyni;
07187 if(nz > 1) {
07188 nzn = nzni;
07189 } else {
07190 nzn = 1;
07191 }
07192 } else {
07193 nyn = 1; nzn = 1;
07194 }
07195 lsd = nx-2 + 2 - nx%2;
07196 lsdn = nxn + 2 - nxn%2;
07197
07198 EMData *temp_ft = this->copy();
07199 EMData *ret = this->copy();
07200 ret->set_size(lsdn, nyn, nzn);
07201 ret->to_zero();
07202 float *fout = ret->get_data();
07203 float *fint = temp_ft->get_data();
07204
07205
07206
07207 float anorm = (float) nxn* (float) nyn* (float) nzn/(float) nx/ (float) ny/ (float) nz;
07208 for (i = 0; i < lsd*ny*nz; i++) fint[i] *= anorm;
07209 inx = nxn-(nx-2); iny = nyn - ny; inz = nzn - nz;
07210 for (j=1; j<=nyn; j++)
07211 for (i=1; i<=lsdn; i++)
07212 fout(i,j,1)=fint((i-1)/2*4+2-i%2,j*2-1,1);
07213 ret->set_complex(true);
07214 ret->set_ri(1);
07215
07216
07217 if (nxn%2 == 1) {ret->set_fftodd(true);} else {ret->set_fftodd(false);}
07218 if(RetReal) {
07219 ret->do_ift_inplace();
07220 ret->depad();
07221 }
07222 ret->update();
07223
07224 delete temp_ft;
07225 temp_ft = 0;
07226 return ret;
07227 }
07228
07229 EMData *EMData::Four_shuf_ds_cen_us(int nxn, int nyni, int, bool RetReal) {
07230
07231 int nyn, nzn, lsd, lsdn, inx, iny, inz;
07232 int i, j;
07233
07234 nyn = nyni;
07235 nzn = 1;
07236 lsd = nx;
07237 lsdn = nxn + 2 - nxn%2;
07238
07239 EMData *temp_ft = this->copy();
07240 EMData *ret = this->copy();
07241 ret->set_size(lsdn, nyn, nzn);
07242 ret->to_zero();
07243 float *fout = ret->get_data();
07244 float *fint = temp_ft->get_data();
07245
07246
07247 float sq2 = 1.0f/std::sqrt(2.0f);
07248
07249 for (size_t i = 0; i < (size_t)lsd*ny*nz; i++) fint[i] *= 4;
07250
07251 inx = nxn-(nx-2); iny = nyn - ny; inz = nzn - nz;
07252 for (j=1; j<=ny/4; j++)
07253 for (i=1; i<=(nx-2)/2+2; i++) {
07254 int g = (i-1)/2+1;
07255 if ((g+j)%2 == 0) {
07256 fout(i,j,1)=fint(g*4-2-i%2,j*2-1+ny/2,1);
07257 } else {
07258 fout(i,j,1)=-fint(g*4-2-i%2,j*2-1+ny/2,1);
07259 }
07260 }
07261
07262 for (j=ny/4+1; j<=ny/4+1; j++)
07263 for (i=1; i<=(nx-2)/2+2; i++) {
07264 int g = (i-1)/2+1;
07265 if ((g+j)%2 == 0) {
07266 fout(i,j,1)=fint(g*4-2-i%2,j*2-1-ny/2,1);
07267 } else {
07268 fout(i,j,1)=-fint(g*4-2-i%2,j*2-1-ny/2,1);
07269 }
07270 }
07271
07272 for (j=ny/4+2; j<=ny/2; j++)
07273 for (i=1; i<=(nx-2)/2+2; i++) {
07274 int g = (i-1)/2+1;
07275 if ((g+j)%2 == 0) {
07276 fout(i,j+ny/2,1)=fint(g*4-2-i%2,j*2-1-ny/2,1);
07277 } else {
07278 fout(i,j+ny/2,1)=-fint(g*4-2-i%2,j*2-1-ny/2,1);
07279 }
07280 }
07281
07282 if (nx%2 == 0) {
07283 for (j=1; j<=nyn; j++) {
07284 fout((nx-2)/2+1,j,1) *= sq2;
07285 fout((nx-2)/2+2,j,1) *= sq2;
07286 }
07287 for (i=1; i<=lsd/2+1; i++) {
07288 fout(i,ny/4+1+ny/2,1) = sq2*fout(i,ny/4+1,1);
07289 fout(i,ny/4+1,1) *= sq2;
07290 }
07291 }
07292
07293 ret->set_complex(true);
07294 ret->set_ri(1);
07295
07296 if (nxn%2 == 1) {ret->set_fftodd(true);} else {ret->set_fftodd(false);}
07297 if(RetReal) {
07298 ret->do_ift_inplace();
07299 ret->depad();
07300 }
07301 ret->update();
07302
07303 delete temp_ft;
07304 temp_ft = 0;
07305 return ret;
07306 }
07307
07308 #undef fint
07309 #undef fout
07310
07311 #define fint(jx,jy,jz) fint[jx + (jy + jz*ny)*(size_t)nox]
07312 EMData *EMData::filter_by_image(EMData* image, bool RetReal) {
07313
07314
07315 bool complex_input = this->is_complex();
07316 nx = this->get_xsize();
07317 ny = this->get_ysize();
07318 nz = this->get_zsize();
07319 int nox;
07320 if (complex_input) nox = (nx - 2 + this->is_fftodd()); else nox = nx;
07321
07322 int lsd2 = (nox + 2 - nox%2) / 2;
07323
07324 EMData* fp = NULL;
07325 if(complex_input) {
07326
07327 fp = this->copy();
07328 } else {
07329 fp = this->norm_pad( false, 1);
07330 fp->do_fft_inplace();
07331 }
07332 fp->set_array_offsets(1,1,1);
07333 int nx2 = nox/2;
07334 int ny2 = ny/2;
07335 int nz2 = nz/2;
07336 float *fint = image->get_data();
07337 for ( int iz = 1; iz <= nz; iz++) {
07338 int jz=nz2-iz+1; if(jz<0) jz += nz;
07339 for ( int iy = 1; iy <= ny; iy++) {
07340 int jy=ny2-iy+1; if(jy<0) jy += ny;
07341 for ( int ix = 1; ix <= lsd2; ix++) {
07342 int jx = nx2-ix+1;
07343 fp->cmplx(ix,iy,iz) *= fint(jx,jy,jz);
07344 }
07345 }
07346 }
07347
07348 fp->set_ri(1);
07349 fp->set_fftpad(true);
07350 fp->set_attr("npad", 1);
07351 if (nx%2 == 1) fp->set_fftodd(true);
07352 else fp->set_fftodd(false);
07353 if(RetReal) {
07354 fp->do_ift_inplace();
07355 fp->depad();
07356 }
07357 fp->set_array_offsets(0,0,0);
07358 fp->update();
07359
07360 return fp;
07361 }
07362 #undef fint
07363 #define fint(jx,jy,jz) fint[jx + (jy + jz*ny)*(size_t)nx]
07364 #define fout(jx,jy,jz) fout[jx + (jy + jz*ny)*(size_t)nx]
07365 EMData *EMData::replace_amplitudes(EMData* image, bool RetReal) {
07366
07367
07368 bool complex_input = this->is_complex();
07369 nx = this->get_xsize();
07370 ny = this->get_ysize();
07371 nz = this->get_zsize();
07372 int nox;
07373 if (complex_input) nox = (nx - 2 + this->is_fftodd()); else nox = nx;
07374
07375 EMData* fp = NULL;
07376 if(complex_input) {
07377
07378 fp = this->copy();
07379 } else {
07380 fp = this->norm_pad( false, 1);
07381 fp->do_fft_inplace();
07382 }
07383 float *fout = fp->get_data();
07384 float *fint = image->get_data();
07385 for ( int iz = 0; iz < nz; iz++) {
07386 for ( int iy = 0; iy < ny; iy++) {
07387 for ( int ix = 0; ix < nx; ix+=2) {
07388 float qt = fint(ix,iy,iz)*fint(ix,iy,iz)+fint(ix+1,iy,iz)*fint(ix+1,iy,iz);
07389 float rt = fout(ix,iy,iz)*fout(ix,iy,iz)+fout(ix+1,iy,iz)*fout(ix+1,iy,iz);
07390 if(rt > 1.0e-20) {
07391 fout(ix,iy,iz) *= (qt/rt);
07392 fout(ix+1,iy,iz) *= (qt/rt);
07393 } else {
07394 qt = std::sqrt(qt/2.0f);
07395 fout(ix,iy,iz) = qt;
07396 fout(ix+1,iy,iz) = qt;
07397 }
07398 }
07399 }
07400 }
07401
07402 fp->set_ri(1);
07403 fp->set_fftpad(true);
07404 fp->set_attr("npad", 1);
07405 if (nx%2 == 1) fp->set_fftodd(true);
07406 else fp->set_fftodd(false);
07407 if(RetReal) {
07408 fp->do_ift_inplace();
07409 fp->depad();
07410 }
07411 fp->set_array_offsets(0,0,0);
07412 fp->update();
07413
07414 return fp;
07415 }
07416 #undef fint
07417 #undef fout
07418
07419
07420 #undef QUADPI
07421 #undef DGR_TO_RAD