Skip to content

GitLab

Commit 8f4d17a7 authored by David Rohr's avatar David Rohr
Browse files

replace 2d loops by 1d loops where possible

parent 82c718f2
Hide whitespace changes
Inline Side-by-side
bioem.cpp
View file @ 8f4d17a7
......@@ -402,13 +402,10 @@ int bioem::compareRefMaps2(int iOrient, int iConv, mycomplex_t* localConvFFT,myf
/////////////NEW ROUTINE ////////////////
inline int bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC,myfloat_t sumsquareC, mycomplex_t* localConvFFT,mycomplex_t* localCCT,myfloat_t* lCC)
{
for(int i=0; i < param.param_device.NumberPixels ; i++ )
for(int i = 0;i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D;i++)
{
for(int j=0; j < param.param_device.NumberFFTPixels1D ; j++ )
{
localCCT[i*param.param_device.NumberFFTPixels1D+j][0]=localConvFFT[i*param.param_device.NumberFFTPixels1D+j][0]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberFFTPixels1D+j][0]+localConvFFT[i*param.param_device.NumberFFTPixels1D+j][1]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberFFTPixels1D+j][1];
localCCT[i*param.param_device.NumberFFTPixels1D+j][1]=localConvFFT[i*param.param_device.NumberFFTPixels1D+j][1]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberFFTPixels1D+j][0]-localConvFFT[i*param.param_device.NumberFFTPixels1D+j][0]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberFFTPixels1D+j][1];
}
localCCT[i][0] = localConvFFT[i][0] * RefMap.RefMapFFT[iRefMap].cpoints[i][0] + localConvFFT[i][1] * RefMap.RefMapFFT[iRefMap].cpoints[i][1];
localCCT[i][1] = localConvFFT[i][1] * RefMap.RefMapFFT[iRefMap].cpoints[i][0] - localConvFFT[i][0] * RefMap.RefMapFFT[iRefMap].cpoints[i][1];
}
myfftw_execute_dft_c2r(param.fft_plan_c2r_backward,localCCT,lCC);
......@@ -592,14 +589,11 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
/**** Multiplying FFTmap with corresponding kernel ****/
for(int i=0; i < param.param_device.NumberPixels ; i++ )
for(int i=0;i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D;i++)
{
for(int j=0; j < param.param_device.NumberFFTPixels1D ; j++ )
{ //Projection*CONJ(KERNEL)
localmultFFT[i*param.param_device.NumberFFTPixels1D+j][0]=lproj[i*param.param_device.NumberFFTPixels1D+j][0]*param.refCTF[iConv].cpoints[i*param.param_device.NumberFFTPixels1D+j][0]+lproj[i*param.param_device.NumberFFTPixels1D+j][1]*param.refCTF[iConv].cpoints[i*param.param_device.NumberFFTPixels1D+j][1];
localmultFFT[i*param.param_device.NumberFFTPixels1D+j][1]=lproj[i*param.param_device.NumberFFTPixels1D+j][1]*param.refCTF[iConv].cpoints[i*param.param_device.NumberFFTPixels1D+j][0]-lproj[i*param.param_device.NumberFFTPixels1D+j][0]*param.refCTF[iConv].cpoints[i*param.param_device.NumberFFTPixels1D+j][1];
// cout << "GG " << i << " " << j << " " << param.refCTF[iConv].cpoints[i*param.param_device.NumberFFTPixels1D+j][0] << " " <<param.refCTF[iConv].cpoints[i*param.param_device.NumberFFTPixels1D+j][1] <<" " <<lproj[i*param.param_device.NumberFFTPixels1D+j][0] <<" " <<lproj[i*param.param_device.NumberFFTPixels1D+j][1] << "\n";
}
localmultFFT[i][0] = lproj[i][0] * param.refCTF[iConv].cpoints[i][0] + lproj[i][1] * param.refCTF[iConv].cpoints[i][1];
localmultFFT[i][1] = lproj[i][1] * param.refCTF[iConv].cpoints[i][0] - lproj[i][0] * param.refCTF[iConv].cpoints[i][1];
// cout << "GG " << i << " " << j << " " << param.refCTF[iConv].cpoints[i][0] << " " <<param.refCTF[iConv].cpoints[i][1] <<" " <<lproj[i][0] <<" " <<lproj[i][1] << "\n";
}
//FFTW_C2R will destroy the input array, so we have to work on a copy here
......@@ -620,18 +614,17 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
/*** Calculating Cross-correlations of cal-convoluted map with its self *****/
sumC=0;
sumsquareC=0;
for(int i=0; i < param.param_device.NumberPixels ; i++ )
for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++)
{
for(int j=0; j < param.param_device.NumberPixels ; j++ )
{
sumC+=localconvFFT[i*param.param_device.NumberPixels+j];
sumsquareC+=localconvFFT[i*param.param_device.NumberPixels+j]*localconvFFT[i*param.param_device.NumberPixels+j];
}
sumC += localconvFFT[i];
sumsquareC += localconvFFT[i] * localconvFFT[i];
}
/*** The DTF gives an unnormalized value so have to divded by the total number of pixels in Fourier ***/
// Normalizing
sumC=sumC/ (myfloat_t) (param.param_device.NumberPixels*param.param_device.NumberPixels);
sumsquareC=sumsquareC / pow((myfloat_t) param.param_device.NumberPixels,4);
myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
myfloat_t norm4 = norm2 * norm2;
sumC = sumC / norm2;
sumsquareC = sumsquareC / norm4;
/**** Freeing fftw_complex created (dont know if omp critical is necessary) ****/
myfftw_free(localconvFFT);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment