Skip to content

Commit

Permalink
PERF: Fill jsh (JacobianOfSpatialJacobian) in-place and remove jsh1
Browse files Browse the repository at this point in the history
Used modern C++ range-based `for` loops to iterate over the data of `jsh`. Reduced dynamic memory usage by removing `jsh1`.
  • Loading branch information
N-Dekker committed May 10, 2023
1 parent 8477fdb commit b9ea3a8
Showing 1 changed file with 8 additions and 10 deletions.
18 changes: 8 additions & 10 deletions Common/Transforms/itkAdvancedCombinationTransform.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -1193,7 +1193,6 @@ AdvancedCombinationTransform<TScalarType, NDimensions>::GetJacobianOfSpatialHess
SpatialJacobianType sj0;
SpatialHessianType sh0;
JacobianOfSpatialJacobianType jsj1;
JacobianOfSpatialHessianType jsh1;

/** Transform the input point. */
// \todo: this has already been computed and it is expensive.
Expand All @@ -1207,19 +1206,19 @@ AdvancedCombinationTransform<TScalarType, NDimensions>::GetJacobianOfSpatialHess
/** Assume/demand that GetJacobianOfSpatialJacobian returns
* the same nonZeroJacobianIndices as the GetJacobianOfSpatialHessian. */
m_CurrentTransform->GetJacobianOfSpatialJacobian(transformedPoint, jsj1, nonZeroJacobianIndices);
m_CurrentTransform->GetJacobianOfSpatialHessian(transformedPoint, jsh1, nonZeroJacobianIndices);
m_CurrentTransform->GetJacobianOfSpatialHessian(transformedPoint, jsh, nonZeroJacobianIndices);

typename SpatialJacobianType::InternalMatrixType sj0tvnl = sj0.GetTranspose();
SpatialJacobianType sj0t(sj0tvnl);

jsh.resize(nonZeroJacobianIndices.size());

/** Combine them in one overall Jacobian of spatial Hessian. */
for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu)
for (auto & spatialHessian : jsh)
{
for (unsigned int dim = 0; dim < SpaceDimension; ++dim)
for (auto & matrix : spatialHessian)
{
jsh[mu][dim] = sj0t * (jsh1[mu][dim] * sj0);
matrix = sj0t * (matrix * sj0);
}
}

Expand Down Expand Up @@ -1256,7 +1255,6 @@ AdvancedCombinationTransform<TScalarType, NDimensions>::GetJacobianOfSpatialHess
SpatialJacobianType sj0, sj1;
SpatialHessianType sh0, sh1;
JacobianOfSpatialJacobianType jsj1;
JacobianOfSpatialHessianType jsh1;

/** Transform the input point. */
// \todo this has already been computed and it is expensive.
Expand All @@ -1272,18 +1270,18 @@ AdvancedCombinationTransform<TScalarType, NDimensions>::GetJacobianOfSpatialHess
* nonZeroJacobianIndices as the GetJacobianOfSpatialHessian.
*/
m_CurrentTransform->GetJacobianOfSpatialJacobian(transformedPoint, sj1, jsj1, nonZeroJacobianIndices);
m_CurrentTransform->GetJacobianOfSpatialHessian(transformedPoint, sh1, jsh1, nonZeroJacobianIndices);
m_CurrentTransform->GetJacobianOfSpatialHessian(transformedPoint, sh1, jsh, nonZeroJacobianIndices);

typename SpatialJacobianType::InternalMatrixType sj0tvnl = sj0.GetTranspose();
SpatialJacobianType sj0t(sj0tvnl);
jsh.resize(nonZeroJacobianIndices.size());

/** Combine them in one overall Jacobian of spatial Hessian. */
for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu)
for (auto & spatialHessian : jsh)
{
for (unsigned int dim = 0; dim < SpaceDimension; ++dim)
for (auto & matrix : spatialHessian)
{
jsh[mu][dim] = sj0t * (jsh1[mu][dim] * sj0);
matrix = sj0t * (matrix * sj0);
}
}

Expand Down

0 comments on commit b9ea3a8

Please sign in to comment.