Merge remote-tracking branch 'origin/main'

README.md

@@ -4,7 +4,7 @@
 
 [![Test](https://github.com/paquiteau/snake-fmri/actions/workflows/test.yml/badge.svg)](https://github.com/paquiteau/snake-fmri/actions/workflows/test.yml)
 [![deploy-docs](https://github.com/paquiteau/snake-fmri/actions/workflows/deploy-docs.yml/badge.svg)](https://paquiteau.github.io/snake-fmri)
-[![arxiv](https://img.shields.io/badge/preprint-2024.6455-darkred?logo=arXiv&logoColor=white)](httsp://arxiv.org)
+[![HAL](https://img.shields.io/badge/preprint-04533862-purple?logo=HAL&logoColor=white)](https://hal.science/hal-04533862)
 
 [![python](https://img.shields.io/badge/python-3.10%2B-blue?logo=python&logoColor=blue)](https://pypi.org/project/snake-fmri)
 ![black](https://img.shields.io/badge/code--style-black-black)
@@ -26,7 +26,7 @@ To install SNAKE-fMRI, you can use pip.
 ```bash
 pip install snake-fmri 
 # Required for the reconstruction 
-pip install git+github.com/pysap-fmri
+pip install git+https://github.com/paquiteau/pysap-fmri
 # Recommended for the nufft acceleration 
 pip install gpunufft # or cufinufft 
 ```
@@ -46,7 +46,7 @@ import snkf
 ```
 
 # Documentation
-The documentation is available at https://paquiteau.github.io/snake-fmri/, our [preprint (TBA)](XXXX) describe also the framework in details.
+The documentation is available at https://paquiteau.github.io/snake-fmri/, our [preprint](https://hal.science/hal-04533862) describe also the framework in details.
 
 Don't hesitate to also check the [examples gallery (TBA)](https://paquiteau.github.io/snake-fmri).
 
@@ -64,16 +64,27 @@ The configuration  are located in `snkf/conf` and articulates over 3 main compon
 - the validation via statistical analysis
 
 
-<!--
 
 # Citing SNAKE-fMRI
 If you use SNAKE-fMRI in your research, please cite the following paper:
 
-> 
-> 
-
-```
+>  Pierre-Antoine Comby, Alexandre Vignaud, Philippe Ciuciu. SNAKE-fMRI: A modular fMRI data simulator from the space-time domain to k-space and back. 2024. ⟨hal-04533862⟩ 
+
+```bibtex
+@unpublished{comby:hal-04533862,
+  TITLE = {{SNAKE-fMRI: A modular fMRI data simulator from the space-time domain to k-space and back}},
+  AUTHOR = {Comby, Pierre-Antoine and Vignaud, Alexandre and Ciuciu, Philippe},
+  URL = {https://hal.science/hal-04533862},
+  NOTE = {working paper or preprint},
+  YEAR = {2024},
+  MONTH = Mar,
+  KEYWORDS = {fMRI ; Brain Imaging ; Accelerated sampling ; Compressed Sensing ; Simulation ; Open Source ; Python ; np ndarray , np newaxis] ; np ; ndarray , np ; newaxis]},
+  PDF = {https://hal.science/hal-04533862/file/main.pdf},
+  HAL_ID = {hal-04533862},
+  HAL_VERSION = {v1},
+}
 ```
--->
+
 # License
+
 SNAKE-fMRI is licensed under the MIT License. See the LICENSE file for more information.

src/snkf/reconstructors/__init__.py

@@ -6,6 +6,7 @@
     SequentialReconstructor,
     ZeroFilledReconstructor,
     LowRankPlusSparseReconstructor,
+    LowRankPlusTVReconstructor,
 )
 
 
@@ -16,4 +17,5 @@
     "SequentialReconstructor",
     "ZeroFilledReconstructor",
     "LowRankPlusSparseReconstructor",
+    "LowRankPlusTVReconstructor",
 ]

src/snkf/reconstructors/pysap.py

@@ -78,7 +78,6 @@ def get_fourier_operator(
     if "stacked" in backend_name:
         kwargs["z_index"] = "auto"
     if "nufft" in backend_name or "stacked" in backend_name:
-
         return LazySpaceFourier(
             backend=backend_name,
             samples=sim.kspace_mask,
@@ -308,9 +307,81 @@ def reconstruct(
 
 
 class LowRankPlusTVReconstructor(LowRankPlusSparseReconstructor):
-    """Low Rank + TV."""
+    """Low Rank + TV Benchmark reconstructors.
 
-    ...
+    Parameters
+    ----------
+    algorithm_tv
+        Algorithm for the TV proximal operator.
+    max_iter_tv
+        maximal number of iteration for the TV proximal operator.
+    """
+
+    __reconstructor_name__ = "lr_tv"
+
+    algorithm_tv: str = "condat"
+    max_iter_tv: int = 100
+
+    def __post_init__(self):
+        super().__post_init__()
+        if self.algorithm == "otazo_raw":
+            self.algorithm = "otazo"
+
+    def __str__(self):
+        if isinstance(self.lambda_s, float):
+            return f"LRTV-{self.lambda_l:.2e}-{self.lambda_s:.2e}"
+        return f"LRTV-{self.lambda_l:.2e}-{self.lambda_s}"
+
+    def setup(self, sim: SimData) -> None:
+        """Set up the reconstructor."""
+        from fmri.reconstructors.time_aware import LowRankPlusSparseReconstructor
+        from fmri.operators.utils import sure_est, sigma_mad
+        from fmri.operators.proximity import ProxTV1d
+        from fmri.operators.time_op import TimeFourier
+        from fmri.operators.svt import FlattenSVT
+
+        if self.fourier_op is None:
+            self.fourier_op = get_fourier_operator(
+                sim,
+                **self.nufft_kwargs,
+            )
+
+        logger.debug("Space Fourier operator initialized")
+        if self.time_linear_op is None:
+            self.time_linear_op = TimeFourier(time_axis=0)
+
+        logger.debug("Time Fourier operator initialized")
+        if self.lambda_s == "sure":
+            adj_data = self.fourier_op.adj_op(sim.kspace_data)
+            sure_thresh = np.zeros(np.prod(adj_data.shape[1:]))
+            tf = self.time_linear_op.op(adj_data).reshape(len(adj_data), -1)
+            for i in range(len(sure_thresh)):
+                sure_thresh[i] = sure_est(tf[:, i]) * sigma_mad(tf[:, i])
+
+            self.lambda_s = np.median(sure_thresh) / 2
+            logger.info(f"SURE threshold: {self.lambda_s:.4e}")
+
+        if self.time_prox_op is None and self.time_linear_op is not None:
+            self.time_prox_op = ProxTV1d(
+                self.lambda_s, method=self.algorithm_tv, max_iter=self.max_iter_tv
+            )
+
+        logger.debug("Prox Time  operator initialized")
+        if self.space_prox_op is None:
+            self.space_prox_op = FlattenSVT(
+                self.lambda_l, initial_rank=10, thresh_type="soft-rel"
+            )
+        logger.debug("Prox Space operator initialized")
+
+        self.reconstructor: LowRankPlusSparseReconstructor = (
+            LowRankPlusSparseReconstructor(
+                self.fourier_op,
+                space_prox_op=self.space_prox_op,
+                time_prox_op=self.time_prox_op,
+                cost="auto",
+            )
+        )
+        logger.debug("Reconstructor initialized")
 
 
 class LowRankPlusWaveletReconstructor(LowRankPlusSparseReconstructor):