Skip to content

Commit

Permalink
add bounds and fix bugs in find_cenerCol and find_tau
Browse files Browse the repository at this point in the history
  • Loading branch information
@kylechampley
kylechampley committed Nov 15, 2024
1 parent c68eead commit 2eaf458
Show file tree
Hide file tree
Showing 7 changed files with 80 additions and 49 deletions.
67 changes: 44 additions & 23 deletions src/find_center_cpu.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -318,7 +318,7 @@ float interpolate2D(float* data, float ind_1, float ind_2, int N_1, int N_2)
d_1 * ((1.0 - d_2) * data[ind_1_hi * N_2 + ind_2_lo] + d_2 * data[ind_1_hi * N_2 + ind_2_hi]);
}

float findCenter_cpu(float* g, parameters* params, int iRow, bool find_tau)
float findCenter_cpu(float* g, parameters* params, int iRow, bool find_tau, float* searchBounds)
{
if (g == NULL || params == NULL)
return 0.0;
Expand All @@ -337,15 +337,15 @@ float findCenter_cpu(float* g, parameters* params, int iRow, bool find_tau)
return 0.0;
}
else
return findCenter_parallel_cpu(g, params, iRow);
return findCenter_parallel_cpu(g, params, iRow, searchBounds);
}
else if (params->geometry == parameters::FAN)
return findCenter_fan_cpu(g, params, iRow, find_tau);
return findCenter_fan_cpu(g, params, iRow, find_tau, searchBounds);
else if (params->geometry == parameters::CONE)
{
if (params->helicalPitch != 0.0)
printf("Warning: find_centerCol/find_tau will likely not work for helical data\n");
return findCenter_cone_cpu(g, params, iRow, find_tau);
return findCenter_cone_cpu(g, params, iRow, find_tau, searchBounds);
}
else
{
Expand All @@ -354,12 +354,12 @@ float findCenter_cpu(float* g, parameters* params, int iRow, bool find_tau)
}
}

float findCenter_fan_cpu(float* g, parameters* params, int iRow, bool find_tau)
float findCenter_fan_cpu(float* g, parameters* params, int iRow, bool find_tau, float* searchBounds)
{
return findCenter_cone_cpu(g, params, iRow, find_tau);
return findCenter_cone_cpu(g, params, iRow, find_tau, searchBounds);
}

float findCenter_parallel_cpu(float* g, parameters* params, int iRow)
float findCenter_parallel_cpu(float* g, parameters* params, int iRow, float* searchBounds)
{
if (iRow < 0 || iRow > params->numRows - 1)
iRow = max(0, min(params->numRows-1, int(floor(0.5 + params->centerRow))));
Expand All @@ -386,7 +386,13 @@ float findCenter_parallel_cpu(float* g, parameters* params, int iRow)

int centerCol_low, centerCol_high;

setDefaultRange_centerCol(params, centerCol_low, centerCol_high);
if (searchBounds != NULL && 0.0 <= searchBounds[0] && searchBounds[1] <= params->numCols-1 && searchBounds[0] <= searchBounds[1])
{
centerCol_low = int(0.5+searchBounds[0]);
centerCol_high = int(0.5+searchBounds[1]);
}
else
setDefaultRange_centerCol(params, centerCol_low, centerCol_high);

double* shiftCosts = (double*)malloc(sizeof(double) * params->numCols);

Expand All @@ -395,8 +401,8 @@ float findCenter_parallel_cpu(float* g, parameters* params, int iRow)
float phi_min = min(phi_0, phi_end);
float phi_max = max(phi_0, phi_end);

float u_0 = params->u(0);
float u_end = params->u(params->numCols - 1);
//float u_0 = params->u(0);
//float u_end = params->u(params->numCols - 1);

omp_set_num_threads(omp_get_num_procs());
#pragma omp parallel for
Expand All @@ -406,6 +412,7 @@ float findCenter_parallel_cpu(float* g, parameters* params, int iRow)
//double denom = 0.0;

float u_0 = -(float(n) + params->colShiftFromFilter) * params->pixelWidth;
float u_end = params->pixelWidth * (params->numCols - 1) + u_0;

double num = 0.0;
double count = 0.0;
Expand Down Expand Up @@ -472,7 +479,7 @@ float findCenter_parallel_cpu(float* g, parameters* params, int iRow)
return retVal;
}

float findCenter_cone_cpu(float* g, parameters* params, int iRow, bool find_tau)
float findCenter_cone_cpu(float* g, parameters* params, int iRow, bool find_tau, float* searchBounds)
{
if (iRow < 0 || iRow > params->numRows - 1)
iRow = max(0, min(params->numRows - 1, int(floor(0.5 + params->centerRow))));
Expand All @@ -482,6 +489,8 @@ float findCenter_cone_cpu(float* g, parameters* params, int iRow, bool find_tau)
rowStart = iRow;
rowEnd = iRow;

//printf("iRow = %d, tiltAngle = %f\n", iRow, params->tiltAngle);

float* sino = get_rotated_sinogram(g, params, iRow);

int conj_ind = 0;
Expand All @@ -490,9 +499,21 @@ float findCenter_cone_cpu(float* g, parameters* params, int iRow, bool find_tau)
else
conj_ind = int(floor(0.5 + params->phi_inv(params->phis[0] - PI)));

float tau_shift = params->pixelWidth * params->sod / params->sdd;
int centerCol_low, centerCol_high;

setDefaultRange_centerCol(params, centerCol_low, centerCol_high);
if (searchBounds != NULL && find_tau == true)
{
searchBounds[0] = searchBounds[0]/tau_shift;
searchBounds[1] = searchBounds[1]/tau_shift;
}
if (searchBounds != NULL && 0.0 <= searchBounds[0] && searchBounds[1] <= params->numCols-1 && searchBounds[0] <= searchBounds[1])
{
centerCol_low = searchBounds[0];
centerCol_high = searchBounds[1];
}
else
setDefaultRange_centerCol(params, centerCol_low, centerCol_high);

double* shiftCosts = (double*)malloc(sizeof(double) * params->numCols);

Expand All @@ -501,18 +522,16 @@ float findCenter_cone_cpu(float* g, parameters* params, int iRow, bool find_tau)
bool normalizeConeAndFanCoordinateFunctions_save = params->normalizeConeAndFanCoordinateFunctions;
params->normalizeConeAndFanCoordinateFunctions = true;

float tau_shift = params->pixelWidth * params->sod / params->sdd;

float phi_0 = params->phis[0];
float phi_end = params->phis[params->numAngles - 1];
float phi_min = min(phi_0, phi_end);
float phi_max = max(phi_0, phi_end);

float u_0 = params->u(0);
float u_end = params->u(params->numCols-1);
//float u_0 = params->u(0);
//float u_end = params->u(params->numCols-1);

float atanTu = atan(params->pixelWidth / params->sdd);
float T_u = params->pixelWidth;
float atanTu = atan(params->pixelWidth / params->sdd); // radians
float T_u = params->pixelWidth / params->sdd;
if (params->detectorType == parameters::CURVED)
T_u = atanTu;

Expand All @@ -537,6 +556,7 @@ float findCenter_cone_cpu(float* g, parameters* params, int iRow, bool find_tau)
{
u_0 = -(float(n) + params->colShiftFromFilter) * T_u;
}
float u_end = T_u*(params->numCols-1) + u_0;

double num = 0.0;
double count = 0.0;
Expand All @@ -552,9 +572,9 @@ float findCenter_cone_cpu(float* g, parameters* params, int iRow, bool find_tau)
for (int k = 0; k < params->numCols; k++)
{
//float u = params->u(k);
float u = (k * params->pixelWidth + u_0) / params->sdd;
float u = k * T_u + u_0;
if (params->detectorType == parameters::CURVED)
u = atan((k * atanTu + u_0));
u = tan((k * atanTu + u_0));

float u_conj = (-u + A) / (1.0 + u*A);
float phi_conj = phi - 2.0 * atan(u) + atanA + PI;
Expand All @@ -564,9 +584,10 @@ float findCenter_cone_cpu(float* g, parameters* params, int iRow, bool find_tau)
{
int phi_conj_ind = int(0.5 + params->phi_inv(phi_conj));

int u_conj_ind = int(0.5 + ((u_conj * params->sdd) - u_0) / params->pixelWidth);
//int u_conj_ind = int(0.5 + ((u_conj * params->sdd) - u_0) / params->pixelWidth);
int u_conj_ind = int(0.5 + (u_conj - u_0) / T_u);
if (params->detectorType == parameters::CURVED)
u_conj_ind = int(0.5 + (tan(u_conj) - u_0) / atanTu);
u_conj_ind = int(0.5 + (atan(u_conj) - u_0) / atanTu);

float val = lineA[k];
//float val_conj = g[uint64(phi_conj_ind) * uint64(params->numRows * params->numCols) + uint64(j * params->numCols + u_conj_ind)];
Expand Down Expand Up @@ -723,7 +744,7 @@ float* get_rotated_sinogram(float* g, parameters* params, int iRow)
if (params->tiltAngle == 0.0)
{
for (int iCol = 0; iCol < params->numCols; iCol++)
sino_line[iCol] = aProj[uint64(iRow)*uint64(params->numRows) + uint64(iCol)];
sino_line[iCol] = aProj[uint64(iRow)*uint64(params->numCols) + uint64(iCol)];
}
else
{
Expand Down
8 changes: 4 additions & 4 deletions src/find_center_cpu.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,11 +22,11 @@
* which is also known as a half-fan or half-cone.
*/

float findCenter_cpu(float* g, parameters* params, int iRow = -1, bool find_tau = false);
float findCenter_cpu(float* g, parameters* params, int iRow = -1, bool find_tau = false, float* searchBounds = NULL);

float findCenter_parallel_cpu(float* g, parameters* params, int iRow = -1);
float findCenter_fan_cpu(float* g, parameters* params, int iRow = -1, bool find_tau = false);
float findCenter_cone_cpu(float* g, parameters* params, int iRow = -1, bool find_tau = false);
float findCenter_parallel_cpu(float* g, parameters* params, int iRow = -1, float* searchBounds = NULL);
float findCenter_fan_cpu(float* g, parameters* params, int iRow = -1, bool find_tau = false, float* searchBounds = NULL);
float findCenter_cone_cpu(float* g, parameters* params, int iRow = -1, bool find_tau = false, float* searchBounds = NULL);

float estimateTilt(float* g, parameters* params);
bool getConjugateDifference(float* g, parameters* params, float alpha, float centerCol, float* diff);
Expand Down
28 changes: 18 additions & 10 deletions src/leapctype.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -843,7 +843,7 @@ def set_centerCol(self, centerCol):
self.libprojectors.set_centerCol.argtypes = [ctypes.c_float]
return self.libprojectors.set_centerCol(centerCol)

def find_centerCol(self, g, iRow=-1):
def find_centerCol(self, g, iRow=-1, searchBounds=None):
r"""Find the centerCol parameter
This function works by minimizing the difference of conjugate rays, by changing the detector column sample locations. The cost functions
Expand Down Expand Up @@ -875,16 +875,20 @@ def find_centerCol(self, g, iRow=-1):
"""
if iRow is None:
iRow = -1
if searchBounds is None:
searchBounds = -1.0*np.ones(2, dtype=np.float32)
else:
searchBounds = np.array(searchBounds, dtype=np.float32)
self.libprojectors.find_centerCol.restype = ctypes.c_float
self.set_model()
if has_torch == True and type(g) is torch.Tensor:
self.libprojectors.find_centerCol.argtypes = [ctypes.c_void_p, ctypes.c_int, ctypes.c_bool]
return self.libprojectors.find_centerCol(g.data_ptr(), iRow, g.is_cuda == False)
self.libprojectors.find_centerCol.argtypes = [ctypes.c_void_p, ctypes.c_int, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_bool]
return self.libprojectors.find_centerCol(g.data_ptr(), iRow, searchBounds, g.is_cuda == False)
else:
self.libprojectors.find_centerCol.argtypes = [ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int, ctypes.c_bool]
return self.libprojectors.find_centerCol(g, iRow, True)
self.libprojectors.find_centerCol.argtypes = [ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_bool]
return self.libprojectors.find_centerCol(g, iRow, searchBounds, True)

def find_tau(self, g, iRow=-1):
def find_tau(self, g, iRow=-1, searchBounds=None):
r"""Find the tau parameter
This function works by minimizing the difference of conjugate rays, by changing the horizontal source position shift (equivalent to rotation stage shifts).
Expand Down Expand Up @@ -915,14 +919,18 @@ def find_tau(self, g, iRow=-1):
"""
if iRow is None:
iRow = -1
if searchBounds is None:
searchBounds = -1.0*np.ones(2, dtype=np.float32)
else:
searchBounds = np.array(searchBounds, dtype=np.float32)
self.libprojectors.find_tau.restype = ctypes.c_float
self.set_model()
if has_torch == True and type(g) is torch.Tensor:
self.libprojectors.find_tau.argtypes = [ctypes.c_void_p, ctypes.c_int, ctypes.c_bool]
return self.libprojectors.find_tau(g.data_ptr(), iRow, g.is_cuda == False)
self.libprojectors.find_tau.argtypes = [ctypes.c_void_p, ctypes.c_int, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_bool]
return self.libprojectors.find_tau(g.data_ptr(), iRow, searchBounds, g.is_cuda == False)
else:
self.libprojectors.find_tau.argtypes = [ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int, ctypes.c_bool]
return self.libprojectors.find_tau(g, iRow, True)
self.libprojectors.find_tau.argtypes = [ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_bool]
return self.libprojectors.find_tau(g, iRow, searchBounds, True)

def estimate_tilt(self, g):
"""Estimates the tilt angle (around the optical axis) of the detector
Expand Down
8 changes: 4 additions & 4 deletions src/tomographic_models.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3914,21 +3914,21 @@ bool tomographicModels::synthesize_symmetry(float* f_radial, float* f)
return symObject.synthesizeSymmetry(f_radial, f);
}

float tomographicModels::find_centerCol(float* g, int iRow, bool data_on_cpu)
float tomographicModels::find_centerCol(float* g, int iRow, float* searchBounds, bool data_on_cpu)
{
if (data_on_cpu)
return findCenter_cpu(g, &params, iRow);
return findCenter_cpu(g, &params, iRow, false, searchBounds);
else
{
printf("Error: find_centerCol not yet implemented for data on the GPU\n");
return 0.0;
}
}

float tomographicModels::find_tau(float* g, int iRow, bool data_on_cpu)
float tomographicModels::find_tau(float* g, int iRow, float* searchBounds, bool data_on_cpu)
{
if (data_on_cpu)
return findCenter_cpu(g, &params, iRow, true);
return findCenter_cpu(g, &params, iRow, true, searchBounds);
else
{
printf("Error: find_tau not yet implemented for data on the GPU\n");
Expand Down
6 changes: 4 additions & 2 deletions src/tomographic_models.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -686,20 +686,22 @@ class tomographicModels
* \brief finds centerCol of parallel-, fan-, or cone-beam data using conjugate rays
* \param[in] g pointer to the projection data
* \param[in] iRow the detector row index to use the estimate centerCol
* \param[in] searchBounds 2-element array specifying the search bounds
* \param[in] data_on_cpu true if data (g) is on the cpu, false if they are on the gpu
* \return true if operation was sucessful, false otherwise
*/
float find_centerCol(float* g, int iRow, bool data_on_cpu);
float find_centerCol(float* g, int iRow, float* searchBounds, bool data_on_cpu);

/**
* \fn find_tau
* \brief finds tau of fan- or cone-beam data using conjugate rays
* \param[in] g pointer to the projection data
* \param[in] iRow the detector row index to use the estimate centerCol
* \param[in] searchBounds 2-element array specifying the search bounds
* \param[in] data_on_cpu true if data (g) is on the cpu, false if they are on the gpu
* \return true if operation was sucessful, false otherwise
*/
float find_tau(float* g, int iRow, bool data_on_cpu);
float find_tau(float* g, int iRow, float* searchBounds, bool data_on_cpu);

/**
* \fn estimate_tilt
Expand Down
8 changes: 4 additions & 4 deletions src/tomographic_models_c_interface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -912,14 +912,14 @@ float get_z0()
return tomo()->params.z_0();
}

float find_centerCol(float* g, int iRow, bool data_on_cpu)
float find_centerCol(float* g, int iRow, float* searchBounds, bool data_on_cpu)
{
return tomo()->find_centerCol(g, iRow, data_on_cpu);
return tomo()->find_centerCol(g, iRow, searchBounds, data_on_cpu);
}

float find_tau(float* g, int iRow, bool data_on_cpu)
float find_tau(float* g, int iRow, float* searchBounds, bool data_on_cpu)
{
return tomo()->find_tau(g, iRow, data_on_cpu);
return tomo()->find_tau(g, iRow, searchBounds, data_on_cpu);
}

float consistency_cost(float* g, float Delta_centerRow, float Delta_centerCol, float Delta_tau, float Delta_tilt, bool data_on_cpu)
Expand Down
4 changes: 2 additions & 2 deletions src/tomographic_models_c_interface.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -215,8 +215,8 @@ extern "C" PROJECTOR_API float get_offsetY();
extern "C" PROJECTOR_API float get_offsetZ();
extern "C" PROJECTOR_API float get_z0();

extern "C" PROJECTOR_API float find_centerCol(float* g, int iRow, bool data_on_cpu);
extern "C" PROJECTOR_API float find_tau(float* g, int iRow, bool data_on_cpu);
extern "C" PROJECTOR_API float find_centerCol(float* g, int iRow, float* searchBounds, bool data_on_cpu);
extern "C" PROJECTOR_API float find_tau(float* g, int iRow, float* searchBounds, bool data_on_cpu);
extern "C" PROJECTOR_API float consistency_cost(float* g, float Delta_centerRow, float Delta_centerCol, float Delta_tau, float Delta_tilt, bool data_on_cpu);
extern "C" PROJECTOR_API float estimate_tilt(float* g, bool data_on_cpu);
extern "C" PROJECTOR_API float conjugate_difference(float* g, float alpha, float centerCol, float* diff, bool data_on_cpu);
Expand Down

0 comments on commit 2eaf458

Please sign in to comment.