Skip to content
New issue

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

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

Already on GitHub? Sign in to your account

TV-regularised reconstruction and TV-based denoising #596

Draft
wants to merge 118 commits into
base: main
Choose a base branch
from
Draft

Conversation

ckolbPTB
Copy link
Collaborator

@ckolbPTB ckolbPTB commented Jan 6, 2025

Requires #426

Copy link
Contributor

github-actions bot commented Jan 6, 2025

Coverage

Coverage Report
FileStmtsMissCoverMissing
src/mrpro/algorithms/csm
   inati.py24196%44
   walsh.py16194%34
src/mrpro/algorithms/dcf
   dcf_voronoi.py53492%15, 48–49, 76
src/mrpro/algorithms/optimizers
   adam.py20195%69
   pdhg.py79396%173–174, 180
src/mrpro/algorithms/reconstruction
   DirectReconstruction.py281643%51–71, 85
   IterativeSENSEReconstruction.py13192%76
   Reconstruction.py502256%42, 54–56, 80–87, 104–113
   RegularizedIterativeSENSEReconstruction.py411759%96–100, 114–139
   TotalVariationDenoising.py331458%54–57, 62, 77–100
   TotalVariationRegularizedReconstruction.py381658%81–83, 97–127
src/mrpro/data
   AcqInfo.py128398%26, 169, 207
   CsmData.py29390%15, 82–84
   DcfData.py45882%18, 66, 78–83
   IData.py67987%119, 125, 129, 159–167
   IHeader.py75791%75, 109, 127–131
   KHeader.py1531789%25, 119–123, 150, 199, 210, 217–218, 221, 228, 260–271
   KNoise.py311552%39–52, 56–61
   KTrajectory.py811285%108–113, 116–118, 203–207
   MoveDataMixin.py1401887%15, 113, 129, 143–145, 207, 323–325, 338, 417, 437–438, 440, 455–456, 458
   QData.py39782%42, 65–73
   Rotation.py6743595%100, 198, 335, 433, 477, 495, 581, 583, 592, 626, 628, 691, 768, 773, 776, 791, 808, 813, 889, 1077, 1082, 1085, 1109, 1113, 1240, 1242, 1250–1251, 1315, 1397, 1690, 1846, 1881, 1885, 1996
   SpatialDimension.py2322191%34, 104, 141, 148, 154, 274–276, 289–291, 325, 343, 356, 369, 382, 395, 404–405, 420, 429
   acq_filters.py12192%47
src/mrpro/data/_kdata
   KData.py1341887%109–110, 125, 132, 142, 150, 204–205, 243, 248–249, 268–279
   KDataRemoveOsMixin.py29293%44, 46
   KDataSelectMixin.py19289%48, 63
   KDataSplitMixin.py48394%53, 84, 93
src/mrpro/data/traj_calculators
   KTrajectoryCalculator.py25292%23, 45
   KTrajectoryIsmrmrd.py13285%41, 50
   KTrajectoryPulseq.py23196%55
src/mrpro/operators
   CartesianSamplingOp.py89397%118, 157, 280
   ConstraintsOp.py60297%46, 48
   EndomorphOperator.py65297%228, 234
   FiniteDifferenceOp.py27293%40, 105
   FourierOp.py158398%263, 381, 386
   Functional.py71593%20–22, 117, 119
   GridSamplingOp.py136993%72–73, 82–83, 90–91, 94, 96, 98
   LinearOperator.py168995%55, 91, 190, 220, 261, 270, 278, 295, 320
   LinearOperatorMatrix.py1581690%82, 119, 152, 161, 166, 175–178, 191–194, 203, 215, 304, 331, 359
   MultiIdentityOp.py13285%43, 48
   Operator.py78297%25, 74
   ProximableFunctionalSeparableSum.py39392%50, 103, 110
   SliceProjectionOp.py173895%44, 61, 63, 69, 206, 227, 260, 300
   WaveletOp.py120596%152, 170, 205, 210, 233
   ZeroPadOp.py16194%30
src/mrpro/utils
   filters.py62297%44, 49
   reshape.py60198%191
   slice_profiles.py46687%20, 36, 113–116, 149
   sliding_window.py34197%34
   split_idx.py10280%43, 47
   summarize_tensorvalues.py11918%20–29
   typing.py181139%8–23
   zero_pad_or_crop.py31681%26, 30, 54, 57, 60, 63
TOTAL508639292% 

Tests Skipped Failures Errors Time
2271 0 💤 0 ❌ 0 🔥 1m 9s ⏱️

Copy link
Contributor

github-actions bot commented Jan 7, 2025

📚 Documentation

📁 Download as zip
🔍 View online

class TotalVariationRegularizedReconstruction(DirectReconstruction):
r"""TV-regularized reconstruction.

This algorithm solves the problem :math:`min_x \frac{1}{2}||(Ax - y)||_2^2 + ||L\nabla x||_1`
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe it would be good to distinguish between \lambda * | \nabla x |_1 (for scalar reg parameter) and | \Lambda \nabla x |_1 for an entire lambda-map. Also, strictly speaking, we would need to introduce the notation to have different regularization parameters for different directions, but this could get too complicated and cumbersome for this here?

:math:`\nabla` is the finite difference operator applied to :math:`x`.
"""

n_iterations: int
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do we want to have n_iterations or max_iterations with a tolerance? Our PDHG allows for the second as well.

noise
KNoise used for prewhitening. If None, no prewhitening is performed
dcf
K-space sampling density compensation. If None, set up based on kdata. The dcf is only used to calculate a
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we want to restrict ourselves to this or also allow the problem min_x 1/2*||W^(1/2)(Ax-y)||_2^2?

kdata = prewhiten_kspace(kdata, self.noise)

# Create the acquisition model A = F S if the CSM S is defined otherwise A = F with the Fourier operator F
acquisition_operator = self.fourier_op @ self.csm.as_operator() if self.csm is not None else self.fourier_op
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

potentially W^(1/2) F C?

r"""TV denoising.

This algorithm solves the problem :math:`min_x \frac{1}{2}||(x - y)||_2^2 + ||L\nabla x||_1`
by using the PDHG-algorithm. :math:`y` is the target image, :math:`L` is the strength of the regularization and
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe better:
:math:y is the given noisy image


# %%
data_weight = 0.5
n_adam_iterations = 4
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

n_admm_iterations instead of n_adam_iterations

#
# by doing
#
# $x_{k+1} = argmin_x \lambda \| \nabla x \|_1 + \frac{\rho}{2}||x - z_k + u_k||_2^2$
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

\argmin

#
# using PDHG.
#
# Because we have 2D dynamic images we can apply the TV-regularization along x,y and time.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

suggestion:

Because we have 2D dynamic images, we can apply the TV-regularization along x,y, and time.


# %% [markdown]
# #### TV-Regularized Reconstruction using ADMM
# In the above example we need to apply the acquisition operator during the PDHG iterations which is computationally
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

suggestion:

In the above example, PDHG repeatedly applies the acquisition operator and its adjoint during the iterations, which is computationally demanding and hence takes a long time. Another option is to use the Alternating Direction Method of Multipliers (ADMM) [S. Boyd et al, 2011], which solves the general problem

# $u_{k+1} = u_k + x_{k+1} - z_{k+1}$
#
# The first step is TV-based denoising of $x$, the second step is a regularized iterative SENSE update of $z$ and the
# final step updates the helper variable $u$.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

updates the dual variable $u$.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants