Skip to content

Commit

Permalink
add 3D OTF support and correct argument error (#30)
Browse files Browse the repository at this point in the history 
* add 3D OTF support and correct argument error

* correct arguments in RichardsonLucy_GPU()

* correct arguments

* Update linearDecon.h

* Update linearDecon.h

---------

Co-authored-by: Talley Lambert <[email protected]>
  • Loading branch information
@zichenzachwang @tlambert03
zichenzachwang and tlambert03 authored Jul 18, 2024
1 parent c84af91 commit e16428c
Show file tree
Hide file tree
Showing 2 changed files with 88 additions and 23 deletions.
2 changes: 1 addition & 1 deletion src/RL-Biggs-Andrews.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -690,7 +690,7 @@ int RL_interface(const unsigned short * const raw_data,
RichardsonLucy_GPU(raw_image, background, d_interpOTF, nIters,
deskewFactor, deskewedXdim, extraShift, napodize, nZblend, rotMatrix,
rfftplanGPU, rfftplanInvGPU, raw_deskewed, &deviceProp,
bFlatStartGuess, my_median, bDoRescale, padVal, bDupRevStack, bSkewedDecon, false);
bFlatStartGuess, my_median, bDoRescale, bSkewedDecon, padVal, bDupRevStack, false, 0);

// Copy deconvolved data, stored in raw_image, to "result" for return:
memcpy(result, raw_image.data(), raw_image.size() * sizeof(float));
Expand Down
109 changes: 87 additions & 22 deletions src/radialft_interface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,9 @@ void cleanup(std::complex<float> *otfkxkz, int nx, int nz, float dkr, float dkz,
void radialft(std::complex<float> *bands, int nx, int ny, int nz, std::complex<float> *avg);

bool fixorigin(std::complex<float> *otfkxkz, int nx, int nz, int kx2);
void fixorigin(CImg<> &otf3d, int kx2);
void rescale(std::complex<float> *otfkxkz, int nx, int nz);
void rescale(CImg<> &otf3d);


#ifdef _WIN32
Expand All @@ -52,7 +54,7 @@ RADIALFT_API int makeOTF(const char *const ifiles, const char *const ofiles,
int lambdanm = 520, float dz = 0.102, int interpkr = 10,
bool bUserBackground = false, float background = 90,
float NA = 1.25, float NIMM = 1.3, float dr = 0.102,
int krmax = 0, bool bDoCleanup = false);
int krmax = 0, bool bDoCleanup = false, bool b3Dout = false);
}


Expand Down Expand Up @@ -92,7 +94,7 @@ fftwf_plan rfftplan3d;

int makeOTF(const char *const ifiles, const char *const ofiles, int lambdanm,
float dz, int interpkr, bool bUserBackground, float background,
float NA, float NIMM, float dr, int krmax, bool bDoCleanup)
float NA, float NIMM, float dr, int krmax, bool bDoCleanup, bool b3Dout)
{
printf("called");

Expand All @@ -109,7 +111,8 @@ int makeOTF(const char *const ifiles, const char *const ofiles, int lambdanm,
dkr = 1/(ny*dr);
dkz = 1/(nz*dz);

nxy=(nx+2)*ny;
int nx2 = (nx/2+1)*2; // in case nx is odd numbered, nx2 is nx+1, not nx+2
nxy=nx2*ny;

floatimage = (float *) malloc(nxy*nz*sizeof(float));
bands = (std::complex<float> *) floatimage;
Expand All @@ -119,7 +122,7 @@ int makeOTF(const char *const ifiles, const char *const ofiles, int lambdanm,
for(z=0; z<nz; z++)
for (i=0; i<ny; i++) {
for (j=0; j<nx; j++)
floatimage[z*nxy+i*(nx+2)+j] = rawtiff(j, i, z);
floatimage[z*nxy+i*nx2+j] = rawtiff(j, i, z);
}

/* Before FFT, estimate bead center position */
Expand All @@ -135,7 +138,7 @@ int makeOTF(const char *const ifiles, const char *const ofiles, int lambdanm,
for(z=0; z<nz; z++) {
for(i=0;i<ny;i++)
for(j=0;j<nx;j++)
floatimage[z*nxy + i*(nx+2) + j] -= background;
floatimage[z*nxy + i*nx2 + j] -= background;
}

rfftplan3d = fftwf_plan_dft_r2c_3d(nz, ny, nx, floatimage,
Expand All @@ -155,6 +158,13 @@ int makeOTF(const char *const ifiles, const char *const ofiles, int lambdanm,
shift_center(bands, nx, ny, nz, xcofm, ycofm, zcofm);

CImg<> output_tiff(nz*2, nx/2+1, 1, 1, 0.f);
if (b3Dout) {
output_tiff.assign(floatimage, nx2, ny, nz, 1, true);
if (interpkr>1)
fixorigin(output_tiff, interpkr);
rescale(output_tiff);
}
else {
avg_output = (std::complex<float> *) output_tiff.data();

radialft(bands, nx, ny, nz, avg_output);
Expand All @@ -176,7 +186,7 @@ int makeOTF(const char *const ifiles, const char *const ofiles, int lambdanm,

// printf("%d\n", interpkr);
rescale(avg_output, nx, nz);

}
/* For side bands, combine bandre's and bandim's into bandplus */
/* Shouldn't this be done later on the averaged bands? */

Expand All @@ -188,26 +198,27 @@ int makeOTF(const char *const ifiles, const char *const ofiles, int lambdanm,
/* locate peak pixel to subpixel accuracy by fitting parabolas */
void determine_center_and_background(float *stack3d, int nx, int ny, int nz, float *xc, float *yc, float *zc, float *background)
{
int i, j, k, maxi, maxj, maxk, ind, nxy2, infocus_sec;
int i, j, k, maxi=0, maxj=0, maxk=0, ind, nxy2, infocus_sec;
int iminus, iplus, jminus, jplus, kminus, kplus;
float maxval, reval, valminus, valplus;
double sum;

//printf("In determine_center_and_background()\n");
nxy2 = (nx+2)*ny;
int nx2 = (nx/2+1)*2; // in case nx is odd numbered, nx2 is nx+1, not nx+2
nxy2 = nx2*ny;

/* Search for the peak pixel */
/* Be aware that stack3d is of dimension (nx+2)xnyxnz */
maxval=0.0;
for(k=0;k<nz;k++)
for(i=0;i<ny;i++)
for(j=0;j<nx;j++) {
ind=k*nxy2+i*(nx+2)+j;
reval=stack3d[ind];
if( reval > maxval ) {
maxval = reval;
maxi=i; maxj=j;
maxk=k;
ind=k*nxy2+i*nx2+j;
reval=stack3d[ind];
if( reval > maxval ) {
maxval = reval;
maxi=i; maxj=j;
maxk=k;
}
}

Expand All @@ -221,25 +232,25 @@ void determine_center_and_background(float *stack3d, int nx, int ny, int nz, flo
if( kminus<0 ) kminus+=nz;
if( kplus>=nz ) kplus-=nz;

valminus = stack3d[kminus*nxy2+maxi*(nx+2)+maxj];
valplus = stack3d[kplus *nxy2+maxi*(nx+2)+maxj];
valminus = stack3d[kminus*nxy2+maxi*nx2+maxj];
valplus = stack3d[kplus *nxy2+maxi*nx2+maxj];
*zc = maxk + fitparabola(valminus, maxval, valplus);

*zc += 0.6;

valminus = stack3d[maxk*nxy2+iminus*(nx+2)+maxj];
valplus = stack3d[maxk*nxy2+iplus *(nx+2)+maxj];
valminus = stack3d[maxk*nxy2+iminus*nx2+maxj];
valplus = stack3d[maxk*nxy2+iplus *nx2+maxj];
*yc = maxi + fitparabola(valminus, maxval, valplus);

valminus = stack3d[maxk*nxy2+maxi*(nx+2)+jminus];
valplus = stack3d[maxk*nxy2+maxi*(nx+2)+jplus];
valminus = stack3d[maxk*nxy2+maxi*nx2+jminus];
valplus = stack3d[maxk*nxy2+maxi*nx2+jplus];
*xc = maxj + fitparabola(valminus, maxval, valplus);

sum = 0;
infocus_sec = floor(*zc);
for (i=0; i<*yc-20; i++)
for (j=0; j<nx; j++)
sum += stack3d[infocus_sec*nxy2 + i*(nx+2) + j];
sum += stack3d[infocus_sec*nxy2 + i*nx2 + j];
*background = sum / ((*yc-20)*nx);
}

Expand Down Expand Up @@ -431,6 +442,49 @@ bool fixorigin(std::complex<float> *otfkxkz, int nx, int nz, int kx2)
return true;
}

void fixorigin(CImg<> &otf3d, int kx2)
{
// Project 3D OTF onto kx-ky plane for the +/- 45-degree radial lines
int ny = otf3d.height();
int nx = otf3d.width()/2; // real and imaginary

CImg<> projected(nx*2, 1, 1, 1, 0.);

for (int z=0; z<otf3d.depth(); z++)
for (int x=1; x<nx; x++) {
int xst = x<<1;
// +45 degree line:
projected(xst) += otf3d(xst, x, z);
projected(xst+1) += otf3d(xst+1, x, z);
// -45 degree line:
projected(xst) += otf3d(xst, ny-x, z);
projected(xst+1) += otf3d(xst+1, ny-x, z);
}

//
std::complex<double> mean_val=0;
std::complex<double> slope_numerator=0;
std::complex<double> slope_denominat=0;
std::complex<float> * projectedC = (std::complex<float> *) projected.data();

for (int x=1; x<=kx2; x++)
mean_val += projectedC[x];
mean_val /= kx2;

double mean_kx = (kx2+1)/2.; // the mean of [1, 2, ..., n] is (n+1)/2
for (int x=1; x<=kx2; x++) {
std::complex<double> complexval = projectedC[x];
slope_numerator += (x - mean_kx) * (complexval - mean_val);
slope_denominat += (x - mean_kx) * (x - mean_kx);
}
std::complex<double> slope = slope_numerator / slope_denominat;
std::complex<double> valOrigin = mean_val - slope * mean_kx; // intercept at kx=0
std::cout << "otf3d(0) = " << otf3d(0, 0, 0) << "+i" << otf3d(1,0,0) << std::endl;
otf3d(0, 0, 0) = valOrigin.real();
otf3d(1, 0, 0) = valOrigin.imag();
std::cout << "otf3d(0) = " << otf3d(0, 0, 0) << "+i" << otf3d(1,0,0) << std::endl;
}

void rescale(std::complex<float> *otfkxkz, int nx, int nz)
{
int nxz, ind;
Expand All @@ -449,5 +503,16 @@ void rescale(std::complex<float> *otfkxkz, int nx, int nz)
}
}

void rescale(CImg<> &otf3d)
{
std::complex<float> * otf3dC = (std::complex<float> *) otf3d.data();
size_t nxyz = otf3d.size() / 2;
float valmax = 0., mag;
for (unsigned i=0; i<nxyz; i++) {
mag = std::abs(otf3dC[i]);
if (mag > valmax)
valmax = mag;
}


otf3d /= valmax;
}

0 comments on commit e16428c

Please sign in to comment.