use gram in iterative reco. remove dcf from examples

examples/iterative_sense_reconstruction.py

@@ -24,18 +24,16 @@
 # Let's assume we have obtained the k-space data $y$ from an image $x$ with an acquisition model (Fourier transforms,
 # coil sensitivity maps...) $A$ then we can formulate the forward problem as:
 #
-# $ y = Ax + n $
+# $ y = Ax + \eta $
 #
-# where $n$ describes complex Gaussian noise. The image $x$ can be obtained by minimizing the functional $F$
+# where $\eta$ describes complex Gaussian noise. The image $x$ can be obtained by minimizing the functional $F$
 #
-# $ F(x) = ||W^{\frac{1}{2}}(Ax - y)||_2^2 $
+# $ F(x) = ||(Ax - y)||_2^2 $
 #
-# where $W^\frac{1}{2}$ is the square root of the density compensation function (which corresponds to a diagonal
-# operator).
 #
 # Setting the derivative of the functional $F$ to zero and rearranging yields
 #
-# $ A^H W A x = A^H W y$
+# $ A^H A x = A^H y$
 #
 # which is a linear system $Hx = b$ that needs to be solved for $x$.
 # %%
@@ -91,7 +89,7 @@
 # Calculate coil maps
 # Note that operators return a tuple of tensors, so we need to unpack it,
 # even though there is only one tensor returned from adjoint operator.
-img_coilwise = mrpro.data.IData.from_tensor_and_kheader(*fourier_operator.H(*dcf_operator(kdata.data)), kdata.header)
+img_coilwise = mrpro.data.IData.from_tensor_and_kheader(*(fourier_operator.H @ dcf_operator)(kdata.data), kdata.header)
 csm_operator = mrpro.data.CsmData.from_idata_walsh(img_coilwise).as_operator()
 
 # Create the acquisition operator A
@@ -101,14 +99,14 @@
 # ##### Calculate the right-hand-side of the linear system $b = A^H W y$
 
 # %%
-(right_hand_side,) = acquisition_operator.H(dcf_operator(kdata.data)[0])
+(right_hand_side,) = acquisition_operator.H(kdata.data)
 
 
 # %% [markdown]
-# ##### Set-up the linear self-adjoint operator $H = A^H W A$
+# ##### Set-up the linear self-adjoint gram operator $H = A^H A$
 
 # %%
-operator = acquisition_operator.H @ dcf_operator @ acquisition_operator
+operator = acquisition_operator.gram
 
 # %% [markdown]
 # ##### Run conjugate gradient

examples/regularized_iterative_sense_reconstruction.py

@@ -28,20 +28,20 @@
 #
 # where $n$ describes complex Gaussian noise. The image $x$ can be obtained by minimizing the functionl $F$
 #
-# $ F(x) = ||W^{\frac{1}{2}}(Ax - y)||_2^2 $
+# $ F(x) = ||(Ax - y)||_2^2 $
 #
-# where $W^\frac{1}{2}$ is the square root of the density compensation function (which corresponds to a diagonal
-# operator). Because this is an ill-posed problem, we can add a regularization term to stabilize the problem and obtain
+# Because this is an ill-posed problem, we can add a regularization term to stabilize the problem and obtain
 # a solution with certain properties:
 #
-# $ F(x) = ||W^{\frac{1}{2}}(Ax - y)||_2^2 + l||Bx - x_{reg}||_2^2$
+# $ F(x) = ||(Ax - y)||_2^2 + \lambda||Bx - x_{reg}||_2^2$
 #
-# where $l$ is the strength of the regularization, $B$ is a linear operator and $x_{reg}$ is a regularization image.
+# where $\lambda$ is the strength of the regularization, $B$ is a linear operator and $x_{reg}$ is a
+# regularization image.
 # With this functional $F$ we obtain a solution which is close to $x_{reg}$ and to the acquired data $y$.
 #
 # Setting the derivative of the functional $F$ to zero and rearranging yields
 #
-# $ (A^H W A + l B) x = A^H W y + l x_{reg}$
+# $ (A^H A + \lambda B) x = A^H y + \lambda x_{reg}$
 #
 # which is a linear system $Hx = b$ that needs to be solved for $x$.
 #
@@ -143,9 +143,7 @@
 # ##### Calculate the right-hand-side of the linear system $b = A^H W y + l x_{reg}$
 
 # %%
-right_hand_side = (
-    acquisition_operator.H(dcf_operator(kdata_us.data)[0])[0] + regularization_weight * img_iterative_sense.data
-)
+right_hand_side = acquisition_operator.H(kdata_us.data) + regularization_weight * img_iterative_sense.data
 
 
 # %% [markdown]
@@ -154,9 +152,7 @@
 # %%
 from mrpro.operators import IdentityOp
 
-operator = acquisition_operator.H @ dcf_operator @ acquisition_operator + IdentityOp() * torch.as_tensor(
-    regularization_weight
-)
+operator = acquisition_operator.gram + IdentityOp() * torch.as_tensor(regularization_weight)
 
 # %% [markdown]
 # ##### Run conjugate gradient

src/mrpro/algorithms/reconstruction/IterativeSENSEReconstruction.py

@@ -17,14 +17,13 @@
 class IterativeSENSEReconstruction(RegularizedIterativeSENSEReconstruction):
     r"""Iterative SENSE reconstruction.
 
-    This algorithm solves the problem :math:`min_x \frac{1}{2}||W^\frac{1}{2} (Ax - y)||_2^2`
+    This algorithm solves the problem :math:`min_x \frac{1}{2}||(Ax - y)||_2^2`
     by using a conjugate gradient algorithm to solve
-    :math:`H x = b` with :math:`H = A^H W A` and :math:`b = A^H W y` where :math:`A` is the acquisition model
-    (coil sensitivity maps, Fourier operator, k-space sampling), :math:`y` is the acquired k-space data and :math:`W`
-    describes the density compensation [PRU2001]_ .
+    :math:`H x = b` with :math:`H = A^H A` and :math:`b = A^H y` where :math:`A` is the acquisition model
+    (coil sensitivity maps, Fourier operator, k-space sampling) and :math:`y` is the acquired k-space data [PRU2001]_ .
 
-    Note: In [PRU2001]_ a k-space filter is applied as a final step to null all k-space values outside the k-space
-    coverage. This is not done here.
+    Note: In [PRU2001]_, the density function is used to to reweight the reconstruction objective and a k-space filter
+    is applied as a final step to null all k-space values outside the k-space coverage. Both is not done here.
 
     .. [PRU2001] Pruessmann K, Weiger M, Boernert P, and Boesiger P (2001), Advances in sensitivity encoding with
        arbitrary k-space trajectories. MRI 46, 638-651. https://doi.org/10.1002/mrm.1241
@@ -64,7 +63,8 @@ def __init__(
         noise
             KNoise used for prewhitening. If None, no prewhitening is performed
         dcf
-            K-space sampling density compensation. If None, set up based on kdata.
+            K-space sampling density compensation.
+            Only used to obtain the starting point of the iterative reconstruction.
         n_iterations
             Number of CG iterations
 

src/mrpro/algorithms/reconstruction/RegularizedIterativeSENSEReconstruction.py

@@ -21,23 +21,24 @@
 class RegularizedIterativeSENSEReconstruction(DirectReconstruction):
     r"""Regularized iterative SENSE reconstruction.
 
-    This algorithm solves the problem :math:`min_x \frac{1}{2}||W^\frac{1}{2} (Ax - y)||_2^2 +
-    \frac{1}{2}L||Bx - x_0||_2^2`
+    This algorithm solves the problem :math:`min_x \frac{1}{2}||(Ax - y)||_2^2 +
+    \frac{1}{2} \lambda ||Bx - x_r||_2^2`
     by using a conjugate gradient algorithm to solve
-    :math:`H x = b` with :math:`H = A^H W A + L B` and :math:`b = A^H W y + L x_0` where :math:`A`
+    :math:`H x = b` with :math:`H = A^H A + \lambda B` and :math:`b = A^H y + \lambda x_r` where :math:`A`
     is the acquisition model (coil sensitivity maps, Fourier operator, k-space sampling), :math:`y` is the acquired
-    k-space data, :math:`W` describes the density compensation, :math:`L` is the strength of the regularization and
-    :math:`x_0` is the regularization image (i.e. the prior). :math:`B` is a linear operator applied to :math:`x`.
+    k-space data, :math:`\lambda` is the regularization weight.
+    :math:`x_0` is the regularization data (i.e. the prior) and :math:`B` is a linear operator applied to :math:`x`
+    in the regularization term.
     """
 
     n_iterations: int
     """Number of CG iterations."""
 
     regularization_data: torch.Tensor
-    """Regularization data (i.e. prior) :math:`x_0`."""
+    """Regularization data (i.e. prior) :math:`x_r`."""
 
     regularization_weight: torch.Tensor
-    """Strength of the regularization :math:`L`."""
+    r"""Strength of the regularization :math:`\lambda`."""
 
     regularization_op: LinearOperator
     """Linear operator :math:`B` applied to the current estimate in the regularization term."""
@@ -76,7 +77,8 @@ def __init__(
         noise
             KNoise used for prewhitening. If None, no prewhitening is performed
         dcf
-            K-space sampling density compensation. If None, set up based on kdata.
+            K-space sampling density compensation.
+            Only used to obtain the starting point of the iterative reconstruction.
         n_iterations
             Number of CG iterations
         regularization_data
@@ -114,24 +116,29 @@ def forward(self, kdata: KData) -> IData:
         if self.noise is not None:
             kdata = prewhiten_kspace(kdata, self.noise)
 
-        # Create the normal operator as H = A^H W A if the DCF is not None otherwise as H = A^H A.
-        # The acquisition model is A = F S if the CSM S is defined otherwise A = F with the Fourier operator F
+        # Create the the acquisition model A = Fourier CSM # and normal operator H = A^H A
         csm_op = self.csm.as_operator() if self.csm is not None else IdentityOp()
-        precondition_op = self.dcf.as_operator() if self.dcf is not None else IdentityOp()
-        operator = (self.fourier_op @ csm_op).H @ precondition_op @ (self.fourier_op @ csm_op)
+        acquisition_model = self.fourier_op @ csm_op
+        operator = acquisition_model.gram
 
-        # Calculate the right-hand-side as b = A^H W y if the DCF is not None otherwise as b = A^H y.
-        (right_hand_side,) = (self.fourier_op @ csm_op).H(precondition_op(kdata.data)[0])
+        # Calculate the right-hand-side as b = A^H y
+        (right_hand_side,) = (acquisition_model.H)(kdata.data)
 
         # Add regularization
         if not torch.all(self.regularization_weight == 0):
-            operator = operator + IdentityOp() @ (self.regularization_weight * self.regularization_op)
+            operator = operator + (self.regularization_weight * self.regularization_op)
             right_hand_side += self.regularization_weight * self.regularization_data
 
+        # if available, use dcf for initial value
+        if self.dcf is not None:
+            (initial_value,) = (acquisition_model.H @ self.dcf.as_operator())(kdata.data)
+        else:
+            initial_value = right_hand_side
+
         img_tensor = cg(
             operator,
             right_hand_side,
-            initial_value=right_hand_side,
+            initial_value=initial_value,
             max_iterations=self.n_iterations,
             tolerance=0.0,
         )