Skills › Scientific Agent Skills › nwb-conversion › convert_session.py
convert_session.py
249 lines · skills/nwb-conversion/scripts/convert_session.py at 5f2bc03
1"""Convert planar TIFF imaging and timestamped position CSV to validated NWB."""
2from __future__ import annotations
4import argparse
5import csv
6from datetime import datetime
7import hashlib
8import importlib.metadata
9import json
10from pathlib import Path
13def digest(path):
14 checksum = hashlib.sha256()
15 with Path(path).open("rb") as handle:
16 for chunk in iter(lambda: handle.read(1024 * 1024), b""):
17 checksum.update(chunk)
18 return checksum.hexdigest()
21def read_numeric_csv(path, columns):
22 import numpy as np
23 with Path(path).open(newline="") as handle:
24 reader = csv.DictReader(handle)
25 if reader.fieldnames != columns:
26 raise ValueError(f"{path}: expected exactly {columns}")
27 try:
28 rows = []
29 for row in reader:
30 if None in row or any(row[k] is None for k in columns):
31 raise ValueError("Row width differs from the header")
32 rows.append([float(row[k]) for k in columns])
33 values = np.asarray(rows, dtype=float)
34 except (TypeError, ValueError) as error:
35 raise ValueError(f"{path}: nonnumeric or missing data") from error
36 if values.ndim != 2 or len(values) < 2 or not np.isfinite(values).all():
37 raise ValueError(f"{path}: need at least two finite rows")
38 if (np.diff(values[:, 0]) <= 0).any():
39 raise ValueError(f"{path}: timestamps must be strictly increasing")
40 return values
43def align_behavior(times, synchronization, base):
44 import numpy as np
45 times = np.asarray(times, dtype=float)
46 if times.ndim != 1 or len(times) < 2 or not np.isfinite(times).all() or (np.diff(times) <= 0).any():
47 raise ValueError("Behavior timestamps must be finite and strictly increasing")
48 shared = synchronization.get("shared_clock_evidence")
49 pulses = synchronization.get("pulse_pairs_csv")
50 if bool(shared) == bool(pulses):
51 raise ValueError("Supply either shared_clock_evidence or pulse_pairs_csv")
52 if shared:
53 if not isinstance(shared, str) or not shared.strip():
54 raise ValueError("Shared-clock evidence must be a nonempty description")
55 return times.copy(), {"method": "shared_clock", "evidence": shared, "slope": 1., "offset_s": 0.}
56 pairs = read_numeric_csv(base / pulses, ["device_time_s", "reference_time_s"])
57 if len(pairs) < 3 or (np.diff(pairs[:, 1]) <= 0).any():
58 raise ValueError("At least three ordered, matched pulse pairs are required")
59 tolerance = float(synchronization["max_residual_s"])
60 if not np.isfinite(tolerance) or tolerance <= 0:
61 raise ValueError("max_residual_s must be finite and positive")
62 # Center the device clock to avoid fitting against a large absolute origin.
63 device_origin, reference_origin = pairs[0]
64 slope, centered_offset = np.polyfit(pairs[:, 0] - device_origin, pairs[:, 1] - reference_origin, 1)
65 offset = reference_origin + centered_offset - slope * device_origin
66 residuals = slope * (pairs[:, 0] - device_origin) + centered_offset - (pairs[:, 1] - reference_origin)
67 residual = float(np.max(np.abs(residuals)))
68 if not np.isfinite([slope, offset, residual]).all() or slope <= 0 or residual > tolerance:
69 raise ValueError(f"Clock fit failed: maximum residual {residual:g} s")
70 if times[0] < pairs[0, 0] or times[-1] > pairs[-1, 0]:
71 raise ValueError("Behavior times exceed pulse support; affine extrapolation is not validated")
72 aligned = slope * (times - device_origin) + centered_offset + reference_origin
73 if not np.isfinite(aligned).all() or (np.diff(aligned) <= 0).any():
74 raise ValueError("Clock mapping produced invalid or indistinguishable timestamps")
75 return aligned, {"method": "matched_pulse_affine", "slope": float(slope),
76 "offset_s": float(offset), "max_residual_s": residual, "tolerance_s": tolerance,
77 "device_origin_s": float(device_origin), "reference_origin_s": float(reference_origin),
78 "centered_offset_s": float(centered_offset), "pulse_residuals_s": residuals.tolist(),
79 "device_support_s": pairs[[0, -1], 0].tolist(), "reference_support_s": pairs[[0, -1], 1].tolist(),
80 "pulse_pairs": len(pairs), "pulse_sha256": digest(base / pulses)}
83def convert(config_path, output):
84 import numpy as np
85 import tifffile
86 from neuroconv.datainterfaces import TiffImagingInterface
87 from pynwb import NWBFile, NWBHDF5IO, validate
88 from pynwb.behavior import Position, SpatialSeries
89 from pynwb.file import Subject
90 from nwbinspector import inspect_nwbfile
92 config_path, output = Path(config_path), Path(output)
93 config = json.loads(config_path.read_text())
94 base = config_path.resolve().parent
95 if output.suffix.lower() != ".nwb":
96 raise ValueError("Output must use the .nwb extension")
97 report_path = output.with_suffix(".validation.json")
98 for path in (output, report_path):
99 if path.exists():
100 raise FileExistsError(f"Refusing to overwrite {path}")
101 start = datetime.fromisoformat(config["session_start_time"])
102 if start.utcoffset() is None:
103 raise ValueError("session_start_time needs an explicit timezone")
104 for field in ("identifier", "session_description", "institution", "lab", "position_reference_frame", "position_description"):
105 if not isinstance(config.get(field), str) or not config[field].strip():
106 raise ValueError(f"Missing {field}")
107 if not isinstance(config.get("experimenter"), list) or not config["experimenter"] or any(
108 not isinstance(name, str) or not name.strip() for name in config["experimenter"]
109 ):
110 raise ValueError("experimenter must be a nonempty list of names")
111 optics = config["imaging"]
112 for key in ("num_channels", "num_planes"):
113 if type(optics.get(key)) is not int or optics[key] != 1:
114 raise ValueError(f"imaging.{key} must explicitly equal 1; interleaved channels/planes are unsupported")
115 tiff_path = base / config["tiff"]
116 frame_times = read_numeric_csv(base / config["frame_times_csv"], ["time_s"])[:, 0]
117 behavior = read_numeric_csv(base / config["position_csv"], ["time_s", "x", "y"])
118 with tifffile.TiffFile(tiff_path) as tiff:
119 if len(tiff.pages) != len(frame_times) or len(tiff.pages[0].shape) != 2:
120 raise ValueError("Expected one grayscale 2D TIFF page per frame timestamp")
121 shape, dtype = tiff.pages[0].shape, tiff.pages[0].dtype
122 if any(p.shape != shape or p.dtype != dtype for p in tiff.pages):
123 raise ValueError("TIFF pages differ in shape or dtype")
124 if any(p.samplesperpixel != 1 or p.photometric not in (0, 1) for p in tiff.pages):
125 raise ValueError("TIFF pages must contain grayscale intensity samples")
126 if len(tiff.series) != 1 or any(
127 size > 1 and axis in "CZS" for axis, size in zip(tiff.series[0].axes, tiff.series[0].shape)
128 ):
129 raise ValueError("TIFF metadata identifies multiple series, channels, samples or depth planes")
130 position_units = {"m": 1., "cm": .01, "mm": .001}
131 if config["position_unit"] not in position_units:
132 raise ValueError("Position unit must be m, cm or mm; pixels require calibrated coordinates")
133 position_m = behavior[:, 1:] * position_units[config["position_unit"]]
134 aligned, alignment = align_behavior(behavior[:, 0], config["synchronization"], base)
135 if min(frame_times[0], aligned[0]) < 0:
136 raise ValueError("This converter requires nonnegative seconds since session_start_time")
137 if optics.get("modality") != "two-photon":
138 raise ValueError("imaging.modality must explicitly be two-photon for this converter")
139 for key in ("device", "device_description", "description", "indicator", "location", "unit", "optical_channel_description"):
140 if not isinstance(optics.get(key), str) or not optics[key].strip():
141 raise ValueError(f"Missing imaging.{key}; do not invent acquisition metadata")
142 for key in ("excitation_nm", "emission_nm"):
143 if not np.isfinite(float(optics[key])) or float(optics[key]) <= 0:
144 raise ValueError(f"imaging.{key} must be positive and finite")
145 subject_metadata = dict(config["subject"])
146 if "date_of_birth" in subject_metadata:
147 birth = datetime.fromisoformat(subject_metadata["date_of_birth"])
148 if birth.utcoffset() is None:
149 raise ValueError("subject.date_of_birth needs an explicit timezone")
150 subject_metadata["date_of_birth"] = birth
151 nwb = NWBFile(session_description=config["session_description"], identifier=config["identifier"],
152 session_start_time=start, timestamps_reference_time=start, experimenter=config["experimenter"],
153 institution=config["institution"], lab=config["lab"], subject=Subject(**subject_metadata))
154 interface = TiffImagingInterface(file_paths=[str(tiff_path)], sampling_frequency=float(1 / np.median(np.diff(frame_times))),
155 num_channels=1, num_planes=1, metadata_key="imaging", verbose=False)
156 interface.set_aligned_timestamps(aligned_timestamps=frame_times)
157 metadata = interface.get_metadata()
158 metadata["Devices"] = {"microscope": {"name": optics["device"], "description": optics["device_description"]}}
159 metadata["Ophys"]["ImagingPlanes"] = {"plane": {
160 "name": "ImagingPlane", "description": optics["description"], "device_metadata_key": "microscope",
161 "excitation_lambda": float(optics["excitation_nm"]), "indicator": optics["indicator"], "location": optics["location"],
162 "optical_channel": [{"name": "channel", "description": optics["optical_channel_description"],
163 "emission_lambda": float(optics["emission_nm"])}]}}
164 metadata["Ophys"]["MicroscopySeries"]["imaging"].update(
165 name="Imaging", description=optics["description"], unit=optics["unit"], imaging_plane_metadata_key="plane")
166 interface.add_to_nwbfile(nwbfile=nwb, metadata=metadata, always_write_timestamps=True)
167 series = SpatialSeries(name="Position", data=position_m, timestamps=aligned,
168 unit="meters", reference_frame=config["position_reference_frame"],
169 description=f"Source coordinates in {config['position_unit']}; converted to meters. {config['position_description']}")
170 nwb.create_processing_module(name="behavior", description="Measured animal position").add(Position(spatial_series=series))
171 provenance = {"source_sha256": {key: digest(base / config[key]) for key in ("tiff", "frame_times_csv", "position_csv")},
172 "config_sha256": digest(config_path), "configuration": config, "alignment": alignment,
173 "converter_sha256": digest(__file__), "timestamps_reference_time": start.isoformat(),
174 "axis_mapping": "TIFF (time,y,x) -> NWB (time,x,y)",
175 "versions": {name: importlib.metadata.version(name) for name in ("neuroconv", "pynwb", "nwbinspector", "roiextractors", "tifffile", "zarr", "hdmf-zarr")}}
176 nwb.add_scratch(json.dumps(provenance), name="conversion_provenance", description="Source metadata, checksums and clock mapping")
177 output.parent.mkdir(parents=True, exist_ok=True)
178 with NWBHDF5IO(str(output), "w-") as io:
179 io.write(nwb)
180 errors = [str(error) for error in validate(path=str(output))]
181 findings = []
182 for message in inspect_nwbfile(output, skip_validate=True):
183 if message is not None:
184 findings.append({"importance": message.importance.name, "check": message.check_function_name,
185 "message": message.message, "location": message.location})
186 with NWBHDF5IO(str(output), "r") as io, tifffile.TiffFile(tiff_path) as tiff:
187 restored = io.read()
188 image = restored.acquisition["Imaging"]
189 np.testing.assert_array_equal(image.timestamps[:], frame_times)
190 for index, page in enumerate(tiff.pages):
191 np.testing.assert_array_equal(image.data[index], page.asarray().T)
192 restored_position = restored.processing["behavior"]["Position"]["Position"]
193 np.testing.assert_array_equal(restored_position.data[:], position_m)
194 np.testing.assert_array_equal(restored_position.timestamps[:], aligned)
195 if image.unit != optics["unit"] or restored_position.unit != "meters":
196 raise ValueError("Units changed during round trip")
197 if image.data.dtype != dtype or any(s.conversion != 1.0 or s.offset != 0.0 for s in (image, restored_position)):
198 raise ValueError("Data dtype or unit scaling changed during round trip")
199 plane = image.imaging_plane
200 if (plane.device.name != optics["device"] or plane.location != optics["location"]
201 or plane.device.description != optics["device_description"] or plane.description != optics["description"]
202 or plane.indicator != optics["indicator"] or plane.excitation_lambda != float(optics["excitation_nm"])
203 or len(plane.optical_channel) != 1
204 or plane.optical_channel[0].description != optics["optical_channel_description"]
205 or plane.optical_channel[0].emission_lambda != float(optics["emission_nm"])):
206 raise ValueError("Optical metadata changed during round trip")
207 if restored_position.reference_frame != config["position_reference_frame"]:
208 raise ValueError("Position reference frame changed during round trip")
209 if (restored.identifier != config["identifier"] or restored.session_start_time != start
210 or restored.timestamps_reference_time != start):
211 raise ValueError("Session metadata changed during round trip")
212 for key in ("session_description", "institution", "lab"):
213 if getattr(restored, key) != config[key]:
214 raise ValueError(f"Session {key} changed during round trip")
215 if list(restored.experimenter) != config["experimenter"]:
216 raise ValueError("Experimenter metadata changed during round trip")
217 for key in subject_metadata:
218 if getattr(restored.subject, key) != getattr(nwb.subject, key):
219 raise ValueError(f"Subject {key} changed during round trip")
220 if json.loads(restored.scratch["conversion_provenance"].data) != provenance:
221 raise ValueError("Conversion provenance changed during round trip")
222 for finding in findings:
223 if finding["check"] == "check_data_orientation" and finding["location"] == "/acquisition/Imaging":
224 finding["review_evidence"] = "Time axis verified against all frame timestamps and every source TIFF page; short recordings can trigger the longest-axis heuristic"
225 report = {**provenance, "schema_errors": errors, "inspector_findings": findings,
226 "inspector_requires_review": any(f["importance"] in {"ERROR", "PYNWB_VALIDATION", "CRITICAL"} for f in findings),
227 "roundtrip": {"all_frames_equal": True, "position_equal": True, "timestamps_equal": True,
228 "metadata_equal": True, "provenance_equal": True, "identity_unit_scaling": True,
229 "frame_count": len(frame_times), "position_samples": len(behavior)}}
230 with report_path.open("x") as handle:
231 handle.write(json.dumps(report, indent=2) + "\n")
232 if errors:
233 raise ValueError("NWB schema validation failed; inspect the validation JSON")
234 return report
237def main():
238 parser = argparse.ArgumentParser(description=__doc__)
239 parser.add_argument("config", type=Path)
240 parser.add_argument("--output", type=Path, required=True)
241 args = parser.parse_args()
242 report = convert(args.config, args.output)
243 print(json.dumps({"output": str(args.output), "schema_errors": report["schema_errors"],
244 "inspector_requires_review": report["inspector_requires_review"],
245 "inspector_findings": len(report["inspector_findings"]), "roundtrip": report["roundtrip"]}, indent=2))
248if __name__ == "__main__":
249 main()