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