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
00223 EMData *EMData::FH2F(int Size, float OverSamplekB, int IntensityFlag)
00224 {
00225 int nx=get_xsize();
00226 int ny=get_ysize();
00227 int nz=get_zsize();
00228 float ScalFactor=4.1f;
00229 int Center = (int) floor((Size+1.0)/2.0 +.1);
00230 int CenterM= Center-1;
00231 int CountMax = (Center+1)*Center/2;
00232
00233 int *PermMatTr = new int[CountMax];
00234 float *RValsSorted = new float[CountMax];
00235 float *weightofkValsSorted = new float[CountMax];
00236 int *SizeReturned = new int[1];
00237 Util::Radialize(PermMatTr, RValsSorted,weightofkValsSorted,Size, SizeReturned);
00238 int RIntMax= SizeReturned[0];
00239
00240
00241 int mMax = (int) floor( ScalFactor*RValsSorted[RIntMax-1]+10.0);
00242
00243 int kIntMax = 2+ (int) floor( RValsSorted[RIntMax-1]*OverSamplekB);
00244 float *kVec2Use = new float[kIntMax];
00245 for (int kk=0; kk<kIntMax; kk++){
00246 kVec2Use[kk]= ((float) kk)/OverSamplekB;}
00247
00248
00249
00250 #ifdef DEBUG
00251 printf("nx=%d, ny=%d, nz=%d Center=%d mMax=%d CountMax=%d kIntMax=%d Centerm1=%d Size=%d\n\n",
00252 nx,ny,nz, Center, mMax, CountMax, kIntMax, CenterM, Size);
00253 #endif
00254
00255 EMData * rhoOfkmB = this;
00256
00257
00258
00259
00260 if ( (nx==2*(mMax+1)) && (ny==kIntMax) &&(nz==1) ) {
00261
00262 EMData *rhoOfkandm = copy();
00263 rhoOfkandm ->set_size(2*(mMax+1),RIntMax);
00264 rhoOfkandm ->to_zero();
00265
00266
00267 for (int mr=0; mr <2*(mMax+1); mr++){
00268 float *Row= new float[kIntMax];
00269 float *RowOut= new float[RIntMax];
00270 for (int ii=0; ii<kIntMax; ii++){ Row[ii]=(*rhoOfkmB)(mr,ii);}
00271 Util::spline_mat(kVec2Use, Row, kIntMax, RValsSorted, RowOut, RIntMax );
00272 for (int ii=0; ii<RIntMax; ii++){
00273 (*rhoOfkandm)(mr,ii) = RowOut[ii];
00274
00275 }
00276
00277
00278 }
00279 rhoOfkandm ->update();
00280
00281
00282
00283 EMData* outCopy = rhoOfkandm ->copy();
00284 outCopy->set_size(2*Size,Size,1);
00285 outCopy->to_zero();
00286
00287
00288 int Count =0, kInt, kIntm1;
00289 std::complex <float> ImfTemp;
00290 float kValue, thetak;
00291
00292 for (int jkx=0; jkx <Center; jkx++) {
00293 for (int jky=0; jky<=jkx; jky++){
00294 kInt = PermMatTr[Count];
00295 kIntm1= kInt-1;
00296 Count++;
00297 float fjkx = float(jkx);
00298 float fjky = float(jky);
00299
00300 kValue = std::sqrt(fjkx*fjkx + fjky*fjky ) ;
00301
00302
00303
00304
00305 thetak = atan2(fjky,fjkx);
00306 ImfTemp = (*rhoOfkandm)(0, kIntm1) ;
00307 for (int mm= 1; mm <mMax;mm++) {
00308 std::complex <float> fact(0,-mm*thetak);
00309 std::complex <float> expfact= exp(fact);
00310 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1),(*rhoOfkandm)(2*mm+1,kIntm1));
00311 float mmFac = float(1-2*(mm%2));
00312 if (IntensityFlag==1){ mmFac=1;}
00313 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00314 }
00315 (*outCopy)(2*(CenterM+jkx),CenterM+jky) = ImfTemp.real();
00316 (*outCopy)(2*(CenterM+jkx)+1,CenterM+jky) = ImfTemp.imag();
00317
00318
00319 if (jky>0) {
00320 thetak = atan2(-fjky,fjkx);
00321 ImfTemp = (*rhoOfkandm)(0,kIntm1);
00322 for (int mm= 1; mm<mMax; mm++) {
00323 std::complex <float> fact(0,-mm*thetak);
00324 std::complex <float> expfact= exp(fact);
00325 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1), (*rhoOfkandm)(2*mm+1,kIntm1));
00326 float mmFac = float(1-2*(mm%2));
00327 if (IntensityFlag==1){ mmFac=1;}
00328 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00329 }
00330 (*outCopy)(2*(CenterM+jkx),CenterM-jky) = ImfTemp.real();
00331
00332 (*outCopy)(2*(CenterM+jkx)+1,CenterM-jky) = ImfTemp.imag();
00333 }
00334
00335 if (jkx>0) {
00336 thetak = atan2(fjky,-fjkx);
00337 ImfTemp = (*rhoOfkandm)(0,kIntm1);
00338 for (int mm= 1; mm<mMax; mm++) {
00339 std::complex <float> fact(0,-mm*thetak);
00340 std::complex <float> expfact= exp(fact);
00341 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1), (*rhoOfkandm)(2*mm+1,kIntm1));
00342 float mmFac = float(1-2*(mm%2));
00343 if (IntensityFlag==1){ mmFac=1;}
00344 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00345 }
00346 (*outCopy)(2*(CenterM-jkx) ,CenterM+jky) = ImfTemp.real();
00347 (*outCopy)(2*(CenterM-jkx)+1,CenterM+jky) = ImfTemp.imag();
00348 }
00349
00350 if (jkx>0 && jky>0) {
00351 thetak = atan2(-fjky,-fjkx);
00352 ImfTemp = (*rhoOfkandm)(0 , kIntm1);
00353 for (int mm= 1; mm<mMax; mm++) {
00354 std::complex <float> fact(0,-mm*thetak);
00355 std::complex <float> expfact= exp(fact);
00356 std::complex <float> tempRho( (*rhoOfkandm)(2*mm,kIntm1),(*rhoOfkandm)(2*mm+1,kIntm1) );
00357 float mmFac = float(1-2*(mm%2));
00358 if (IntensityFlag==1){ mmFac=1;}
00359 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00360 }
00361 (*outCopy)(2*(CenterM-jkx) ,CenterM-jky) = ImfTemp.real();
00362 (*outCopy)(2*(CenterM-jkx)+1,CenterM-jky) = ImfTemp.imag();
00363 }
00364
00365 if (jky< jkx) {
00366 thetak = atan2(fjkx,fjky);
00367 ImfTemp = (*rhoOfkandm)(0,kIntm1);
00368 for (int mm= 1; mm<mMax; mm++){
00369 std::complex <float> fact(0,-mm*thetak);
00370 std::complex <float> expfact= exp(fact);
00371 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1),(*rhoOfkandm)(2*mm+1,kIntm1));
00372 float mmFac = float(1-2*(mm%2));
00373 if (IntensityFlag==1){ mmFac=1;}
00374 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00375 }
00376 (*outCopy)(2*(CenterM+jky) ,CenterM+jkx) = ImfTemp.real();
00377 (*outCopy)(2*(CenterM+jky)+1,CenterM+jkx) = ImfTemp.imag();
00378
00379 if (jky>0){
00380 thetak = atan2(fjkx,-fjky);
00381 ImfTemp = (*rhoOfkandm)(0, kIntm1);
00382 for (int mm= 1; mm <mMax; mm++) {
00383 std::complex <float> fact(0,-mm*thetak);
00384 std::complex <float> expfact= exp(fact);
00385 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1),(*rhoOfkandm)(2*mm+1,kIntm1));
00386 float mmFac = float(1-2*(mm%2));
00387 if (IntensityFlag==1){ mmFac=1;}
00388 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00389 }
00390 (*outCopy)(2*(CenterM-jky) ,CenterM+jkx) = ImfTemp.real();
00391 (*outCopy)(2*(CenterM-jky)+1,CenterM+jkx) = ImfTemp.imag();
00392 }
00393
00394 if (jkx>0) {
00395 thetak = atan2(-fjkx,fjky);
00396 ImfTemp = (*rhoOfkandm)(0,kIntm1);
00397 for (int mm= 1; mm <mMax; mm++) {
00398 std::complex <float> fact(0,-mm*thetak);
00399 std::complex <float> expfact= exp(fact);
00400 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1),(*rhoOfkandm)(2*mm+1,kIntm1));
00401 float mmFac = float(1-2*(mm%2));
00402 if (IntensityFlag==1){ mmFac=1;}
00403 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00404 }
00405 (*outCopy)(2*(CenterM+jky) ,CenterM-jkx) = ImfTemp.real();
00406 (*outCopy)(2*(CenterM+jky)+1,CenterM-jkx) = ImfTemp.imag();
00407 }
00408
00409 if (jkx>0 && jky>0) {
00410 thetak = atan2(-fjkx,-fjky);
00411 ImfTemp = (*rhoOfkandm)(0,kIntm1) ;
00412 for (int mm= 1; mm <mMax; mm++) {
00413 std::complex <float> fact(0,-mm*thetak);
00414 std::complex <float> expfact= exp(fact);
00415 std::complex <float> tempRho((*rhoOfkandm)(2*mm,kIntm1) ,(*rhoOfkandm)(2*mm+1,kIntm1) );
00416 float mmFac = float(1-2*(mm%2));
00417 if (IntensityFlag==1){ mmFac=1;}
00418 ImfTemp += expfact * tempRho + mmFac *conj(expfact*tempRho);
00419 }
00420 (*outCopy)(2*(CenterM-jky) ,CenterM-jkx) = ImfTemp.real();
00421 (*outCopy)(2*(CenterM-jky)+1,CenterM-jkx) = ImfTemp.imag();
00422 }
00423 }
00424
00425
00426 }
00427 }
00428 outCopy->update();
00429 outCopy->set_complex(true);
00430 if(outCopy->get_ysize()==1 && outCopy->get_zsize()==1) outCopy->set_complex_x(true);
00431 outCopy->set_ri(true);
00432 outCopy->set_FH(false);
00433 outCopy->set_fftodd(true);
00434 outCopy->set_shuffled(true);
00435 return outCopy;
00436 } else {
00437 LOGERR("can't be an FH image not this size");
00438 throw ImageFormatException("something strange about this image: not a FH");
00439
00440 }
00441 }
00442
00443
00444 EMData *EMData::FH2Real(int Size, float OverSamplekB, int)
00445 {
00446 EMData* FFT= FH2F(Size,OverSamplekB,0);
00447 FFT->process_inplace("xform.fourierorigin.tocorner");
00448 EMData* eguess= FFT ->do_ift();
00449 return eguess;
00450 }
00451
00452 float dist(int lnlen, const float* line_1, const float* line_2)
00453 {
00454 float dis2=0.0;
00455 for( int i=0; i < lnlen; ++i) {
00456 float tmp = line_1[i] - line_2[i];
00457 dis2 += tmp*tmp;
00458 }
00459
00460 return dis2;
00461 }
00462
00463 float dist_r(int lnlen, const float* line_1, const float* line_2)
00464 {
00465 double dis2 = 0.0;
00466 for( int i=0; i < lnlen; ++i ) {
00467 float tmp = line_1[lnlen-1-i] - line_2[i];
00468 dis2 += tmp*tmp;
00469 }
00470 return static_cast<float>(std::sqrt(dis2));
00471 }
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503 float EMData::cm_euc(EMData* sinoj, int n1, int n2)
00504 {
00505 int lnlen = get_xsize();
00506 float* line_1 = get_data() + n1 * lnlen;
00507 float* line_2 = sinoj->get_data() + n2 * lnlen;
00508 return dist(lnlen, line_1, line_2);
00509 }
00510
00511 EMData* EMData::rotavg() {
00512
00513 int rmax;
00514
00515 ENTERFUNC;
00516
00517 if (ny<2 && nz <2) {
00518 LOGERR("No 1D images.");
00519 throw ImageDimensionException("No 1D images!");
00520 }
00521 vector<int> saved_offsets = get_array_offsets();
00522 set_array_offsets(-nx/2,-ny/2,-nz/2);
00523 #ifdef _WIN32
00524
00525 if ( nz == 1 ) {
00526 rmax = _MIN(nx/2 + nx%2, ny/2 + ny%2);
00527 } else {
00528 rmax = _MIN(nx/2 + nx%2, _MIN(ny/2 + ny%2, nz/2 + nz%2));
00529 }
00530 #else
00531
00532 if ( nz == 1 ) {
00533 rmax = std::min(nx/2 + nx%2, ny/2 + ny%2);
00534 } else {
00535 rmax = std::min(nx/2 + nx%2, std::min(ny/2 + ny%2, nz/2 + nz%2));
00536 }
00537 #endif //_WIN32
00538 EMData* ret = new EMData();
00539 ret->set_size(rmax+1, 1, 1);
00540 ret->to_zero();
00541 vector<float> count(rmax+1);
00542 for (int k = -nz/2; k < nz/2 + nz%2; k++) {
00543 if (abs(k) > rmax) continue;
00544 for (int j = -ny/2; j < ny/2 + ny%2; j++) {
00545 if (abs(j) > rmax) continue;
00546 for (int i = -nx/2; i < nx/2 + nx%2; i++) {
00547 float r = std::sqrt(float(k*k) + float(j*j) + float(i*i));
00548 int ir = int(r);
00549 if (ir >= rmax) continue;
00550 float frac = r - float(ir);
00551 (*ret)(ir) += (*this)(i,j,k)*(1.0f - frac);
00552 (*ret)(ir+1) += (*this)(i,j,k)*frac;
00553 count[ir] += 1.0f - frac;
00554 count[ir+1] += frac;
00555 }
00556 }
00557 }
00558 for (int ir = 0; ir <= rmax; ir++) {
00559 #ifdef _WIN32
00560 (*ret)(ir) /= _MAX(count[ir],1.0f);
00561 #else
00562 (*ret)(ir) /= std::max(count[ir],1.0f);
00563 #endif //_WIN32
00564 }
00565
00566 set_array_offsets(saved_offsets);
00567 ret->update();
00568 EXITFUNC;
00569 return ret;
00570 }
00571
00572 EMData* EMData::rotavg_i() {
00573
00574 int rmax;
00575 ENTERFUNC;
00576 if ( ny == 1 && nz == 1 ) {
00577 LOGERR("Input image must be 2-D or 3-D!");
00578 throw ImageDimensionException("Input image must be 2-D or 3-D!");
00579 }
00580
00581 EMData* avg1D = new EMData();
00582 EMData* result = new EMData();
00583
00584 result->set_size(nx,ny,nz);
00585 result->to_zero();
00586 result->set_array_offsets(-nx/2, -ny/2, -nz/2);
00587
00588 if ( nz == 1 ) {
00589 #ifdef _WIN32
00590 rmax = _cpp_min(nx/2 + nx%2, ny/2 + ny%2);
00591 } else {
00592 rmax = _cpp_min(nx/2 + nx%2, _cpp_min(ny/2 + ny%2, nz/2 + nz%2));
00593 #else
00594 rmax = std::min(nx/2 + nx%2, ny/2 + ny%2);
00595 } else {
00596 rmax = std::min(nx/2 + nx%2, std::min(ny/2 + ny%2, nz/2 + nz%2));
00597 #endif //_WIN32
00598 }
00599
00600 avg1D = rotavg();
00601 float padded_value = 0.0, r;
00602 int i, j, k, ir;
00603 long int number_of_pixels = 0;
00604 for ( k = -nz/2; k < nz/2 + nz%2; k++) {
00605 if (abs(k) > rmax) continue;
00606 for ( j = -ny/2; j < ny/2 + ny%2; j++) {
00607 if (abs(j) > rmax) continue;
00608 for (i = -nx/2; i < nx/2 + nx%2; i++) {
00609 r = std::sqrt(float(k*k) + float(j*j) + float(i*i));
00610 ir = int(r);
00611 if (ir > rmax || ir < rmax-2 ) continue ;
00612 else {
00613 padded_value += (*avg1D)(ir) ;
00614 number_of_pixels++ ;
00615 }
00616 }
00617 }
00618 }
00619 padded_value /= number_of_pixels;
00620 for ( k = -nz/2; k < nz/2 + nz%2; k++) {
00621 for ( j = -ny/2; j < ny/2 + ny%2; j++) {
00622 for ( i = -nx/2; i < nx/2 + nx%2; i++) {
00623 r = std::sqrt(float(k*k) + float(j*j) + float(i*i));
00624 ir = int(r);
00625 if (ir >= rmax) (*result)(i,j,k) = padded_value ;
00626 else (*result)(i,j,k) = (*avg1D)(ir)+((*avg1D)(ir+1)-(*avg1D)(ir))*(r - float(ir));
00627
00628 }
00629 }
00630 }
00631 result->update();
00632 result->set_array_offsets(0,0,0);
00633 EXITFUNC;
00634 return result;
00635 }
00636
00637
00638 EMData* EMData::mult_radial(EMData* radial) {
00639
00640 ENTERFUNC;
00641 if ( ny == 1 && nz == 1 ) {
00642 LOGERR("Input image must be 2-D or 3-D!");
00643 throw ImageDimensionException("Input image must be 2-D or 3-D!");
00644 }
00645
00646 EMData* result = this->copy_head();
00647
00648 result->to_zero();
00649 result->set_array_offsets(-nx/2, -ny/2, -nz/2);
00650 this->set_array_offsets(-nx/2, -ny/2, -nz/2);
00651 int rmax = radial->get_xsize();
00652 int i, j, k, ir;
00653 float r;
00654 for ( k = -nz/2; k < nz/2+nz%2; k++) {
00655 for ( j = -ny/2; j < ny/2+ny%2; j++) {
00656 for ( i = -nx/2; i < nx/2+nx%2; i++) {
00657 r = std::sqrt(float(k*k) + float(j*j) + float(i*i));
00658 ir = int(r);
00659 if(ir < rmax-1) (*result)(i,j,k) = (*this)(i,j,k) * ((*radial)(ir)+((*radial)(ir+1)-(*radial)(ir))*(r - float(ir)));
00660 }
00661 }
00662 }
00663 result->update();
00664 result->set_array_offsets(0,0,0);
00665 this->set_array_offsets(0,0,0);
00666 EXITFUNC;
00667 return result;
00668 }
00669
00670 #define rdata(i,j,k) rdata[(i-1)+((j-1)+(k-1)*ny)*nx]
00671 #define square(x) ((x)*(x))
00672 vector<float> EMData::cog() {
00673
00674 vector<float> cntog;
00675 int ndim = get_ndim();
00676 int i=1,j=1,k=1;
00677 float val,sum1=0.f,MX=0.f,RG=0.f,MY=0.f,MZ=0.f,r=0.f;
00678
00679 if (ndim == 1) {
00680 for ( i = 1;i <= nx; i++) {
00681 val = rdata(i,j,k);
00682 sum1 += val;
00683 MX += ((i-1)*val);
00684 }
00685 MX=(MX/sum1);
00686 for ( i = 1;i <= nx; i++) {
00687 val = rdata(i,j,k);
00688 sum1 += val;
00689 RG += val*(square(MX - (i-1)));
00690 }
00691 RG=std::sqrt(RG/sum1);
00692 MX=MX-(nx/2);
00693 cntog.push_back(MX);
00694 cntog.push_back(RG);
00695 #ifdef _WIN32
00696 cntog.push_back((float)Util::round(MX));
00697 #else
00698 cntog.push_back(round(MX));
00699 #endif //_WIN32
00700 } else if (ndim == 2) {
00701 for (j=1;j<=ny;j++) {
00702 for (i=1;i<=nx;i++) {
00703 val = rdata(i,j,k);
00704 sum1 += val;
00705 MX += ((i-1)*val);
00706 MY += ((j-1)*val);
00707 }
00708 }
00709 MX=(MX/sum1);
00710 MY=(MY/sum1);
00711 sum1=0.f;
00712 RG=0.f;
00713 for (j=1;j<=ny;j++) {
00714 r = (square(MY-(j-1)));
00715 for (i=1;i<=nx;i++) {
00716 val = rdata(i,j,k);
00717 sum1 += val;
00718 RG += val*(square(MX - (i-1)) + r);
00719 }
00720 }
00721 RG = std::sqrt(RG/sum1);
00722 MX = MX - nx/2;
00723 MY = MY - ny/2;
00724 cntog.push_back(MX);
00725 cntog.push_back(MY);
00726 cntog.push_back(RG);
00727 #ifdef _WIN32
00728 cntog.push_back((float)Util::round(MX));cntog.push_back((float)Util::round(MY));
00729 #else
00730 cntog.push_back(round(MX));cntog.push_back(round(MY));
00731 #endif //_WIN32
00732 } else {
00733 for (k = 1;k <= nz;k++) {
00734 for (j=1;j<=ny;j++) {
00735 for (i=1;i<=nx;i++) {
00736 val = rdata(i,j,k);
00737 sum1 += val;
00738 MX += ((i-1)*val);
00739 MY += ((j-1)*val);
00740 MZ += ((k-1)*val);
00741 }
00742 }
00743 }
00744 MX = MX/sum1;
00745 MY = MY/sum1;
00746 MZ = MZ/sum1;
00747 sum1=0.f;
00748 RG=0.f;
00749 for (k = 1;k <= nz;k++) {
00750 for (j=1;j<=ny;j++) {
00751 float r = (square(MY-(j-1)) + square(MZ - (k-1)));
00752 for (i=1;i<=nx;i++) {
00753 val = rdata(i,j,k);
00754 sum1 += val;
00755 RG += val*(square(MX - (i-1)) + r);
00756 }
00757 }
00758 }
00759 RG = std::sqrt(RG/sum1);
00760 MX = MX - nx/2;
00761 MY = MY - ny/2;
00762 MZ = MZ - nz/2;
00763 cntog.push_back(MX);
00764 cntog.push_back(MY);
00765 cntog.push_back(MZ);
00766 cntog.push_back(RG);
00767 #ifdef _WIN32
00768 cntog.push_back((float)Util::round(MX));cntog.push_back((float)Util::round(MY));cntog.push_back((float)Util::round(MZ));
00769 #else
00770 cntog.push_back(round(MX));cntog.push_back(round(MY));cntog.push_back(round(MZ));
00771 #endif //_WIN32
00772 }
00773 return cntog;
00774 }
00775 #undef square
00776 #undef rdata
00777
00778
00779
00780
00781
00782 vector < float >EMData::calc_fourier_shell_correlation(EMData * with, float w)
00783 {
00784 ENTERFUNC;
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809 int needfree=0, nx, ny, nz, nx2, ny2, nz2, ix, iy, iz, kz, ky, ii;
00810 float dx2, dy2, dz2, argx, argy, argz;
00811
00812 if (!with) {
00813 throw NullPointerException("NULL input image");
00814 }
00815
00816
00817 EMData *f = this;
00818 EMData *g = with;
00819
00820 nx = f->get_xsize();
00821 ny = f->get_ysize();
00822 nz = f->get_zsize();
00823
00824 if (ny==0 && nz==0) {
00825 throw ImageFormatException( "Cannot calculate FSC for 1D images");
00826 }
00827
00828 if (f->is_complex()) nx = (nx - 2 + f->is_fftodd());
00829 int lsd2 = (nx + 2 - nx%2) ;
00830
00831
00832 EMData* fpimage = NULL;
00833 if(f->is_complex()) fpimage = f;
00834 else {fpimage= f->norm_pad(false, 1); fpimage->do_fft_inplace(); needfree|=1;}
00835
00836
00837
00838 EMData* gpimage = NULL;
00839 if(g->is_complex()) gpimage = g;
00840 else {gpimage= g->norm_pad(false, 1); gpimage->do_fft_inplace(); needfree|=2;}
00841
00842
00843 float *d1 = fpimage->get_data();
00844 float *d2 = gpimage->get_data();
00845
00846 nx2=nx/2; ny2 = ny/2; nz2 = nz/2;
00847 dx2 = 1.0f/float(nx2)/float(nx2);
00848 dy2 = 1.0f/float(ny2)/float(ny2);
00849
00850 #ifdef _WIN32
00851 dz2 = 1.0f / _MAX(float(nz2),1.0f)/_MAX(float(nz2),1.0f);
00852
00853 int inc = Util::round(float( _MAX( _MAX(nx2,ny2),nz2) )/w );
00854 #else
00855 dz2 = 1.0f/std::max(float(nz2),1.0f)/std::max(float(nz2),1.0f);
00856
00857 int inc = Util::round(float(std::max(std::max(nx2,ny2),nz2))/w);
00858 #endif //_WIN32
00859
00860 double *ret = new double[inc+1];
00861 double *n1 = new double[inc+1];
00862 double *n2 = new double[inc+1];
00863 float *lr = new float[inc+1];
00864 for (int i = 0; i <= inc; i++) {
00865 ret[i] = 0; n1[i] = 0; n2[i] = 0; lr[i]=0;
00866 }
00867
00868 for ( iz = 0; iz <= nz-1; iz++) {
00869 if(iz>nz2) kz=iz-nz; else kz=iz; argz = float(kz*kz)*dz2;
00870 for ( iy = 0; iy <= ny-1; iy++) {
00871 if(iy>ny2) ky=iy-ny; else ky=iy; argy = argz + float(ky*ky)*dy2;
00872 for ( ix = 0; ix <= lsd2-1; ix+=2) {
00873
00874 if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00875 argx = 0.5f*std::sqrt(argy + float(ix*ix)*0.25f*dx2);
00876 int r = Util::round(inc*2*argx);
00877 if(r <= inc) {
00878 ii = ix + (iy + iz * ny)* lsd2;
00879 ret[r] += d1[ii] * double(d2[ii]) + d1[ii + 1] * double(d2[ii + 1]);
00880 n1[r] += d1[ii] * double(d1[ii]) + d1[ii + 1] * double(d1[ii + 1]);
00881 n2[r] += d2[ii] * double(d2[ii]) + d2[ii + 1] * double(d2[ii + 1]);
00882 lr[r] += 2;
00883 }
00884 }
00885 }
00886 }
00887 }
00888
00889 int linc = 0;
00890 for (int i = 0; i <= inc; i++) if(lr[i]>0) linc++;
00891
00892 vector < float >result(linc*3);
00893
00894 ii = -1;
00895 for (int i = 0; i <= inc; i++) {
00896 if(lr[i]>0) {
00897 ii++;
00898 result[ii] = float(i)/float(2*inc);
00899 result[ii+linc] = float(ret[i] / (std::sqrt(n1[i] * n2[i])));
00900 result[ii+2*linc] = lr[i] ;
00901 }
00902
00903
00904
00905
00906 }
00907
00908 if( ret ) {
00909 delete[]ret;
00910 ret = 0;
00911 }
00912
00913 if( n1 ) {
00914 delete[]n1;
00915 n1 = 0;
00916 }
00917 if( n2 ) {
00918 delete[]n2;
00919 n2 = 0;
00920 }
00921
00922 if (needfree&1) {
00923 if( fpimage ) {
00924 delete fpimage;
00925 fpimage = 0;
00926 }
00927 }
00928 if (needfree&2) {
00929 if( gpimage ) {
00930 delete gpimage;
00931 gpimage = 0;
00932 }
00933 }
00934
00935 EXITFUNC;
00936 return result;
00937 }
00938
00939
00940 EMData* EMData::symvol(string symString) {
00941 ENTERFUNC;
00942 int nsym = Transform::get_nsym(symString);
00943 Transform sym;
00944
00945 EMData *svol = new EMData;
00946 svol->set_size(nx, ny, nz);
00947 svol->to_zero();
00948
00949 EMData* symcopy = new EMData;
00950 symcopy->set_size(nx, ny, nz);
00951
00952
00953 for (int isym = 0; isym < nsym; isym++) {
00954 Transform rm = sym.get_sym(symString, isym);
00955 symcopy = this -> rot_scale_trans(rm);
00956 *svol += (*symcopy);
00957 }
00958 *svol /= ((float) nsym);
00959 svol->update();
00960 EXITFUNC;
00961 return svol;
00962 }
00963
00964 #define proj(ix,iy,iz) proj[ix-1+(iy-1+(iz-1)*ny)*nx]
00965 #define pnewimg(ix,iy,iz) pnewimg[ix-1+(iy-1+(iz-1)*ny)*nx]
00966 EMData* EMData::average_circ_sub() const
00967 {
00968
00969
00970 ENTERFUNC;
00971 const EMData* const image = this;
00972 EMData* newimg = copy_head();
00973 float *proj = image->get_data();
00974 float *pnewimg = newimg->get_data();
00975
00976 float r2 = static_cast<float>( (nx/2)*(nx/2) );
00977 float qs=0.0f;
00978 int m=0;
00979 int ncz = nz/2 + 1;
00980 int ncy = ny/2 + 1;
00981 int ncx = nx/2 + 1;
00982 for (int iz = 1; iz <= nz; iz++) {
00983 float yy = static_cast<float>( (iz-ncz)*(iz-ncz) );
00984 for (int iy = 1; iy <=ny; iy++) { float xx = yy + (iy-ncy)*(iy-ncy);
00985 for (int ix = 1; ix <= nx; ix++) {
00986 if ( xx+float((ix-ncx)*(ix-ncx)) > r2 ) {
00987 qs += proj(ix,iy,iz);
00988 m++;
00989 }
00990 }
00991 }
00992 }
00993
00994
00995 if( m > 0 ) qs /= m;
00996
00997 for (int iz = 1; iz <= nz; iz++)
00998 for (int iy = 1; iy <= ny; iy++)
00999 for (int ix = 1; ix <= nx; ix++)
01000 pnewimg(ix,iy,iz) = proj(ix,iy,iz) - qs;
01001 newimg->update();
01002 return newimg;
01003 EXITFUNC;
01004 }
01005
01006
01007
01008
01009
01010 void EMData::onelinenn(int j, int n, int n2, EMData* wptr, EMData* bi, const Transform& tf)
01011 {
01012
01013 int jp = (j >= 0) ? j+1 : n+j+1;
01014
01015
01016 for (int i = 0; i <= n2; i++) {
01017 if (((i*i+j*j) < n*n/4) && !((0 == i) && (j < 0))) {
01018
01019 float xnew = i*tf[0][0] + j*tf[1][0];
01020 float ynew = i*tf[0][1] + j*tf[1][1];
01021 float znew = i*tf[0][2] + j*tf[1][2];
01022 std::complex<float> btq;
01023 if (xnew < 0.) {
01024 xnew = -xnew;
01025 ynew = -ynew;
01026 znew = -znew;
01027 btq = conj(bi->cmplx(i,jp));
01028 } else {
01029 btq = bi->cmplx(i,jp);
01030 }
01031 int ixn = int(xnew + 0.5 + n) - n;
01032 int iyn = int(ynew + 0.5 + n) - n;
01033 int izn = int(znew + 0.5 + n) - n;
01034 if ((ixn <= n2) && (iyn >= -n2) && (iyn <= n2) && (izn >= -n2) && (izn <= n2)) {
01035 if (ixn >= 0) {
01036 int iza, iya;
01037 if (izn >= 0) iza = izn + 1;
01038 else iza = n + izn + 1;
01039
01040 if (iyn >= 0) iya = iyn + 1;
01041 else iya = n + iyn + 1;
01042
01043 cmplx(ixn,iya,iza) += btq;
01044
01045 (*wptr)(ixn,iya,iza)++;
01046 } else {
01047 int izt, iyt;
01048 if (izn > 0) izt = n - izn + 1;
01049 else izt = -izn + 1;
01050
01051 if (iyn > 0) iyt = n - iyn + 1;
01052 else iyt = -iyn + 1;
01053
01054 cmplx(-ixn,iyt,izt) += conj(btq);
01055
01056 (*wptr)(-ixn,iyt,izt)++;
01057 }
01058 }
01059 }
01060 }
01061 }
01062
01063
01064 void EMData::onelinenn_mult(int j, int n, int n2, EMData* wptr, EMData* bi, const Transform& tf, int mult)
01065 {
01066
01067 int jp = (j >= 0) ? j+1 : n+j+1;
01068
01069
01070 for (int i = 0; i <= n2; i++) {
01071 if (((i*i+j*j) < n*n/4) && !((0 == i) && (j < 0))) {
01072
01073 float xnew = i*tf[0][0] + j*tf[1][0];
01074 float ynew = i*tf[0][1] + j*tf[1][1];
01075 float znew = i*tf[0][2] + j*tf[1][2];
01076 std::complex<float> btq;
01077 if (xnew < 0.) {
01078 xnew = -xnew;
01079 ynew = -ynew;
01080 znew = -znew;
01081 btq = conj(bi->cmplx(i,jp));
01082 } else {
01083 btq = bi->cmplx(i,jp);
01084 }
01085 int ixn = int(xnew + 0.5 + n) - n;
01086 int iyn = int(ynew + 0.5 + n) - n;
01087 int izn = int(znew + 0.5 + n) - n;
01088 if ((ixn <= n2) && (iyn >= -n2) && (iyn <= n2) && (izn >= -n2) && (izn <= n2)) {
01089 if (ixn >= 0) {
01090 int iza, iya;
01091 if (izn >= 0) iza = izn + 1;
01092 else iza = n + izn + 1;
01093
01094 if (iyn >= 0) iya = iyn + 1;
01095 else iya = n + iyn + 1;
01096
01097 cmplx(ixn,iya,iza) += btq*float(mult);
01098 (*wptr)(ixn,iya,iza)+=float(mult);
01099 } else {
01100 int izt, iyt;
01101 if (izn > 0) izt = n - izn + 1;
01102 else izt = -izn + 1;
01103
01104 if (iyn > 0) iyt = n - iyn + 1;
01105 else iyt = -iyn + 1;
01106
01107 cmplx(-ixn,iyt,izt) += conj(btq)*float(mult);
01108 (*wptr)(-ixn,iyt,izt)+=float(mult);
01109 }
01110 }
01111 }
01112 }
01113 }
01114
01115 void EMData::nn(EMData* wptr, EMData* myfft, const Transform& tf, int mult)
01116 {
01117 ENTERFUNC;
01118 int nxc = attr_dict["nxc"];
01119
01120 vector<int> saved_offsets = get_array_offsets();
01121 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01122 set_array_offsets(0,1,1);
01123 myfft->set_array_offsets(0,1);
01124
01125 if( mult == 1 ) {
01126 for (int iy = -ny/2 + 1; iy <= ny/2; iy++) onelinenn(iy, ny, nxc, wptr, myfft, tf);
01127 } else {
01128 for (int iy = -ny/2 + 1; iy <= ny/2; iy++) onelinenn_mult(iy, ny, nxc, wptr, myfft, tf, mult);
01129 }
01130
01131 set_array_offsets(saved_offsets);
01132 myfft->set_array_offsets(myfft_saved_offsets);
01133 EXITFUNC;
01134 }
01135
01136 void EMData::nn_SSNR(EMData* wptr, EMData* wptr2, EMData* myfft, const Transform& tf, int)
01137 {
01138 ENTERFUNC;
01139 int nxc = attr_dict["nxc"];
01140
01141 vector<int> saved_offsets = get_array_offsets();
01142 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01143
01144 set_array_offsets(0,1,1);
01145 myfft->set_array_offsets(0,1);
01146
01147 int iymin = is_fftodd() ? -ny/2 : -ny/2 + 1 ;
01148 int iymax = ny/2;
01149 int izmin = is_fftodd() ? -nz/2 : -nz/2 + 1 ;
01150 int izmax = nz/2;
01151
01152 for (int iy = iymin; iy <= iymax; iy++) {
01153 int jp = iy >= 0 ? iy+1 : ny+iy+1;
01154 for (int ix = 0; ix <= nxc; ix++) {
01155 if (( 4*(ix*ix+iy*iy) < ny*ny ) && !( ix == 0 && iy < 0 ) ) {
01156 float xnew = ix*tf[0][0] + iy*tf[1][0];
01157 float ynew = ix*tf[0][1] + iy*tf[1][1];
01158 float znew = ix*tf[0][2] + iy*tf[1][2];
01159 std::complex<float> btq;
01160 if (xnew < 0.0) {
01161 xnew = -xnew;
01162 ynew = -ynew;
01163 znew = -znew;
01164 btq = conj(myfft->cmplx(ix,jp));
01165 } else {
01166 btq = myfft->cmplx(ix,jp);
01167 }
01168 int ixn = int(xnew + 0.5 + nx) - nx;
01169 int iyn = int(ynew + 0.5 + ny) - ny;
01170 int izn = int(znew + 0.5 + nz) - nz;
01171 if ((ixn <= nxc) && (iyn >= iymin) && (iyn <= iymax) && (izn >= izmin) && (izn <= izmax)) {
01172 if (ixn >= 0) {
01173 int iza, iya;
01174 if (izn >= 0) iza = izn + 1;
01175 else iza = nz + izn + 1;
01176
01177 if (iyn >= 0) iya = iyn + 1;
01178 else iya = ny + iyn + 1;
01179
01180 cmplx(ixn,iya,iza) += btq;
01181 (*wptr)(ixn,iya,iza)++;
01182 (*wptr2)(ixn,iya,iza) += norm(btq);
01183 } else {
01184 int izt, iyt;
01185 if (izn > 0) izt = nz - izn + 1;
01186 else izt = -izn + 1;
01187
01188 if (iyn > 0) iyt = ny - iyn + 1;
01189 else iyt = -iyn + 1;
01190
01191 cmplx(-ixn,iyt,izt) += conj(btq);
01192 (*wptr)(-ixn,iyt,izt)++;
01193 (*wptr2)(-ixn,iyt,izt) += norm(btq);
01194 }
01195 }
01196 }
01197 }
01198 }
01199 set_array_offsets(saved_offsets);
01200 myfft->set_array_offsets(myfft_saved_offsets);
01201 EXITFUNC;
01202 }
01203
01204
01205
01206 void EMData::symplane0(EMData* wptr) {
01207 ENTERFUNC;
01208 int nxc = attr_dict["nxc"];
01209 int n = nxc*2;
01210
01211 vector<int> saved_offsets = get_array_offsets();
01212 set_array_offsets(0,1,1);
01213 for (int iza = 2; iza <= nxc; iza++) {
01214 for (int iya = 2; iya <= nxc; iya++) {
01215 cmplx(0,iya,iza) += conj(cmplx(0,n-iya+2,n-iza+2));
01216 (*wptr)(0,iya,iza) += (*wptr)(0,n-iya+2,n-iza+2);
01217 cmplx(0,n-iya+2,n-iza+2) = conj(cmplx(0,iya,iza));
01218 (*wptr)(0,n-iya+2,n-iza+2) = (*wptr)(0,iya,iza);
01219 cmplx(0,n-iya+2,iza) += conj(cmplx(0,iya,n-iza+2));
01220 (*wptr)(0,n-iya+2,iza) += (*wptr)(0,iya,n-iza+2);
01221 cmplx(0,iya,n-iza+2) = conj(cmplx(0,n-iya+2,iza));
01222 (*wptr)(0,iya,n-iza+2) = (*wptr)(0,n-iya+2,iza);
01223 }
01224 }
01225 for (int iya = 2; iya <= nxc; iya++) {
01226 cmplx(0,iya,1) += conj(cmplx(0,n-iya+2,1));
01227 (*wptr)(0,iya,1) += (*wptr)(0,n-iya+2,1);
01228 cmplx(0,n-iya+2,1) = conj(cmplx(0,iya,1));
01229 (*wptr)(0,n-iya+2,1) = (*wptr)(0,iya,1);
01230 }
01231 for (int iza = 2; iza <= nxc; iza++) {
01232 cmplx(0,1,iza) += conj(cmplx(0,1,n-iza+2));
01233 (*wptr)(0,1,iza) += (*wptr)(0,1,n-iza+2);
01234 cmplx(0,1,n-iza+2) = conj(cmplx(0,1,iza));
01235 (*wptr)(0,1,n-iza+2) = (*wptr)(0,1,iza);
01236 }
01237 EXITFUNC;
01238 }
01239
01240 void EMData::symplane1(EMData* wptr, EMData* wptr2) {
01241 ENTERFUNC;
01242 int nxc = attr_dict["nxc"];
01243 int n = nxc*2;
01244 vector<int> saved_offsets = get_array_offsets();
01245 set_array_offsets(0,1,1);
01246 for (int iza = 2; iza <= nxc; iza++) {
01247 for (int iya = 2; iya <= nxc; iya++) {
01248 cmplx(0,iya,iza) += conj(cmplx(0,n-iya+2,n-iza+2));
01249 (*wptr)(0,iya,iza) += (*wptr)(0,n-iya+2,n-iza+2);
01250 (*wptr2)(0,iya,iza) += (*wptr2)(0,n-iya+2,n-iza+2);
01251 cmplx(0,n-iya+2,n-iza+2) = conj(cmplx(0,iya,iza));
01252 (*wptr)(0,n-iya+2,n-iza+2) = (*wptr)(0,iya,iza);
01253 (*wptr2)(0,n-iya+2,n-iza+2) = (*wptr2)(0,iya,iza);
01254 cmplx(0,n-iya+2,iza) += conj(cmplx(0,iya,n-iza+2));
01255 (*wptr)(0,n-iya+2,iza) += (*wptr)(0,iya,n-iza+2);
01256 (*wptr2)(0,n-iya+2,iza) += (*wptr2)(0,iya,n-iza+2);
01257 cmplx(0,iya,n-iza+2) = conj(cmplx(0,n-iya+2,iza));
01258 (*wptr)(0,iya,n-iza+2) = (*wptr)(0,n-iya+2,iza);
01259 (*wptr2)(0,iya,n-iza+2) = (*wptr2)(0,n-iya+2,iza);
01260 }
01261 }
01262 for (int iya = 2; iya <= nxc; iya++) {
01263 cmplx(0,iya,1) += conj(cmplx(0,n-iya+2,1));
01264 (*wptr)(0,iya,1) += (*wptr)(0,n-iya+2,1);
01265 (*wptr2)(0,iya,1) += (*wptr2)(0,n-iya+2,1);
01266 cmplx(0,n-iya+2,1) = conj(cmplx(0,iya,1));
01267 (*wptr)(0,n-iya+2,1) = (*wptr)(0,iya,1);
01268 (*wptr2)(0,n-iya+2,1) = (*wptr2)(0,iya,1);
01269 }
01270 for (int iza = 2; iza <= nxc; iza++) {
01271 cmplx(0,1,iza) += conj(cmplx(0,1,n-iza+2));
01272 (*wptr)(0,1,iza) += (*wptr)(0,1,n-iza+2);
01273 (*wptr2)(0,1,iza) += (*wptr2)(0,1,n-iza+2);
01274 cmplx(0,1,n-iza+2) = conj(cmplx(0,1,iza));
01275 (*wptr)(0,1,n-iza+2) = (*wptr)(0,1,iza);
01276 (*wptr2)(0,1,n-iza+2) = (*wptr2)(0,1,iza);
01277 }
01278 EXITFUNC;
01279 }
01280
01281 void EMData::symplane2(EMData* wptr, EMData* wptr2, EMData* wptr3) {
01282 ENTERFUNC;
01283 int nxc = attr_dict["nxc"];
01284 int n = nxc*2;
01285 vector<int> saved_offsets = get_array_offsets();
01286 set_array_offsets(0,1,1);
01287 for (int iza = 2; iza <= nxc; iza++) {
01288 for (int iya = 2; iya <= nxc; iya++) {
01289 cmplx(0,iya,iza) += conj(cmplx(0,n-iya+2,n-iza+2));
01290 (*wptr)(0,iya,iza) += (*wptr)(0,n-iya+2,n-iza+2);
01291 (*wptr2)(0,iya,iza) += (*wptr2)(0,n-iya+2,n-iza+2);
01292 (*wptr3)(0,iya,iza) += (*wptr3)(0,n-iya+2,n-iza+2);
01293
01294 cmplx(0,n-iya+2,n-iza+2) = conj(cmplx(0,iya,iza));
01295 (*wptr)(0,n-iya+2,n-iza+2) = (*wptr)(0,iya,iza);
01296 (*wptr2)(0,n-iya+2,n-iza+2) = (*wptr2)(0,iya,iza);
01297 (*wptr3)(0,n-iya+2,n-iza+2) = (*wptr3)(0,iya,iza);
01298
01299 cmplx(0,n-iya+2,iza) += conj(cmplx(0,iya,n-iza+2));
01300 (*wptr)(0,n-iya+2,iza) += (*wptr)(0,iya,n-iza+2);
01301 (*wptr2)(0,n-iya+2,iza) += (*wptr2)(0,iya,n-iza+2);
01302 (*wptr3)(0,n-iya+2,iza) += (*wptr3)(0,iya,n-iza+2);
01303
01304 cmplx(0,iya,n-iza+2) = conj(cmplx(0,n-iya+2,iza));
01305 (*wptr)(0,iya,n-iza+2) = (*wptr)(0,n-iya+2,iza);
01306 (*wptr2)(0,iya,n-iza+2) = (*wptr2)(0,n-iya+2,iza);
01307 (*wptr3)(0,iya,n-iza+2) = (*wptr3)(0,n-iya+2,iza);
01308 }
01309 }
01310 for (int iya = 2; iya <= nxc; iya++) {
01311 cmplx(0,iya,1) += conj(cmplx(0,n-iya+2,1));
01312 (*wptr)(0,iya,1) += (*wptr)(0,n-iya+2,1);
01313 (*wptr2)(0,iya,1) += (*wptr2)(0,n-iya+2,1);
01314 (*wptr3)(0,iya,1) += (*wptr3)(0,n-iya+2,1);
01315
01316 cmplx(0,n-iya+2,1) = conj(cmplx(0,iya,1));
01317 (*wptr)(0,n-iya+2,1) = (*wptr)(0,iya,1);
01318 (*wptr2)(0,n-iya+2,1) = (*wptr2)(0,iya,1);
01319 (*wptr3)(0,n-iya+2,1) = (*wptr3)(0,iya,1);
01320 }
01321 for (int iza = 2; iza <= nxc; iza++) {
01322 cmplx(0,1,iza) += conj(cmplx(0,1,n-iza+2));
01323 (*wptr)(0,1,iza) += (*wptr)(0,1,n-iza+2);
01324 (*wptr2)(0,1,iza) += (*wptr2)(0,1,n-iza+2);
01325 (*wptr3)(0,1,iza) += (*wptr3)(0,1,n-iza+2);
01326
01327 cmplx(0,1,n-iza+2) = conj(cmplx(0,1,iza));
01328 (*wptr)(0,1,n-iza+2) = (*wptr)(0,1,iza);
01329 (*wptr2)(0,1,n-iza+2) = (*wptr2)(0,1,iza);
01330 (*wptr3)(0,1,n-iza+2) = (*wptr3)(0,1,iza);
01331 }
01332 EXITFUNC;
01333 }
01334
01335
01336 class ctf_store
01337 {
01338 public:
01339
01340 static void init( int winsize, const Ctf* ctf )
01341 {
01342 Dict params = ctf->to_dict();
01343
01344 m_winsize = winsize;
01345
01346 m_voltage = params["voltage"];
01347 m_pixel = params["apix"];
01348 m_cs = params["cs"];
01349 m_ampcont = params["ampcont"];
01350 m_bfactor = params["bfactor"];
01351 m_defocus = params["defocus"];
01352 m_winsize2= m_winsize*m_winsize;
01353 m_vecsize = m_winsize2/4;
01354 }
01355
01356 static float get_ctf( int r2 )
01357 {
01358 float ak = std::sqrt( r2/float(m_winsize2) )/m_pixel;
01359 return Util::tf( m_defocus, ak, m_voltage, m_cs, m_ampcont, m_bfactor, 1 );
01360 }
01361
01362 private:
01363
01364 static int m_winsize, m_winsize2, m_vecsize;
01365 static float m_cs;
01366 static float m_voltage;
01367 static float m_pixel;
01368 static float m_ampcont;
01369 static float m_bfactor;
01370 static float m_defocus;
01371 };
01372
01373
01374 int ctf_store::m_winsize, ctf_store::m_winsize2, ctf_store::m_vecsize;
01375
01376 float ctf_store::m_cs, ctf_store::m_voltage, ctf_store::m_pixel;
01377 float ctf_store::m_ampcont, ctf_store::m_bfactor;
01378 float ctf_store::m_defocus;
01379
01380
01381
01382 void EMData::onelinenn_ctf(int j, int n, int n2,
01383 EMData* w, EMData* bi, const Transform& tf, int mult) {
01384
01385 int remove = bi->get_attr_default( "remove", 0 );
01386
01387 int jp = (j >= 0) ? j+1 : n+j+1;
01388
01389 for (int i = 0; i <= n2; i++) {
01390 int r2 = i*i+j*j;
01391 if ( (r2<n*n/4) && !( (0==i) && (j<0) ) ) {
01392 float ctf = ctf_store::get_ctf( r2 );
01393 float xnew = i*tf[0][0] + j*tf[1][0];
01394 float ynew = i*tf[0][1] + j*tf[1][1];
01395 float znew = i*tf[0][2] + j*tf[1][2];
01396 std::complex<float> btq;
01397 if (xnew < 0.) {
01398 xnew = -xnew;
01399 ynew = -ynew;
01400 znew = -znew;
01401 btq = conj(bi->cmplx(i,jp));
01402 } else btq = bi->cmplx(i,jp);
01403 int ixn = int(xnew + 0.5 + n) - n;
01404 int iyn = int(ynew + 0.5 + n) - n;
01405 int izn = int(znew + 0.5 + n) - n;
01406 if ((ixn <= n2) && (iyn >= -n2) && (iyn <= n2) && (izn >= -n2) && (izn <= n2)) {
01407 if (ixn >= 0) {
01408 int iza, iya;
01409 if (izn >= 0) iza = izn + 1;
01410 else iza = n + izn + 1;
01411
01412 if (iyn >= 0) iya = iyn + 1;
01413 else iya = n + iyn + 1;
01414
01415 if(remove > 0 ) {
01416 cmplx(ixn,iya,iza) -= btq*ctf*float(mult);
01417 (*w)(ixn,iya,iza) -= ctf*ctf*mult;
01418 } else {
01419 cmplx(ixn,iya,iza) += btq*ctf*float(mult);
01420 (*w)(ixn,iya,iza) += ctf*ctf*mult;
01421 }
01422
01423
01424 } else {
01425 int izt, iyt;
01426 if (izn > 0) izt = n - izn + 1;
01427 else izt = -izn + 1;
01428
01429 if (iyn > 0) iyt = n - iyn + 1;
01430 else iyt = -iyn + 1;
01431
01432 if( remove > 0 ) {
01433 cmplx(-ixn,iyt,izt) -= conj(btq)*ctf*float(mult);
01434 (*w)(-ixn,iyt,izt) -= ctf*ctf*float(mult);
01435 } else {
01436 cmplx(-ixn,iyt,izt) += conj(btq)*ctf*float(mult);
01437 (*w)(-ixn,iyt,izt) += ctf*ctf*float(mult);
01438 }
01439
01440
01441 }
01442 }
01443 }
01444 }
01445 }
01446
01447 void EMData::onelinenn_ctf_applied(int j, int n, int n2,
01448 EMData* w, EMData* bi, const Transform& tf, int mult) {
01449
01450 int remove = bi->get_attr_default( "remove", 0 );
01451
01452 int jp = (j >= 0) ? j+1 : n+j+1;
01453
01454 for (int i = 0; i <= n2; i++) {
01455 int r2 = i*i + j*j;
01456 if ( (r2< n*n/4) && !((0==i) && (j< 0)) ) {
01457 float ctf = ctf_store::get_ctf(r2);
01458
01459
01460 float xnew = i*tf[0][0] + j*tf[1][0];
01461 float ynew = i*tf[0][1] + j*tf[1][1];
01462 float znew = i*tf[0][2] + j*tf[1][2];
01463 std::complex<float> btq;
01464 if (xnew < 0.) {
01465 xnew = -xnew;
01466 ynew = -ynew;
01467 znew = -znew;
01468 btq = conj(bi->cmplx(i,jp));
01469 } else btq = bi->cmplx(i,jp);
01470 int ixn = int(xnew + 0.5 + n) - n;
01471 int iyn = int(ynew + 0.5 + n) - n;
01472 int izn = int(znew + 0.5 + n) - n;
01473
01474 if ((ixn <= n2) && (iyn >= -n2) && (iyn <= n2) && (izn >= -n2) && (izn <= n2)) {
01475 if (ixn >= 0) {
01476 int iza, iya;
01477 if (izn >= 0) iza = izn + 1;
01478 else iza = n + izn + 1;
01479
01480 if (iyn >= 0) iya = iyn + 1;
01481 else iya = n + iyn + 1;
01482
01483 if( remove > 0 ) {
01484 cmplx(ixn,iya,iza) -= btq*float(mult);
01485 (*w)(ixn,iya,iza) -= mult*ctf*ctf;
01486 } else {
01487 cmplx(ixn,iya,iza) += btq*float(mult);
01488 (*w)(ixn,iya,iza) += mult*ctf*ctf;
01489 }
01490
01491 } else {
01492 int izt, iyt;
01493 if (izn > 0) izt = n - izn + 1;
01494 else izt = -izn + 1;
01495
01496 if (iyn > 0) iyt = n - iyn + 1;
01497 else iyt = -iyn + 1;
01498
01499
01500 if( remove > 0 ) {
01501 cmplx(-ixn,iyt,izt) -= conj(btq)*float(mult);
01502 (*w)(-ixn,iyt,izt) -= mult*ctf*ctf;
01503 } else {
01504 cmplx(-ixn,iyt,izt) += conj(btq)*float(mult);
01505 (*w)(-ixn,iyt,izt) += mult*ctf*ctf;
01506 }
01507
01508 }
01509 }
01510 }
01511 }
01512 }
01513
01514 void
01515 EMData::nn_ctf(EMData* w, EMData* myfft, const Transform& tf, int mult) {
01516 ENTERFUNC;
01517 int nxc = attr_dict["nxc"];
01518
01519 vector<int> saved_offsets = get_array_offsets();
01520 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01521 set_array_offsets(0,1,1);
01522 myfft->set_array_offsets(0,1);
01523
01524
01525 Ctf* ctf = myfft->get_attr("ctf");
01526 ctf_store::init( ny, ctf );
01527 if(ctf) {delete ctf; ctf=0;}
01528
01529
01530 for (int iy = -ny/2 + 1; iy <= ny/2; iy++) onelinenn_ctf(iy, ny, nxc, w, myfft, tf, mult);
01531 set_array_offsets(saved_offsets);
01532 myfft->set_array_offsets(myfft_saved_offsets);
01533 EXITFUNC;
01534 }
01535
01536 void
01537 EMData::nn_ctf_applied(EMData* w, EMData* myfft, const Transform& tf, int mult) {
01538 ENTERFUNC;
01539 int nxc = attr_dict["nxc"];
01540
01541 vector<int> saved_offsets = get_array_offsets();
01542 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01543 set_array_offsets(0,1,1);
01544 myfft->set_array_offsets(0,1);
01545
01546 Ctf* ctf = myfft->get_attr( "ctf" );
01547 ctf_store::init( ny, ctf );
01548 if(ctf) {delete ctf; ctf=0;}
01549
01550
01551
01552 for (int iy = -ny/2 + 1; iy <= ny/2; iy++) onelinenn_ctf_applied(iy, ny, nxc, w, myfft, tf, mult);
01553 set_array_offsets(saved_offsets);
01554 myfft->set_array_offsets(myfft_saved_offsets);
01555 EXITFUNC;
01556 }
01557
01558
01559
01560 void EMData::nn_SSNR_ctf(EMData* wptr, EMData* wptr2, EMData* wptr3, EMData* myfft, const Transform& tf, int)
01561 {
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571 ENTERFUNC;
01572 int nxc = attr_dict["nxc"];
01573 vector<int> saved_offsets = get_array_offsets();
01574 vector<int> myfft_saved_offsets = myfft->get_array_offsets();
01575 set_array_offsets(0,1,1);
01576 myfft->set_array_offsets(0,1);
01577
01578 Ctf* ctf = myfft->get_attr("ctf");
01579 ctf_store::init( ny, ctf );
01580 int iymin = is_fftodd() ? -ny/2 : -ny/2 + 1;
01581 int iymax = ny/2;
01582 int izmin = is_fftodd() ? -nz/2 : -nz/2 + 1;
01583 int izmax = nz/2;
01584
01585 for (int iy = iymin; iy <= iymax; iy++) {
01586 int jp = iy >= 0 ? iy+1 : ny+iy+1;
01587 for (int ix = 0; ix <= nxc; ix++) {
01588 int r2 = ix*ix+iy*iy;
01589 if (( 4*r2 < ny*ny ) && !( ix == 0 && iy < 0 ) ) {
01590 float ctf = ctf_store::get_ctf( r2 )*10.f;
01591 float xnew = ix*tf[0][0] + iy*tf[1][0];
01592 float ynew = ix*tf[0][1] + iy*tf[1][1];
01593 float znew = ix*tf[0][2] + iy*tf[1][2];
01594 std::complex<float> btq;
01595 if (xnew < 0.0) {
01596 xnew = -xnew;
01597 ynew = -ynew;
01598 znew = -znew;
01599 btq = conj(myfft->cmplx(ix,jp));
01600 } else {
01601 btq = myfft->cmplx(ix,jp);
01602 }
01603 int ixn = int(xnew + 0.5 + nx) - nx;
01604 int iyn = int(ynew + 0.5 + ny) - ny;
01605 int izn = int(znew + 0.5 + nz) - nz;
01606 if ((ixn <= nxc) && (iyn >= iymin) && (iyn <= iymax) && (izn >= izmin) && (izn <= izmax)) {
01607 if (ixn >= 0) {
01608 int iza, iya;
01609 if (izn >= 0) iza = izn + 1;
01610 else iza = nz + izn + 1;
01611
01612 if (iyn >= 0) iya = iyn + 1;
01613 else iya = ny + iyn + 1;
01614
01615 cmplx(ixn,iya,iza) += btq*ctf;
01616 (*wptr)(ixn,iya,iza) += ctf*ctf;
01617 (*wptr2)(ixn,iya,iza) += std::norm(btq);
01618 (*wptr3)(ixn,iya,iza) += 1;
01619 } else {
01620 int izt, iyt;
01621 if (izn > 0) izt = nz - izn + 1;
01622 else izt = -izn + 1;
01623
01624 if (iyn > 0) iyt = ny - iyn + 1;
01625 else iyt = -iyn + 1;
01626
01627 cmplx(-ixn,iyt,izt) += std::conj(btq)*ctf;
01628 (*wptr) (-ixn,iyt,izt) += ctf*ctf;
01629 (*wptr2)(-ixn,iyt,izt) += std::norm(btq);
01630 (*wptr3)(-ixn,iyt,izt) += 1;
01631 }
01632 }
01633 }
01634 }
01635 }
01636 set_array_offsets(saved_offsets);
01637 myfft->set_array_offsets(myfft_saved_offsets);
01638 if(ctf) {delete ctf; ctf=0;}
01639 EXITFUNC;
01640 }
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740 void EMData::symplane0_ctf(EMData* w) {
01741 ENTERFUNC;
01742 int nxc = attr_dict["nxc"];
01743 int n = nxc*2;
01744
01745 vector<int> saved_offsets = get_array_offsets();
01746 set_array_offsets(0,1,1);
01747 for (int iza = 2; iza <= nxc; iza++) {
01748 for (int iya = 2; iya <= nxc; iya++) {
01749 cmplx(0,iya,iza) += conj(cmplx(0,n-iya+2,n-iza+2));
01750 (*w)(0,iya,iza) += (*w)(0,n-iya+2,n-iza+2);
01751 cmplx(0,n-iya+2,n-iza+2) = conj(cmplx(0,iya,iza));
01752 (*w)(0,n-iya+2,n-iza+2) = (*w)(0,iya,iza);
01753 cmplx(0,n-iya+2,iza) += conj(cmplx(0,iya,n-iza+2));
01754 (*w)(0,n-iya+2,iza) += (*w)(0,iya,n-iza+2);
01755 cmplx(0,iya,n-iza+2) = conj(cmplx(0,n-iya+2,iza));
01756 (*w)(0,iya,n-iza+2) = (*w)(0,n-iya+2,iza);
01757 }
01758 }
01759 for (int iya = 2; iya <= nxc; iya++) {
01760 cmplx(0,iya,1) += conj(cmplx(0,n-iya+2,1));
01761 (*w)(0,iya,1) += (*w)(0,n-iya+2,1);
01762 cmplx(0,n-iya+2,1) = conj(cmplx(0,iya,1));
01763 (*w)(0,n-iya+2,1) = (*w)(0,iya,1);
01764 }
01765 for (int iza = 2; iza <= nxc; iza++) {
01766 cmplx(0,1,iza) += conj(cmplx(0,1,n-iza+2));
01767 (*w)(0,1,iza) += (*w)(0,1,n-iza+2);
01768 cmplx(0,1,n-iza+2) = conj(cmplx(0,1,iza));
01769 (*w)(0,1,n-iza+2) = (*w)(0,1,iza);
01770 }
01771 EXITFUNC;
01772 }
01773
01774 EMData* EMData::rot_scale_trans2D(float angDeg, float delx, float dely, float scale) {
01775 float ang=angDeg*M_PI/180.0f;
01776 if (1 >= ny)
01777 throw ImageDimensionException("Can't rotate 1D image");
01778 if (nz<2) {
01779 vector<int> saved_offsets = get_array_offsets();
01780 set_array_offsets(0,0,0);
01781 if (0.0f == scale) scale = 1.0f;
01782 EMData* ret = copy_head();
01783 delx = restrict2(delx, nx);
01784 dely = restrict2(dely, ny);
01785
01786 int xc = nx/2;
01787 int yc = ny/2;
01788
01789 float shiftxc = xc + delx;
01790 float shiftyc = yc + dely;
01791
01792 float cang = cos(ang);
01793 float sang = sin(ang);
01794 for (int iy = 0; iy < ny; iy++) {
01795 float y = float(iy) - shiftyc;
01796 float ycang = y*cang/scale + yc;
01797 float ysang = -y*sang/scale + xc;
01798 for (int ix = 0; ix < nx; ix++) {
01799 float x = float(ix) - shiftxc;
01800 float xold = x*cang/scale + ysang ;
01801 float yold = x*sang/scale + ycang ;
01802
01803 (*ret)(ix,iy) = Util::quadri(xold+1.0f, yold+1.0f, nx, ny, get_data());
01804
01805 }
01806 }
01807 set_array_offsets(saved_offsets);
01808 return ret;
01809 } else {
01810 throw ImageDimensionException("Volume not currently supported");
01811 }
01812 }
01813
01814 EMData* EMData::rot_scale_trans2D_background(float angDeg, float delx, float dely, float scale) {
01815 float ang=angDeg*M_PI/180.0f;
01816 if (1 >= ny)
01817 throw ImageDimensionException("Can't rotate 1D image");
01818 if (nz<2) {
01819 vector<int> saved_offsets = get_array_offsets();
01820 set_array_offsets(0,0,0);
01821 if (0.0f == scale) scale = 1.0f;
01822 EMData* ret = copy_head();
01823 delx = restrict2(delx, nx);
01824 dely = restrict2(dely, ny);
01825
01826 int xc = nx/2;
01827 int yc = ny/2;
01828
01829 float shiftxc = xc + delx;
01830 float shiftyc = yc + dely;
01831
01832 float cang = cos(ang);
01833 float sang = sin(ang);
01834 for (int iy = 0; iy < ny; iy++) {
01835 float y = float(iy) - shiftyc;
01836 float ycang = y*cang/scale + yc;
01837 float ysang = -y*sang/scale + xc;
01838 for (int ix = 0; ix < nx; ix++) {
01839 float x = float(ix) - shiftxc;
01840 float xold = x*cang/scale + ysang ;
01841 float yold = x*sang/scale + ycang ;
01842
01843 (*ret)(ix,iy) = Util::quadri_background(xold+1.0f, yold+1.0f, nx, ny, get_data(),ix+1,iy+1);
01844
01845 }
01846 }
01847 set_array_offsets(saved_offsets);
01848 return ret;
01849 } else {
01850 throw ImageDimensionException("Volume not currently supported");
01851 }
01852 }
01853
01854 #define in(i,j,k) in[i+(j+(k*ny))*nx]
01855 EMData*
01856 EMData::rot_scale_trans(const Transform &RA) {
01857
01858 EMData* ret = copy_head();
01859 float *in = this->get_data();
01860 vector<int> saved_offsets = get_array_offsets();
01861 set_array_offsets(0,0,0);
01862 Vec3f translations = RA.get_trans();
01863 Transform RAinv = RA.inverse();
01864
01865 if (1 >= ny) throw ImageDimensionException("Can't rotate 1D image");
01866 if (nz < 2) {
01867 float p1, p2, p3, p4;
01868 float delx = translations.at(0);
01869 float dely = translations.at(1);
01870 delx = restrict2(delx, nx);
01871 dely = restrict2(dely, ny);
01872 int xc = nx/2;
01873 int yc = ny/2;
01874
01875 float shiftxc = xc + delx;
01876 float shiftyc = yc + dely;
01877 for (int iy = 0; iy < ny; iy++) {
01878 float y = float(iy) - shiftyc;
01879 float ysang = y*RAinv[0][1]+xc;
01880 float ycang = y*RAinv[1][1]+yc;
01881 for (int ix = 0; ix < nx; ix++) {
01882 float x = float(ix) - shiftxc;
01883 float xold = x*RAinv[0][0] + ysang;
01884 float yold = x*RAinv[1][0] + ycang;
01885
01886 xold = restrict1(xold, nx);
01887 yold = restrict1(yold, ny);
01888
01889 int xfloor = int(xold);
01890 int yfloor = int(yold);
01891 float t = xold-xfloor;
01892 float u = yold-yfloor;
01893 if(xfloor == nx -1 && yfloor == ny -1) {
01894
01895 p1 =in[xfloor + yfloor*ny];
01896 p2 =in[ yfloor*ny];
01897 p3 =in[0];
01898 p4 =in[xfloor];
01899 } else if(xfloor == nx - 1) {
01900
01901 p1 =in[xfloor + yfloor*ny];
01902 p2 =in[ yfloor*ny];
01903 p3 =in[ (yfloor+1)*ny];
01904 p4 =in[xfloor + (yfloor+1)*ny];
01905 } else if(yfloor == ny - 1) {
01906
01907 p1 =in[xfloor + yfloor*ny];
01908 p2 =in[xfloor+1 + yfloor*ny];
01909 p3 =in[xfloor+1 ];
01910 p4 =in[xfloor ];
01911 } else {
01912 p1 =in[xfloor + yfloor*ny];
01913 p2 =in[xfloor+1 + yfloor*ny];
01914 p3 =in[xfloor+1 + (yfloor+1)*ny];
01915 p4 =in[xfloor + (yfloor+1)*ny];
01916 }
01917 (*ret)(ix,iy) = p1 + u * ( p4 - p1) + t * ( p2 - p1 + u *(p3-p2-p4+p1));
01918 }
01919 }
01920 set_array_offsets(saved_offsets);
01921 return ret;
01922 } else {
01923
01924
01925 float delx = translations.at(0);
01926 float dely = translations.at(1);
01927 float delz = translations.at(2);
01928 delx = restrict2(delx, nx);
01929 dely = restrict2(dely, ny);
01930 delz = restrict2(delz, nz);
01931 int xc = nx/2;
01932 int yc = ny/2;
01933 int zc = nz/2;
01934
01935 float shiftxc = xc + delx;
01936 float shiftyc = yc + dely;
01937 float shiftzc = zc + delz;
01938
01939 for (int iz = 0; iz < nz; iz++) {
01940 float z = float(iz) - shiftzc;
01941 float xoldz = z*RAinv[0][2]+xc;
01942 float yoldz = z*RAinv[1][2]+yc;
01943 float zoldz = z*RAinv[2][2]+zc;
01944 for (int iy = 0; iy < ny; iy++) {
01945 float y = float(iy) - shiftyc;
01946 float xoldzy = xoldz + y*RAinv[0][1] ;
01947 float yoldzy = yoldz + y*RAinv[1][1] ;
01948 float zoldzy = zoldz + y*RAinv[2][1] ;
01949 for (int ix = 0; ix < nx; ix++) {
01950 float x = float(ix) - shiftxc;
01951 float xold = xoldzy + x*RAinv[0][0] ;
01952 float yold = yoldzy + x*RAinv[1][0] ;
01953 float zold = zoldzy + x*RAinv[2][0] ;
01954
01955 xold = restrict1(xold, nx);
01956 yold = restrict1(yold, ny);
01957 zold = restrict1(zold, nz);
01958
01959
01960 int IOX = int(xold);
01961 int IOY = int(yold);
01962 int IOZ = int(zold);
01963
01964 #ifdef _WIN32
01965 int IOXp1 = _MIN( nx-1 ,IOX+1);
01966 #else
01967 int IOXp1 = std::min( nx-1 ,IOX+1);
01968 #endif //_WIN32
01969
01970 #ifdef _WIN32
01971 int IOYp1 = _MIN( ny-1 ,IOY+1);
01972 #else
01973 int IOYp1 = std::min( ny-1 ,IOY+1);
01974 #endif //_WIN32
01975
01976 #ifdef _WIN32
01977 int IOZp1 = _MIN( nz-1 ,IOZ+1);
01978 #else
01979 int IOZp1 = std::min( nz-1 ,IOZ+1);
01980 #endif //_WIN32
01981
01982 float dx = xold-IOX;
01983 float dy = yold-IOY;
01984 float dz = zold-IOZ;
01985
01986 float a1 = in(IOX,IOY,IOZ);
01987 float a2 = in(IOXp1,IOY,IOZ) - in(IOX,IOY,IOZ);
01988 float a3 = in(IOX,IOYp1,IOZ) - in(IOX,IOY,IOZ);
01989 float a4 = in(IOX,IOY,IOZp1) - in(IOX,IOY,IOZ);
01990 float a5 = in(IOX,IOY,IOZ) - in(IOXp1,IOY,IOZ) - in(IOX,IOYp1,IOZ) + in(IOXp1,IOYp1,IOZ);
01991 float a6 = in(IOX,IOY,IOZ) - in(IOXp1,IOY,IOZ) - in(IOX,IOY,IOZp1) + in(IOXp1,IOY,IOZp1);
01992 float a7 = in(IOX,IOY,IOZ) - in(IOX,IOYp1,IOZ) - in(IOX,IOY,IOZp1) + in(IOX,IOYp1,IOZp1);
01993 float a8 = in(IOXp1,IOY,IOZ) + in(IOX,IOYp1,IOZ)+ in(IOX,IOY,IOZp1)
01994 - in(IOX,IOY,IOZ)- in(IOXp1,IOYp1,IOZ) - in(IOXp1,IOY,IOZp1)
01995 - in(IOX,IOYp1,IOZp1) + in(IOXp1,IOYp1,IOZp1);
01996 (*ret)(ix,iy,iz) = a1 + dz*(a4 + a6*dx + (a7 + a8*dx)*dy) + a3*dy + dx*(a2 + a5*dy);
01997 }
01998 }
01999 }
02000
02001 set_array_offsets(saved_offsets);
02002 return ret;
02003
02004
02005
02006
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032
02033
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
02107
02108
02109
02110
02111
02112
02113
02114
02115
02116
02117
02118 }
02119 }
02120 #undef in
02121
02122
02123 #define in(i,j,k) in[i+(j+(k*ny))*nx]
02124 EMData*
02125 EMData::rot_scale_trans_background(const Transform &RA) {
02126 EMData* ret = copy_head();
02127 float *in = this->get_data();
02128 vector<int> saved_offsets = get_array_offsets();
02129 set_array_offsets(0,0,0);
02130 Vec3f translations = RA.get_trans();
02131 Transform RAinv = RA.inverse();
02132
02133 if (1 >= ny) throw ImageDimensionException("Can't rotate 1D image");
02134 if (nz < 2) {
02135 float p1, p2, p3, p4;
02136 float delx = translations.at(0);
02137 float dely = translations.at(1);
02138 delx = restrict2(delx, nx);
02139 dely = restrict2(dely, ny);
02140 int xc = nx/2;
02141 int yc = ny/2;
02142
02143 float shiftxc = xc + delx;
02144 float shiftyc = yc + dely;
02145 for (int iy = 0; iy < ny; iy++) {
02146 float y = float(iy) - shiftyc;
02147 float ysang = y*RAinv[0][1]+xc;
02148 float ycang = y*RAinv[1][1]+yc;
02149 for (int ix = 0; ix < nx; ix++) {
02150 float x = float(ix) - shiftxc;
02151 float xold = x*RAinv[0][0] + ysang;
02152 float yold = x*RAinv[1][0] + ycang;
02153
02154
02155
02156 if ( (xold < 0.0f) || (xold >= (float)(nx)) || (yold < 0.0f) || (yold >= (float)(ny)) ){
02157 xold = (float)ix;
02158 yold = (float)iy;
02159 }
02160
02161 int xfloor = int(xold);
02162 int yfloor = int(yold);
02163 float t = xold-xfloor;
02164 float u = yold-yfloor;
02165 if(xfloor == nx -1 && yfloor == ny -1) {
02166
02167 p1 =in[xfloor + yfloor*ny];
02168 p2 =in[ yfloor*ny];
02169 p3 =in[0];
02170 p4 =in[xfloor];
02171 } else if(xfloor == nx - 1) {
02172
02173 p1 =in[xfloor + yfloor*ny];
02174 p2 =in[ yfloor*ny];
02175 p3 =in[ (yfloor+1)*ny];
02176 p4 =in[xfloor + (yfloor+1)*ny];
02177 } else if(yfloor == ny - 1) {
02178
02179 p1 =in[xfloor + yfloor*ny];
02180 p2 =in[xfloor+1 + yfloor*ny];
02181 p3 =in[xfloor+1 ];
02182 p4 =in[xfloor ];
02183 } else {
02184
02185 p1 =in[xfloor + yfloor*ny];
02186 p2 =in[xfloor+1 + yfloor*ny];
02187 p3 =in[xfloor+1 + (yfloor+1)*ny];
02188 p4 =in[xfloor + (yfloor+1)*ny];
02189 }
02190 (*ret)(ix,iy) = p1 + u * ( p4 - p1) + t * ( p2 - p1 + u *(p3-p2-p4+p1));
02191 }
02192 }
02193 set_array_offsets(saved_offsets);
02194 return ret;
02195 } else {
02196
02197
02198 float delx = translations.at(0);
02199 float dely = translations.at(1);
02200 float delz = translations.at(2);
02201 delx = restrict2(delx, nx);
02202 dely = restrict2(dely, ny);
02203 delz = restrict2(delz, nz);
02204 int xc = nx/2;
02205 int yc = ny/2;
02206 int zc = nz/2;
02207
02208 float shiftxc = xc + delx;
02209 float shiftyc = yc + dely;
02210 float shiftzc = zc + delz;
02211
02212 for (int iz = 0; iz < nz; iz++) {
02213 float z = float(iz) - shiftzc;
02214 float xoldz = z*RAinv[0][2]+xc;
02215 float yoldz = z*RAinv[1][2]+yc;
02216 float zoldz = z*RAinv[2][2]+zc;
02217 for (int iy = 0; iy < ny; iy++) {
02218 float y = float(iy) - shiftyc;
02219 float xoldzy = xoldz + y*RAinv[0][1] ;
02220 float yoldzy = yoldz + y*RAinv[1][1] ;
02221 float zoldzy = zoldz + y*RAinv[2][1] ;
02222 for (int ix = 0; ix < nx; ix++) {
02223 float x = float(ix) - shiftxc;
02224 float xold = xoldzy + x*RAinv[0][0] ;
02225 float yold = yoldzy + x*RAinv[1][0] ;
02226 float zold = zoldzy + x*RAinv[2][0] ;
02227
02228
02229
02230 if ( (xold < 0.0f) || (xold >= (float)(nx)) || (yold < 0.0f) || (yold >= (float)(ny)) || (zold < 0.0f) || (zold >= (float)(nz)) ){
02231 xold = (float)ix;
02232 yold = (float)iy;
02233 zold = (float)iz;
02234 }
02235
02236 int IOX = int(xold);
02237 int IOY = int(yold);
02238 int IOZ = int(zold);
02239
02240 #ifdef _WIN32
02241 int IOXp1 = _MIN( nx-1 ,IOX+1);
02242 #else
02243 int IOXp1 = std::min( nx-1 ,IOX+1);
02244 #endif //_WIN32
02245
02246 #ifdef _WIN32
02247 int IOYp1 = _MIN( ny-1 ,IOY+1);
02248 #else
02249 int IOYp1 = std::min( ny-1 ,IOY+1);
02250 #endif //_WIN32
02251
02252 #ifdef _WIN32
02253 int IOZp1 = _MIN( nz-1 ,IOZ+1);
02254 #else
02255 int IOZp1 = std::min( nz-1 ,IOZ+1);
02256 #endif //_WIN32
02257
02258 float dx = xold-IOX;
02259 float dy = yold-IOY;
02260 float dz = zold-IOZ;
02261
02262 float a1 = in(IOX,IOY,IOZ);
02263 float a2 = in(IOXp1,IOY,IOZ) - in(IOX,IOY,IOZ);
02264 float a3 = in(IOX,IOYp1,IOZ) - in(IOX,IOY,IOZ);
02265 float a4 = in(IOX,IOY,IOZp1) - in(IOX,IOY,IOZ);
02266 float a5 = in(IOX,IOY,IOZ) - in(IOXp1,IOY,IOZ) - in(IOX,IOYp1,IOZ) + in(IOXp1,IOYp1,IOZ);
02267 float a6 = in(IOX,IOY,IOZ) - in(IOXp1,IOY,IOZ) - in(IOX,IOY,IOZp1) + in(IOXp1,IOY,IOZp1);
02268 float a7 = in(IOX,IOY,IOZ) - in(IOX,IOYp1,IOZ) - in(IOX,IOY,IOZp1) + in(IOX,IOYp1,IOZp1);
02269 float a8 = in(IOXp1,IOY,IOZ) + in(IOX,IOYp1,IOZ)+ in(IOX,IOY,IOZp1)
02270 - in(IOX,IOY,IOZ)- in(IOXp1,IOYp1,IOZ) - in(IOXp1,IOY,IOZp1)
02271 - in(IOX,IOYp1,IOZp1) + in(IOXp1,IOYp1,IOZp1);
02272 (*ret)(ix,iy,iz) = a1 + dz*(a4 + a6*dx + (a7 + a8*dx)*dy) + a3*dy + dx*(a2 + a5*dy);
02273 }
02274 }
02275 }
02276
02277 set_array_offsets(saved_offsets);
02278 return ret;
02279
02280 }
02281 }
02282 #undef in
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314
02315
02316
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339
02340
02341
02342
02343
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353
02354
02355
02356
02357
02358
02359
02360 EMData* EMData::rot_scale_conv(float ang, float delx, float dely, Util::KaiserBessel& kb, float scale_input) {
02361 int nxn, nyn, nzn;
02362 if(scale_input == 0.0f) scale_input = 1.0f;
02363
02364 float scale = 0.5f*scale_input;
02365 float sum, w;
02366 if (1 >= ny)
02367 throw ImageDimensionException("Can't rotate 1D image");
02368 if (1 < nz)
02369 throw ImageDimensionException("Volume not currently supported");
02370 nxn=nx/2;nyn=ny/2;nzn=nz/2;
02371
02372 int K = kb.get_window_size();
02373 int kbmin = -K/2;
02374 int kbmax = -kbmin;
02375 int kbc = kbmax+1;
02376 vector<int> saved_offsets = get_array_offsets();
02377 set_array_offsets(0,0,0);
02378 EMData* ret = this->copy_head();
02379 #ifdef _WIN32
02380 ret->set_size(nxn, _MAX(nyn,1), _MAX(nzn,1));
02381 #else
02382 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
02383 #endif //_WIN32
02384
02385 delx = restrict2(delx, nx);
02386 dely = restrict2(dely, ny);
02387
02388 int xc = nxn;
02389 int ixs = nxn%2;
02390 int yc = nyn;
02391 int iys = nyn%2;
02392
02393 int xcn = nxn/2;
02394 int ycn = nyn/2;
02395
02396 float shiftxc = xcn + delx;
02397 float shiftyc = ycn + dely;
02398
02399 float ymin = -ny/2.0f;
02400 float xmin = -nx/2.0f;
02401 float ymax = -ymin;
02402 float xmax = -xmin;
02403 if (0 == nx%2) xmax--;
02404 if (0 == ny%2) ymax--;
02405
02406 float *t = (float*)calloc(kbmax-kbmin+1, sizeof(float));
02407
02408
02409 float cang = cos(ang);
02410 float sang = sin(ang);
02411 for (int iy = 0; iy < nyn; iy++) {
02412 float y = float(iy) - shiftyc;
02413 float ycang = y*cang/scale + yc;
02414 float ysang = -y*sang/scale + xc;
02415 for (int ix = 0; ix < nxn; ix++) {
02416 float x = float(ix) - shiftxc;
02417 float xold = x*cang/scale + ysang-ixs;
02418 float yold = x*sang/scale + ycang-iys;
02419
02420 xold = restrict1(xold, nx);
02421 yold = restrict1(yold, ny);
02422
02423 int inxold = int(Util::round(xold)); int inyold = int(Util::round(yold));
02424 sum=0.0f; w=0.0f;
02425 for (int m1 =kbmin; m1 <=kbmax; m1++) t[m1-kbmin] = kb.i0win_tab(xold - inxold-m1);
02426 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 ) {
02427 for (int m2 =kbmin; m2 <=kbmax; m2++) {
02428 float qt = kb.i0win_tab(yold - inyold-m2);
02429 for (int m1 =kbmin; m1 <=kbmax; m1++) {
02430 float q = t[m1-kbmin]*qt;
02431 sum += (*this)((inxold+m1+nx)%nx,(inyold+m2+ny)%ny)*q; w+=q;
02432 }
02433 }
02434 } else {
02435 for (int m2 =kbmin; m2 <=kbmax; m2++) {
02436 float qt = kb.i0win_tab(yold - inyold-m2);
02437 for (int m1 =kbmin; m1 <=kbmax; m1++) {
02438 float q = t[m1-kbmin]*qt;
02439 sum += (*this)(inxold+m1,inyold+m2)*q; w+=q;}
02440 }
02441 }
02442 (*ret)(ix,iy)=sum/w;
02443 }
02444 }
02445 if (t) free(t);
02446 set_array_offsets(saved_offsets);
02447 return ret;
02448 }
02449
02450
02451
02452 EMData* EMData::rot_scale_conv7(float ang, float delx, float dely, Util::KaiserBessel& kb, float scale_input) {
02453 int nxn, nyn, nzn;
02454 float scale = 0.5f*scale_input;
02455 float sum, w;
02456 if (1 >= ny)
02457 throw ImageDimensionException("Can't rotate 1D image");
02458 if (1 < nz)
02459 throw ImageDimensionException("Volume not currently supported");
02460 nxn = nx/2; nyn=ny/2; nzn=nz/2;
02461
02462 int K = kb.get_window_size();
02463 int kbmin = -K/2;
02464 int kbmax = -kbmin;
02465 int kbc = kbmax+1;
02466 vector<int> saved_offsets = get_array_offsets();
02467 set_array_offsets(0,0,0);
02468 EMData* ret = this->copy_head();
02469 #ifdef _WIN32
02470 ret->set_size(nxn, _MAX(nyn,1), _MAX(nzn,1));
02471 #else
02472 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
02473 #endif //_WIN32
02474
02475 delx = restrict2(delx, nx);
02476 dely = restrict2(dely, ny);
02477
02478 int xc = nxn;
02479 int ixs = nxn%2;
02480 int yc = nyn;
02481 int iys = nyn%2;
02482
02483 int xcn = nxn/2;
02484 int ycn = nyn/2;
02485
02486 float shiftxc = xcn + delx;
02487 float shiftyc = ycn + dely;
02488
02489 float ymin = -ny/2.0f;
02490 float xmin = -nx/2.0f;
02491 float ymax = -ymin;
02492 float xmax = -xmin;
02493 if (0 == nx%2) xmax--;
02494 if (0 == ny%2) ymax--;
02495
02496 float *t = (float*)calloc(kbmax-kbmin+1, sizeof(float));
02497
02498
02499 float cang = cos(ang);
02500 float sang = sin(ang);
02501 for (int iy = 0; iy < nyn; iy++) {
02502 float y = float(iy) - shiftyc;
02503 float ycang = y*cang/scale + yc;
02504 float ysang = -y*sang/scale + xc;
02505 for (int ix = 0; ix < nxn; ix++) {
02506 float x = float(ix) - shiftxc;
02507 float xold = x*cang/scale + ysang-ixs;
02508 float yold = x*sang/scale + ycang-iys;
02509
02510 xold = restrict1(xold, nx);
02511 yold = restrict1(yold, ny);
02512
02513 int inxold = int(Util::round(xold)); int inyold = int(Util::round(yold));
02514 sum=0.0f; w=0.0f;
02515
02516 float tablex1 = kb.i0win_tab(xold-inxold+3);
02517 float tablex2 = kb.i0win_tab(xold-inxold+2);
02518 float tablex3 = kb.i0win_tab(xold-inxold+1);
02519 float tablex4 = kb.i0win_tab(xold-inxold);
02520 float tablex5 = kb.i0win_tab(xold-inxold-1);
02521 float tablex6 = kb.i0win_tab(xold-inxold-2);
02522 float tablex7 = kb.i0win_tab(xold-inxold-3);
02523
02524 float tabley1 = kb.i0win_tab(yold-inyold+3);
02525 float tabley2 = kb.i0win_tab(yold-inyold+2);
02526 float tabley3 = kb.i0win_tab(yold-inyold+1);
02527 float tabley4 = kb.i0win_tab(yold-inyold);
02528 float tabley5 = kb.i0win_tab(yold-inyold-1);
02529 float tabley6 = kb.i0win_tab(yold-inyold-2);
02530 float tabley7 = kb.i0win_tab(yold-inyold-3);
02531
02532 int x1, x2, x3, x4, x5, x6, x7, y1, y2, y3, y4, y5, y6, y7;
02533
02534 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 ) {
02535 x1 = (inxold-3+nx)%nx;
02536 x2 = (inxold-2+nx)%nx;
02537 x3 = (inxold-1+nx)%nx;
02538 x4 = (inxold +nx)%nx;
02539 x5 = (inxold+1+nx)%nx;
02540 x6 = (inxold+2+nx)%nx;
02541 x7 = (inxold+3+nx)%nx;
02542
02543 y1 = (inyold-3+ny)%ny;
02544 y2 = (inyold-2+ny)%ny;
02545 y3 = (inyold-1+ny)%ny;
02546 y4 = (inyold +ny)%ny;
02547 y5 = (inyold+1+ny)%ny;
02548 y6 = (inyold+2+ny)%ny;
02549 y7 = (inyold+3+ny)%ny;
02550 } else {
02551 x1 = inxold-3;
02552 x2 = inxold-2;
02553 x3 = inxold-1;
02554 x4 = inxold;
02555 x5 = inxold+1;
02556 x6 = inxold+2;
02557 x7 = inxold+3;
02558
02559 y1 = inyold-3;
02560 y2 = inyold-2;
02561 y3 = inyold-1;
02562 y4 = inyold;
02563 y5 = inyold+1;
02564 y6 = inyold+2;
02565 y7 = inyold+3;
02566 }
02567 sum = ( (*this)(x1,y1)*tablex1 + (*this)(x2,y1)*tablex2 + (*this)(x3,y1)*tablex3 +
02568 (*this)(x4,y1)*tablex4 + (*this)(x5,y1)*tablex5 + (*this)(x6,y1)*tablex6 +
02569 (*this)(x7,y1)*tablex7 ) * tabley1 +
02570 ( (*this)(x1,y2)*tablex1 + (*this)(x2,y2)*tablex2 + (*this)(x3,y2)*tablex3 +
02571 (*this)(x4,y2)*tablex4 + (*this)(x5,y2)*tablex5 + (*this)(x6,y2)*tablex6 +
02572 (*this)(x7,y2)*tablex7 ) * tabley2 +
02573 ( (*this)(x1,y3)*tablex1 + (*this)(x2,y3)*tablex2 + (*this)(x3,y3)*tablex3 +
02574 (*this)(x4,y3)*tablex4 + (*this)(x5,y3)*tablex5 + (*this)(x6,y3)*tablex6 +
02575 (*this)(x7,y3)*tablex7 ) * tabley3 +
02576 ( (*this)(x1,y4)*tablex1 + (*this)(x2,y4)*tablex2 + (*this)(x3,y4)*tablex3 +
02577 (*this)(x4,y4)*tablex4 + (*this)(x5,y4)*tablex5 + (*this)(x6,y4)*tablex6 +
02578 (*this)(x7,y4)*tablex7 ) * tabley4 +
02579 ( (*this)(x1,y5)*tablex1 + (*this)(x2,y5)*tablex2 + (*this)(x3,y5)*tablex3 +
02580 (*this)(x4,y5)*tablex4 + (*this)(x5,y5)*tablex5 + (*this)(x6,y5)*tablex6 +
02581 (*this)(x7,y5)*tablex7 ) * tabley5 +
02582 ( (*this)(x1,y6)*tablex1 + (*this)(x2,y6)*tablex2 + (*this)(x3,y6)*tablex3 +
02583 (*this)(x4,y6)*tablex4 + (*this)(x5,y6)*tablex5 + (*this)(x6,y6)*tablex6 +
02584 (*this)(x7,y6)*tablex7 ) * tabley6 +
02585 ( (*this)(x1,y7)*tablex1 + (*this)(x2,y7)*tablex2 + (*this)(x3,y7)*tablex3 +
02586 (*this)(x4,y7)*tablex4 + (*this)(x5,y7)*tablex5 + (*this)(x6,y7)*tablex6 +
02587 (*this)(x7,y7)*tablex7 ) * tabley7;
02588
02589 w = (tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7) *
02590 (tabley1+tabley2+tabley3+tabley4+tabley5+tabley6+tabley7);
02591
02592 (*ret)(ix,iy)=sum/w;
02593 }
02594 }
02595 if (t) free(t);
02596 set_array_offsets(saved_offsets);
02597 return ret;
02598 }
02599
02600 EMData* EMData::downsample(Util::sincBlackman& kb, float scale) {
02601
02602
02603
02604
02605
02606 int nxn, nyn, nzn;
02607 nxn = (int)(nx*scale); nyn = (int)(ny*scale); nzn = (int)(nz*scale);
02608
02609 vector<int> saved_offsets = get_array_offsets();
02610 set_array_offsets(0,0,0);
02611 EMData* ret = this->copy_head();
02612 #ifdef _WIN32
02613 ret->set_size(nxn, _MAX(nyn,1), _MAX(nzn,1));
02614 #else
02615 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
02616 #endif //_WIN32
02617 ret->to_zero();
02618
02619
02620 for (int iy =0; iy < nyn; iy++) {
02621 float y = float(iy)/scale;
02622 for (int ix = 0; ix < nxn; ix++) {
02623 float x = float(ix)/scale;
02624 (*ret)(ix,iy) = this->get_pixel_filtered(x, y, 1.0f, kb);
02625 }
02626 }
02627 set_array_offsets(saved_offsets);
02628 return ret;
02629 }
02630
02631
02632 EMData* EMData::rot_scale_conv_new(float ang, float delx, float dely, Util::KaiserBessel& kb, float scale_input) {
02633
02634 int nxn, nyn, nzn;
02635
02636 if (scale_input == 0.0f) scale_input = 1.0f;
02637 float scale = 0.5f*scale_input;
02638
02639 if (1 >= ny)
02640 throw ImageDimensionException("Can't rotate 1D image");
02641 if (1 < nz)
02642 throw ImageDimensionException("Volume not currently supported");
02643 nxn = nx/2; nyn = ny/2; nzn = nz/2;
02644
02645 int K = kb.get_window_size();
02646 int kbmin = -K/2;
02647 int kbmax = -kbmin;
02648
02649 vector<int> saved_offsets = get_array_offsets();
02650 set_array_offsets(0,0,0);
02651 EMData* ret = this->copy_head();
02652 #ifdef _WIN32
02653 ret->set_size(nxn, _MAX(nyn,1), _MAX(nzn,1));
02654 #else
02655 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
02656 #endif //_WIN32
02657
02658 delx = restrict2(delx, nx);
02659 dely = restrict2(dely, ny);
02660
02661 int xc = nxn;
02662 int ixs = nxn%2;
02663 int yc = nyn;
02664 int iys = nyn%2;
02665
02666 int xcn = nxn/2;
02667 int ycn = nyn/2;
02668
02669 float shiftxc = xcn + delx;
02670 float shiftyc = ycn + dely;
02671
02672 float ymin = -ny/2.0f;
02673 float xmin = -nx/2.0f;
02674 float ymax = -ymin;
02675 float xmax = -xmin;
02676 if (0 == nx%2) xmax--;
02677 if (0 == ny%2) ymax--;
02678
02679 float *t = (float*)calloc(kbmax-kbmin+1, sizeof(float));
02680
02681 float* data = this->get_data();
02682
02683
02684 float cang = cos(ang);
02685 float sang = sin(ang);
02686 for (int iy = 0; iy < nyn; iy++) {
02687 float y = float(iy) - shiftyc;
02688 float ycang = y*cang/scale + yc;
02689 float ysang = -y*sang/scale + xc;
02690 for (int ix = 0; ix < nxn; ix++) {
02691 float x = float(ix) - shiftxc;
02692 float xold = x*cang/scale + ysang-ixs;
02693 float yold = x*sang/scale + ycang-iys;
02694
02695 (*ret)(ix,iy) = Util::get_pixel_conv_new(nx, ny, 1, xold, yold, 1, data, kb);
02696 }
02697 }
02698 if (t) free(t);
02699 set_array_offsets(saved_offsets);
02700 return ret;
02701 }
02702
02703 EMData* EMData::rot_scale_conv_new_background(float ang, float delx, float dely, Util::KaiserBessel& kb, float scale_input) {
02704
02705 int nxn, nyn, nzn;
02706
02707 if (scale_input == 0.0f) scale_input = 1.0f;
02708 float scale = 0.5f*scale_input;
02709
02710 if (1 >= ny)
02711 throw ImageDimensionException("Can't rotate 1D image");
02712 if (1 < nz)
02713 throw ImageDimensionException("Volume not currently supported");
02714 nxn = nx/2; nyn = ny/2; nzn = nz/2;
02715
02716 int K = kb.get_window_size();
02717 int kbmin = -K/2;
02718 int kbmax = -kbmin;
02719
02720 vector<int> saved_offsets = get_array_offsets();
02721 set_array_offsets(0,0,0);
02722 EMData* ret = this->copy_head();
02723 #ifdef _WIN32
02724 ret->set_size(nxn, _MAX(nyn,1), _MAX(nzn,1));
02725 #else
02726 ret->set_size(nxn, std::max(nyn,1), std::max(nzn,1));
02727 #endif //_WIN32
02728
02729 delx = restrict2(delx, nx);
02730 dely = restrict2(dely, ny);
02731
02732 int xc = nxn;
02733 int ixs = nxn%2;
02734 int yc = nyn;
02735 int iys = nyn%2;
02736
02737 int xcn = nxn/2;
02738 int ycn = nyn/2;
02739
02740 float shiftxc = xcn + delx;
02741 float shiftyc = ycn + dely;
02742
02743 float ymin = -ny/2.0f;
02744 float xmin = -nx/2.0f;
02745 float ymax = -ymin;
02746 float xmax = -xmin;
02747 if (0 == nx%2) xmax--;
02748 if (0 == ny%2) ymax--;
02749
02750 float *t = (float*)calloc(kbmax-kbmin+1, sizeof(float));
02751
02752 float* data = this->get_data();
02753
02754
02755 float cang = cos(ang);
02756 float sang = sin(ang);
02757 for (int iy = 0; iy < nyn; iy++) {
02758 float y = float(iy) - shiftyc;
02759 float ycang = y*cang/scale + yc;
02760 float ysang = -y*sang/scale + xc;
02761 for (int ix = 0; ix < nxn; ix++) {
02762 float x = float(ix) - shiftxc;
02763 float xold = x*cang/scale + ysang-ixs;
02764 float yold = x*sang/scale + ycang-iys;
02765
02766 (*ret)(ix,iy) = Util::get_pixel_conv_new_background(nx, ny, 1, xold, yold, 1, data, kb, ix,iy);
02767 }
02768 }
02769 if (t) free(t);
02770 set_array_offsets(saved_offsets);
02771 return ret;
02772 }
02773
02774 float EMData::get_pixel_conv(float delx, float dely, float delz, Util::KaiserBessel& kb) {
02775
02776
02777 int K = kb.get_window_size();
02778 int kbmin = -K/2;
02779 int kbmax = -kbmin;
02780 int kbc = kbmax+1;
02781
02782 float pixel =0.0f;
02783 float w=0.0f;
02784
02785 delx = restrict2(delx, nx);
02786 int inxold = int(Util::round(delx));
02787 if(ny<2) {
02788 if(inxold <= kbc || inxold >=nx-kbc-2 ) {
02789
02790 for (int m1 =kbmin; m1 <=kbmax; m1++) {
02791 float q = kb.i0win_tab(delx - inxold-m1);
02792 pixel += (*this)((inxold+m1+nx)%nx)*q; w+=q;
02793 }
02794 } else {
02795 for (int m1 =kbmin; m1 <=kbmax; m1++) {
02796 float q = kb.i0win_tab(delx - inxold-m1);
02797 pixel += (*this)(inxold+m1)*q; w+=q;
02798 }
02799 }
02800
02801 } else if(nz<2) {
02802 dely = restrict2(dely, ny);
02803 int inyold = int(Util::round(dely));
02804 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 ) {
02805
02806 for (int m2 =kbmin; m2 <=kbmax; m2++){ for (int m1 =kbmin; m1 <=kbmax; m1++) {
02807 float q = kb.i0win_tab(delx - inxold-m1)*kb.i0win_tab(dely - inyold-m2);
02808 pixel += (*this)((inxold+m1+nx)%nx,(inyold+m2+ny)%ny)*q; w+=q;}
02809 }
02810 } else {
02811 for (int m2 =kbmin; m2 <=kbmax; m2++){ for (int m1 =kbmin; m1 <=kbmax; m1++) {
02812 float q = kb.i0win_tab(delx - inxold-m1)*kb.i0win_tab(dely - inyold-m2);
02813 pixel += (*this)(inxold+m1,inyold+m2)*q; w+=q;}
02814 }
02815 }
02816 } else {
02817 dely = restrict2(dely, ny);
02818 int inyold = int(Util::round(dely));
02819 delz = restrict2(delz, nz);
02820 int inzold = int(Util::round(delz));
02821
02822 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 || inzold <= kbc || inzold >=nz-kbc-2 ) {
02823
02824 for (int m3 =kbmin; m3 <=kbmax; m3++){ for (int m2 =kbmin; m2 <=kbmax; m2++){ for (int m1 =kbmin; m1 <=kbmax; m1++) {
02825 float q = kb.i0win_tab(delx - inxold-m1)*kb.i0win_tab(dely - inyold-m2)*kb.i0win_tab(delz - inzold-m3);
02826
02827 pixel += (*this)((inxold+m1+nx)%nx,(inyold+m2+ny)%ny,(inzold+m3+nz)%nz)*q ;w+=q;}}
02828 }
02829 } else {
02830 for (int m3 =kbmin; m3 <=kbmax; m3++){ for (int m2 =kbmin; m2 <=kbmax; m2++){ for (int m1 =kbmin; m1 <=kbmax; m1++) {
02831 float q = kb.i0win_tab(delx - inxold-m1)*kb.i0win_tab(dely - inyold-m2)*kb.i0win_tab(delz - inzold-m3);
02832
02833 pixel += (*this)(inxold+m1,inyold+m2,inzold+m3)*q; w+=q;}}
02834 }
02835 }
02836 }
02837 return pixel/w;
02838 }
02839
02840
02841 float EMData::get_pixel_filtered(float delx, float dely, float delz, Util::sincBlackman& kb) {
02842
02843
02844 int K = kb.get_sB_size();
02845 int kbmin = -K/2;
02846 int kbmax = -kbmin;
02847 int kbc = kbmax+1;
02848
02849 float pixel =0.0f;
02850 float w=0.0f;
02851
02852
02853 int inxold = int(Util::round(delx));
02854
02855
02856
02857
02858
02859
02860
02861
02862
02863
02864
02865
02866
02867
02868
02869
02870 int inyold = int(Util::round(dely));
02871 if(inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 ) {
02872
02873 for (int m2 =kbmin; m2 <=kbmax; m2++){
02874 float t = kb.sBwin_tab(dely - inyold-m2);
02875 for (int m1 =kbmin; m1 <=kbmax; m1++) {
02876 float q = kb.sBwin_tab(delx - inxold-m1)*t;
02877 pixel += (*this)((inxold+m1+nx)%nx,(inyold+m2+ny)%ny)*q;
02878 w += q;
02879 }
02880 }
02881 } else {
02882 for (int m2 =kbmin; m2 <=kbmax; m2++){
02883 float t = kb.sBwin_tab(dely - inyold-m2);
02884 for (int m1 =kbmin; m1 <=kbmax; m1++) {
02885 float q = kb.sBwin_tab(delx - inxold-m1)*t;
02886 pixel += (*this)(inxold+m1,inyold+m2)*q;
02887 w += q;
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 return pixel/w;
02913 }
02914
02915
02916
02917
02918 float EMData::get_pixel_conv7(float delx, float dely, float delz, Util::KaiserBessel& kb) {
02919
02920
02921 float *image=(this->get_data());
02922 int nx = this->get_xsize();
02923 int ny = this->get_ysize();
02924 int nz = this->get_zsize();
02925
02926 float result;
02927
02928 result = Util::get_pixel_conv_new(nx,ny,nz,delx,dely,delz,image,kb);
02929 return result;
02930 }
02931
02932 float EMData::getconvpt2d_kbi0(float x, float y, Util::KaiserBessel::kbi0_win win, int size) {
02933 const int nxhalf = nx/2;
02934 const int nyhalf = ny/2;
02935 const int bd = size/2;
02936 float* wxarr = new float[size];
02937 float* wyarr = new float[size];
02938 float* wx = wxarr + bd;
02939 float* wy = wyarr + bd;
02940 int ixc = int(x + 0.5f*Util::sgn(x));
02941 int iyc = int(y + 0.5f*Util::sgn(y));
02942 if (abs(ixc) > nxhalf)
02943 throw InvalidValueException(ixc, "getconv: X value out of range");
02944 if (abs(iyc) > nyhalf)
02945 throw InvalidValueException(ixc, "getconv: Y value out of range");
02946 for (int i = -bd; i <= bd; i++) {
02947 int iyp = iyc + i;
02948 wy[i] = win(y - iyp);
02949 int ixp = ixc + i;
02950 wx[i] = win(x - ixp);
02951 }
02952 vector<int> saved_offsets = get_array_offsets();
02953 set_array_offsets(-nxhalf, -nyhalf);
02954 float conv = 0.f, wsum = 0.f;
02955 for (int iy = -bd; iy <= bd; iy++) {
02956 int iyp = iyc + iy;
02957 for (int ix = -bd; ix <= bd; ix++) {
02958 int ixp = ixc + ix;
02959 float wg = wx[ix]*wy[iy];
02960 conv += (*this)(ixp,iyp)*wg;
02961 wsum += wg;
02962 }
02963 }
02964 set_array_offsets(saved_offsets);
02965 delete [] wxarr;
02966 delete [] wyarr;
02967
02968 return conv;
02969 }
02970
02971 std::complex<float> EMData::extractpoint(float nuxnew, float nuynew, Util::KaiserBessel& kb) {
02972 if (2 != get_ndim())
02973 throw ImageDimensionException("extractpoint needs a 2-D image.");
02974 if (!is_complex())
02975 throw ImageFormatException("extractpoint requires a fourier image");
02976 int nxreal = nx - 2;
02977 if (nxreal != ny)
02978 throw ImageDimensionException("extractpoint requires ny == nx");
02979 int nhalf = nxreal/2;
02980 int kbsize = kb.get_window_size();
02981 int kbmin = -kbsize/2;
02982 int kbmax = -kbmin;
02983 bool flip = (nuxnew < 0.f);
02984 if (flip) {
02985 nuxnew *= -1;
02986 nuynew *= -1;
02987 }
02988
02989
02990
02991 int ixn = int(Util::round(nuxnew));
02992 int iyn = int(Util::round(nuynew));
02993
02994 float* wy0 = new float[kbmax - kbmin + 1];
02995 float* wy = wy0 - kbmin;
02996 float* wx0 = new float[kbmax - kbmin + 1];
02997 float* wx = wx0 - kbmin;
02998 for (int i = kbmin; i <= kbmax; i++) {
02999 int iyp = iyn + i;
03000 wy[i] = kb.i0win_tab(nuynew - iyp);
03001 int ixp = ixn + i;
03002 wx[i] = kb.i0win_tab(nuxnew - ixp);
03003 }
03004
03005 int iymin = 0;
03006 for (int iy = kbmin; iy <= -1; iy++) {
03007 if (wy[iy] != 0.f) {
03008 iymin = iy;
03009 break;
03010 }
03011 }
03012 int iymax = 0;
03013 for (int iy = kbmax; iy >= 1; iy--) {
03014 if (wy[iy] != 0.f) {
03015 iymax = iy;
03016 break;
03017 }
03018 }
03019 int ixmin = 0;
03020 for (int ix = kbmin; ix <= -1; ix++) {
03021 if (wx[ix] != 0.f) {
03022 ixmin = ix;
03023 break;
03024 }
03025 }
03026 int ixmax = 0;
03027 for (int ix = kbmax; ix >= 1; ix--) {
03028 if (wx[ix] != 0.f) {
03029 ixmax = ix;
03030 break;
03031 }
03032 }
03033 float wsum = 0.0f;
03034 for (int iy = iymin; iy <= iymax; iy++)
03035 for (int ix = ixmin; ix <= ixmax; ix++)
03036 wsum += wx[ix]*wy[iy];
03037 std::complex<float> result(0.f,0.f);
03038 if ((ixn >= -kbmin) && (ixn <= nhalf-1-kbmax) && (iyn >= -nhalf-kbmin) && (iyn <= nhalf-1-kbmax)) {
03039
03040 for (int iy = iymin; iy <= iymax; iy++) {
03041 int iyp = iyn + iy;
03042 for (int ix = ixmin; ix <= ixmax; ix++) {
03043 int ixp = ixn + ix;
03044 float w = wx[ix]*wy[iy];
03045 std::complex<float> val = cmplx(ixp,iyp);
03046 result += val*w;
03047 }
03048 }
03049 } else {
03050
03051 for (int iy = iymin; iy <= iymax; iy++) {
03052 int iyp = iyn + iy;
03053 for (int ix = ixmin; ix <= ixmax; ix++) {
03054 int ixp = ixn + ix;
03055 bool mirror = false;
03056 int ixt= ixp, iyt= iyp;
03057 if (ixt < 0) {
03058 ixt = -ixt;
03059 iyt = -iyt;
03060 mirror = !mirror;
03061 }
03062 if (ixt > nhalf) {
03063 ixt = nxreal - ixt;
03064 iyt = -iyt;
03065 mirror = !mirror;
03066 }
03067 if (iyt > nhalf-1) iyt -= nxreal;
03068 if (iyt < -nhalf) iyt += nxreal;
03069 float w = wx[ix]*wy[iy];
03070 std::complex<float> val = this->cmplx(ixt,iyt);
03071 if (mirror) result += conj(val)*w;
03072 else result += val*w;
03073 }
03074 }
03075 }
03076 if (flip) result = conj(result)/wsum;
03077 else result /= wsum;
03078 delete [] wx0;
03079 delete [] wy0;
03080 return result;
03081 }
03082
03083 EMData* EMData::extractline(Util::KaiserBessel& kb, float nuxnew, float nuynew)
03084 {
03085 if (!is_complex())
03086 throw ImageFormatException("extractline requires a fourier image");
03087 if (nx%2 != 0)
03088 throw ImageDimensionException("extractline requires nx to be even");
03089 int nxreal = nx - 2;
03090 if (nxreal != ny)
03091 throw ImageDimensionException("extractline requires ny == nx");
03092
03093 EMData* res = new EMData();
03094 res->set_size(nx,1,1);
03095 res->to_zero();
03096 res->set_complex(true);
03097 res->set_fftodd(false);
03098 res->set_fftpad(true);
03099 res->set_ri(true);
03100
03101 int n = nxreal;
03102 int nhalf = n/2;
03103 vector<int> saved_offsets = get_array_offsets();
03104 set_array_offsets(0,-nhalf,-nhalf);
03105
03106
03107 int kbsize = kb.get_window_size();
03108 int kbmin = -kbsize/2;
03109 int kbmax = -kbmin;
03110 float* wy0 = new float[kbmax - kbmin + 1];
03111 float* wy = wy0 - kbmin;
03112 float* wx0 = new float[kbmax - kbmin + 1];
03113 float* wx = wx0 - kbmin;
03114
03115 int count = 0;
03116 float wsum = 0.f;
03117 bool flip = (nuxnew < 0.f);
03118
03119 for (int jx = 0; jx <= nhalf; jx++) {
03120 float xnew = jx*nuxnew, ynew = jx*nuynew;
03121 count++;
03122 std::complex<float> btq(0.f,0.f);
03123 if (flip) {
03124 xnew = -xnew;
03125 ynew = -ynew;
03126 }
03127 int ixn = int(Util::round(xnew));
03128 int iyn = int(Util::round(ynew));
03129
03130 for (int i=kbmin; i <= kbmax; i++) {
03131 int iyp = iyn + i;
03132 wy[i] = kb.i0win_tab(ynew - iyp);
03133 int ixp = ixn + i;
03134 wx[i] = kb.i0win_tab(xnew - ixp);
03135 }
03136
03137
03138 int lnby = 0;
03139 for (int iy = kbmin; iy <= -1; iy++) {
03140 if (wy[iy] != 0.f) {
03141 lnby = iy;
03142 break;
03143 }
03144 }
03145 int lney = 0;
03146 for (int iy = kbmax; iy >= 1; iy--) {
03147 if (wy[iy] != 0.f) {
03148 lney = iy;
03149 break;
03150 }
03151 }
03152 int lnbx = 0;
03153 for (int ix = kbmin; ix <= -1; ix++) {
03154 if (wx[ix] != 0.f) {
03155 lnbx = ix;
03156 break;
03157 }
03158 }
03159 int lnex = 0;
03160 for (int ix = kbmax; ix >= 1; ix--) {
03161 if (wx[ix] != 0.f) {
03162 lnex = ix;
03163 break;
03164 }
03165 }
03166 if (ixn >= -kbmin && ixn <= nhalf-1-kbmax
03167 && iyn >= -nhalf-kbmin && iyn <= nhalf-1-kbmax) {
03168
03169 for (int ly=lnby; ly<=lney; ly++) {
03170 int iyp = iyn + ly;
03171 for (int lx=lnbx; lx<=lnex; lx++) {
03172 int ixp = ixn + lx;
03173 float wg = wx[lx]*wy[ly];
03174 btq += cmplx(ixp,iyp)*wg;
03175 wsum += wg;
03176 }
03177 }
03178 } else {
03179
03180 for (int ly=lnby; ly<=lney; ly++) {
03181 int iyp = iyn + ly;
03182 for (int lx=lnbx; lx<=lnex; lx++) {
03183 int ixp = ixn + lx;
03184 float wg = wx[lx]*wy[ly];
03185 bool mirror = false;
03186 int ixt(ixp), iyt(iyp);
03187 if (ixt > nhalf || ixt < -nhalf) {
03188 ixt = Util::sgn(ixt)*(n - abs(ixt));
03189 iyt = -iyt;
03190 mirror = !mirror;
03191 }
03192 if (iyt >= nhalf || iyt < -nhalf) {
03193 if (ixt != 0) {
03194 ixt = -ixt;
03195 iyt = Util::sgn(iyt)*(n - abs(iyt));
03196 mirror = !mirror;
03197 } else {
03198 iyt -= n*Util::sgn(iyt);
03199 }
03200 }
03201 if (ixt < 0) {
03202 ixt = -ixt;
03203 iyt = -iyt;
03204 mirror = !mirror;
03205 }
03206 if (iyt == nhalf) iyt = -nhalf;
03207 if (mirror) btq += conj(cmplx(ixt,iyt))*wg;
03208 else btq += cmplx(ixt,iyt)*wg;
03209 wsum += wg;
03210 }
03211 }
03212 }
03213 if (flip) res->cmplx(jx) = conj(btq);
03214 else res->cmplx(jx) = btq;
03215 }
03216 for (int jx = 0; jx <= nhalf; jx++) res->cmplx(jx) *= count/wsum;
03217
03218 delete[] wx0; delete[] wy0;
03219 set_array_offsets(saved_offsets);
03220 res->set_array_offsets(0,0,0);
03221 return res;
03222 }
03223
03224
03226 inline void swapx(float* a, float* b, float* temp, size_t nbytes) {
03227 memcpy(temp, a, nbytes);
03228 memcpy(a, b, nbytes);
03229 memcpy(b, temp, nbytes);
03230 }
03231
03232 void EMData::fft_shuffle() {
03233 if (!is_complex())
03234 throw ImageFormatException("fft_shuffle requires a fourier image");
03235 vector<int> offsets = get_array_offsets();
03236 set_array_offsets();
03237 EMData& self = *this;
03238 int nyhalf = ny/2;
03239 int nzhalf = nz/2;
03240 int nbytes = nx*sizeof(float);
03241 float* temp = new float[nx];
03242 for (int iz=0; iz < nz; iz++)
03243 for (int iy=0; iy < nyhalf; iy++)
03244 swapx(&self(0,iy,iz),&self(0,iy+nyhalf,iz),temp,nbytes);
03245 if (nz > 1) {
03246 for (int iy=0; iy < ny; iy++)
03247 for (int iz=0; iz < nzhalf; iz++)
03248 swapx(&self(0,iy,iz),&self(0,iy,iz+nzhalf),temp,nbytes);
03249 }
03250 set_shuffled(!is_shuffled());
03251 set_array_offsets(offsets);
03252 update();
03253 delete[] temp;
03254 }
03255
03256 void EMData::pad_corner(float *pad_image) {
03257 int nbytes = nx*sizeof(float);
03258 for (int iy=0; iy<ny; iy++)
03259 memcpy(&(*this)(0,iy), pad_image+3+(iy+3)*nx, nbytes);
03260 }
03261
03262 void EMData::shuffle_pad_corner(float *pad_image) {
03263 int nyhalf = ny/2;
03264 int nbytes = nx*sizeof(float);
03265 for (int iy = 0; iy < nyhalf; iy++)
03266 memcpy(&(*this)(0,iy), pad_image+6+(iy+nyhalf+3)*nx, nbytes);
03267 for (int iy = nyhalf; iy < ny; iy++)
03268 memcpy(&(*this)(0,iy), pad_image+6+(iy-nyhalf+3)*nx, nbytes);
03269 }
03270
03271 #define QUADPI 3.141592653589793238462643383279502884197
03272 #define DGR_TO_RAD QUADPI/180
03273
03274
03275
03276
03277
03278
03279
03280
03281
03282
03283
03284
03285
03286
03287
03288
03289
03290
03291
03292
03293
03294
03295
03296
03297
03298
03299
03300
03301
03302
03303
03304
03305
03306
03307
03308
03309
03310
03311
03312
03313
03314
03315
03316
03317
03318
03319
03320
03321
03322
03323
03324
03325
03326
03327 EMData* EMData::fouriergridrot2d(float ang, float scale, Util::KaiserBessel& kb) {
03328 if (2 != get_ndim())
03329 throw ImageDimensionException("fouriergridrot2d needs a 2-D image.");
03330 if (!is_complex())
03331 throw ImageFormatException("fouriergridrot2d requires a fourier image");
03332 int nxreal = nx - 2 + int(is_fftodd());
03333 if (nxreal != ny)
03334 throw ImageDimensionException("fouriergridrot2d requires ny == nx(real)");
03335 if (0 != nxreal%2)
03336 throw ImageDimensionException("fouriergridrot2d needs an even image.");
03337 if (scale == 0.0f) scale = 1.0f;
03338 int nxhalf = nxreal/2;
03339 int nyhalf = ny/2;
03340 float cir = (nxhalf-1)*(nxhalf-1);
03341
03342 if (!is_shuffled()) fft_shuffle();
03343
03344 EMData* result = copy_head();
03345 set_array_offsets(0,-nyhalf);
03346 result->set_array_offsets(0,-nyhalf);
03347
03348
03349
03350 ang = ang*(float)DGR_TO_RAD;
03351 float cang = cos(ang);
03352 float sang = sin(ang);
03353 for (int iy = -nyhalf; iy < nyhalf; iy++) {
03354 float ycang = iy*cang;
03355 float ysang = iy*sang;
03356 for (int ix = 0; ix <= nxhalf; ix++) {
03357 float nuxold = (ix*cang - ysang)*scale;
03358 float nuyold = (ix*sang + ycang)*scale;
03359 if (nuxold*nuxold+nuyold*nuyold<cir) result->cmplx(ix,iy) = Util::extractpoint2(nx, ny, nuxold, nuyold, this, kb);
03360
03361 }
03362 }
03363 result->set_array_offsets();
03364 result->fft_shuffle();
03365 result->update();
03366 set_array_offsets();
03367 fft_shuffle();
03368 return result;
03369 }
03370
03371 EMData* EMData::fouriergridrot_shift2d(float ang, float sx, float sy, Util::KaiserBessel& kb) {
03372 if (2 != get_ndim())
03373 throw ImageDimensionException("fouriergridrot_shift2d needs a 2-D image.");
03374 if (!is_complex())
03375 throw ImageFormatException("fouriergridrot_shift2d requires a fourier image");
03376 int nxreal = nx - 2 + int(is_fftodd());
03377 if (nxreal != ny)
03378 throw ImageDimensionException("fouriergridrot_shift2d requires ny == nx(real)");
03379 if (0 != nxreal%2)
03380 throw ImageDimensionException("fouriergridrot_shift2d needs an even image.");
03381 int nxhalf = nxreal/2;
03382 int nyhalf = ny/2;
03383
03384 if (!is_shuffled()) fft_shuffle();
03385
03386 EMData* result = copy_head();
03387 set_array_offsets(0, -nyhalf);
03388 result->set_array_offsets(0, -nyhalf);
03389
03390 ang = ang*(float)DGR_TO_RAD;
03391 float cang = cos(ang);
03392 float sang = sin(ang);
03393 float temp = -2.0f*M_PI/nxreal;
03394 for (int iy = -nyhalf; iy < nyhalf; iy++) {
03395 float ycang = iy*cang;
03396 float ysang = iy*sang;
03397 for (int ix = 0; ix <= nxhalf; ix++) {
03398 float nuxold = ix*cang - ysang;
03399 float nuyold = ix*sang + ycang;
03400 result->cmplx(ix,iy) = Util::extractpoint2(nx, ny, nuxold, nuyold, this, kb);
03401
03402 float phase_ang = temp*(sx*ix+sy*iy);
03403 result->cmplx(ix,iy) *= complex<float>(cos(phase_ang), sin(phase_ang));
03404 }
03405 }
03406 result->set_array_offsets();
03407 result->fft_shuffle();
03408 result->update();
03409 set_array_offsets();
03410 fft_shuffle();
03411 return result;
03412 }
03413
03414 #undef QUADPI
03415 #undef DGR_TO_RAD
03416
03417 void EMData::divkbsinh(const Util::KaiserBessel& kb) {
03418 if (is_complex())
03419 throw ImageFormatException("divkbsinh requires a real image.");
03420 vector<int> saved_offsets = get_array_offsets();
03421 set_array_offsets(0,0,0);
03422
03423
03424
03425
03426 for (int iz=0; iz < nz; iz++) {
03427 float wz = kb.sinhwin(static_cast<float>(iz-nz/2));
03428 for (int iy=0; iy < ny; iy++) {
03429 float wy = kb.sinhwin(static_cast<float>(iy-ny/2));
03430 for (int ix=0; ix < nx; ix++) {
03431 float wx = kb.sinhwin(static_cast<float>(ix-nx/2));
03432 float w = wx*wy*wz;
03433 (*this)(ix,iy,iz) /= w;
03434 }
03435 }
03436 }
03437 set_array_offsets(saved_offsets);
03438 }
03439
03440
03441
03442
03443
03444
03445
03446
03447
03448
03449
03450
03451
03452
03453
03454
03455
03456
03457
03458
03459
03460
03461
03462
03463
03464
03465 EMData* EMData::extract_plane(const Transform& tf, Util::KaiserBessel& kb) {
03466 if (!is_complex())
03467 throw ImageFormatException("extractplane requires a complex image");
03468 if (nx%2 != 0)
03469 throw ImageDimensionException("extractplane requires nx to be even");
03470 int nxreal = nx - 2;
03471 if (nxreal != ny || nxreal != nz)
03472 throw ImageDimensionException("extractplane requires ny == nx == nz");
03473
03474 EMData* res = new EMData();
03475 res->set_size(nx,ny,1);
03476 res->to_zero();
03477 res->set_complex(true);
03478 res->set_fftodd(false);
03479 res->set_fftpad(true);
03480 res->set_ri(true);
03481
03482 int n = nxreal;
03483 int nhalf = n/2;
03484 vector<int> saved_offsets = get_array_offsets();
03485 set_array_offsets(0,-nhalf,-nhalf);
03486 res->set_array_offsets(0,-nhalf,0);
03487
03488 int kbsize = kb.get_window_size();
03489 int kbmin = -kbsize/2;
03490 int kbmax = -kbmin;
03491 float* wy0 = new float[kbmax - kbmin + 1];
03492 float* wy = wy0 - kbmin;
03493 float* wx0 = new float[kbmax - kbmin + 1];
03494 float* wx = wx0 - kbmin;
03495 float* wz0 = new float[kbmax - kbmin + 1];
03496 float* wz = wz0 - kbmin;
03497 float rim = nhalf*float(nhalf);
03498 int count = 0;
03499 float wsum = 0.f;
03500 Transform tftrans = tf;
03501 tftrans.invert();
03502 for (int jy = -nhalf; jy < nhalf; jy++) {
03503 for (int jx = 0; jx <= nhalf; jx++) {
03504 Vec3f nucur((float)jx, (float)jy, 0.f);
03505 Vec3f nunew = tftrans*nucur;
03506 float xnew = nunew[0], ynew = nunew[1], znew = nunew[2];
03507 if (xnew*xnew+ynew*ynew+znew*znew <= rim) {
03508 count++;
03509 std::complex<float> btq(0.f,0.f);
03510 bool flip = false;
03511 if (xnew < 0.f) {
03512 flip = true;
03513 xnew = -xnew;
03514 ynew = -ynew;
03515 znew = -znew;
03516 }
03517 int ixn = int(Util::round(xnew));
03518 int iyn = int(Util::round(ynew));
03519 int izn = int(Util::round(znew));
03520
03521 for (int i=kbmin; i <= kbmax; i++) {
03522 int izp = izn + i;
03523 wz[i] = kb.i0win_tab(znew - izp);
03524 int iyp = iyn + i;
03525 wy[i] = kb.i0win_tab(ynew - iyp);
03526 int ixp = ixn + i;
03527 wx[i] = kb.i0win_tab(xnew - ixp);
03528
03529 }
03530
03531 int lnbz = 0;
03532 for (int iz = kbmin; iz <= -1; iz++) {
03533 if (wz[iz] != 0.f) {
03534 lnbz = iz;
03535 break;
03536 }
03537 }
03538 int lnez = 0;
03539 for (int iz = kbmax; iz >= 1; iz--) {
03540 if (wz[iz] != 0.f) {
03541 lnez = iz;
03542 break;
03543 }
03544 }
03545 int lnby = 0;
03546 for (int iy = kbmin; iy <= -1; iy++) {
03547 if (wy[iy] != 0.f) {
03548 lnby = iy;
03549 break;
03550 }
03551 }
03552 int lney = 0;
03553 for (int iy = kbmax; iy >= 1; iy--) {
03554 if (wy[iy] != 0.f) {
03555 lney = iy;
03556 break;
03557 }
03558 }
03559 int lnbx = 0;
03560 for (int ix = kbmin; ix <= -1; ix++) {
03561 if (wx[ix] != 0.f) {
03562 lnbx = ix;
03563 break;
03564 }
03565 }
03566 int lnex = 0;
03567 for (int ix = kbmax; ix >= 1; ix--) {
03568 if (wx[ix] != 0.f) {
03569 lnex = ix;
03570 break;
03571 }
03572 }
03573 if (ixn >= -kbmin && ixn <= nhalf-1-kbmax
03574 && iyn >= -nhalf-kbmin && iyn <= nhalf-1-kbmax
03575 && izn >= -nhalf-kbmin && izn <= nhalf-1-kbmax) {
03576
03577 for (int lz = lnbz; lz <= lnez; lz++) {
03578 int izp = izn + lz;
03579 for (int ly=lnby; ly<=lney; ly++) {
03580 int iyp = iyn + ly;
03581 float ty = wz[lz]*wy[ly];
03582 for (int lx=lnbx; lx<=lnex; lx++) {
03583 int ixp = ixn + lx;
03584 float wg = wx[lx]*ty;
03585 btq += cmplx(ixp,iyp,izp)*wg;
03586 wsum += wg;
03587 }
03588 }
03589 }
03590 } else {
03591
03592 for (int lz = lnbz; lz <= lnez; lz++) {
03593 int izp = izn + lz;
03594 for (int ly=lnby; ly<=lney; ly++) {
03595 int iyp = iyn + ly;
03596 float ty = wz[lz]*wy[ly];
03597 for (int lx=lnbx; lx<=lnex; lx++) {
03598 int ixp = ixn + lx;
03599 float wg = wx[lx]*ty;
03600 bool mirror = false;
03601 int ixt(ixp), iyt(iyp), izt(izp);
03602 if (ixt > nhalf || ixt < -nhalf) {
03603 ixt = Util::sgn(ixt)
03604 *(n - abs(ixt));
03605 iyt = -iyt;
03606 izt = -izt;
03607 mirror = !mirror;
03608 }
03609 if (iyt >= nhalf || iyt < -nhalf) {
03610 if (ixt != 0) {
03611 ixt = -ixt;
03612 iyt = Util::sgn(iyt)
03613 *(n - abs(iyt));
03614 izt = -izt;
03615 mirror = !mirror;
03616 } else {
03617 iyt -= n*Util::sgn(iyt);
03618 }
03619 }
03620 if (izt >= nhalf || izt < -nhalf) {
03621 if (ixt != 0) {
03622 ixt = -ixt;
03623 iyt = -iyt;
03624 izt = Util::sgn(izt)
03625 *(n - abs(izt));
03626 mirror = !mirror;
03627 } else {
03628 izt -= Util::sgn(izt)*n;
03629 }
03630 }
03631 if (ixt < 0) {
03632 ixt = -ixt;
03633 iyt = -iyt;
03634 izt = -izt;
03635 mirror = !mirror;
03636 }
03637 if (iyt == nhalf) iyt = -nhalf;
03638 if (izt == nhalf) izt = -nhalf;
03639 if (mirror) btq += conj(cmplx(ixt,iyt,izt))*wg;
03640 else btq += cmplx(ixt,iyt,izt)*wg;
03641 wsum += wg;
03642 }
03643 }
03644 }
03645 }
03646 if (flip) res->cmplx(jx,jy) = conj(btq);
03647 else res->cmplx(jx,jy) = btq;
03648 }
03649 }
03650 }
03651 for (int jy = -nhalf; jy < nhalf; jy++)
03652 for (int jx = 0; jx <= nhalf; jx++)
03653 res->cmplx(jx,jy) *= count/wsum;
03654 delete[] wx0; delete[] wy0; delete[] wz0;
03655 set_array_offsets(saved_offsets);
03656 res->set_array_offsets(0,0,0);
03657 res->set_shuffled(true);
03658 return res;
03659 }
03660
03661
03662 bool EMData::peakcmp(const Pixel& p1, const Pixel& p2) {
03663 return (p1.value > p2.value);
03664 }
03665
03666 ostream& operator<< (ostream& os, const Pixel& peak) {
03667 os << peak.x << peak.y << peak.z << peak.value;
03668 return os;
03669 }
03670
03671 vector<float> EMData::max_search() {
03672
03673 EMData& buf = *this;
03674
03675 int nx = buf.get_xsize();
03676 int ny = buf.get_ysize();
03677 int nz = buf.get_zsize();
03678
03679 int dim = buf.get_ndim();
03680
03681 vector<float> result;
03682
03683 if (dim == 1) {
03684 float max = -1e20f;
03685 int index = -1;
03686 for (int i=0; i<nx; i++) {
03687 if (buf(i)>max) {
03688 max = buf(i);
03689 index = i;
03690 }
03691 }
03692 result.push_back((float)index);
03693 result.push_back(max);
03694 } else if (dim == 2) {
03695 float max = -1e20f;
03696 int index1 = -1;
03697 int index2 = -1;
03698 for (int i=0; i<nx; i++) {
03699 for (int j=0; j<ny; j++) {
03700 if (buf(i, j)>max) {
03701 max = buf(i, j);
03702 index1 = i;
03703 index2 = j;
03704 }
03705 }
03706 }
03707 result.push_back((float)index1);
03708 result.push_back((float)index2);
03709 result.push_back(max);
03710 } else {
03711 float max = -1e20f;
03712 int index1 = -1;
03713 int index2 = -1;
03714 int index3 = -1;
03715 for (int i=0; i<nx; i++) {
03716 for (int j=0; j<ny; j++) {
03717 for (int k=0; k<nz; k++) {
03718 if (buf(i, j, k)>max) {
03719 max = buf(i, j, k);
03720 index1 = i;
03721 index2 = j;
03722 index3 = k;
03723 }
03724 }
03725 }
03726 }
03727 result.push_back((float)index1);
03728 result.push_back((float)index2);
03729 result.push_back((float)index3);
03730 result.push_back(max);
03731 }
03732 return result;
03733 }
03734
03735 vector<float> EMData::peak_search(int ml, float invert) {
03736
03737 EMData& buf = *this;
03738 vector<Pixel> peaks;
03739 int img_dim;
03740 int i, j, k;
03741 int i__1, i__2;
03742 int j__1, j__2;
03743
03744 bool peak_check;
03745 img_dim=buf.get_ndim();
03746 vector<int> ix, jy, kz;
03747 vector<float>res;
03748 int nx = buf.get_xsize();
03749 int ny = buf.get_ysize();
03750 int nz = buf.get_zsize();
03751 if(invert <= 0.0f) invert=-1.0f;
03752 else invert=1.0f ;
03753 int count = 0;
03754 switch (img_dim) {
03755 case(1):
03756 for(i=0;i<=nx-1;++i) {
03757 i__1 = (i-1+nx)%nx;
03758 i__2 = (i+1)%nx;
03759
03760
03761
03762 float qbf = buf(i)*invert;
03763 peak_check = qbf > buf(i__1)*invert && qbf > buf(i__2)*invert;
03764 if(peak_check) {
03765 if(count < ml) {
03766 count++;
03767 peaks.push_back( Pixel(i, 0, 0, qbf) );
03768 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
03769 } else {
03770 if( qbf > (peaks.back()).value ) {
03771
03772 peaks.pop_back();
03773 peaks.push_back( Pixel(i, 0, 0, qbf) );
03774 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
03775 }
03776 }
03777 }
03778 }
03779 break;
03780 case(2):
03781
03782
03783
03784
03785
03786
03787
03788
03789 for(j=1;j<=ny-2;++j) {
03790 j__1 = j-1;
03791 j__2 = j+1;
03792 for(i=1;i<=nx-2;++i) {
03793 i__1 = i-1;
03794 i__2 = i+1;
03795 float qbf = buf(i,j)*invert;
03796 peak_check = (qbf > buf(i,j__1)*invert) && (qbf > buf(i,j__2)*invert);
03797 if(peak_check) {
03798 peak_check = (qbf > buf(i__1,j)*invert) && (qbf > buf(i__2,j)*invert);
03799 if(peak_check) {
03800 peak_check = (qbf > buf(i__1,j__1)*invert) && (qbf > buf(i__1,j__2)*invert);
03801 if(peak_check) {
03802 peak_check = (qbf > buf(i__2,j__1)*invert) && (qbf > buf(i__2,j__2)*invert);
03803 if(peak_check) {
03804 if(count < ml) {
03805 count++;
03806 peaks.push_back( Pixel(i, j, 0, qbf) );
03807 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
03808 } else {
03809 if( qbf > (peaks.back()).value ) {
03810
03811 peaks.pop_back();
03812 peaks.push_back( Pixel(i, j, 0, qbf) );
03813 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
03814 }
03815 }
03816 }
03817 }
03818 }
03819 }
03820 }
03821 }
03822 break;
03823 case(3):
03824
03825
03826
03827
03828
03829
03830
03831
03832
03833
03834
03835
03836
03837
03838
03839
03840
03841
03842
03843
03844 for(k=1; k<=nz-2; ++k) {
03845 kz.clear();
03846 kz.push_back(k-1);
03847 kz.push_back(k);
03848 kz.push_back(k+1);
03849 for(j=1; j<=ny-2; ++j) {
03850 jy.clear();
03851 jy.push_back(j-1);
03852 jy.push_back(j);
03853 jy.push_back(j+1);
03854 for(i=1; i<=nx-2; ++i) {
03855 ix.clear();
03856 ix.push_back(i-1);
03857 ix.push_back(i);
03858 ix.push_back(i+1);
03859 float qbf = buf(i,j,k)*invert;
03860 peak_check = qbf > buf(ix[0],jy[0],kz[0])*invert;
03861 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[0])*invert;
03862 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[0])*invert;
03863 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[0])*invert;
03864 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[0])*invert;
03865 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[0])*invert;
03866 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[0])*invert;
03867 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[0])*invert;
03868 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[0])*invert;
03869 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[1])*invert;
03870 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[1])*invert;
03871 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[1])*invert;
03872 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[1])*invert;
03873
03874 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[1])*invert;
03875 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[1])*invert;
03876 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[1])*invert;
03877 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[1])*invert;
03878 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[2])*invert;
03879 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[2])*invert;
03880 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[2])*invert;
03881 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[2])*invert;
03882 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[2])*invert;
03883 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[2])*invert;
03884 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[2])*invert;
03885 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[2])*invert;
03886 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[2])*invert;
03887 if(peak_check) {
03888 if(count < ml) {
03889 count++;
03890 peaks.push_back( Pixel(i, j, k, qbf) );
03891 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
03892 } else {
03893 if( qbf > (peaks.back()).value ) {
03894
03895
03896 peaks.pop_back();
03897 peaks.push_back( Pixel(i, j, k, qbf) );
03898 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
03899 }
03900 }
03901 }
03902 }}}}}}}}}}}}}}}}}}}}}}}}}
03903 }
03904 }
03905 }
03906
03907
03908 for(k=1; k<=nz-2; ++k) {
03909 kz.clear();
03910 kz.push_back(k-1);
03911 kz.push_back(k);
03912 kz.push_back(k+1);
03913 for(j=1; j<=ny-2; ++j) {
03914 jy.clear();
03915 jy.push_back(j-1);
03916 jy.push_back(j);
03917 jy.push_back(j+1);
03918 for(i=0; i<=0; ++i) {
03919 ix.clear();
03920 ix.push_back(nx-1);
03921 ix.push_back(i);
03922 ix.push_back(i+1);
03923 float qbf = buf(i,j,k)*invert;
03924 peak_check = qbf > buf(ix[0],jy[0],kz[0])*invert;
03925 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[0])*invert;
03926 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[0])*invert;
03927 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[0])*invert;
03928 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[0])*invert;
03929 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[0])*invert;
03930 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[0])*invert;
03931 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[0])*invert;
03932 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[0])*invert;
03933 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[1])*invert;
03934 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[1])*invert;
03935 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[1])*invert;
03936 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[1])*invert;
03937
03938 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[1])*invert;
03939 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[1])*invert;
03940 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[1])*invert;
03941 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[1])*invert;
03942 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[2])*invert;
03943 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[2])*invert;
03944 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[2])*invert;
03945 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[2])*invert;
03946 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[2])*invert;
03947 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[2])*invert;
03948 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[2])*invert;
03949 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[2])*invert;
03950 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[2])*invert;
03951 if(peak_check) {
03952 if(count < ml) {
03953 count++;
03954 peaks.push_back( Pixel(i, j, k, qbf) );
03955 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
03956 } else {
03957 if( qbf > (peaks.back()).value ) {
03958
03959
03960 peaks.pop_back();
03961 peaks.push_back( Pixel(i, j, k, qbf) );
03962 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
03963 }
03964 }
03965 }
03966 }}}}}}}}}}}}}}}}}}}}}}}}}
03967 }
03968 for(i=nx-1; i<=nx-1; ++i) {
03969 ix.clear();
03970 ix.push_back(i-1);
03971 ix.push_back(i);
03972 ix.push_back(0);
03973 float qbf = buf(i,j,k)*invert;
03974 peak_check = qbf > buf(ix[0],jy[0],kz[0])*invert;
03975 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[0])*invert;
03976 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[0])*invert;
03977 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[0])*invert;
03978 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[0])*invert;
03979 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[0])*invert;
03980 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[0])*invert;
03981 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[0])*invert;
03982 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[0])*invert;
03983 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[1])*invert;
03984 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[1])*invert;
03985 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[1])*invert;
03986 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[1])*invert;
03987
03988 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[1])*invert;
03989 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[1])*invert;
03990 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[1])*invert;
03991 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[1])*invert;
03992 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[0],kz[2])*invert;
03993 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[0],kz[2])*invert;
03994 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[0],kz[2])*invert;
03995 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[1],kz[2])*invert;
03996 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[1],kz[2])*invert;
03997 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[1],kz[2])*invert;
03998 if( peak_check ) {peak_check = qbf > buf(ix[0],jy[2],kz[2])*invert;
03999 if( peak_check ) {peak_check = qbf > buf(ix[1],jy[2],kz[2])*invert;
04000 if( peak_check ) {peak_check = qbf > buf(ix[2],jy[2],kz[2])*invert;
04001 if(peak_check) {
04002 if(count < ml) {
04003 count++;
04004 peaks.push_back( Pixel(i, j, k, qbf) );
04005 if(count == ml-1) sort(peaks.begin(), peaks.end(), peakcmp);
04006 } else {
04007 if( qbf > (peaks.back()).value ) {
04008
04009
04010 peaks.pop_back();
04011 peaks.push_back( Pixel(i, j, k, qbf) );
04012 if(ml > 1) sort(peaks.begin(), peaks.end(), peakcmp);
04013 }
04014 }
04015 }
04016 }}}}}}}}}}}}}}}}}}}}}}}}}
04017 }
04018 }
04019 }
04020 break;
04021
04022
04023
04024
04025
04026
04027
04028
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
04084
04085
04086
04087
04088
04089
04090
04091
04092
04093
04094
04095
04096
04097
04098
04099
04100
04101
04102
04103
04104
04105
04106
04107
04108
04109
04110
04111
04112
04113
04114
04115
04116
04117
04118
04119
04120
04121
04122
04123
04124
04125
04126
04127
04128
04129
04130
04131
04132
04133
04134
04135
04136
04137
04138
04139
04140
04141
04142
04143
04144
04145
04146
04147
04148
04149
04150
04151
04152
04153
04154
04155
04156
04157
04158
04159
04160
04161
04162
04163
04164
04165
04166
04167
04168
04169
04170
04171
04172
04173
04174
04175
04176
04177
04178
04179
04180
04181
04182
04183
04184
04185
04186
04187
04188
04189
04190
04191
04192
04193
04194
04195
04196
04197
04198
04199
04200
04201
04202
04203
04204
04205
04206
04207
04208
04209
04210
04211
04212
04213
04214
04215
04216
04217
04218
04219
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
04249
04250
04251
04252
04253
04254
04255
04256
04257
04258
04259
04260
04261
04262
04263
04264
04265
04266
04267
04268
04269
04270
04271
04272
04273
04274
04275
04276
04277
04278
04279
04280
04281
04282
04283
04284
04285
04286
04287
04288
04289
04290
04291
04292
04293
04294
04295
04296
04297
04298
04299
04300
04301
04302
04303
04304
04305
04306
04307
04308
04309
04310
04311
04312
04313
04314
04315
04316
04317
04318
04319
04320
04321
04322
04323
04324
04325
04326
04327
04328
04329
04330
04331
04332
04333
04334
04335
04336
04337
04338
04339
04340
04341
04342
04343
04344
04345
04346
04347
04348
04349
04350
04351
04352
04353
04354
04355
04356
04357
04358
04359
04360
04361
04362
04363
04364
04365
04366
04367
04368
04369
04370
04371
04372
04373
04374
04375
04376
04377
04378
04379
04380
04381
04382
04383
04384
04385
04386
04387
04388
04389
04390
04391
04392
04393
04394
04395
04396
04397
04398
04399
04400
04401
04402
04403
04404
04405
04406
04407
04408
04409
04410
04411
04412
04413
04414
04415
04416
04417
04418
04419
04420
04421
04422
04423
04424
04425
04426
04427
04428
04429
04430
04431
04432
04433
04434
04435
04436
04437
04438
04439
04440
04441
04442
04443
04444
04445
04446
04447
04448
04449
04450
04451
04452
04453
04454
04455
04456
04457
04458
04459
04460 }
04461
04462 if (peaks.begin() != peaks.end()) {
04463
04464 sort(peaks.begin(), peaks.end(), peakcmp);
04465
04466 int count = 0;
04467
04468 float xval = (*peaks.begin()).value;
04469
04470 for (vector<Pixel>::iterator it = peaks.begin(); it != peaks.end(); it++) {
04471
04472 count++;
04473
04474 if(count <= ml) {
04475
04476 res.push_back((*it).value);
04477 res.push_back(static_cast<float>((*it).x));
04478
04479 if(img_dim > 1) {
04480 res.push_back(static_cast<float>((*it).y));
04481 if(nz > 1) res.push_back(static_cast<float>((*it).z));
04482 }
04483
04484 if(xval != 0.0) res.push_back((*it).value/xval);
04485 else res.push_back((*it).value);
04486 res.push_back((*it).x-float(int(nx/2)));
04487 if(img_dim >1) {
04488 res.push_back((*it).y-float(int(ny/2)));
04489 if(nz>1) res.push_back((*it).z-float(nz/2));
04490 }
04491 }
04492 }
04493 res.insert(res.begin(),1,img_dim);
04494 } else {
04495
04496 res.push_back(buf(0,0,0));
04497 res.insert(res.begin(),1,0.0);
04498 }
04499
04500
04501 return res;
04502 }
04503
04504 #define rdata(i,j,k) rdata[(i-1)+((j-1)+(k-1)*ny)*nx]
04505 #define X(i) X[i-1]
04506 #define Y(j) Y[j-1]
04507 #define Z(k) Z[k-1]
04508 vector<float> EMData::phase_cog()
04509 {
04510 vector<float> ph_cntog;
04511 int i=1,j=1,k=1;
04512 float C=0.f,S=0.f,P=0.f,F1=0.f,SNX;
04513 if (get_ndim()==1) {
04514 P = 8*atan(1.0f)/nx;
04515 for (i=1;i<=nx;i++) {
04516 C += cos(P * (i-1)) * rdata(i,j,k);
04517 S += sin(P * (i-1)) * rdata(i,j,k);
04518 }
04519 F1 = atan2(S,C);
04520 if (F1 < 0.0) F1 += 8*atan(1.0f);
04521 SNX = F1/P +1.0f;
04522 SNX = SNX - ((nx/2)+1);
04523 ph_cntog.push_back(SNX);
04524 #ifdef _WIN32
04525 ph_cntog.push_back((float)Util::round(SNX));
04526 #else
04527 ph_cntog.push_back(round(SNX));
04528 #endif //_WIN32
04529 } else if (get_ndim()==2) {
04530 #ifdef _WIN32
04531 float SNY;
04532 float T=0.0f;
04533 vector<float> X;
04534 X.resize(nx);
04535 #else
04536 float SNY,X[nx],T=0.f;
04537 #endif //_WIN32
04538 for ( i=1;i<=nx;i++) X(i)=0.0;
04539 P = 8*atan(1.0f)/ny;
04540 for(j=1;j<=ny;j++) {
04541 T=0.f;
04542 for(i=1;i<=nx;i++) {
04543 T += rdata(i,j,k);
04544 X(i)+=rdata(i,j,k);
04545 }
04546 C += cos(P*(j-1))*T;
04547 S += sin(P*(j-1))*T;
04548 }
04549 F1=atan2(S,C);
04550 if(F1<0.0) F1 += 8*atan(1.0f);
04551 SNY = F1/P +1.0f;
04552 C=0.f; S=0.f;
04553 P = 8*atan(1.0f)/nx;
04554 for(i=1;i<=nx;i++) {
04555 C += cos(P*(i-1))*X(i);
04556 S += sin(P*(i-1))*X(i);
04557 }
04558 F1=atan2(S,C);
04559 if(F1<0.0) F1 += 8*atan(1.0f);
04560 SNX = F1/P +1.0f;
04561 SNX = SNX - ((nx/2)+1);
04562 SNY = SNY - ((ny/2)+1);
04563 ph_cntog.push_back(SNX); ph_cntog.push_back(SNY);
04564 #ifdef _WIN32
04565 ph_cntog.push_back((float)Util::round(SNX)); ph_cntog.push_back((float)Util::round(SNY));
04566 #else
04567 ph_cntog.push_back(round(SNX)); ph_cntog.push_back(round(SNY));
04568 #endif //_WIN32
04569 } else {
04570 #ifdef _WIN32
04571 float val=0.f,sum1=0.f, SNY,SNZ;
04572 vector<float> X;
04573 X.resize(nx);
04574 vector<float> Y;
04575 Y.resize(ny);
04576 vector<float> Z;
04577 Z.resize(nz);
04578 #else
04579 float val=0.f, sum1=0.f, X[nx], Y[ny], Z[nz], SNY, SNZ;
04580 #endif //_WIN32
04581 for (i=1;i<=nx;i++) X(i)=0.0;
04582 for (j=1;j<=ny;j++) Y(j)=0.0;
04583 for (k=1;k<=nz;k++) Z(k)=0.0;
04584 for(k=1;k<=nz;k++) {
04585 for(j=1;j<=ny;j++) {
04586 sum1=0.f;
04587 for(i=1;i<=nx;i++) {
04588 val = rdata(i,j,k);
04589 sum1 += val;
04590 X(i) += val;
04591 }
04592 Y(j) += sum1;
04593 Z(k) += sum1;
04594 }
04595 }
04596 P = 8*atan(1.0f)/nx;
04597 for (i=1;i<=nx;i++) {
04598 C += cos(P*(i-1))*X(i);
04599 S += sin(P*(i-1))*X(i);
04600 }
04601 F1=atan2(S,C);
04602 if(F1<0.0) F1 += 8*atan(1.0f);
04603 SNX = F1/P +1.0f;
04604 C=0.f; S=0.f;
04605 P = 8*atan(1.0f)/ny;
04606 for(j=1;j<=ny;j++) {
04607 C += cos(P*(j-1))*Y(j);
04608 S += sin(P*(j-1))*Y(j);
04609 }
04610 F1=atan2(S,C);
04611 if(F1<0.0) F1 += 8*atan(1.0f);
04612 SNY = F1/P +1.0f;
04613 C=0.f; S=0.f;
04614 P = 8*atan(1.0f)/nz;
04615 for(k=1;k<=nz;k++) {
04616 C += cos(P*(k-1))*Z(k);
04617 S += sin(P*(k-1))*Z(k);
04618 }
04619 F1=atan2(S,C);
04620 if(F1<0.0) F1 += 8*atan(1.0f);
04621 SNZ = F1/P +1.0f;
04622 SNX = SNX - ((nx/2)+1);
04623 SNY = SNY - ((ny/2)+1);
04624 SNZ = SNZ - ((nz/2)+1);
04625 ph_cntog.push_back(SNX); ph_cntog.push_back(SNY); ph_cntog.push_back(SNZ);
04626 #ifdef _WIN32
04627 ph_cntog.push_back((float)Util::round(SNX)); ph_cntog.push_back((float)Util::round(SNY)); ph_cntog.push_back((float)Util::round(SNZ));
04628 #else
04629 ph_cntog.push_back(round(SNX)); ph_cntog.push_back(round(SNY));ph_cntog.push_back(round(SNZ));
04630 #endif
04631 }
04632 return ph_cntog;
04633 }
04634 #undef rdata
04635 #undef X
04636 #undef Y
04637 #undef Z
04638
04639 #define avagadro (6.023*(double)pow(10.0,23.0))
04640 #define density_protein (1.36)
04641 #define R (0.61803399f)
04642 #define C (1.f-R)
04643 float EMData::find_3d_threshold(float mass, float pixel_size)
04644 {
04645
04646 if(get_ndim()!=3)
04647 throw ImageDimensionException("The image should be 3D");
04648
04649
04650
04651 float density_1_mole, vol_1_mole, vol_angstrom;
04652 int vol_voxels;
04653 density_1_mole = static_cast<float>( (mass*1000.0f)/avagadro );
04654 vol_1_mole = static_cast<float>( density_1_mole/density_protein );
04655 vol_angstrom = static_cast<float>( vol_1_mole*(double)pow((double)pow(10.0,8),3) );
04656 vol_voxels = static_cast<int> (vol_angstrom/(double)pow(pixel_size,3));
04657
04658
04659
04660 float thr1 = get_attr("maximum");
04661 float thr3 = get_attr("minimum");
04662 float thr2 = (thr1-thr3)/2 + thr3;
04663 int size = nx*ny*nz;
04664 float x0 = thr1,x3 = thr3,x1,x2,THR=0;
04665
04666 #ifdef _WIN32
04667 int ILE = _MIN(nx*ny*nx,_MAX(1,vol_voxels));
04668 #else
04669 int ILE = std::min(nx*ny*nx,std::max(1,vol_voxels));
04670 #endif //_WIN32
04671
04672 if (abs(thr3-thr2)>abs(thr2-thr1)) {
04673 x1=thr2;
04674 x2=thr2+C*(thr3-thr2);
04675 } else {
04676 x2=thr2;
04677 x1=thr2-C*(thr2-thr1);
04678 }
04679
04680 int cnt1=0,cnt2=0;
04681 for (int i=0;i<size;i++) {
04682 if(rdata[i]>=x1) cnt1++;
04683 if(rdata[i]>=x2) cnt2++;
04684 }
04685 float LF1 = static_cast<float>( cnt1 - ILE );
04686 float F1 = LF1*LF1;
04687 float LF2 = static_cast<float>( cnt2 - ILE );
04688 float F2 = LF2*LF2;
04689
04690 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)))
04691 {
04692 if(F2 < F1) {
04693 x0=x1;
04694 x1=x2;
04695 x2 = R*x1 + C*x3;
04696 F1=F2;
04697 int cnt=0;
04698 for(int i=0;i<size;i++)
04699 if(rdata[i]>=x2)
04700 cnt++;
04701 LF2 = static_cast<float>( cnt - ILE );
04702 F2 = LF2*LF2;
04703 } else {
04704 x3=x2;
04705 x2=x1;
04706 x1=R*x2 + C*x0;
04707 F2=F1;
04708 int cnt=0;
04709 for(int i=0;i<size;i++)
04710 if(rdata[i]>=x1)
04711 cnt++;
04712 LF1 = static_cast<float>( cnt - ILE );
04713 F1 = LF1*LF1;
04714 }
04715 }
04716
04717 if(F1 < F2) {
04718 ILE = static_cast<int> (LF1 + ILE);
04719 THR = x1;
04720 } else {
04721 ILE = static_cast<int> (LF2 + ILE);
04722 THR = x2;
04723 }
04724 return THR;
04725
04726 }
04727 #undef avagadro
04728 #undef density_protein
04729 #undef R
04730 #undef C
04731
04732
04733
04734
04735
04736
04737
04738 vector<float> EMData::peak_ccf(float hf_p)
04739 {
04740
04741
04742
04743 EMData & buf = *this;
04744 vector<Pixel> peaks;
04745 int half=int(hf_p);
04746 float hf_p2 = hf_p*hf_p;
04747 int i,j;
04748 int i__1,i__2;
04749 int j__1,j__2;
04750 vector<float>res;
04751 int nx = buf.get_xsize()-half;
04752 int ny = buf.get_ysize()-half;
04753
04754 for(i=half; i<=nx; ++i) {
04755
04756 i__1 = i-1;
04757 i__2 = i+1;
04758 for (j=half;j<=ny;++j) {
04759 j__1 = j-1;
04760 j__2 = j+1;
04761
04762 if((buf(i,j)>0.0f)&&buf(i,j)>buf(i,j__1)) {
04763 if(buf(i,j)>buf(i,j__2)) {
04764 if(buf(i,j)>buf(i__1,j)) {
04765 if(buf(i,j)>buf(i__2,j)) {
04766 if(buf(i,j)>buf(i__1,j__1)) {
04767 if((buf(i,j))> buf(i__1,j__2)) {
04768 if(buf(i,j)>buf(i__2,j__1)) {
04769 if(buf(i,j)> buf(i__2,j__2)) {
04770
04771
04772
04773 if (peaks.size()==0) {
04774
04775 peaks.push_back(Pixel(i,j,0,buf(i,j)));
04776
04777 } else {
04778
04779
04780 bool overlap = false;
04781
04782
04783
04784
04785
04786 std::stack <vector<Pixel>::iterator> delete_stack;
04787
04788
04789 for (vector<Pixel>::iterator it=peaks.begin();it!=peaks.end();++it) {
04790
04791
04792
04793
04794 float radius=((*it).x-float(i))*((*it).x-float(i))+((*it).y-float(j))*((*it).y-float(j));
04795 if (radius <= hf_p2 ) {
04796
04797 if( buf(i,j) > (*it).value) {
04798
04799
04800
04801
04802
04803
04804 delete_stack.push(it);
04805
04806
04807
04808
04809 } else {
04810 overlap = true;
04811
04812
04813 break;
04814 }
04815 }
04816 }
04817
04818
04819
04820 if (false == overlap) {
04821 vector<Pixel>::iterator delete_iterator;
04822 while (!delete_stack.empty()) {
04823
04824
04825 delete_iterator = delete_stack.top();
04826 peaks.erase(delete_iterator);
04827 delete_stack.pop();
04828 }
04829
04830
04831
04832 if (! (peaks.size() < 2000 )) {
04833
04834
04835
04836
04837 sort(peaks.begin(), peaks.end(), peakcmp);
04838
04839
04840 peaks.pop_back();
04841 }
04842
04843
04844 peaks.push_back(Pixel(i,j,0,buf(i,j)));
04845
04846
04847 } else {
04848
04849
04850 while (!delete_stack.empty()) delete_stack.pop();
04851 }
04852 }
04853 }
04854 }}}}}}}
04855 }
04856 }
04857
04858
04859 if(peaks.size()>0) {
04860
04861 sort(peaks.begin(),peaks.end(), peakcmp);
04862
04863 for (vector<Pixel>::iterator it = peaks.begin(); it != peaks.end(); it++) {
04864
04865 res.push_back((*it).value);
04866 res.push_back(static_cast<float>((*it).x));
04867 res.push_back(static_cast<float>((*it).y));
04868 }
04869 } else {
04870
04871 res.push_back(buf(0,0,0));
04872 res.insert(res.begin(),1,0.0);
04873 }
04874 return res;
04875 }
04876
04877 EMData* EMData::get_pow(float n_pow)
04878 {
04879 EMData* buf_new = this->copy_head();
04880 float *in = this->get_data();
04881 float *out = buf_new->get_data();
04882 for(int i=0; i<nx*ny*nz; i++) out[i] = pow(in[i],n_pow);
04883 return buf_new;
04884 }
04885
04886 EMData* EMData::conjg()
04887 {
04888 if(this->is_complex()) {
04889 EMData* buf_new = this->copy_head();
04890 float *in = this->get_data();
04891 float *out = buf_new->get_data();
04892 for(int i=0; i<nx*ny*nz; i+=2) {out[i] = in[i]; out[i+1] = -in[i+1];}
04893 return buf_new;
04894 } else throw ImageFormatException("image has to be complex");
04895 }
04896
04897 EMData* EMData::delete_disconnected_regions(int ix, int iy, int iz) {
04898 if (3 != get_ndim())
04899 throw ImageDimensionException("delete_disconnected_regions needs a 3-D image.");
04900 if (is_complex())
04901 throw ImageFormatException("delete_disconnected_regions requires a real image");
04902 if ((*this)(ix+nx/2,iy+ny/2,iz+nz/2) == 0)
04903 throw ImageDimensionException("delete_disconnected_regions starting point is zero.");
04904
04905 EMData* result = this->copy_head();
04906 result->to_zero();
04907 (*result)(ix+nx/2,iy+ny/2,iz+nz/2) = (*this)(ix+nx/2,iy+ny/2,iz+nz/2);
04908 bool kpt = true;
04909
04910 while(kpt) {
04911 kpt = false;
04912 for (int cz = 1; cz < nz-1; cz++) {
04913 for (int cy = 1; cy < ny-1; cy++) {
04914 for (int cx = 1; cx < nx-1; cx++) {
04915 if((*result)(cx,cy,cz) == 1) {
04916 for (int lz = -1; lz <= 1; lz++) {
04917 for (int ly = -1; ly <= 1; ly++) {
04918 for (int lx = -1; lx <= 1; lx++) {
04919 if(((*this)(cx+lx,cy+ly,cz+lz) == 1) && ((*result)(cx+lx,cy+ly,cz+lz) == 0)) {
04920 (*result)(cx+lx,cy+ly,cz+lz) = 1;
04921 kpt = true;
04922 }
04923 }
04924 }
04925 }
04926 }
04927 }
04928 }
04929 }
04930 }
04931 result->update();
04932 return result;
04933 }
04934
04935 #define QUADPI 3.141592653589793238462643383279502884197
04936 #define DGR_TO_RAD QUADPI/180
04937 EMData* EMData::helicise(float pixel_size, float dp, float dphi, float section_use, float radius, float minrad) {
04938 if (3 != get_ndim())
04939 throw ImageDimensionException("helicise needs a 3-D image.");
04940 if (is_complex())
04941 throw ImageFormatException("helicise requires a real image");
04942
04943 EMData* result = this->copy_head();
04944 result->to_zero();
04945 int nyc = ny/2;
04946 int nxc = nx/2;
04947 int nb = int(nx*(1.0f - section_use)/2.);
04948 int ne = nx - nb -1;
04949 int numst = int((ne - nb)/dp*pixel_size + 0.5);
04950
04951 int nst = int(nz*pixel_size/dp+0.5);
04952 float r2, ir;
04953 if(radius < 0.0f) r2 = (float)((nxc-1)*(nxc-1));
04954 else r2 = radius*radius;
04955 if(minrad < 0.0f) ir = 0.0f;
04956 else ir = minrad*minrad;
04957 for (int k = 0; k<nz; k++) {
04958 for (int j = 0; j<ny; j++) {
04959 int jy = j - nyc;
04960 int jj = jy*jy;
04961 for (int i = 0; i<nx; i++) {
04962 int ix = i - nxc;
04963 float d2 = (float)(ix*ix + jj);
04964 if(d2 <= r2 && d2>=ir) {
04965 int nq = 1;
04966 for ( int ist = -nst; ist <= nst; ist++) {
04967 float zold = (k*pixel_size + ist*dp)/pixel_size;
04968 int IOZ = int(zold);
04969 if(IOZ >= nb && IOZ <= ne) {
04970
04971 float cphi = ist*dphi*(float)DGR_TO_RAD;
04972 float ca = cos(cphi);
04973 float sa = sin(cphi);
04974 float xold = ix*ca - jy*sa + nxc;
04975 float yold = ix*sa + jy*ca + nyc;
04976 nq++;
04977
04978
04979
04980 int IOX = int(xold);
04981 int IOY = int(yold);
04982
04983
04984 #ifdef _WIN32
04985 int IOXp1 = _MIN( nx-1 ,IOX+1);
04986 #else
04987 int IOXp1 = std::min( nx-1 ,IOX+1);
04988 #endif //_WIN32
04989
04990 #ifdef _WIN32
04991 int IOYp1 = _MIN( ny-1 ,IOY+1);
04992 #else
04993 int IOYp1 = std::min( ny-1 ,IOY+1);
04994 #endif //_WIN32
04995
04996 #ifdef _WIN32
04997 int IOZp1 = _MIN( nz-1 ,IOZ+1);
04998 #else
04999 int IOZp1 = std::min( nz-1 ,IOZ+1);
05000 #endif //_WIN32
05001
05002 float dx = xold-IOX;
05003 float dy = yold-IOY;
05004 float dz = zold-IOZ;
05005
05006 float a1 = (*this)(IOX,IOY,IOZ);
05007 float a2 = (*this)(IOXp1,IOY,IOZ) - (*this)(IOX,IOY,IOZ);
05008 float a3 = (*this)(IOX,IOYp1,IOZ) - (*this)(IOX,IOY,IOZ);
05009 float a4 = (*this)(IOX,IOY,IOZp1) - (*this)(IOX,IOY,IOZ);
05010 float a5 = (*this)(IOX,IOY,IOZ) - (*this)(IOXp1,IOY,IOZ) - (*this)(IOX,IOYp1,IOZ) + (*this)(IOXp1,IOYp1,IOZ);
05011 float a6 = (*this)(IOX,IOY,IOZ) - (*this)(IOXp1,IOY,IOZ) - (*this)(IOX,IOY,IOZp1) + (*this)(IOXp1,IOY,IOZp1);
05012 float a7 = (*this)(IOX,IOY,IOZ) - (*this)(IOX,IOYp1,IOZ) - (*this)(IOX,IOY,IOZp1) + (*this)(IOX,IOYp1,IOZp1);
05013 float a8 = (*this)(IOXp1,IOY,IOZ) + (*this)(IOX,IOYp1,IOZ)+ (*this)(IOX,IOY,IOZp1)
05014 - (*this)(IOX,IOY,IOZ)- (*this)(IOXp1,IOYp1,IOZ) - (*this)(IOXp1,IOY,IOZp1)
05015 - (*this)(IOX,IOYp1,IOZp1) + (*this)(IOXp1,IOYp1,IOZp1);
05016 (*result)(i,j,k) += a1 + dz*(a4 + a6*dx + (a7 + a8*dx)*dy) + a3*dy + dx*(a2 + a5*dy);
05017
05018
05019 if(nq == numst) break;
05020 }
05021 }
05022 if(nq != numst)
05023 throw InvalidValueException(nq, "incorrect number of repeats encoutered.");
05024 }
05025 }
05026 }
05027 }
05028 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 ;
05029
05030 result->update();
05031 return result;
05032 }
05033
05034
05035
05036
05037
05038
05039
05040
05041
05042
05043
05044 void EMData::depad() {
05045 if (is_complex())
05046 throw ImageFormatException("Depadding of complex images not supported");
05047 vector<int> saved_offsets = get_array_offsets();
05048 set_array_offsets(0,0,0);
05049 int npad = attr_dict["npad"];
05050 if (0 == npad) npad = 1;
05051 int offset = is_fftodd() ? 1 : 2;
05052 int nxold = (nx - offset)/npad;
05053 #ifdef _WIN32
05054 int nyold = _cpp_max(ny/npad, 1);
05055 int nzold = _cpp_max(nz/npad, 1);
05056 #else
05057 int nyold = std::max<int>(ny/npad, 1);
05058 int nzold = std::max<int>(nz/npad, 1);
05059 #endif //_WIN32
05060 int xstart = 0, ystart = 0, zstart = 0;
05061 if( npad > 1) {
05062 xstart = (nx - offset - nxold)/2 + nxold%2;
05063 if(ny > 1) {
05064 ystart = (ny - nyold)/2 + nyold%2;
05065 if(nz > 1) {
05066 zstart = (nz - nzold)/2 + nzold%2;
05067 }
05068 }
05069 }
05070 int bytes = nxold*sizeof(float);
05071 float* dest = get_data();
05072 for (int iz=0; iz < nzold; iz++) {
05073 for (int iy = 0; iy < nyold; iy++) {
05074 memmove(dest, &(*this)(xstart,iy+ystart,iz+zstart), bytes);
05075 dest += nxold;
05076 }
05077 }
05078 set_size(nxold, nyold, nzold);
05079 set_attr("npad", 1);
05080 set_fftpad(false);
05081 set_fftodd(false);
05082 set_complex(false);
05083 if(ny==1 && nz==1) set_complex_x(false);
05084 set_array_offsets(saved_offsets);
05085 update();
05086 EXITFUNC;
05087 }
05088
05089
05090
05091
05092
05093
05094
05095
05096 void EMData::depad_corner() {
05097 if(is_complex())
05098 throw ImageFormatException("Depadding of complex images not allowed");
05099 vector<int> saved_offsets = get_array_offsets();
05100 set_array_offsets(0,0,0);
05101 int npad = attr_dict["npad"];
05102 if(0 == npad) npad = 1;
05103 int offset = is_fftodd() ? 1 : 2;
05104 int nxold = (nx - offset)/npad;
05105 #ifdef _WIN32
05106 int nyold = _cpp_max(ny/npad, 1);
05107 int nzold = _cpp_max(nz/npad, 1);
05108 #else
05109 int nyold = std::max<int>(ny/npad, 1);
05110 int nzold = std::max<int>(nz/npad, 1);
05111 #endif //_WIN32
05112 int bytes = nxold*sizeof(float);
05113 float* dest = get_data();
05114 for (int iz=0; iz < nzold; iz++) {
05115 for (int iy = 0; iy < nyold; iy++) {
05116 memmove(dest, &(*this)(0,iy,iz), bytes);
05117 dest += nxold;
05118 }
05119 }
05120 set_size(nxold, nyold, nzold);
05121 set_attr("npad", 1);
05122 set_fftpad(false);
05123 set_fftodd(false);
05124 set_complex(false);
05125 if(ny==1 && nz==1) set_complex_x(false);
05126 set_array_offsets(saved_offsets);
05127 update();
05128 EXITFUNC;
05129 }
05130
05131
05132
05133 float circumference( EMData* emdata, int npixel )
05134 {
05135 int nx = emdata->get_xsize();
05136 int ny = emdata->get_ysize();
05137 int nz = emdata->get_zsize();
05138
05139 float* data = emdata->get_data();
05140 if( ny==1 && nz==1 )
05141 {
05142
05143 float sumf=0.0;
05144 int sumn=0;
05145 for( int i=0; i < npixel; ++i )
05146 {
05147 sumf += data[i];
05148 sumf += data[nx-1-i];
05149 sumn += 2;
05150 }
05151 return sumf/sumn;
05152 }
05153
05154 if( nz==1 )
05155 {
05156 float sumf=0.0;
05157 int sumn=0;
05158 int id=0;
05159 for( int iy=0; iy < ny; ++iy )
05160 {
05161 for( int ix=0; ix < nx; ++ix )
05162 {
05163 if( iy<npixel || iy>ny-1-npixel || ix<npixel || ix>nx-1-npixel )
05164 {
05165 sumf += data[id];
05166 sumn += 1;
05167 }
05168 id++;
05169 }
05170 }
05171
05172 Assert( id==nx*ny );
05173 Assert( sumn == nx*ny - (nx-2*npixel)*(ny-2*npixel) );
05174 return sumf/sumn;
05175 }
05176
05177
05178
05179 float sumf = 0.0;
05180 int sumn = 0;
05181 int id = 0;
05182 for( int iz=0; iz < nz; ++iz)
05183 {
05184 for( int iy=0; iy < ny; ++iy)
05185 {
05186 for( int ix=0; ix < nx; ++ix )
05187 {
05188 if( iz<npixel||iz>nz-1-npixel||iy<npixel||iy>ny-1-npixel||ix<npixel||ix>nx-1-npixel)
05189 {
05190 sumf += data[id];
05191 sumn += 1;
05192 }
05193 id++;
05194 }
05195 }
05196 }
05197
05198
05199 Assert( id==nx*ny*nz);
05200 Assert( sumn==nx*ny*nz-(nx-2*npixel)*(ny-2*npixel)*(nz-2*npixel) );
05201 return sumf/sumn;
05202 }
05203
05204
05205
05206
05207
05208
05209
05210
05211 EMData* EMData::norm_pad(bool donorm, int npad, int valtype) {
05212 if (this->is_complex())
05213 throw ImageFormatException("Padding of complex images not supported");
05214 int nx = this->get_xsize();
05215 int ny = this->get_ysize();
05216 int nz = this->get_zsize();
05217 float mean = 0., stddev = 1.;
05218 if(donorm) {
05219 mean = this->get_attr("mean");
05220 stddev = this->get_attr("sigma");
05221 }
05222
05223 if (npad < 1) npad = 1;
05224 int nxpad = npad*nx;
05225 int nypad = npad*ny;
05226 int nzpad = npad*nz;
05227 if (1 == ny) {
05228
05229
05230 nypad = ny;
05231 nzpad = nz;
05232 } else if (nz == 1) {
05233
05234 nzpad = nz;
05235 }
05236 size_t bytes;
05237 size_t offset;
05238
05239 offset = 2 - nxpad%2;
05240 bytes = nx*sizeof(float);
05241 EMData* fpimage = copy_head();
05242 fpimage->set_size(nxpad+offset, nypad, nzpad);
05243 if( npad > 1) {
05244 if( valtype==0 ) {
05245 fpimage->to_zero();
05246 } else {
05247 float val = circumference(this, 1);
05248 float* data = fpimage->get_data();
05249 int nxyz = (nxpad+offset)*nypad*nzpad;
05250 for( int i=0; i < nxyz; ++i ) data[i] = val;
05251 }
05252 }
05253
05254
05255
05256 int xstart = 0, ystart = 0, zstart = 0;
05257 if (npad > 1) {
05258 xstart = (nxpad - nx)/2 + nx%2;
05259 if(ny > 1) {
05260 ystart = (nypad - ny)/2 + ny%2;
05261 if(nz > 1) {
05262 zstart = (nzpad - nz)/2 + nz%2;
05263 }
05264 }
05265 }
05266
05267
05268 vector<int> saved_offsets = this->get_array_offsets();
05269 this->set_array_offsets( 0, 0, 0 );
05270 for (int iz = 0; iz < nz; iz++) {
05271 for (int iy = 0; iy < ny; iy++) {
05272 memcpy(&(*fpimage)(xstart,iy+ystart,iz+zstart), &(*this)(0,iy,iz), bytes);
05273 }
05274 }
05275 this->set_array_offsets( saved_offsets );
05276
05277
05278
05279
05280 if (donorm) {
05281 for (int iz = zstart; iz < nz+zstart; iz++)
05282 for (int iy = ystart; iy < ny+ystart; iy++)
05283 for (int ix = xstart; ix < nx+xstart; ix++)
05284 (*fpimage)(ix,iy,iz) = ((*fpimage)(ix,iy,iz)-mean)/stddev;
05285 }
05286
05287 fpimage->set_fftpad(true);
05288 fpimage->set_attr("npad", npad);
05289 if (offset == 1) fpimage->set_fftodd(true);
05290 else fpimage->set_fftodd(false);
05291 return fpimage;
05292 }
05293
05294 void EMData::center_origin()
05295 {
05296 ENTERFUNC;
05297 if (is_complex()) {
05298 LOGERR("Real image expected. Input image is complex.");
05299 throw ImageFormatException("Real image expected. Input image is complex.");
05300 }
05301 for (int iz = 0; iz < nz; iz++) {
05302 for (int iy = 0; iy < ny; iy++) {
05303 for (int ix = 0; ix < nx; ix++) {
05304
05305 (*this)(ix,iy,iz) *= -2*((ix+iy+iz)%2) + 1;
05306 }
05307 }
05308 }
05309 update();
05310 EXITFUNC;
05311 }
05312
05313 void EMData::center_origin_yz()
05314 {
05315 ENTERFUNC;
05316 if (is_complex()) {
05317 LOGERR("Real image expected. Input image is complex.");
05318 throw ImageFormatException("Real image expected. Input image is complex.");
05319 }
05320 for (int iz = 0; iz < nz; iz++) {
05321 for (int iy = (iz+1)%2; iy < ny; iy+=2) {
05322 for (int ix = 0; ix < nx; ix++) {
05323 (*this)(ix,iy,iz) *= -1;
05324 }
05325 }
05326 }
05327 update();
05328 EXITFUNC;
05329 }
05330
05331 void EMData::center_origin_fft()
05332 {
05333 ENTERFUNC;
05334 if (!is_complex()) {
05335 LOGERR("complex image expected. Input image is real image.");
05336 throw ImageFormatException("complex image expected. Input image is real image.");
05337 }
05338
05339 if (!is_ri()) {
05340 LOGWARN("Only RI should be used. ");
05341 }
05342 vector<int> saved_offsets = get_array_offsets();
05343
05344
05345
05346 set_array_offsets(0,1,1);
05347 int nxc = nx/2;
05348
05349 if (is_fftodd()) {
05350 for (int iz = 1; iz <= nz; iz++) {
05351 for (int iy = 1; iy <= ny; iy++) {
05352 for (int ix = 0; ix < nxc; ix++) {
05353 cmplx(ix,iy,iz) *= float(-2*((ix+iy+iz)%2) + 1);
05354 float temp = float(iz-1+iy-1+ix)/float(ny)*M_PI;
05355 complex<float> temp2 = complex<float>(cos(temp), -sin(temp));
05356 cmplx(ix,iy,iz) *= temp2;
05357 }
05358 }
05359 }
05360 } else {
05361 for (int iz = 1; iz <= nz; iz++) {
05362 for (int iy = 1; iy <= ny; iy++) {
05363 for (int ix = 0; ix < nxc; ix++) {
05364
05365 cmplx(ix,iy,iz) *= float(-2*((ix+iy+iz)%2) + 1);
05366 }
05367 }
05368 }
05369 }
05370 set_array_offsets(saved_offsets);
05371 update();
05372 EXITFUNC;
05373 }
05374
05375
05376 #define fint(i,j,k) fint[(i-1) + ((j-1) + (k-1)*ny)*lsd]
05377 #define fout(i,j,k) fout[(i-1) + ((j-1) + (k-1)*nyn)*lsdn]
05378 EMData *EMData::FourInterpol(int nxn, int nyni, int nzni, bool RetReal) {
05379
05380 int nyn, nzn, lsd, lsdn, inx, iny, inz;
05381 int i, j, k;
05382 if (is_complex())
05383 throw ImageFormatException("Input image has to be real");
05384
05385 if(ny > 1) {
05386 nyn = nyni;
05387 if(nz > 1) {
05388 nzn = nzni;
05389 } else {
05390 nzn = 1;
05391 }
05392 } else {
05393 nyn = 1; nzn = 1;
05394 }
05395 if(nxn<nx || nyn<ny || nzn<nz) throw ImageDimensionException("Cannot reduce the image size");
05396 lsd = nx + 2 - nx%2;
05397 lsdn = nxn + 2 - nxn%2;
05398
05399 EMData *temp_ft = do_fft();
05400 EMData *ret = this->copy();
05401 ret->set_size(lsdn, nyn, nzn);
05402 ret->to_zero();
05403 float *fout = ret->get_data();
05404 float *fint = temp_ft->get_data();
05405
05406
05407 float sq2 = 1.0f/std::sqrt(2.0f);
05408 float anorm = (float) nxn* (float) nyn* (float) nzn/(float) nx/ (float) ny/ (float) nz;
05409 for (i = 0; i < lsd*ny*nz; i++) fint[i] *= anorm;
05410 inx = nxn-nx; iny = nyn - ny; inz = nzn - nz;
05411 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);
05412 if(nyn>1) {
05413
05414 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);
05415 if(nzn>1) {
05416 for (k=nz/2+2+inz; k<=nzn; k++) {
05417 for (j=1; j<=ny/2+1; j++) {
05418 for (i=1; i<=lsd; i++) {
05419 fout(i,j,k)=fint(i,j,k-inz);
05420 }
05421 }
05422 for (j=ny/2+2+iny; j<=nyn; j++) {
05423 for (i=1; i<=lsd; i++) {
05424 fout(i,j,k)=fint(i,j-iny,k-inz);
05425 }
05426 }
05427 }
05428 }
05429 }
05430
05431
05432
05433 if(nx%2 == 0 && inx !=0) {
05434 for (k=1; k<=nzn; k++) {
05435 for (j=1; j<=nyn; j++) {
05436 fout(nx+1,j,k) *= sq2;
05437 fout(nx+2,j,k) *= sq2;
05438 }
05439 }
05440 if(nyn>1) {
05441 for (k=1; k<=nzn; k++) {
05442 for (i=1; i<=lsd; i++) {
05443 fout(i,ny/2+1+iny,k) = sq2*fout(i,ny/2+1,k);
05444 fout(i,ny/2+1,k) *= sq2;
05445 }
05446 }
05447 if(nzn>1) {
05448 for (j=1; j<=nyn; j++) {
05449 for (i=1; i<=lsd; i++) {
05450 fout(i,j,nz/2+1+inz) = sq2*fout(i,j,nz/2+1);
05451 fout(i,j,nz/2+1) *= sq2;
05452 }
05453 }
05454 }
05455 }
05456 }
05457 ret->set_complex(true);
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 ret->set_ri(1);
05492 ret->set_fftpad(true);
05493 ret->set_attr("npad", 1);
05494 if (nxn%2 == 1) {ret->set_fftodd(true);} else {ret->set_fftodd(false);}
05495 if(RetReal) {
05496 ret->do_ift_inplace();
05497 ret->depad();
05498 }
05499 ret->update();
05500
05501
05502
05503
05504
05505
05506
05507 delete temp_ft;
05508 temp_ft = 0;
05509 return ret;
05510 }
05511
05512 EMData *EMData::FourTruncate(int nxn, int nyni, int nzni, bool RetReal) {
05513
05514 int nyn, nzn, lsd, lsdn, inx, iny, inz;
05515 int i, j, k;
05516 float *fint;
05517 EMData *temp_ft = NULL;
05518
05519
05520
05521 if(ny > 1) {
05522 nyn = nyni;
05523 if(nz > 1) {
05524 nzn = nzni;
05525 } else {
05526 nzn = 1;
05527 }
05528 } else {
05529 nyn = 1; nzn = 1;
05530 }
05531 if (is_complex()) {
05532 nx = nx - 2 + nx%2;
05533 fint = get_data();
05534 } else {
05535
05536 temp_ft = do_fft();
05537 fint = temp_ft->get_data();
05538 }
05539 if(nxn>nx || nyn>ny || nzn>nz) throw ImageDimensionException("Cannot increase the image size");
05540 lsd = nx + 2 - nx%2;
05541 lsdn = nxn + 2 - nxn%2;
05542 EMData *ret = this->copy_head();
05543 ret->set_size(lsdn, nyn, nzn);
05544 float *fout = ret->get_data();
05545
05546
05547
05548 float anorm = (float) nxn* (float) nyn* (float) nzn/(float) nx/ (float) ny/ (float) nz;
05549 for (i = 0; i < lsd*ny*nz; i++) fint[i] *= anorm;
05550 inx = nx - nxn; iny = ny - nyn; inz = nz - nzn;
05551 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);
05552 if(nyn>1) {
05553 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);
05554 if(nzn>1) {
05555 for (k=nzn/2+2; k<=nzn; k++) {
05556 for (j=1; j<=nyn/2+1; j++) {
05557 for (i=1; i<=lsdn; i++) {
05558 fout(i,j,k)=fint(i,j,k+inz);
05559 }
05560 }
05561 for (j=nyn/2+2; j<=nyn; j++) {
05562 for (i=1; i<=lsdn; i++) {
05563 fout(i,j,k)=fint(i,j+iny,k+inz);
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 ret->set_complex(true);
05598 ret->set_ri(1);
05599 ret->set_fftpad(true);
05600 ret->set_attr("npad", 1);
05601 if (nxn%2 == 1) {ret->set_fftodd(true);} else {ret->set_fftodd(false);}
05602 if(RetReal) {
05603 ret->do_ift_inplace();
05604 ret->depad();
05605 }
05606 ret->update();
05607
05608
05609
05610
05611
05612
05613
05614 if (!is_complex()) {
05615 delete temp_ft;
05616 temp_ft = 0;
05617 }
05618 return ret;
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 EMData *EMData::Four_ds(int nxn, int nyni, int nzni, bool RetReal) {
05716
05717 int nyn, nzn, lsd, lsdn, inx, iny, inz;
05718 int i, j;
05719
05720 if(ny > 1) {
05721 nyn = nyni;
05722 if(nz > 1) {
05723 nzn = nzni;
05724 } else {
05725 nzn = 1;
05726 }
05727 } else {
05728 nyn = 1; nzn = 1;
05729 }
05730 lsd = nx-2 + 2 - nx%2;
05731 lsdn = nxn + 2 - nxn%2;
05732
05733 EMData *temp_ft = this->copy();
05734 EMData *ret = this->copy();
05735 ret->set_size(lsdn, nyn, nzn);
05736 ret->to_zero();
05737 float *fout = ret->get_data();
05738 float *fint = temp_ft->get_data();
05739
05740
05741
05742 float anorm = (float) nxn* (float) nyn* (float) nzn/(float) nx/ (float) ny/ (float) nz;
05743 for (i = 0; i < lsd*ny*nz; i++) fint[i] *= anorm;
05744 inx = nxn-(nx-2); iny = nyn - ny; inz = nzn - nz;
05745 for (j=1; j<=nyn; j++)
05746 for (i=1; i<=lsdn; i++)
05747 fout(i,j,1)=fint((i-1)/2*4+2-i%2,j*2-1,1);
05748 ret->set_complex(true);
05749 ret->set_ri(1);
05750
05751
05752 if (nxn%2 == 1) {ret->set_fftodd(true);} else {ret->set_fftodd(false);}
05753 if(RetReal) {
05754 ret->do_ift_inplace();
05755 ret->depad();
05756 }
05757 ret->update();
05758
05759 delete temp_ft;
05760 temp_ft = 0;
05761 return ret;
05762 }
05763
05764 EMData *EMData::Four_shuf_ds_cen_us(int nxn, int nyni, int, bool RetReal) {
05765
05766 int nyn, nzn, lsd, lsdn, inx, iny, inz;
05767 int i, j;
05768
05769 nyn = nyni;
05770 nzn = 1;
05771 lsd = nx;
05772 lsdn = nxn + 2 - nxn%2;
05773
05774 EMData *temp_ft = this->copy();
05775 EMData *ret = this->copy();
05776 ret->set_size(lsdn, nyn, nzn);
05777 ret->to_zero();
05778 float *fout = ret->get_data();
05779 float *fint = temp_ft->get_data();
05780
05781
05782 float sq2 = 1.0f/std::sqrt(2.0f);
05783
05784 for (i = 0; i < lsd*ny*nz; i++) fint[i] *= 4;
05785
05786 inx = nxn-(nx-2); iny = nyn - ny; inz = nzn - nz;
05787 for (j=1; j<=ny/4; j++)
05788 for (i=1; i<=(nx-2)/2+2; i++) {
05789 int g = (i-1)/2+1;
05790 if ((g+j)%2 == 0) {
05791 fout(i,j,1)=fint(g*4-2-i%2,j*2-1+ny/2,1);
05792 } else {
05793 fout(i,j,1)=-fint(g*4-2-i%2,j*2-1+ny/2,1);
05794 }
05795 }
05796
05797 for (j=ny/4+1; j<=ny/4+1; j++)
05798 for (i=1; i<=(nx-2)/2+2; i++) {
05799 int g = (i-1)/2+1;
05800 if ((g+j)%2 == 0) {
05801 fout(i,j,1)=fint(g*4-2-i%2,j*2-1-ny/2,1);
05802 } else {
05803 fout(i,j,1)=-fint(g*4-2-i%2,j*2-1-ny/2,1);
05804 }
05805 }
05806
05807 for (j=ny/4+2; j<=ny/2; j++)
05808 for (i=1; i<=(nx-2)/2+2; i++) {
05809 int g = (i-1)/2+1;
05810 if ((g+j)%2 == 0) {
05811 fout(i,j+ny/2,1)=fint(g*4-2-i%2,j*2-1-ny/2,1);
05812 } else {
05813 fout(i,j+ny/2,1)=-fint(g*4-2-i%2,j*2-1-ny/2,1);
05814 }
05815 }
05816
05817 if (nx%2 == 0) {
05818 for (j=1; j<=nyn; j++) {
05819 fout((nx-2)/2+1,j,1) *= sq2;
05820 fout((nx-2)/2+2,j,1) *= sq2;
05821 }
05822 for (i=1; i<=lsd/2+1; i++) {
05823 fout(i,ny/4+1+ny/2,1) = sq2*fout(i,ny/4+1,1);
05824 fout(i,ny/4+1,1) *= sq2;
05825 }
05826 }
05827
05828 ret->set_complex(true);
05829 ret->set_ri(1);
05830
05831 if (nxn%2 == 1) {ret->set_fftodd(true);} else {ret->set_fftodd(false);}
05832 if(RetReal) {
05833 ret->do_ift_inplace();
05834 ret->depad();
05835 }
05836 ret->update();
05837
05838 delete temp_ft;
05839 temp_ft = 0;
05840 return ret;
05841 }
05842
05843 #undef fint
05844 #undef fout
05845
05846 #define fint(jx,jy,jz) fint[jx + (jy + jz*ny)*nox]
05847 EMData *EMData::filter_by_image(EMData* image, bool RetReal) {
05848
05849
05850 bool complex_input = this->is_complex();
05851 nx = this->get_xsize();
05852 ny = this->get_ysize();
05853 nz = this->get_zsize();
05854 int nox;
05855 if (complex_input) nox = (nx - 2 + this->is_fftodd()); else nox = nx;
05856
05857 int lsd2 = (nox + 2 - nox%2) / 2;
05858
05859 EMData* fp = NULL;
05860 if(complex_input) {
05861
05862 fp = this->copy();
05863 } else {
05864 fp = this->norm_pad( false, 1);
05865 fp->do_fft_inplace();
05866 }
05867 fp->set_array_offsets(1,1,1);
05868 int nx2 = nox/2;
05869 int ny2 = ny/2;
05870 int nz2 = nz/2;
05871 float *fint = image->get_data();
05872 for ( int iz = 1; iz <= nz; iz++) {
05873 int jz=nz2-iz+1; if(jz<0) jz += nz;
05874 for ( int iy = 1; iy <= ny; iy++) {
05875 int jy=ny2-iy+1; if(jy<0) jy += ny;
05876 for ( int ix = 1; ix <= lsd2; ix++) {
05877 int jx = nx2-ix+1;
05878 fp->cmplx(ix,iy,iz) *= fint(jx,jy,jz);
05879 }
05880 }
05881 }
05882
05883 fp->set_ri(1);
05884 fp->set_fftpad(true);
05885 fp->set_attr("npad", 1);
05886 if (nx%2 == 1) fp->set_fftodd(true);
05887 else fp->set_fftodd(false);
05888 if(RetReal) {
05889 fp->do_ift_inplace();
05890 fp->depad();
05891 }
05892 fp->set_array_offsets(0,0,0);
05893 fp->update();
05894
05895 return fp;
05896 }
05897 #undef fint
05898 #define fint(jx,jy,jz) fint[jx + (jy + jz*ny)*nx]
05899 #define fout(jx,jy,jz) fout[jx + (jy + jz*ny)*nx]
05900 EMData *EMData::replace_amplitudes(EMData* image, bool RetReal) {
05901
05902
05903 bool complex_input = this->is_complex();
05904 nx = this->get_xsize();
05905 ny = this->get_ysize();
05906 nz = this->get_zsize();
05907 int nox;
05908 if (complex_input) nox = (nx - 2 + this->is_fftodd()); else nox = nx;
05909
05910 EMData* fp = NULL;
05911 if(complex_input) {
05912
05913 fp = this->copy();
05914 } else {
05915 fp = this->norm_pad( false, 1);
05916 fp->do_fft_inplace();
05917 }
05918 float *fout = fp->get_data();
05919 float *fint = image->get_data();
05920 for ( int iz = 0; iz < nz; iz++) {
05921 for ( int iy = 0; iy < ny; iy++) {
05922 for ( int ix = 0; ix < nx; ix+=2) {
05923 float qt = fint(ix,iy,iz)*fint(ix,iy,iz)+fint(ix+1,iy,iz)*fint(ix+1,iy,iz);
05924 float rt = fout(ix,iy,iz)*fout(ix,iy,iz)+fout(ix+1,iy,iz)*fout(ix+1,iy,iz);
05925 if(rt > 1.0e-20) {
05926 fout(ix,iy,iz) *= (qt/rt);
05927 fout(ix+1,iy,iz) *= (qt/rt);
05928 } else {
05929 qt = std::sqrt(qt/2.0f);
05930 fout(ix,iy,iz) = qt;
05931 fout(ix+1,iy,iz) = qt;
05932 }
05933 }
05934 }
05935 }
05936
05937 fp->set_ri(1);
05938 fp->set_fftpad(true);
05939 fp->set_attr("npad", 1);
05940 if (nx%2 == 1) fp->set_fftodd(true);
05941 else fp->set_fftodd(false);
05942 if(RetReal) {
05943 fp->do_ift_inplace();
05944 fp->depad();
05945 }
05946 fp->set_array_offsets(0,0,0);
05947 fp->update();
05948
05949 return fp;
05950 }
05951 #undef fint
05952 #undef fout
05953
05954
05955 #undef QUADPI
05956 #undef DGR_TO_RAD