Connection size affine

CHANGELOG.md

@@ -1,3 +1,5 @@
+- 1.0.6
+  - Replace `VoxelGridPhantom.size` attribute with a property that sets `affine` correctly
 - 1.0.5
   - fix .seq soft-delays: default to 0.0 for missing values, add documentation
 - 1.0.4

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "mrzero_core"
-version = "1.0.5"
+version = "1.0.6"
 edition = "2021"
 
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

documentation/src/api-phantom.md

@@ -42,7 +42,7 @@ Contains the data used for simulation (together with the [`Sequence`](api-sequen
 
 ```python
 # Usually not constructed directly but via the `.build()` function on a phantom
-SimData(PD, T1, T2, T2dash, D, B0, B1, coil_sens, size, voxel_pos, nyquist, dephasing_func, recover_func=None, phantom_motion=None, voxel_motion=None, tissue_masks=None)
+SimData(PD, T1, T2, T2dash, D, B0, B1, coil_sens, voxel_pos, nyquist, dephasing_func, recover_func=None, phantom_motion=None, voxel_motion=None, tissue_masks=None)
 ```
 
 | Parameter | Description |

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "maturin"
 
 [project]
 name = "mrzerocore"
-version = "1.0.5"
+version = "1.0.6"
 description = "Core functionality of MRzero"
 authors = [
     {name = "Jonathan Endres", email = "jonathan.endres@uk-erlangen.de"},

python/MRzeroCore/phantom/custom_voxel_phantom.py

@@ -220,7 +220,7 @@ def build_dephasing_func(shape: str, size: torch.Tensor,
     elif shape == "gauss":
         return lambda t, _: gauss(t, size)
     else:
-        raise ValueError("shape not implemented:", self.voxel_shape)
+        raise ValueError("shape not implemented:", shape)
 
 
 def recover(voxel_size: torch.Tensor, voxel_shape: str, sim_data: SimData

python/MRzeroCore/phantom/tissue_dict.py

@@ -248,7 +248,6 @@ def to_full(sparse):
                                     to_full(sim_data.B0[tissue_begin:tissue_end]),
                                     to_full(sim_data.B1[:, tissue_begin:tissue_end]),
                                     to_full(sim_data.coil_sens[:, tissue_begin:tissue_end]),
-                                    sim_data.size,
                                     sim_data.affine,
                                 )
                             )
@@ -322,15 +321,13 @@ def lp(cfg):
 
     def rs(arr):
         return _resample_nifti(arr, nifti_affine, target_shape, aff_mm)
-    size = target_shape * np.linalg.norm(aff_mm[:3, :3], axis=0) /1000 # np.abs(target_shape @ aff_mm[:3, :3]) / 1000
     density = rs(density)
     T1, T2, T2dash, ADC, B0 = rs(T1), rs(T2), rs(T2dash), rs(ADC), rs(B0)
     B1   = [rs(b) for b in B1]
     coil = [rs(c) for c in coil]
 
     return VoxelGridPhantom(
         PD=torch.as_tensor(density),
-        size=torch.as_tensor(size),
         T1=torch.as_tensor(T1),
         T2=torch.as_tensor(T2),
         T2dash=torch.as_tensor(T2dash),

python/MRzeroCore/phantom/voxel_grid_phantom.py

@@ -99,8 +99,6 @@ class VoxelGridPhantom:
         (coil_count, sx, sy, sz) tensor of RF coil profiles
     coil_sens : torch.Tensor
         (coil_count, sx, sy, sz) tensor of coil sensitivities
-    size : torch.Tensor
-        Size of the data, in meters.
     affine : torch.Tensor
         Affine matrix of the phantom data, in millimeters.
     tissue_masks : Dict[str, torch.Tensor] | None
@@ -117,7 +115,6 @@ def __init__(
         B0: torch.Tensor,
         B1: torch.Tensor,
         coil_sens: torch.Tensor,
-        size: torch.Tensor,
         affine: torch.Tensor,
         phantom_motion=None,
         voxel_motion=None,
@@ -139,11 +136,37 @@ def __init__(
         if self.tissue_masks is None:
             self.tissue_masks = {}
         self.coil_sens = torch.as_tensor(coil_sens, dtype=torch.complex64)
-        self.size = torch.as_tensor(size, dtype=torch.float32)
         self.affine = torch.as_tensor(affine, dtype=torch.float32)
 
         self.phantom_motion = phantom_motion
         self.voxel_motion = voxel_motion
+
+    @property
+    def size(self):
+        shape = torch.tensor(self.PD.shape, dtype=torch.float32)
+        return shape * torch.linalg.norm(self.affine[:3, :3], dim=0) / 1000
+
+    @size.setter
+    def size(self, value):
+        value = torch.as_tensor(value, dtype=torch.float32)
+        shape = torch.tensor(self.PD.shape, dtype=torch.float32)
+
+        old_affine = self.affine
+        directions = old_affine[:3, :3]
+        old_norms = torch.linalg.norm(directions, dim=0)
+        unit_dirs = directions / old_norms
+
+        voxel_size = value / shape * 1000
+        new_affine = old_affine.clone()
+        new_affine[:3, :3] = unit_dirs * voxel_size
+
+        # Corner-referenced translation: for even resolutions the shift is exactly
+        # FOV/2, for odd resolutions it's FOV/2 - voxel_size/2
+        half_index = torch.floor(shape / 2)
+        shift = half_index * voxel_size
+        new_affine[:3, 3] = -unit_dirs @ shift
+
+        self.affine = new_affine
 
     def build(self, PD_threshold: float = 1e-6,
               voxel_shape: Literal["sinc", "box", "point"] = "sinc"
@@ -251,7 +274,7 @@ def load(cls, file_name: str) -> VoxelGridPhantom:
 
         return cls(
             PD, T1, T2, T2dash, D, B0, B1,
-            torch.ones(1, *PD.shape), size, affine,
+            torch.ones(1, *PD.shape), affine,
             tissue_masks=tissue_masks
         )
 
@@ -328,7 +351,6 @@ def load_mat(
             data[..., 3],  # B0
             data[..., 4][None, ...],  # B1
             coil_sens=torch.ones(1, *data.shape[:-1]),
-            size=torch.as_tensor(size),
             affine=affine,
         )
 
@@ -363,8 +385,6 @@ def select_multicoil(tensor: torch.Tensor):
         # reindexed to start at 0, so shift the origin to the first selected
         # slice's world position (assumes contiguous slices; for a sparse list
         # slices[0] is used as the reference).
-        new_size = self.size.clone()
-        new_size[2] = self.size[2] * len(slices) / self.PD.shape[2]
         new_affine = self.affine.clone()
         new_affine[:3, 3] = self.affine[:3, 3] + self.affine[:3, 2] * slices[0]
 
@@ -377,7 +397,6 @@ def select_multicoil(tensor: torch.Tensor):
             select(self.B0),
             select_multicoil(self.B1),
             select_multicoil(self.coil_sens),
-            new_size,
             new_affine,
             tissue_masks={
                 key: mask[..., slices] for key, mask in self.tissue_masks.items()
@@ -420,7 +439,6 @@ def scale(map: torch.Tensor) -> torch.Tensor:
             scale(self.B0),
             scale(self.B1.squeeze()).unsqueeze(0),
             scale(self.coil_sens.squeeze()).unsqueeze(0),
-            self.size.clone(),
             rescale_affine(self.affine, self.PD.shape, (x, y, z)),
             tissue_masks={
                 key: scale(mask) for key, mask in self.tissue_masks.items()
@@ -503,7 +521,6 @@ def resample_masks(tensors: Dict) -> Optional[Dict]:
             resample(self.B0),
             resample_multicoil(self.B1),
             resample_multicoil(self.coil_sens),
-            self.size.clone(),
             rescale_affine(self.affine, self.PD.shape, (x, y, z)),
             tissue_masks=resample_masks(self.tissue_masks)
         )
@@ -639,7 +656,6 @@ def to_full(sparse):
         to_full(sim_data.B0),
         to_full(sim_data.B1),
         to_full(sim_data.coil_sens),
-        sim_data.size,
         sim_data.affine,
     )