Skip to content

Commit

Permalink
increase numerical precision of priors
Browse files Browse the repository at this point in the history 
use double in neighbourhood calculations
  • Loading branch information
@KrisThielemans
KrisThielemans committed Apr 28, 2024
1 parent 86c4f7f commit d30729f
Show file tree
Hide file tree
Showing 5 changed files with 25 additions and 69 deletions.
5 changes: 3 additions & 2 deletions documentation/release_6.1.htm
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,8 +59,9 @@ <h3>New functionality</h3>

<h3>Changed functionality</h3>
<ul>
<li>
</li>
<li>
Accumulation in computation of priors now uses doubles, which could result in slightly better precision. <a href=https://github.com/UCL/STIR/pull/1410>Part of PR #1410</a>.
</li>
</ul>


Expand Down
14 changes: 7 additions & 7 deletions src/recon_buildblock/LogcoshPrior.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -289,12 +289,12 @@ LogcoshPrior<elemT>::compute_value(const DiscretisedDensity<3, elemT>& current_i
for (int dx = min_dx; dx <= max_dx; ++dx)
{
// 1/scalar^2 * log(cosh(x * scalar))
elemT voxel_diff = current_image_estimate[z][y][x] - current_image_estimate[z + dz][y + dy][x + dx];
elemT current
double voxel_diff = current_image_estimate[z][y][x] - current_image_estimate[z + dz][y + dy][x + dx];
double current
= weights[dz][dy][dx] * 1 / (this->scalar * this->scalar) * logcosh(this->scalar * voxel_diff);
if (do_kappa)
current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z + dz][y + dy][x + dx];
result += static_cast<double>(current);
result += current;
}
}
}
Expand Down Expand Up @@ -349,21 +349,21 @@ LogcoshPrior<elemT>::compute_gradient(DiscretisedDensity<3, elemT>& prior_gradie
const int min_dx = max(weights[0][0].get_min_index(), min_x - x);
const int max_dx = min(weights[0][0].get_max_index(), max_x - x);

elemT gradient = 0;
double gradient = 0;
for (int dz = min_dz; dz <= max_dz; ++dz)
for (int dy = min_dy; dy <= max_dy; ++dy)
for (int dx = min_dx; dx <= max_dx; ++dx)
{
// 1/scalar * tanh(x * scalar)
elemT voxel_diff = current_image_estimate[z][y][x] - current_image_estimate[z + dz][y + dy][x + dx];
elemT current = weights[dz][dy][dx] * (1 / this->scalar) * tanh(this->scalar * voxel_diff);
double voxel_diff = current_image_estimate[z][y][x] - current_image_estimate[z + dz][y + dy][x + dx];
double current = weights[dz][dy][dx] * (1 / this->scalar) * tanh(this->scalar * voxel_diff);

if (do_kappa)
current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z + dz][y + dy][x + dx];

gradient += current;
}
prior_gradient[z][y][x] = gradient * this->penalisation_factor;
prior_gradient[z][y][x] = static_cast<elemT>(gradient * this->penalisation_factor);
}
}
}
Expand Down
4 changes: 2 additions & 2 deletions src/recon_buildblock/PLSPrior.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -498,12 +498,12 @@ PLSPrior<elemT>::compute_value(const DiscretisedDensity<3, elemT>& current_image
(penalty[z][y][x]) * (*kappa_ptr)[z][y][x];
*/

elemT current = (*penalty_sptr)[z][y][x];
double current = (*penalty_sptr)[z][y][x];

if (do_kappa)
current *= (*kappa_ptr)[z][y][x];

result += static_cast<double>(current);
result += current;
}
}
}
Expand Down
61 changes: 8 additions & 53 deletions src/recon_buildblock/QuadraticPrior.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -278,14 +278,14 @@ QuadraticPrior<elemT>::compute_value(const DiscretisedDensity<3, elemT>& current
for (int dy = min_dy; dy <= max_dy; ++dy)
for (int dx = min_dx; dx <= max_dx; ++dx)
{
elemT current = weights[dz][dy][dx]
* square(current_image_estimate[z][y][x] - current_image_estimate[z + dz][y + dy][x + dx])
/ 4;
double current = weights[dz][dy][dx]
* square(current_image_estimate[z][y][x] - current_image_estimate[z + dz][y + dy][x + dx])
/ 4;

if (do_kappa)
current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z + dz][y + dy][x + dx];

result += static_cast<double>(current);
result += current;
}
}
}
Expand Down Expand Up @@ -348,65 +348,20 @@ QuadraticPrior<elemT>::compute_gradient(DiscretisedDensity<3, elemT>& prior_grad
(current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]) *
(*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx];
*/
#if 1
elemT gradient = 0;
double gradient = 0.;
for (int dz = min_dz; dz <= max_dz; ++dz)
for (int dy = min_dy; dy <= max_dy; ++dy)
for (int dx = min_dx; dx <= max_dx; ++dx)
{
elemT current = weights[dz][dy][dx]
* (current_image_estimate[z][y][x] - current_image_estimate[z + dz][y + dy][x + dx]);
double current = weights[dz][dy][dx]
* (current_image_estimate[z][y][x] - current_image_estimate[z + dz][y + dy][x + dx]);

if (do_kappa)
current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z + dz][y + dy][x + dx];

gradient += current;
}
#else
// attempt to speed up by precomputing the sum of weights.
// The current code gives identical results but is actually slower
// than the above, at least when kappas are present.
// precompute sum of weights
// TODO without kappas, this is just weights.sum() most of the time,
// but not near edges
float sum_of_weights = 0;
{
if (do_kappa)
{
for (int dz = min_dz; dz <= max_dz; ++dz)
for (int dy = min_dy; dy <= max_dy; ++dy)
for (int dx = min_dx; dx <= max_dx; ++dx)
sum_of_weights += weights[dz][dy][dx] * (*kappa_ptr)[z + dz][y + dy][x + dx];
}
else
{
for (int dz = min_dz; dz <= max_dz; ++dz)
for (int dy = min_dy; dy <= max_dy; ++dy)
for (int dx = min_dx; dx <= max_dx; ++dx)
sum_of_weights += weights[dz][dy][dx];
}
}
// now compute contribution of central term
elemT gradient = sum_of_weights * current_image_estimate[z][y][x];
// subtract the rest
for (int dz = min_dz; dz <= max_dz; ++dz)
for (int dy = min_dy; dy <= max_dy; ++dy)
for (int dx = min_dx; dx <= max_dx; ++dx)
{
elemT current = weights[dz][dy][dx] * current_image_estimate[z + dz][y + dy][x + dx];
if (do_kappa)
current *= (*kappa_ptr)[z + dz][y + dy][x + dx];
gradient -= current;
}
// multiply with central kappa
if (do_kappa)
gradient *= (*kappa_ptr)[z][y][x];
#endif
prior_gradient[z][y][x] = gradient * this->penalisation_factor;
prior_gradient[z][y][x] = static_cast<elemT>(gradient * this->penalisation_factor);
}
}
}
Expand Down
10 changes: 5 additions & 5 deletions src/recon_buildblock/RelativeDifferencePrior.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -340,7 +340,7 @@ RelativeDifferencePrior<elemT>::compute_value(const DiscretisedDensity<3, elemT>
for (int dy = min_dy; dy <= max_dy; ++dy)
for (int dx = min_dx; dx <= max_dx; ++dx)
{
elemT current;
double current;
if (this->epsilon == 0.0 && current_image_estimate[z][y][x] == 0.0
&& current_image_estimate[z + dz][y + dy][x + dx] == 0.0)
{
Expand All @@ -355,7 +355,7 @@ RelativeDifferencePrior<elemT>::compute_value(const DiscretisedDensity<3, elemT>
if (do_kappa)
current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z + dz][y + dy][x + dx];

result += static_cast<double>(current);
result += current;
}
}
}
Expand Down Expand Up @@ -412,12 +412,12 @@ RelativeDifferencePrior<elemT>::compute_gradient(DiscretisedDensity<3, elemT>& p
const int min_dx = max(weights[0][0].get_min_index(), min_x - x);
const int max_dx = min(weights[0][0].get_max_index(), max_x - x);

elemT gradient = 0;
double gradient = 0;
for (int dz = min_dz; dz <= max_dz; ++dz)
for (int dy = min_dy; dy <= max_dy; ++dy)
for (int dx = min_dx; dx <= max_dx; ++dx)
{
elemT current
double current
= weights[dz][dy][dx]
* derivative_10(current_image_estimate[z][y][x], current_image_estimate[z + dz][y + dy][x + dx]);
if (do_kappa)
Expand All @@ -426,7 +426,7 @@ RelativeDifferencePrior<elemT>::compute_gradient(DiscretisedDensity<3, elemT>& p
gradient += current;
}

prior_gradient[z][y][x] = gradient * this->penalisation_factor;
prior_gradient[z][y][x] = static_cast<elemT>(gradient * this->penalisation_factor);
}
}
}
Expand Down

0 comments on commit d30729f

Please sign in to comment.