Add zone plate test image

Author: fepegar

src/torchio/datasets/__init__.py

@@ -20,6 +20,7 @@
 from .rsna_miccai import RSNAMICCAI
 from .rsna_spine_fracture import RSNACervicalSpineFracture
 from .slicer import Slicer
+from .zone_plate import ZonePlate
 
 __all__ = [
     'FPG',
@@ -44,4 +45,5 @@
     'FractureMNIST3D',
     'VesselMNIST3D',
     'SynapseMNIST3D',
+    'ZonePlate',
 ]

src/torchio/datasets/zone_plate.py

@@ -0,0 +1,52 @@
+import numpy as np
+
+from ..data import ScalarImage
+from ..data.subject import Subject
+
+
+class ZonePlate(Subject):
+    """Synthetic data generated from a zone plate.
+
+    The zone plate is a circular diffraction grating that produces concentric
+    rings of light and dark bands. This dataset is useful for testing image
+    processing algorithms, particularly those related to frequency analysis and
+    interpolation.
+
+    See equation 10.63 in `Practical Handbook on Image Processing for
+    Scientific Applications <https://www.routledge.com/Practical-Handbook-on-Image-Processing-for-Scientific-and-Technical-Applications/Jahne/p/book/9780849319006?srsltid=AfmBOoptrtzILIlMx9FYqvx6UrGbevfD66x2k242iprFdn_CfyOWXjjH>`_
+    by Bernd Jähne.
+
+    Args:
+        size: The size of the generated image along all dimensions.
+    """
+
+    def __init__(self, size: int = 501):
+        if size < 3:
+            raise ValueError('Size must be at least 3.')
+        self.size = size
+        image = self._generate_image(size)
+        super().__init__(image=image)
+
+    @staticmethod
+    def _generate_image(size: int) -> ScalarImage:
+        if size % 2 == 1:
+            fin = (size - 1) // 2
+            ini = -fin
+        else:
+            fin = size // 2
+            ini = -fin + 1
+        x = np.arange(ini, fin)
+        y = np.arange(ini, fin)
+        z = np.arange(ini, fin)
+        X, Y, Z = np.meshgrid(x, y, z)
+        r = np.sqrt(X**2 + Y**2 + Z**2)
+        km = 0.8 * np.pi
+        rm = ini
+        w = rm / 10
+        term1 = np.sin((km * r**2) / (2 * rm))
+        term2 = 0.5 * np.tanh((rm - r) / w) + 0.5
+        g = term1 * term2
+        affine = np.eye(4)
+        origin = np.array([ini, ini, ini])
+        affine[:3, 3] = origin
+        return ScalarImage(tensor=g[np.newaxis], affine=affine)