diff --git a/README.md b/README.md index ba98c79..e39212b 100644 --- a/README.md +++ b/README.md @@ -91,7 +91,7 @@ PACFISH API. We have examples for both `Python` and `MATLAB`. ## Use case: Implement a conversion adapter - impot pacfish as pf + import pacfish as pf class DeviceSpecificAdapter(pf.BaseAdapter): diff --git a/examples/python/compare_python_matlab_hdf5_file.py b/examples/python/compare_python_matlab_hdf5_file.py deleted file mode 100644 index 003af73..0000000 --- a/examples/python/compare_python_matlab_hdf5_file.py +++ /dev/null @@ -1,22 +0,0 @@ -""" -SPDX-FileCopyrightText: 2021 International Photoacoustics Standardisation Consortium (IPASC) -SPDX-License-Identifier: CC0 -""" -import pacfish as pf -import numpy as np - -from testing.unit_tests.utils import assert_equal_dicts - -FILE_PATH_MATLAB_WRITTEN_FILE = "C:/PACFISH/pacfish_matlab/+pacfish/test.hdf5" -FILE_PATH_PYTHON_WRITTEN_FILE = "C:/standardised-image-reconstruction/1BdSLl4BSxpxXDwWcBKKVV4nHULPe7IS8_ipasc.hdf5" - -pa_data = pf.load_data(FILE_PATH_MATLAB_WRITTEN_FILE) -pa_data_2 = pf.load_data(FILE_PATH_PYTHON_WRITTEN_FILE) - -pf.visualize_device(pa_data.meta_data_device) -pf.visualize_device(pa_data_2.meta_data_device) - -# Testing if the two read files are actually the same -assert np.equal(pa_data.binary_time_series_data, pa_data_2.binary_time_series_data).all() -assert_equal_dicts(pa_data.meta_data_acquisition, pa_data_2.meta_data_acquisition) -assert_equal_dicts(pa_data.meta_data_device, pa_data_2.meta_data_device) diff --git a/examples/python/input/2017_12_11-10_54_55_0.lom b/examples/python/input/2017_12_11-10_54_55_0.lom new file mode 100644 index 0000000..7cdcfe7 Binary files /dev/null and b/examples/python/input/2017_12_11-10_54_55_0.lom differ diff --git a/examples/python/input/20211102134837219776.dcm.lom b/examples/python/input/20211102134837219776.dcm.lom new file mode 100644 index 0000000..e285b89 Binary files /dev/null and b/examples/python/input/20211102134837219776.dcm.lom differ diff --git a/examples/python/input/I0000001.dcm.lom b/examples/python/input/I0000001.dcm.lom new file mode 100644 index 0000000..acabddf Binary files /dev/null and b/examples/python/input/I0000001.dcm.lom differ diff --git a/examples/python/read_imagio_file.py b/examples/python/read_imagio_file.py new file mode 100644 index 0000000..2723d06 --- /dev/null +++ b/examples/python/read_imagio_file.py @@ -0,0 +1,51 @@ +# +# This file reads an example Laser Optic Movie (LOM) produced by the Seno (https://senomedical.com/) +# Imagio System and produces an IPASC-compatible HDF5 file. In addition, if enabled, graphical +# representations of the data are written as .png files. +# + +import os +import numpy as np +import cv2 +import pathlib +from pacfish.api.adapters.Imagio_File_Converter import ImagioFileConverter +from pacfish import write_data +from pacfish import quality_check_pa_data + +# input file is a Laser Optic Movie (LOM) from a scan of a phantom using the Seno Imagio system +# input_file = '../../examples/python/input/I0000001.dcm.lom' # Gen2 (4 OA frames, 1 US image, hair target) + +# other examples +#input_file = '../../examples/python/input/2017_12_11-10_54_55_0.lom' # Gen1a (many OA frames, many US images, lesion) +input_file = '../../examples/python/input/20211102134837219776.dcm.lom' # Gen2 (many OA frames, many US images, hair target) + +converter = ImagioFileConverter(input_file) + +pa_data = converter.generate_pa_data() + +quality_check_pa_data(pa_data, verbose=True, log_file_path="") + +folder = "output/" + os.path.basename(input_file) +pathlib.Path(folder).mkdir(parents=True, exist_ok=True) + +output_png = True # disable if needed +if (output_png): + timestamps = pa_data.get_measurement_time_stamps() + wavelengths = pa_data.get_acquisition_wavelengths() + for idx, oa_frame in enumerate(pa_data.binary_time_series_data[0, 0, :]): + iTick = str(int(timestamps[idx] * 1E3)) + "_msec" + wavelength = str(int(wavelengths[idx] * 1E9)) + "_nm" + file = folder + "/oa_" + str(iTick) + "_" + wavelength + ".png" + cv2.imwrite(file, (oa_frame+32768)/256) # conversion is to ensure image is written with 16-bit values + print(f"INFO: Wrote file '{file}' for OA frame") + + for idx, us_image in enumerate(pa_data.meta_data_acquisition['ultrasound_image_data']): + iTick = str(int(pa_data.meta_data_acquisition['ultrasound_image_timestamps'][idx]*1E3)) + "_msec" + file = folder + "/us_" + str(iTick) + ".png" + cv2.imwrite(file, us_image) + print(f"INFO: Wrote file '{file}' for US frame") + +file = folder + "/" + os.path.basename(input_file) + "_imagio_ipasc.hdf5" +write_data(file, pa_data) +print(f"INFO: Wrote file '{file}' with IPASC HDF5 data") + diff --git a/examples/python/read_ipasc_hdf5_file.py b/examples/python/read_ipasc_hdf5_file.py index ced1d32..6f128a6 100644 --- a/examples/python/read_ipasc_hdf5_file.py +++ b/examples/python/read_ipasc_hdf5_file.py @@ -7,7 +7,7 @@ import numpy as np name = "kwave_2Dsim_random_array_new" -FILE_PATH = r"..\..\data\ipasc_compatible_V1.hdf5" +FILE_PATH = r"..\..\data\sample_ipasc_kwave_2Dsim_circular_array.hdf5" #FILE_PATH = f"C:/kWave-PACFISH-export/{name}.hdf5" # FILE_PATH = "C:/standardised-image-reconstruction/15PPMPX__ZJQLvYSdWe5CxumvueVrizFy_ipasc.hdf5" @@ -21,7 +21,8 @@ shape = np.shape(pa_data.binary_time_series_data) plt.figure() plt.title("Time Series Data") -plt.imshow(np.squeeze(pa_data.binary_time_series_data[:, :, 0].T)) +plt.imshow(np.squeeze(pa_data.binary_time_series_data[:, :, 0].T), aspect=(pa_data.binary_time_series_data.shape[0]/ + pa_data.binary_time_series_data.shape[1])) plt.xlabel("Num detector elements") plt.ylabel("Num time samples") cb = plt.colorbar() diff --git a/examples/python/read_nrrd_file.py b/examples/python/read_nrrd_file.py index bf91555..32fbd44 100644 --- a/examples/python/read_nrrd_file.py +++ b/examples/python/read_nrrd_file.py @@ -19,11 +19,11 @@ pa_data = converter.generate_pa_data() -quality_check_pa_data(pa_data, verbose=True, log_file_path="") +quality_check_pa_data(pa_data, verbose=True) write_data("demodata_ipasc.hdf5", pa_data) -binary = np.rot90(pa_data.binary_time_series_data[:, 500:-2500, 0, 0], -1) +binary = np.rot90(pa_data.binary_time_series_data[:, 500:-2500, 0], -1) binary = binary - np.min(binary) + 1 binary = np.log10(binary) plt.imshow(binary, aspect=0.08, vmin=np.percentile(binary, 1), vmax=np.percentile(binary, 99)) diff --git a/examples/python/visualise_zenodo_example_data.py b/examples/python/visualise_zenodo_example_data.py index 3f73344..faefc16 100644 --- a/examples/python/visualise_zenodo_example_data.py +++ b/examples/python/visualise_zenodo_example_data.py @@ -19,10 +19,11 @@ #data_url = "https://zenodo.org/record/5938838/files/sample_ipasc_kwave_2Dsim_semicircular_array.hdf5" # Downloading file -filename = download_file(data_url) +#filename = download_file(data_url) # Loading file from disk as a pa_data instance -pa_data = pf.load_data(filename) +# "D:/PACFISH/examples/python/output/2017_12_11-10_54_55_0.lom/2017_12_11-10_54_55_0.lom_imagio_ipasc.hdf5" +pa_data = pf.load_data("D:/PACFISH/examples/python/output/20211102134837219776.dcm.lom/20211102134837219776.dcm.lom_imagio_ipasc.hdf5") # perform completeness and consistency checking on the pa_data pf.quality_check_pa_data(pa_data, verbose=True) @@ -35,7 +36,7 @@ shape = np.shape(pa_data.binary_time_series_data) plt.figure(figsize=(4.75, 4)) plt.title("Time Series Data") -plt.imshow(np.squeeze(pa_data.binary_time_series_data).T, aspect=shape[0]/shape[1]) +plt.imshow(np.squeeze(pa_data.binary_time_series_data[:, :, 0]).T, aspect=shape[0]/shape[1]) plt.xlabel("Num detector elements") plt.ylabel("Num time samples") cb = plt.colorbar() diff --git a/pacfish/api/BaseAdapter.py b/pacfish/api/BaseAdapter.py index 9945c31..588d4fd 100644 --- a/pacfish/api/BaseAdapter.py +++ b/pacfish/api/BaseAdapter.py @@ -47,6 +47,9 @@ def __init__(self): meta_data_device = self.generate_device_meta_data() self.pa_data.meta_data_device = meta_data_device + # We hard-set the version number when exporting. + self.pa_data.meta_data_acquisition[MetadataAcquisitionTags.VERSION.tag] = "V2" + @abstractmethod def generate_binary_data(self) -> np.ndarray: diff --git a/pacfish/api/adapters/Imagio_File_Converter.py b/pacfish/api/adapters/Imagio_File_Converter.py new file mode 100644 index 0000000..8f9f499 --- /dev/null +++ b/pacfish/api/adapters/Imagio_File_Converter.py @@ -0,0 +1,273 @@ +# +# This file extends the base PACFISH device class for the Seno (https://senomedical.com/) +# Imagio System. The class initializer reads in a Laser Optic Movie (LOM), which contains raw +# Photoacoustic ("Optoacoustic" in Seno terminology) and Ultrasound frames. +# + +import numpy as np +import os +import glob +import struct +import cv2 +import math + +from pacfish import BaseAdapter, MetaDatum +from pacfish import MetadataAcquisitionTags +from pacfish import DeviceMetaDataCreator, DetectionElementCreator, IlluminationElementCreator + +class ImagioFileConverter(BaseAdapter): + + # from OAFrameHeader.h in Seno Imagio Sw + OAFRAMETYPE_OA = 1 + OAFRAMETYPE_US = 2 + OAFRAMETYPE_CONFIG = 3 + OAFRAMETYPE_ANNOTATION = 4 + OAFRAMETYPE_AUDIO = 5 + OAFRAMETYPE_VERSION = 6 + OAFRAMETYPE_PROBE = 7 + OAFRAMETYPE_TGC = 8 + OAFRAMETYPE_PROBE_POSITION = 9 + + OAFRAME_MAGIC = 0xbee5 + OAFRAME_HEADER_SIZE = 1024 # bytes + OAFRAME_SUBHEADER_SIZE = 20 # bytes + OAFRAME_CRC_SIZE = 4 # bytes + OAFRAME_META_LASER_INFO_OFFSET = 72 # bytes + OAFRAME_META_SAMPLE_INFO_OFFSET = 90 # bytes + OAFRAME_META_WAVELENGTH_OFFSET = 124 # bytes + OAFRAME_META_LASER_ENERGY_OFFSET = 190 # bytes + OAFRAME_DATATYPE_SHORT = 2 + + # see AcquisitionInterface.h in Seno Imagio SW + OAFRAME_DEFAULT_SAMPLES_PER_CHANNEL = 2048 + + # see ObjectBufferMetaDataDefinitions.h in Seno Imagio SW + USFRAME_WIDTH_OFFSET = 28 # bytes + USFRAME_MICRONSX_OFFSET = 80 # bytes + OAFRAME_WAVELENGTH_ALEXANDRITE = 1 + OAFRAME_WAVELENGTH_YAG = 2 + wavelengths_nm = {OAFRAME_WAVELENGTH_ALEXANDRITE : 750, OAFRAME_WAVELENGTH_YAG : 1064} # nanometers + pulsewidth_nsec = {OAFRAME_WAVELENGTH_ALEXANDRITE : 90, OAFRAME_WAVELENGTH_YAG : 7} # nanoseconds + + # generated via https://www.uuidgenerator.net/version4 + uuid = "a522bad9-f9a4-43b3-a8c3-80cde9e21d2e" + + OAFRAME_NOMINAL_ENERGY = 85 # mJ + + fov = [] # field of view. populated during parsing, but used used when creating device metadata + + meta = {} + data = [] + + """ + For the Seno Imagio system. + """ + def __init__(self, filename): + + # parse through the entire Laser Optic Movie (.lom) file. see OAFrameHeader.h in the Imagio SW for the binary format. + try: + f = open(filename, "rb") + except OSError: + print(f"ERROR: Could not open {filename = }") + return + with f: + self.meta[MetadataAcquisitionTags.ACQUISITION_WAVELENGTHS] = [] + self.meta[MetadataAcquisitionTags.PULSE_ENERGY] = [] + self.meta[MetadataAcquisitionTags.MEASUREMENT_TIMESTAMPS] = [] + self.meta[MetadataAcquisitionTags.ULTRASOUND_IMAGE_DATA] = [] + self.meta[MetadataAcquisitionTags.ULTRASOUND_IMAGE_TIMESTAMPS] = [] + self.data = [] + + print(f"INFO: Reading in Seno Imagio Optoacoustic data file '{filename}'") + iSampleRate = iSoundVelocity = iMicronsX = iMicronsY = 0 + fGain = 0.0 + while True: + data = f.read(self.OAFRAME_HEADER_SIZE) + + if not data: + break + + # extracted variables (i.e. "sMagic", "iTick", etc..) line up with those defined in ObjectBufferMetaDataDefinitions.h in the Seno Imagio SW + metaData = (sMagic, sVersion, iTick, lSize, lFrameCounter, sType, sDummy) = struct.unpack(" J + self.meta[MetadataAcquisitionTags.MEASUREMENT_TIMESTAMPS].append(iTick / 1E3) # msec -> sec + if (cWavelength == self.OAFRAME_WAVELENGTH_ALEXANDRITE or cWavelength == self.OAFRAME_WAVELENGTH_YAG): + self.meta[MetadataAcquisitionTags.ACQUISITION_WAVELENGTHS].append(self.wavelengths_nm[cWavelength] * 1E-9) # nanometers -> meters + else: + print(f"WARNING: Wavelength ({cWavelength = }) not as expected for frame. Skipping to next") + continue + + elif (sType == self.OAFRAMETYPE_US): # Ultrasound frame + + (w, h, ss) = struct.unpack(" sec + + if (len(self.fov) == 0): # populate with first set of values we get + self.fov = np.asarray([0, iMicronsX * w * 1E-6, 0, 0, 0, iMicronsY * h * 1E-6]) + + + # required for conversion to HDF5 + for key in self.meta: + self.meta[key] = np.asarray(self.meta[key]) + + # fill out remainder of metadata + self.meta[MetadataAcquisitionTags.ENCODING] = "raw" + self.meta[MetadataAcquisitionTags.COMPRESSION] = "none" + self.meta[MetadataAcquisitionTags.DATA_TYPE] = "unsigned short" + self.meta[MetadataAcquisitionTags.DIMENSIONALITY] = "time" + self.meta[MetadataAcquisitionTags.SIZES] = np.asarray([sNumChans, self.OAFRAME_DEFAULT_SAMPLES_PER_CHANNEL]) + self.meta[MetadataAcquisitionTags.MEASUREMENTS_PER_IMAGE] = self.OAFRAME_DEFAULT_SAMPLES_PER_CHANNEL + self.meta[MetadataAcquisitionTags.PHOTOACOUSTIC_IMAGING_DEVICE_REFERENCE] = self.uuid + self.meta[MetadataAcquisitionTags.UUID] = self.uuid + self.meta[MetadataAcquisitionTags.ACOUSTIC_COUPLING_AGENT] = "gel" + self.meta[MetadataAcquisitionTags.SCANNING_METHOD] = "handheld" + self.meta[MetadataAcquisitionTags.AD_SAMPLING_RATE] = float(iSampleRate) # use last parsed value + self.meta[MetadataAcquisitionTags.SPEED_OF_SOUND] = float(iSoundVelocity) # use last parsed value + self.meta[MetadataAcquisitionTags.OVERALL_GAIN] = fGain + + # these are intentionally not populated but necessary for quality check + self.meta[MetadataAcquisitionTags.ELEMENT_DEPENDENT_GAIN] = np.asarray([]) + self.meta[MetadataAcquisitionTags.FREQUENCY_DOMAIN_FILTER] = np.asarray([-1, -1]) + self.meta[MetadataAcquisitionTags.TIME_GAIN_COMPENSATION] = np.asarray([]) + self.meta[MetadataAcquisitionTags.MEASUREMENT_SPATIAL_POSES] = np.asarray([[0],[0]]) + self.meta[MetadataAcquisitionTags.TEMPERATURE_CONTROL] = np.asarray([]) + self.meta[MetadataAcquisitionTags.REGIONS_OF_INTEREST] = {} + + self.data = np.asarray(self.data) + + super().__init__() + + + def generate_binary_data(self) -> np.ndarray: + """ + The binary data is the raw time series data. + It is internally stored as an N-dimensional numpy array. + The binary data must be formatted the following way: + + [detectors, samples, wavelengths, frames] + e.g (128, 2048) + + Return + ------ + np.ndarray + A numpy array containing the binary data + """ + # self.data is formatted this way: [wavelengths, detectors, samples] + self.data = np.moveaxis(self.data, 0, 2) # swap to [detectors, samples, wavelengths] + if len(np.shape(self.data)) == 3: # in case the frames are missing + data_shape = np.shape(self.data) + self.data = np.reshape(self.data, (data_shape[0], data_shape[1], data_shape[2], 1)) + + return self.data + + def generate_device_meta_data(self) -> dict: + """ + Must be implemented to define a digital twin of the photoacoustic imaging device. + This method can be implemented using the DeviceMetaDataCreator. + + Return + ------ + dict + A dictionary containing all key-value pair necessary to describe a digital twin of a + photoacoustic device. + """ + + device_creator = DeviceMetaDataCreator() + device_creator.set_general_information(uuid=self.uuid, fov=self.fov) + num_channels = self.pa_data.get_sizes()[0] + for element_idx in range(num_channels): + detection_element_creator = DetectionElementCreator() + detection_element_creator.set_detector_position(np.asarray([(5.12 / 100 / num_channels) * element_idx, 0, 0])) # 5.12 cm probe + detection_element_creator.set_detector_geometry_type("CUBOID") + detection_element_creator.set_detector_orientation(np.asarray([0, 0, 1])) + + # intentionally not populated but required for quality check + detection_element_creator.set_detector_geometry(np.asarray([0.0000, 0.0000, 0.0000])) # N/A + detection_element_creator.set_frequency_response(np.asarray([[0], [1]])) # N/A + detection_element_creator.set_angular_response(np.asarray([[0], [1]])) # N/A + + device_creator.add_detection_element(detection_element_creator.get_dictionary()) + + for wavelength in self.wavelengths_nm.keys(): + illumination_element_creator = IlluminationElementCreator() + illumination_element_creator.set_wavelength_range(np.asarray([self.wavelengths_nm[wavelength] * 1E-9, self.wavelengths_nm[wavelength]* 1E-9, 1])) + illumination_element_creator.set_pulse_width(float(self.pulsewidth_nsec[wavelength])) + illumination_element_creator.set_illuminator_geometry_type("CUBOID") + illumination_element_creator.set_illuminator_orientation(np.asarray([0, 0, 1])) + illumination_element_creator.set_beam_energy_profile(np.asarray([ + [self.wavelengths_nm[self.OAFRAME_WAVELENGTH_ALEXANDRITE], self.OAFRAME_NOMINAL_ENERGY * 1e-3], + [self.wavelengths_nm[self.OAFRAME_WAVELENGTH_YAG], self.OAFRAME_NOMINAL_ENERGY * 1e-3]])) + illumination_element_creator.set_beam_divergence_angles(33.0 * math.pi / 180) # radians + illumination_element_creator.set_illuminator_geometry(np.asarray([57.25 * 1E-3, 0, 28.63 * 1e-3])) # position of light bars, see SPEC-4702100100 + illumination_element_creator.set_illuminator_position(np.asarray([28.63 * 1e-3, 0, 0])) # see SPEC-4702100100 + + # intentionally not populated but required for quality check + illumination_element_creator.set_beam_stability_profile(np.asarray([[0], [1]])) # N/A + illumination_element_creator.set_beam_intensity_profile(np.asarray([[0], [1]])) # N/A + + device_creator.add_illumination_element(illumination_element_creator.get_dictionary()) + + return device_creator.finalize_device_meta_data() + + def set_metadata_value(self, metadatum: MetaDatum) -> object: + """ + + This method must be implemented to yield appropriate data for all MetaDatum elements in the + MetadataTags class. + + You are given a certain meta datum nd have to return the appropriate information. + + Parameters + ---------- + metadatum: MetaDatum + The MetaDatum for which to return the correct data. + + Return + ------ + object + The data corresponding to the given MetaDatum + """ + if metadatum in self.meta: + return self.meta[metadatum] + else: + return None + + diff --git a/pacfish/api/adapters/Nrrd_File_Converter.py b/pacfish/api/adapters/Nrrd_File_Converter.py index ceda7aa..1f7f712 100644 --- a/pacfish/api/adapters/Nrrd_File_Converter.py +++ b/pacfish/api/adapters/Nrrd_File_Converter.py @@ -28,7 +28,7 @@ def __init__(self, nrrd_file_path): super().__init__() def generate_binary_data(self) -> np.ndarray: - data = np.reshape(self.data, (self.meta['sizes'][0], self.meta['sizes'][1], 1, self.meta['sizes'][2])) + data = np.reshape(self.data, (self.meta['sizes'][0], self.meta['sizes'][1], self.meta['sizes'][2])) return data def generate_device_meta_data(self) -> dict: @@ -69,6 +69,8 @@ def generate_device_meta_data(self) -> dict: np.ones(100)])) illumination_element_creator.set_beam_stability_profile(np.asarray([np.linspace(700, 900, 100), np.ones(100)])) + illumination_element_creator.set_beam_intensity_profile(np.asarray([np.linspace(0, 1, 100), + np.ones(100)])) illumination_element_creator.set_pulse_width(7e-9) device_creator.add_illumination_element(illumination_element_creator.get_dictionary()) diff --git a/pacfish/api/adapters/__init__.py b/pacfish/api/adapters/__init__.py index fd8de29..4ad1da3 100644 --- a/pacfish/api/adapters/__init__.py +++ b/pacfish/api/adapters/__init__.py @@ -6,3 +6,4 @@ """ from .Nrrd_File_Converter import NrrdFileConverter +from .Imagio_File_Converter import ImagioFileConverter diff --git a/pacfish/core/Metadata.py b/pacfish/core/Metadata.py index 11c0604..253ad3e 100644 --- a/pacfish/core/Metadata.py +++ b/pacfish/core/Metadata.py @@ -353,6 +353,7 @@ class MetadataAcquisitionTags: UUID = UnconstrainedMetaDatum("uuid", True, str) ENCODING = UnconstrainedMetaDatum("encoding", True, str) COMPRESSION = UnconstrainedMetaDatum("compression", True, str) + VERSION = UnconstrainedMetaDatum("version", True, str) DATA_TYPE = UnconstrainedMetaDatum("data_type", True, str) DIMENSIONALITY = EnumeratedString("dimensionality", True, str, permissible_strings=DIMENSIONALITY_STRINGS) @@ -380,12 +381,14 @@ class MetadataAcquisitionTags: AD_SAMPLING_RATE = NonNegativeNumber("ad_sampling_rate", True, float, Units.HERTZ) FREQUENCY_DOMAIN_FILTER = UnconstrainedMetaDatum("frequency_domain_filter", False, np.ndarray) MEASUREMENTS_PER_IMAGE = NonNegativeWholeNumber("measurements_per_image", False, int) + ULTRASOUND_IMAGE_DATA = UnconstrainedMetaDatum("ultrasound_image_data", False, np.ndarray) + ULTRASOUND_IMAGE_TIMESTAMPS = UnconstrainedMetaDatum("ultrasound_image_timestamps", False, np.ndarray) - TAGS_BINARY = [DATA_TYPE, DIMENSIONALITY, SIZES] + TAGS_BINARY = [DATA_TYPE, DIMENSIONALITY, SIZES, VERSION] TAGS_CONTAINER = [UUID, ENCODING, COMPRESSION] TAGS_ACQUISITION = [PHOTOACOUSTIC_IMAGING_DEVICE_REFERENCE, PULSE_ENERGY, ACQUISITION_WAVELENGTHS, TIME_GAIN_COMPENSATION, OVERALL_GAIN, ELEMENT_DEPENDENT_GAIN, TEMPERATURE_CONTROL, ACOUSTIC_COUPLING_AGENT, SCANNING_METHOD, AD_SAMPLING_RATE, FREQUENCY_DOMAIN_FILTER, SPEED_OF_SOUND, MEASUREMENTS_PER_IMAGE, REGIONS_OF_INTEREST, MEASUREMENT_TIMESTAMPS, - MEASUREMENT_SPATIAL_POSES] + MEASUREMENT_SPATIAL_POSES, ULTRASOUND_IMAGE_DATA, ULTRASOUND_IMAGE_TIMESTAMPS] TAGS = TAGS_BINARY + TAGS_ACQUISITION + TAGS_CONTAINER diff --git a/pacfish/iohandler/file_reader.py b/pacfish/iohandler/file_reader.py index 52f60e9..22eb7c6 100644 --- a/pacfish/iohandler/file_reader.py +++ b/pacfish/iohandler/file_reader.py @@ -5,6 +5,7 @@ import h5py from pacfish import PAData +from pacfish.core.Metadata import MetadataAcquisitionTags as Tags import numpy as np @@ -59,4 +60,35 @@ def recursively_load_dictionaries(file, path): pa_data = PAData(binary_data) pa_data.meta_data_acquisition = recursively_load_dictionaries(h5file, "/meta_data/") pa_data.meta_data_device = recursively_load_dictionaries(h5file, "/meta_data_device/") + + # Check for the file version to ensure backwards compatibility of the data format. + if ((Tags.VERSION.tag not in pa_data.meta_data_acquisition) or + (pa_data.meta_data_acquisition[Tags.VERSION.tag] == "V1")): + print("INFO: Importing IPASC V1 data file.") + + if len(np.shape(pa_data.binary_time_series_data)) < 4: + # Even though the data is labelled as being V1, it did only contain one (wavelength/measurement) field. + # In this case, we can simply take the data as is and save it into V2. + if len(np.shape(pa_data.binary_time_series_data)) == 3: + pa_data.binary_time_series_data = pa_data.binary_time_series_data[:, :, :, np.newaxis] + elif len(np.shape(pa_data.binary_time_series_data)) == 2: + pa_data.binary_time_series_data = pa_data.binary_time_series_data[:, :, np.newaxis, np.newaxis] + else: + raise AssertionError("The binary time series data was not compatible with the IPASC standard.") + + # we're in version 1, so we ned to reshape the wavelengths and measurements into the same dimension and + # expand the wavelengths field to be as long as the number of measurements. + (n_detectors, n_timesteps, n_wavelengths, n_measurements) = np.shape(pa_data.binary_time_series_data) + pa_data.binary_time_series_data = pa_data.binary_time_series_data.reshape((n_detectors, n_timesteps, + n_wavelengths * n_measurements)) + + if len(pa_data.get_acquisition_wavelengths().shape) == 0: + wavelengths = [pa_data.get_acquisition_wavelengths().item()] + else: + wavelengths = list(pa_data.get_acquisition_wavelengths()) + + wavelengths = wavelengths * n_measurements + + pa_data.meta_data_acquisition[Tags.ACQUISITION_WAVELENGTHS.tag] = np.asarray(wavelengths).reshape((-1, )) + return pa_data diff --git a/pacfish/iohandler/file_writer.py b/pacfish/iohandler/file_writer.py index 51a2f09..da61b33 100644 --- a/pacfish/iohandler/file_writer.py +++ b/pacfish/iohandler/file_writer.py @@ -6,6 +6,7 @@ import h5py from pacfish import PAData import numpy as np +from pacfish.core.Metadata import MetadataAcquisitionTags as Tags def write_data(file_path: str, pa_data: PAData, file_compression: str = None): @@ -28,6 +29,10 @@ def write_data(file_path: str, pa_data: PAData, file_compression: str = None): This method does not return anything """ + if ((Tags.VERSION.tag not in pa_data.meta_data_acquisition) or + (pa_data.meta_data_acquisition[Tags.VERSION.tag] == "V1")): + pa_data.meta_data_acquisition[Tags.VERSION.tag] = "V2" + def recursively_save_dictionaries(file, path, data_dictionary, compression: str = None): for key, item in data_dictionary.items(): key = str(key) diff --git a/pacfish/qualitycontrol/ConsistencyChecker.py b/pacfish/qualitycontrol/ConsistencyChecker.py index 485806f..2d4427c 100644 --- a/pacfish/qualitycontrol/ConsistencyChecker.py +++ b/pacfish/qualitycontrol/ConsistencyChecker.py @@ -41,12 +41,33 @@ def check_binary_data(self, binary_data) -> bool: The binary data to check """ is_consistent = True + + log_message = "" + log_message += "#Consistency report for binary data\n\n" + if not isinstance(binary_data, np.ndarray): + log_message += "binary data was not a numpy array\n" + is_consistent = False + + if len(np.shape(binary_data)) != 3: + log_message += f"binary data shape was not as expected (expected 3 dimensions but was: {len(np.shape(binary_data))})\n" is_consistent = False for number in np.reshape(binary_data, (-1, )): if not isinstance(number, numbers.Number): is_consistent = False + log_message += "binary data contained elements not classifiable as a number\n" + break + + log_message += "\n" + + if self.verbose: + print(log_message) + + if self.log_file_path is not None: + with open(self.log_file_path + self.save_file_name, "a") as log_file_handle: + log_file_handle.writelines(log_message) + return is_consistent def check_acquisition_meta_data(self, acquisition_meta_data: dict) -> bool: diff --git a/pacfish/qualitycontrol/PADataIntegrityCheck.py b/pacfish/qualitycontrol/PADataIntegrityCheck.py index 2e9d5b1..d470a27 100644 --- a/pacfish/qualitycontrol/PADataIntegrityCheck.py +++ b/pacfish/qualitycontrol/PADataIntegrityCheck.py @@ -31,10 +31,8 @@ def quality_check_pa_data(pa_data: PAData, verbose: bool = False, log_file_path: b1 = completeness.check_acquisition_meta_data(pa_data.meta_data_acquisition) b2 = consistency.check_acquisition_meta_data(pa_data.meta_data_acquisition) - b3 = completeness.check_device_meta_data(pa_data.meta_data_device) b4 = consistency.check_device_meta_data(pa_data.meta_data_device) - b5 = consistency.check_binary_data(pa_data.binary_time_series_data) return b1 and b2 and b3 and b4 and b5 \ No newline at end of file diff --git a/pacfish_matlab/@kwave_adapter/kwave_adapter.m b/pacfish_matlab/@kwave_adapter/kwave_adapter.m index ef9144f..3a4641a 100644 --- a/pacfish_matlab/@kwave_adapter/kwave_adapter.m +++ b/pacfish_matlab/@kwave_adapter/kwave_adapter.m @@ -36,11 +36,11 @@ obj.kgrid = varargin{4}; % Changed order of model and fov so it is easier to not define % fov when calling the function - if (nargin==5) % Checks whether the number of inputs from function call is 5 or not - obj.model = varargin{5}; % Used to determine what type of simulation is being done + if (nargin==5) % Checks whether the number of inputs from function call is 6 or above + obj.fov = varargin{5}; % Takes the fov definition end - if (nargin>=6) % Checks whether the number of inputs from function call is 6 or above - obj.fov = varargin{6}; % Takes the fov definition + if (nargin>=6) % Checks whether the number of inputs from function call is 5 or not + obj.model = varargin{6}; % Used to determine what type of simulation is being done end obj.lineheight = 0.1e-3; % We assume the height of line elements is this end diff --git a/requirements.txt b/requirements.txt index 2fe27cd..c48e160 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,4 +6,5 @@ myst-parser coverage pytest-cov imageio -scipy \ No newline at end of file +scipy +opencv-python \ No newline at end of file diff --git a/testing/unit_tests/test_backwards_compatibility.py b/testing/unit_tests/test_backwards_compatibility.py new file mode 100644 index 0000000..b80ff70 --- /dev/null +++ b/testing/unit_tests/test_backwards_compatibility.py @@ -0,0 +1,42 @@ +# SPDX-FileCopyrightText: 2021 International Photoacoustics Standardisation Consortium (IPASC) +# SPDX-License-Identifier: BSD 3-Clause License + +import os +import inspect +from unittest.case import TestCase +import pacfish as pf +from testing.unit_tests.utils import assert_equal_dicts + + +class TestBackwardsCompatibility(TestCase): + + def setUp(self): + print("setUp") + + def tearDown(self): + print("tearDown") + + def test_read_V1_file(self): + path_sep = "../" + file_path = "data/ipasc_compatible_V1.hdf5" + + path = os.path.abspath(os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) + "/../../" + file_path) + pa_data = pf.load_data(path) + + self.assertTrue(len(pa_data.binary_time_series_data.shape) == 3) + + try: + pf.write_data("compatible_ipasc_test.hdf5", pa_data) + test_data = pf.load_data("compatible_ipasc_test.hdf5") + except Exception as e: + raise e + finally: + # clean up after ipasc_test + if os.path.exists("compatible_ipasc_test.hdf5"): + os.remove("compatible_ipasc_test.hdf5") + + pa_data.meta_data_acquisition[pf.MetadataAcquisitionTags.VERSION.tag] = "V2" + + assert_equal_dicts(pa_data.meta_data_acquisition, test_data.meta_data_acquisition) + assert_equal_dicts(pa_data.meta_data_device, test_data.meta_data_device) + self.assertTrue((pa_data.binary_time_series_data == test_data.binary_time_series_data).all()) diff --git a/testing/unit_tests/test_completeness_and_consistency_checker.py b/testing/unit_tests/test_completeness_and_consistency_checker.py index 19948c9..8e10462 100644 --- a/testing/unit_tests/test_completeness_and_consistency_checker.py +++ b/testing/unit_tests/test_completeness_and_consistency_checker.py @@ -132,7 +132,7 @@ def test_pa_data_check(self): device_dict = create_complete_device_metadata_dictionary() acquisition_dict = create_complete_acquisition_meta_data_dictionary() - pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048]), + pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048, 1]), meta_data_acquisition=acquisition_dict, meta_data_device=device_dict) @@ -142,7 +142,7 @@ def test_check_a_complete_and_consistent_pa_data_instance(self): device_dict = create_complete_device_metadata_dictionary() acquisition_dict = create_complete_acquisition_meta_data_dictionary() - pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048]), + pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048, 1]), meta_data_acquisition=acquisition_dict, meta_data_device=device_dict) @@ -164,7 +164,7 @@ def test_check_a_complete_but_inconsistent_pa_data_instance(self): acquisition_dict = create_complete_acquisition_meta_data_dictionary() acquisition_dict[pf.MetadataAcquisitionTags.DIMENSIONALITY.tag] = "Wrong string" - pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048]), + pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048, 1]), meta_data_acquisition=acquisition_dict, meta_data_device=device_dict) diff --git a/testing/unit_tests/test_io_handling.py b/testing/unit_tests/test_io_handling.py index af96e86..e5be243 100644 --- a/testing/unit_tests/test_io_handling.py +++ b/testing/unit_tests/test_io_handling.py @@ -49,7 +49,7 @@ def test_write_and_read_random_dictionary_with_None_values(self): acquisition_dict[pf.MetadataAcquisitionTags.MEASUREMENT_TIMESTAMPS.tag] = None acquisition_dict[pf.MetadataAcquisitionTags.SCANNING_METHOD.tag] = None - pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048]), + pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048, 1]), meta_data_acquisition=acquisition_dict, meta_data_device=device_dict) @@ -72,7 +72,7 @@ def test_overwrite_hdf5_file(self): device_dict = create_complete_device_metadata_dictionary() acquisition_dict = create_complete_acquisition_meta_data_dictionary() - pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048]), + pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048, 1]), meta_data_acquisition=acquisition_dict, meta_data_device=device_dict) @@ -103,7 +103,7 @@ def test_rewrite_and_read_still_consistent(self): device_dict = create_complete_device_metadata_dictionary() acquisition_dict = create_complete_acquisition_meta_data_dictionary() - pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048]), + pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048, 1]), meta_data_acquisition=acquisition_dict, meta_data_device=device_dict) @@ -130,7 +130,7 @@ def test_double_rewrite_consistency(self): device_dict = create_complete_device_metadata_dictionary() acquisition_dict = create_complete_acquisition_meta_data_dictionary() - pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048]), + pa_data = pf.PAData(binary_time_series_data=np.zeros([256, 2048, 1]), meta_data_acquisition=acquisition_dict, meta_data_device=device_dict) diff --git a/testing/unit_tests/utils.py b/testing/unit_tests/utils.py index f5c1b5a..a181554 100644 --- a/testing/unit_tests/utils.py +++ b/testing/unit_tests/utils.py @@ -8,6 +8,7 @@ def create_complete_acquisition_meta_data_dictionary(): dictionary = dict() + dictionary[MetadataAcquisitionTags.VERSION.tag] = "V2" dictionary[MetadataAcquisitionTags.UUID.tag] = create_random_testing_parameters()['test_string'] dictionary[MetadataAcquisitionTags.ENCODING.tag] = create_random_testing_parameters()['test_string'] dictionary[MetadataAcquisitionTags.COMPRESSION.tag] = create_random_testing_parameters()['test_string'] @@ -30,6 +31,8 @@ def create_complete_acquisition_meta_data_dictionary(): dictionary[MetadataAcquisitionTags.SPEED_OF_SOUND.tag] = np.asarray(1540.0) dictionary[MetadataAcquisitionTags.MEASUREMENT_SPATIAL_POSES.tag] = create_random_testing_parameters()['test_array'] dictionary[MetadataAcquisitionTags.MEASUREMENTS_PER_IMAGE.tag] = 1 + dictionary[MetadataAcquisitionTags.ULTRASOUND_IMAGE_DATA.tag] = np.asarray([2, 2]) + dictionary[MetadataAcquisitionTags.ULTRASOUND_IMAGE_TIMESTAMPS.tag] = np.asarray([2, 2]) return dictionary