Pixel Size Attributes

gallery/tutorials/aspire_introduction.py

@@ -571,7 +571,7 @@ def noise_function(x, y):
 
 # Generate several CTFs.
 ctf_filters = [
-    RadialCTFFilter(pixel_size=5, defocus=d)
+    RadialCTFFilter(pixel_size=vol_ds.pixel_size, defocus=d)
     for d in np.linspace(defocus_min, defocus_max, defocus_ct)
 ]
 

gallery/tutorials/pipeline_demo.py

@@ -63,7 +63,7 @@
 defocus_ct = 7
 
 ctf_filters = [
-    RadialCTFFilter(pixel_size=5, defocus=d)
+    RadialCTFFilter(pixel_size=vol.pixel_size, defocus=d)
     for d in np.linspace(defocus_min, defocus_max, defocus_ct)
 ]
 

gallery/tutorials/tutorials/cov3d_simulation.py

@@ -47,7 +47,9 @@
     L=img_size,
     n=num_imgs,
     vols=vols,
-    unique_filters=[RadialCTFFilter(defocus=d) for d in np.linspace(1.5e4, 2.5e4, 7)],
+    unique_filters=[
+        RadialCTFFilter(pixel_size=10, defocus=d) for d in np.linspace(1.5e4, 2.5e4, 7)
+    ],
     dtype=dtype,
 )
 

src/aspire/image/image.py

@@ -80,7 +80,7 @@ def load_mrc(filepath):
     Load raw data from `.mrc` into an array.
 
     :param filepath: File path (string).
-    :return: numpy array of image data.
+    :return: (numpy array of image data, pixel_size)
     """
 
     # mrcfile tends to yield many warnings about EMPIAR datasets being corrupt
@@ -92,6 +92,7 @@ def load_mrc(filepath):
 
         with mrcfile.open(filepath, mode="r", permissive=True) as mrc:
             im = mrc.data
+            pixel_size = Image._vx_array_to_size(mrc.voxel_size)
 
         # Log each mrcfile warning to debug log, noting the associated file
         for w in ws:
@@ -110,19 +111,29 @@ def load_mrc(filepath):
                 f" Will attempt to continue processing {filepath}"
             )
 
-    return im
+    return im, pixel_size
 
 
 def load_tiff(filepath):
     """
     Load raw data from `.tiff` into an array.
 
+    Note, TIFF does not natively provide equivalent to pixel/voxel_size,
+    so users of TIFF files may need to manually assign `pixel_size` to
+    `Image` instances when required. Defaults to `pixel_size=None`.
+
     :param filepath: File path (string).
-    :return: numpy array of image data.
+    :return: (numpy array of image data, pixel_size=None)
     """
+    # Use PIL to open `filepath`
+    img = PILImage.open(filepath)
+
+    # Future todo, extract `voxel_size` if available in TIFF tags (custom tag?)
+    # For now, default to `None`.
+    pixel_size = None
 
-    # Use PIL to open `filepath` and cast to numpy array.
-    return np.array(PILImage.open(filepath))
+    # Cast image data as numpy array
+    return np.array(img), pixel_size
 
 
 class Image:
@@ -133,7 +144,7 @@ class Image:
         ".tiff": load_tiff,
     }
 
-    def __init__(self, data, dtype=None):
+    def __init__(self, data, pixel_size=None, dtype=None):
         """
         A stack of one or more images.
 
@@ -149,6 +160,10 @@ def __init__(self, data, dtype=None):
 
         :param data: Numpy array containing image data with shape
             `(..., resolution, resolution)`.
+        :param pixel_size: Optional pixel size in angstroms.
+            When provided will be saved with `mrc` metadata.
+            Default of `None` will not write to file,
+            but will be considered unit pixels (1) for FSC.
         :param dtype: Optionally cast `data` to this dtype.
             Defaults to `data.dtype`.
 
@@ -180,6 +195,9 @@ def __init__(self, data, dtype=None):
         self.stack_shape = self._data.shape[:-2]
         self.n_images = np.prod(self.stack_shape)
         self.resolution = self._data.shape[-1]
+        self.pixel_size = None
+        if pixel_size is not None:
+            self.pixel_size = float(pixel_size)
 
         # Numpy interop
         # https://numpy.org/devdocs/user/basics.interoperability.html#the-array-interface-protocol
@@ -238,7 +256,7 @@ def _check_key_dims(self, key):
 
     def __getitem__(self, key):
         self._check_key_dims(key)
-        return self.__class__(self._data[key])
+        return self.__class__(self._data[key], pixel_size=self.pixel_size)
 
     def __setitem__(self, key, value):
         self._check_key_dims(key)
@@ -266,31 +284,34 @@ def stack_reshape(self, *args):
                 f"Number of images {self.n_images} cannot be reshaped to {shape}."
             )
 
-        return self.__class__(self._data.reshape(*shape, *self._data.shape[-2:]))
+        return self.__class__(
+            self._data.reshape(*shape, *self._data.shape[-2:]),
+            pixel_size=self.pixel_size,
+        )
 
     def __add__(self, other):
         if isinstance(other, Image):
             other = other._data
 
-        return self.__class__(self._data + other)
+        return self.__class__(self._data + other, pixel_size=self.pixel_size)
 
     def __sub__(self, other):
         if isinstance(other, Image):
             other = other._data
 
-        return self.__class__(self._data - other)
+        return self.__class__(self._data - other, pixel_size=self.pixel_size)
 
     def __mul__(self, other):
         if isinstance(other, Image):
             other = other._data
 
-        return self.__class__(self._data * other)
+        return self.__class__(self._data * other, pixel_size=self.pixel_size)
 
     def __neg__(self):
-        return self.__class__(-self._data)
+        return self.__class__(-self._data, pixel_size=self.pixel_size)
 
     def sqrt(self):
-        return self.__class__(np.sqrt(self._data))
+        return self.__class__(np.sqrt(self._data), pixel_size=self.pixel_size)
 
     @property
     def T(self):
@@ -312,7 +333,9 @@ def transpose(self):
         im = self.stack_reshape(-1)
         imt = np.transpose(im._data, (0, -1, -2))
 
-        return self.__class__(imt).stack_reshape(original_stack_shape)
+        return self.__class__(imt, pixel_size=self.pixel_size).stack_reshape(
+            original_stack_shape
+        )
 
     def flip(self, axis=-2):
         """
@@ -335,11 +358,15 @@ def flip(self, axis=-2):
                     f"Cannot flip axis {ax}: stack axis. Did you mean {ax-3}?"
                 )
 
-        return self.__class__(np.flip(self._data, axis))
+        return self.__class__(np.flip(self._data, axis), pixel_size=self.pixel_size)
 
     def __repr__(self):
+        px_msg = "."
+        if self.pixel_size is not None:
+            px_msg = f" with pixel_size={self.pixel_size} angstroms."
+
         msg = f"{self.n_images} {self.dtype} images arranged as a {self.stack_shape} stack"
-        msg += f" each of size {self.resolution}x{self.resolution}."
+        msg += f" each of size {self.resolution}x{self.resolution}{px_msg}"
         return msg
 
     def asnumpy(self):
@@ -355,7 +382,7 @@ def asnumpy(self):
         return view
 
     def copy(self):
-        return self.__class__(self._data.copy())
+        return self.__class__(self._data.copy(), pixel_size=self.pixel_size)
 
     def shift(self, shifts):
         """
@@ -412,7 +439,14 @@ def downsample(self, ds_res, zero_nyquist=True):
         out = fft.centered_ifft2(crop_fx).real * (ds_res**2 / self.resolution**2)
         out = xp.asnumpy(out)
 
-        return self.__class__(out).stack_reshape(original_stack_shape)
+        # Optionally scale pixel size
+        ds_pixel_size = self.pixel_size
+        if ds_pixel_size is not None:
+            ds_pixel_size *= self.resolution / ds_res
+
+        return self.__class__(out, pixel_size=ds_pixel_size).stack_reshape(
+            original_stack_shape
+        )
 
     def filter(self, filter):
         """
@@ -441,7 +475,9 @@ def filter(self, filter):
 
         im = xp.asnumpy(im.real)
 
-        return self.__class__(im).stack_reshape(original_stack_shape)
+        return self.__class__(im, pixel_size=self.pixel_size).stack_reshape(
+            original_stack_shape
+        )
 
     def rotate(self):
         raise NotImplementedError
@@ -453,6 +489,9 @@ def save(self, mrcs_filepath, overwrite=False):
         with mrcfile.new(mrcs_filepath, overwrite=overwrite) as mrc:
             # original input format (the image index first)
             mrc.set_data(self._data.astype(np.float32))
+            # Note assigning voxel_size must come after `set_data`
+            if self.pixel_size is not None:
+                mrc.voxel_size = self.pixel_size
 
     @staticmethod
     def load(filepath, dtype=None):
@@ -477,14 +516,14 @@ def load(filepath, dtype=None):
             )
 
         # Call the appropriate file reader
-        im = Image.extensions[ext](filepath)
+        im, pixel_size = Image.extensions[ext](filepath)
 
         # Attempt casting when user provides dtype
         if dtype is not None:
             im = im.astype(dtype, copy=False)
 
         # Return as Image instance
-        return Image(im)
+        return Image(im, pixel_size=pixel_size)
 
     def _im_translate(self, shifts):
         """
@@ -535,7 +574,9 @@ def _im_translate(self, shifts):
         im_translated = xp.asnumpy(im_translated.real)
 
         # Reshape to stack shape
-        return self.__class__(im_translated).stack_reshape(stack_shape)
+        return self.__class__(im_translated, pixel_size=self.pixel_size).stack_reshape(
+            stack_shape
+        )
 
     def norm(self):
         return anorm(self._data)
@@ -602,7 +643,9 @@ def backproject(self, rot_matrices, symmetry_group=None, zero_nyquist=True):
 
         vol /= L
 
-        return aspire.volume.Volume(vol, symmetry_group=symmetry_group)
+        return aspire.volume.Volume(
+            vol, pixel_size=self.pixel_size, symmetry_group=symmetry_group
+        )
 
     def show(self, columns=5, figsize=(20, 10), colorbar=True):
         """
@@ -645,7 +688,7 @@ def show(self, columns=5, figsize=(20, 10), colorbar=True):
 
             plt.show()
 
-    def frc(self, other, cutoff=None, pixel_size=None, method="fft", plot=False):
+    def frc(self, other, cutoff=None, method="fft", plot=False):
         r"""
         Compute the Fourier ring correlation between two images.
 
@@ -663,8 +706,6 @@ def frc(self, other, cutoff=None, pixel_size=None, method="fft", plot=False):
             Default `None` implies `cutoff=1` and excludes
             plotting cutoff line.
 
-        :param pixel_size: Pixel size in angstrom.  Default `None`
-            implies unit in pixels, equivalent to pixel_size=1.
         :param method: Selects either 'fft' (on cartesian grid),
             or 'nufft' (on polar grid). Defaults to 'fft'.
         :param plot: Optionally plot to screen or file.
@@ -684,7 +725,7 @@ def frc(self, other, cutoff=None, pixel_size=None, method="fft", plot=False):
         frc = FourierRingCorrelation(
             a=self.asnumpy(),
             b=other.asnumpy(),
-            pixel_size=pixel_size,
+            pixel_size=self.pixel_size,
             method=method,
         )
 
@@ -695,6 +736,32 @@ def frc(self, other, cutoff=None, pixel_size=None, method="fft", plot=False):
 
         return frc.analyze_correlations(cutoff), frc.correlations
 
+    @staticmethod
+    def _vx_array_to_size(vx):
+        """
+        Utility to convert from several possible `mrcfile.voxel_size`
+        representations to a single (float) value or None.
+        """
+
+        # Convert from recarray to single values,
+        #   checks uniformity.
+        if isinstance(vx, np.recarray):
+            if vx.x != vx.y:
+                raise ValueError(f"Voxel sizes are not uniform: {vx}")
+            vx = vx.x
+
+        # Convert `0` to `None`
+        if (
+            isinstance(vx, int) or isinstance(vx, float) or isinstance(vx, np.ndarray)
+        ) and vx == 0:
+            vx = None
+
+        # Consistently return a `float` when not None
+        if vx is not None:
+            vx = float(vx)
+
+        return vx
+
 
 class CartesianImage(Image):
     def expand(self, basis):

src/aspire/operators/filters.py

@@ -403,7 +403,7 @@ def __init__(self, dim=None):
 class CTFFilter(Filter):
     def __init__(
         self,
-        pixel_size=10,
+        pixel_size=1,
         voltage=200,
         defocus_u=15000,
         defocus_v=15000,
@@ -415,7 +415,7 @@ def __init__(
         """
         A CTF (Contrast Transfer Function) Filter
 
-        :param pixel_size:  Pixel size in angstrom
+        :param pixel_size:  Pixel size in angstrom, default 1.
         :param voltage:     Electron voltage in kV
         :param defocus_u:   Defocus depth along the u-axis in angstrom
         :param defocus_v:   Defocus depth along the v-axis in angstrom
@@ -425,7 +425,7 @@ def __init__(
         :param B:           Envelope decay in inverse square angstrom (default 0)
         """
         super().__init__(dim=2, radial=defocus_u == defocus_v)
-        self.pixel_size = pixel_size
+        self.pixel_size = float(pixel_size)
         self.voltage = voltage
         self.wavelength = voltage_to_wavelength(self.voltage)
         self.defocus_u = defocus_u
@@ -482,7 +482,7 @@ def scale(self, c=1):
 
 class RadialCTFFilter(CTFFilter):
     def __init__(
-        self, pixel_size=10, voltage=200, defocus=15000, Cs=2.26, alpha=0.07, B=0
+        self, pixel_size=1, voltage=200, defocus=15000, Cs=2.26, alpha=0.07, B=0
     ):
         super().__init__(
             pixel_size=pixel_size,

src/aspire/source/coordinates.py

@@ -490,7 +490,9 @@ def _images(self, indices):
                     cropped = self._crop_micrograph(arr, next(coord))
                     im[i] = cropped
         # Finally, apply transforms to resulting Image
-        return self.generation_pipeline.forward(Image(im), indices)
+        return self.generation_pipeline.forward(
+            Image(im, pixel_size=self.pixel_size), indices
+        )
 
     @staticmethod
     def _is_number(text):