Merge pull request #6 from catalystneuro/add_motion_correction

Add motion correction from TIF

src/howe_lab_to_nwb/vu2024/metadata/vu2024_ophys_metadata.yaml

@@ -0,0 +1,19 @@
+Ophys:
+  Device:
+    - name: HamamatsuMicroscope
+      manufacturer: Hamamatsu Photonics
+  ImagingPlane:
+    - name: ImagingPlane
+      description: Imaging plane for the two-photon microscope.
+      device: HamamatsuMicroscope
+      # excitation_lambda
+      # indicatior
+      # location
+  TwoPhotonSeries:
+    - name: TwoPhotonSeries
+      description: Two-photon imaging acquired at 30 Hz with Hamamatsu microscope.
+      imaging_plane: ImagingPlane
+    - name: TwoPhotonSeriesMotionCorrected
+      description: The motion corrected two-photon imaging data.
+      imaging_plane: ImagingPlane
+      unit: n.a.

src/howe_lab_to_nwb/vu2024/vu2024_convert_single_wavelength_session.py

@@ -7,6 +7,7 @@
 from neuroconv.utils import load_dict_from_file, dict_deep_update
 
 from howe_lab_to_nwb.vu2024 import Vu2024NWBConverter
+from howe_lab_to_nwb.vu2024.extractors.bioformats_utils import extract_ome_metadata, parse_ome_metadata
 from howe_lab_to_nwb.vu2024.utils import get_fiber_locations
 from howe_lab_to_nwb.vu2024.utils.add_fiber_photometry import update_fiber_photometry_metadata
 
@@ -18,7 +19,9 @@ def single_wavelength_session_to_nwb(
     excitation_wavelength_in_nm: int,
     indicator: str,
     ttl_file_path: Union[str, Path],
+    motion_corrected_imaging_file_path: Union[str, Path],
     nwbfile_path: Union[str, Path],
+    sampling_frequency: float = None,
     stub_test: bool = False,
 ):
     """
@@ -32,6 +35,9 @@ def single_wavelength_session_to_nwb(
         The path to the .xlsx file containing the fiber locations.
     ttl_file_path : Union[str, Path]
         The path to the .mat file containing the TTL signals.
+    sampling_frequency : float, optional
+        The sampling frequency of the data. If None, the sampling frequency will be read from the .cxd file.
+        If missing from the file, the sampling frequency must be provided.
     """
 
     raw_fiber_photometry_file_path = Path(raw_fiber_photometry_file_path)
@@ -40,8 +46,10 @@ def single_wavelength_session_to_nwb(
     conversion_options = dict()
 
     # Add raw imaging data
-    source_data.update(dict(Imaging=dict(file_path=str(raw_imaging_file_path))))
-    conversion_options.update(dict(Imaging=dict(stub_test=stub_test)))
+    imaging_source_data = dict(file_path=str(raw_imaging_file_path))
+    if sampling_frequency is not None:
+        imaging_source_data.update(sampling_frequency=sampling_frequency)
+    source_data.update(dict(Imaging=imaging_source_data))
 
     # Add raw fiber photometry
     source_data.update(
@@ -52,6 +60,26 @@ def single_wavelength_session_to_nwb(
             )
         )
     )
+    conversion_options.update(dict(FiberPhotometry=dict(stub_test=stub_test)))
+    conversion_options.update(dict(Imaging=dict(stub_test=stub_test, photon_series_index=0)))
+
+    # We need the sampling frequency from the raw imaging data
+    if sampling_frequency is None:
+        ome_metadata = extract_ome_metadata(file_path=raw_imaging_file_path)
+        parsed_metadata = parse_ome_metadata(metadata=ome_metadata)
+        sampling_frequency = parsed_metadata["sampling_frequency"]
+
+    # Add motion corrected imaging data
+    source_data.update(
+        dict(
+            ProcessedImaging=dict(
+                file_path=str(motion_corrected_imaging_file_path), sampling_frequency=sampling_frequency
+            )
+        )
+    )
+    conversion_options.update(
+        dict(ProcessedImaging=dict(stub_test=stub_test, photon_series_index=1, parent_container="processing/ophys"))
+    )
 
     # Add fiber locations
     fiber_locations_metadata = get_fiber_locations(fiber_locations_file_path)
@@ -83,6 +111,9 @@ def single_wavelength_session_to_nwb(
 
     metadata = dict_deep_update(metadata, fiber_photometry_metadata)
 
+    ophys_metadata = load_dict_from_file(Path(__file__).parent / "metadata" / "vu2024_ophys_metadata.yaml")
+    metadata = dict_deep_update(metadata, ophys_metadata)
+
     # Run conversion
     converter.run_conversion(
         metadata=metadata, nwbfile_path=nwbfile_path, conversion_options=conversion_options, overwrite=True
@@ -96,6 +127,10 @@ def single_wavelength_session_to_nwb(
     raw_fiber_photometry_file_path = Path("/Volumes/t7-ssd/Howe/DL18/211110/Data00217_crop_MC_ROIs.mat")
     ttl_file_path = Path("/Volumes/t7-ssd/Howe/DL18/211110/GridDL-18_2021.11.10_16.12.31.mat")
     fiber_locations_file_path = Path("/Volumes/t7-ssd/Howe/DL18/DL18_fiber_locations.xlsx")
+    motion_corrected_imaging_file_path = Path("/Volumes/t7-ssd/Howe/DL18/211110/Data00217_crop_MC.tif")
+
+    # The sampling frequency of the raw imaging data must be provided when it cannot be extracted from the .cxd file
+    sampling_frequency = None
 
     excitation_wavelength_in_nm = 470
     indicator = "dLight1.3b"
@@ -110,6 +145,8 @@ def single_wavelength_session_to_nwb(
         fiber_locations_file_path=fiber_locations_file_path,
         excitation_wavelength_in_nm=excitation_wavelength_in_nm,
         indicator=indicator,
+        motion_corrected_imaging_file_path=motion_corrected_imaging_file_path,
         nwbfile_path=nwbfile_path,
+        sampling_frequency=sampling_frequency,
         stub_test=stub_test,
     )

src/howe_lab_to_nwb/vu2024/vu2024nwbconverter.py

@@ -1,6 +1,8 @@
 """Primary NWBConverter class for this dataset."""
 
 from neuroconv import NWBConverter
+from neuroconv.datainterfaces import TiffImagingInterface
+from neuroconv.utils import DeepDict
 
 from howe_lab_to_nwb.vu2024.interfaces import Vu2024FiberPhotometryInterface, CxdImagingInterface
 
@@ -10,6 +12,7 @@ class Vu2024NWBConverter(NWBConverter):
 
     data_interface_classes = dict(
         Imaging=CxdImagingInterface,
+        ProcessedImaging=TiffImagingInterface,
         FiberPhotometry=Vu2024FiberPhotometryInterface,
     )
 
@@ -23,6 +26,17 @@ def get_metadata_schema(self) -> dict:
 
         return metadata_schema
 
+    def get_metadata(self) -> DeepDict:
+        metadata = super().get_metadata()
+
+        # Overwrite the Ophys imaging metadata with the metadata from the CxdImagingInterface
+        imaging_interface_metadata = self.data_interface_objects["Imaging"].get_metadata()
+        metadata["Ophys"]["Device"] = imaging_interface_metadata["Ophys"]["Device"]
+        metadata["Ophys"]["ImagingPlane"] = imaging_interface_metadata["Ophys"]["ImagingPlane"]
+        metadata["Ophys"]["TwoPhotonSeries"] = imaging_interface_metadata["Ophys"]["TwoPhotonSeries"]
+
+        return metadata
+
     def temporally_align_data_interfaces(self):
         imaging = self.data_interface_objects["Imaging"]
         fiber_photometry = self.data_interface_objects["FiberPhotometry"]