Commit 2de9b5dc authored by David Rohr's avatar David Rohr
Browse files

fix sum and sumsquare calculation: for the non FFT algorithm sum was...

fix sum and sumsquare calculation: for the non FFT algorithm sum was calculated twice before (stupid mistake), for the FFT algorithm there is actually a dependency for the sumsquare calculation on whether NumberPixels is odd or even. Although, the effect is very small for normal maps since high order fourier coefficients are small then. Anyhow, should be correct now
parent f1e3f10a
...@@ -761,15 +761,22 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj ...@@ -761,15 +761,22 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj
sumsquareC = 0; sumsquareC = 0;
if (FFTAlgo) if (FFTAlgo)
{ {
int jloopend = param.param_device.NumberFFTPixels1D;
if ((param.param_device.NumberPixels & 1) == 0) jloopend--;
for(int i = 0; i < param.param_device.NumberPixels; i++) for(int i = 0; i < param.param_device.NumberPixels; i++)
{ {
for (int j = 1;j < param.param_device.NumberFFTPixels1D;j++) for (int j = 1;j < jloopend;j++)
{ {
int k = i * param.param_device.NumberFFTPixels1D + j; int k = i * param.param_device.NumberFFTPixels1D + j;
sumsquareC += (localmultFFT[k][0] * localmultFFT[k][0] + localmultFFT[k][1] * localmultFFT[k][1]) * 2; sumsquareC += (localmultFFT[k][0] * localmultFFT[k][0] + localmultFFT[k][1] * localmultFFT[k][1]) * 2;
} }
int k = i * param.param_device.NumberFFTPixels1D; int k = i * param.param_device.NumberFFTPixels1D;
sumsquareC += localmultFFT[k][0] * localmultFFT[k][0] + localmultFFT[k][1] * localmultFFT[k][1]; sumsquareC += localmultFFT[k][0] * localmultFFT[k][0] + localmultFFT[k][1] * localmultFFT[k][1];
if ((param.param_device.NumberPixels & 1) == 0)
{
k += param.param_device.NumberFFTPixels1D - 2;
sumsquareC += localmultFFT[k][0] * localmultFFT[k][0] + localmultFFT[k][1] * localmultFFT[k][1];
}
} }
myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels); myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
...@@ -786,7 +793,6 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj ...@@ -786,7 +793,6 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj
for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++) for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++)
{ {
sumsquareC += Mapconv[i] * Mapconv[i]; sumsquareC += Mapconv[i] * Mapconv[i];
sumC += Mapconv[i];
} }
myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels); myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
......
Markdown is supported
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