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