Skip to content
New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

#

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

PERF: Speed up ComputeJacobianTerms access to vnl_sparse_matrix #959

Merged
merged 1 commit into from
Sep 1, 2023
Merged

PERF: Speed up ComputeJacobianTerms access to vnl_sparse_matrix #959

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 40 additions & 6 deletions Common/itkComputeJacobianTerms.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -187,7 +187,41 @@ ComputeJacobianTerms<TFixedImage, TTransform>::Compute(double & TrC, double & Tr

/** Initialize covariance matrix. Sparse, diagonal, and band form. */
vnl_sparse_matrix<CovarianceValueType> cov(numberOfParameters, numberOfParameters);
DiagCovarianceMatrixType diagcov(numberOfParameters, 0.0);

// getCovariance(r, c) is faster than cov(r, c), at least when using the vnl_sparse_matrix implementation included
// with ITK 5.3.0 (released on December 20, 2022). The speed improvement is especially large on Debug builds.
const auto getCovariance = [&cov](const unsigned int r, const unsigned int c) -> CovarianceValueType & {
auto & row = cov.get_row(r);

if (row.empty())
{
row.push_back({ c, 0.0 });
return row.back().second;
}

if (c < row.back().first)
{
// Because the last column number in the row is greater than `c`, the following iteration will stop before the end
// of the row.
auto it = row.begin();

while (it->first < c)
{
++it;
}
return (it->first == c ? it : row.insert(it, { c, 0.0 }))->second;
}

// At this point, the row is non-empty and c <= row.back().first.

if (c > row.back().first)
{
row.push_back({ c, 0.0 });
}
return row.back().second;
};

DiagCovarianceMatrixType diagcov(numberOfParameters, 0.0);

/** For temporary storage of J'J. */
CovarianceMatrixType jactjac(sizejacind, sizejacind);
Expand Down Expand Up @@ -257,7 +291,7 @@ ComputeJacobianTerms<TFixedImage, TTransform>::Compute(double & TrC, double & Tr
}
else
{
cov(p, q) += tempval;
getCovariance(p, q) += tempval;
}
}
}
Expand Down Expand Up @@ -295,7 +329,7 @@ ComputeJacobianTerms<TFixedImage, TTransform>::Compute(double & TrC, double & Tr
}
else
{
cov(p, q) += tempval;
getCovariance(p, q) += tempval;
}
}
}
Expand All @@ -313,7 +347,7 @@ ComputeJacobianTerms<TFixedImage, TTransform>::Compute(double & TrC, double & Tr
if (std::abs(tempval) > 1e-14)
{
const unsigned int q = p + bandcovMap2[b];
cov(p, q) = tempval;
getCovariance(p, q) = tempval;
}
}
}
Expand All @@ -332,7 +366,7 @@ ComputeJacobianTerms<TFixedImage, TTransform>::Compute(double & TrC, double & Tr
while (notfinished)
{
const int col = cov.getcolumn();
cov(cov.getrow(), col) /= scales[col];
getCovariance(cov.getrow(), col) /= scales[col];
notfinished = cov.next();
}
}
Expand All @@ -343,7 +377,7 @@ ComputeJacobianTerms<TFixedImage, TTransform>::Compute(double & TrC, double & Tr
if (!cov.empty_row(p))
{
// avoid creation of element if the row is empty
CovarianceValueType & covpp = cov(p, p);
CovarianceValueType & covpp = getCovariance(p, p);
TrC += covpp;
diagcov[p] = covpp;
}
Expand Down