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
03341 if (!is_shuffled()) fft_shuffle();
03342
03343 EMData* result = copy_head();
03344 set_array_offsets(0,-nyhalf);
03345 result->set_array_offsets(0,-nyhalf);
03346
03347 ang = ang*(float)DGR_TO_RAD;
03348 float cang = cos(ang);
03349 float sang = sin(ang);
03350 for (int iy = -nyhalf; iy < nyhalf; iy++) {
03351 float ycang = iy*cang;
03352 float ysang = iy*sang;
03353 for (int ix = 0; ix <= nxhalf; ix++) {
03354 float nuxold = (ix*cang - ysang)*scale;
03355 float nuyold = (ix*sang + ycang)*scale;
03356 result->cmplx(ix,iy) = Util::extractpoint2(nx, ny, nuxold, nuyold, this, kb);
03357
03358 }
03359 }
03360 result->set_array_offsets();
03361 result->fft_shuffle();
03362 result->update();
03363 set_array_offsets();
03364 fft_shuffle();
03365 return result;
03366 }
03367
03368 EMData* EMData::fouriergridrot_shift2d(float ang, float sx, float sy, Util::KaiserBessel& kb) {
03369 if (2 != get_ndim())
03370 throw ImageDimensionException("fouriergridrot_shift2d needs a 2-D image.");
03371 if (!is_complex())
03372 throw ImageFormatException("fouriergridrot_shift2d requires a fourier image");
03373 int nxreal = nx - 2 + int(is_fftodd());
03374 if (nxreal != ny)
03375 throw ImageDimensionException("fouriergridrot_shift2d requires ny == nx(real)");
03376 if (0 != nxreal%2)
03377 throw ImageDimensionException("fouriergridrot_shift2d needs an even image.");
03378 int nxhalf = nxreal/2;
03379 int nyhalf = ny/2;
03380
03381 if (!is_shuffled()) fft_shuffle();
03382
03383 EMData* result = copy_head();
03384 set_array_offsets(0, -nyhalf);
03385 result->set_array_offsets(0, -nyhalf);
03386
03387 ang = ang*(float)DGR_TO_RAD;
03388 float cang = cos(ang);
03389 float sang = sin(ang);
03390 float temp = -2.0f*M_PI/nxreal;
03391 for (int iy = -nyhalf; iy < nyhalf; iy++) {
03392 float ycang = iy*cang;
03393 float ysang = iy*sang;
03394 for (int ix = 0; ix <= nxhalf; ix++) {
03395 float nuxold = ix*cang - ysang;
03396 float nuyold = ix*sang + ycang;
03397 result->cmplx(ix,iy) = Util::extractpoint2(nx, ny, nuxold, nuyold, this, kb);
03398
03399 float phase_ang = temp*(sx*ix+sy*iy);
03400 result->cmplx(ix,iy) *= complex<float>(cos(phase_ang), sin(phase_ang));
03401 }
03402 }
03403 result->set_array_offsets();
03404 result->fft_shuffle();
03405 result->update();
03406 set_array_offsets();
03407 fft_shuffle();
03408 return result;
03409 }
03410
03411 #undef QUADPI
03412 #undef DGR_TO_RAD
03413
03414 void EMData::divkbsinh(const Util::KaiserBessel& kb) {
03415 if (is_complex())
03416 throw ImageFormatException("divkbsinh requires a real image.");
03417 vector<int> saved_offsets = get_array_offsets();
03418 set_array_offsets(0,0,0);
03419
03420
03421
03422
03423 for (int iz=0; iz < nz; iz++) {
03424 float wz = kb.sinhwin(static_cast<float>(iz-nz/2));
03425 for (int iy=0; iy < ny; iy++) {
03426 float wy = kb.sinhwin(static_cast<float>(iy-ny/2));
03427 for (int ix=0; ix < nx; ix++) {
03428 float wx = kb.sinhwin(static_cast<float>(ix-nx/2));
03429 float w = wx*wy*wz;
03430 (*this)(ix,iy,iz) /= w;
03431 }
03432 }
03433 }
03434 set_array_offsets(saved_offsets);
03435 }
03436
03437
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 EMData* EMData::extract_plane(const Transform& tf, Util::KaiserBessel& kb) {
03463 if (!is_complex())
03464 throw ImageFormatException("extractplane requires a complex image");
03465 if (nx%2 != 0)
03466 throw ImageDimensionException("extractplane requires nx to be even");
03467 int nxreal = nx - 2;
03468 if (nxreal != ny || nxreal != nz)
03469 throw ImageDimensionException("extractplane requires ny == nx == nz");
03470
03471 EMData* res = new EMData();
03472 res->set_size(nx,ny,1);
03473 res->to_zero();
03474 res->set_complex(true);
03475 res->set_fftodd(false);
03476 res->set_fftpad(true);
03477 res->set_ri(true);
03478
03479 int n = nxreal;
03480 int nhalf = n/2;
03481 vector<int> saved_offsets = get_array_offsets();
03482 set_array_offsets(0,-nhalf,-nhalf);
03483 res->set_array_offsets(0,-nhalf,0);
03484
03485 int kbsize = kb.get_window_size();
03486 int kbmin = -kbsize/2;
03487 int kbmax = -kbmin;
03488 float* wy0 = new float[kbmax - kbmin + 1];
03489 float* wy = wy0 - kbmin;
03490 float* wx0 = new float[kbmax - kbmin + 1];
03491 float* wx = wx0 - kbmin;
03492 float* wz0 = new float[kbmax - kbmin + 1];
03493 float* wz = wz0 - kbmin;
03494 float rim = nhalf*float(nhalf);
03495 int count = 0;
03496 float wsum = 0.f;
03497 Transform tftrans = tf;
03498 tftrans.invert();
03499 for (int jy = -nhalf; jy < nhalf; jy++) {
03500 for (int jx = 0; jx <= nhalf; jx++) {
03501 Vec3f nucur((float)jx, (float)jy, 0.f);
03502 Vec3f nunew = tftrans*nucur;
03503 float xnew = nunew[0], ynew = nunew[1], znew = nunew[2];
03504 if (xnew*xnew+ynew*ynew+znew*znew <= rim) {
03505 count++;
03506 std::complex<float> btq(0.f,0.f);
03507 bool flip = false;
03508 if (xnew < 0.f) {
03509 flip = true;
03510 xnew = -xnew;
03511 ynew = -ynew;
03512 znew = -znew;
03513 }
03514 int ixn = int(Util::round(xnew));
03515 int iyn = int(Util::round(ynew));
03516 int izn = int(Util::round(znew));
03517
03518 for (int i=kbmin; i <= kbmax; i++) {
03519 int izp = izn + i;
03520 wz[i] = kb.i0win_tab(znew - izp);
03521 int iyp = iyn + i;
03522 wy[i] = kb.i0win_tab(ynew - iyp);
03523 int ixp = ixn + i;
03524 wx[i] = kb.i0win_tab(xnew - ixp);
03525
03526 }
03527
03528 int lnbz = 0;
03529 for (int iz = kbmin; iz <= -1; iz++) {
03530 if (wz[iz] != 0.f) {
03531 lnbz = iz;
03532 break;
03533 }
03534 }
03535 int lnez = 0;
03536 for (int iz = kbmax; iz >= 1; iz--) {
03537 if (wz[iz] != 0.f) {
03538 lnez = iz;
03539 break;
03540 }
03541 }
03542 int lnby = 0;
03543 for (int iy = kbmin; iy <= -1; iy++) {
03544 if (wy[iy] != 0.f) {
03545 lnby = iy;
03546 break;
03547 }
03548 }
03549 int lney = 0;
03550 for (int iy = kbmax; iy >= 1; iy--) {
03551 if (wy[iy] != 0.f) {
03552 lney = iy;
03553 break;
03554 }
03555 }
03556 int lnbx = 0;
03557 for (int ix = kbmin; ix <= -1; ix++) {
03558 if (wx[ix] != 0.f) {
03559 lnbx = ix;
03560 break;
03561 }
03562 }
03563 int lnex = 0;
03564 for (int ix = kbmax; ix >= 1; ix--) {
03565 if (wx[ix] != 0.f) {
03566 lnex = ix;
03567 break;
03568 }
03569 }
03570 if (ixn >= -kbmin && ixn <= nhalf-1-kbmax
03571 && iyn >= -nhalf-kbmin && iyn <= nhalf-1-kbmax
03572 && izn >= -nhalf-kbmin && izn <= nhalf-1-kbmax) {
03573
03574 for (int lz = lnbz; lz <= lnez; lz++) {
03575 int izp = izn + lz;
03576 for (int ly=lnby; ly<=lney; ly++) {
03577 int iyp = iyn + ly;
03578 float ty = wz[lz]*wy[ly];
03579 for (int lx=lnbx; lx<=lnex; lx++) {
03580 int ixp = ixn + lx;
03581 float wg = wx[lx]*ty;
03582 btq += cmplx(ixp,iyp,izp)*wg;
03583 wsum += wg;
03584 }
03585 }
03586 }
03587 } else {
03588
03589 for (int lz = lnbz; lz <= lnez; lz++) {
03590 int izp = izn + lz;
03591 for (int ly=lnby; ly<=lney; ly++) {
03592 int iyp = iyn + ly;
03593 float ty = wz[lz]*wy[ly];
03594 for (int lx=lnbx; lx<=lnex; lx++) {
03595 int ixp = ixn + lx;
03596 float wg = wx[lx]*ty;
03597 bool mirror = false;
03598 int ixt(ixp), iyt(iyp), izt(izp);
03599 if (ixt > nhalf || ixt < -nhalf) {
03600 ixt = Util::sgn(ixt)
03601 *(n - abs(ixt));
03602 iyt = -iyt;
03603 izt = -izt;
03604 mirror = !mirror;
03605 }
03606 if (iyt >= nhalf || iyt < -nhalf) {
03607 if (ixt != 0) {
03608 ixt = -ixt;
03609 iyt = Util::sgn(iyt)
03610 *(n - abs(iyt));
03611 izt = -izt;
03612 mirror = !mirror;
03613 } else {
03614 iyt -= n*Util::sgn(iyt);
03615 }
03616 }
03617 if (izt >= nhalf || izt < -nhalf) {
03618 if (ixt != 0) {
03619 ixt = -ixt;
03620 iyt = -iyt;
03621 izt = Util::sgn(izt)
03622 *(n - abs(izt));
03623 mirror = !mirror;
03624 } else {
03625 izt -= Util::sgn(izt)*n;
03626 }
03627 }
03628 if (ixt < 0) {
03629 ixt = -ixt;
03630 iyt = -iyt;
03631 izt = -izt;
03632 mirror = !mirror;
03633 }
03634 if (iyt == nhalf) iyt = -nhalf;
03635 if (izt == nhalf) izt = -nhalf;
03636 if (mirror) btq += conj(cmplx(ixt,iyt,izt))*wg;
03637 else btq += cmplx(ixt,iyt,izt)*wg;
03638 wsum += wg;
03639 }
03640 }
03641 }
03642 }
03643