diff --git a/exact/util/cellvizio.py b/exact/util/cellvizio.py index 92a424eb..2e4e9662 100644 --- a/exact/util/cellvizio.py +++ b/exact/util/cellvizio.py @@ -1,204 +1,209 @@ """ - This is SlideRunner - An Open Source Annotation Tool + This is SlideRunner - An Open Source Annotation Tool for Digital Histology Slides. - Marc Aubreville, Pattern Recognition Lab, - Friedrich-Alexander University Erlangen-Nuremberg + Marc Aubreville, Pattern Recognition Lab, + Friedrich-Alexander University Erlangen-Nuremberg marc.aubreville@fau.de If you use this software in research, please citer our paper: M. Aubreville, C. Bertram, R. Klopfleisch and A. Maier: - SlideRunner - A Tool for Massive Cell Annotations in Whole Slide Images. - In: Bildverarbeitung für die Medizin 2018. + SlideRunner - A Tool for Massive Cell Annotations in Whole Slide Images. + In: Bildverarbeitung für die Medizin 2018. Springer Vieweg, Berlin, Heidelberg, 2018. pp. 309-314. - This file: Support for CellVizio MKT files. + This file: Support for CellVizio MKT files. """ import numpy as np -from pydicom.encaps import decode_data_sequence from PIL import Image -import io -import os import struct -from os import stat import openslide from util.enums import FrameType + +_CHUNK_HEADER_SIZE = 16 +_CSTA_MAGIC = b'\n tuple[list[tuple[int, int, int]], dict[str, str]]: + """Parse a type-5 MKT chunk body. + + Returns (index_entries, kv) where index_entries is a list of + (chunk_type, file_offset, data_size) and kv is all key=value metadata. + """ + n = struct.unpack(' tuple[list[tuple[int, int, int]], dict[str, str]]: + """Locate and parse the type-5 index/metadata chunk in an MKT file. + + Scans the last 64 KB (sufficient for all trailing metadata chunks) to find + the type-5 chunk, which always appears near the end of the file. + Returns (index_entries, kv_dict). """ + with open(filename, 'rb') as f: + f.seek(0, 2) + file_size = f.tell() + scan_size = min(65536, file_size) + f.seek(file_size - scan_size) + tail = f.read() + + pos = 0 + while pos < len(tail): + idx = tail.find(_CSTA_MAGIC, pos) + if idx == -1: + break + hdr = tail[idx : idx + _CHUNK_HEADER_SIZE] + if len(hdr) < _CHUNK_HEADER_SIZE: + break + chunk_type = struct.unpack('>H', hdr[8:10])[0] + chunk_size = struct.unpack('>I', hdr[10:14])[0] + if chunk_type == 5: + body = tail[idx + _CHUNK_HEADER_SIZE : idx + _CHUNK_HEADER_SIZE + chunk_size] + return _parse_type5_body(body) + pos = idx + 1 + + return [], {} - def __init__(self,filename): - #print('Opening:',filename) - self.fileName = filename; + +class ReadableCellVizioMKTDataset(): + + def __init__(self, filename): + self.fileName = filename self.fi = fileinfo() - self.fileHandle = open(filename, 'rb'); - self.fileHandle.seek(5) # we find the FPS at position 05 - self.fileHandle.seek(10) # we find the image size at position 10 - fSizeByte = self.fileHandle.read(4) - self.fi.size = int.from_bytes(fSizeByte, byteorder='big', signed=True) - self.filestats = stat(self.fileName) - self.fi.nImages=1000 - self.fi.nImages = int((self.filestats.st_size-self.fi.offset) / (self.fi.size+self.fi.gapBetweenImages)) - self.numberOfFrames = self.fi.nImages - - meta = self.getMostRelevantMetaInfo() - self.fi.width = int(meta.get('width', 0)) + + # Read image data size from the first chunk header (bytes 10–13, big-endian uint32) + with open(filename, 'rb') as f: + f.seek(10) + self.fi.size = struct.unpack('>I', f.read(4))[0] + + # Parse the type-5 chunk for the chunk index and all key=value metadata + chunk_index, self._kv = _load_chunk_index(filename) + + # Collect data offsets for every type-2 (image frame) chunk + self._frame_offsets = [ + entry[1] + _CHUNK_HEADER_SIZE + for entry in chunk_index if entry[0] == 2 + ] + self.numberOfFrames = len(self._frame_offsets) + self.fi.nImages = self.numberOfFrames + + meta = self._kv + self.fi.width = int(meta.get('width', 0)) self.fi.height = int(meta.get('height', 0)) self.fps = float(meta.get('framerate', 0)) - - self.geometry_imsize = [self.fi.height, self.fi.width] - self.imsize = [self.fi.height, self.fi.width] + + self.geometry_imsize = [self.fi.height, self.fi.width] + self.imsize = [self.fi.height, self.fi.width] self.geometry_tilesize = [(self.fi.height, self.fi.width)] - self.geometry_rows = [1] - self.geometry_columns = [1] - self.levels = [1] - self.channels = 1 - - self.fovx = float(meta.get('fovx', 250)) - self.fovy = float(meta.get('fovy', 250)) + self.geometry_rows = [1] + self.geometry_columns = [1] + self.levels = [1] + self.channels = 1 + + self.fovx = float(meta.get('fovx', 250)) + self.fovy = float(meta.get('fovy', 250)) print(f"fovx: {self.fovx}, fovy: {self.fovy}") - self.mpp_x = self.fovx/self.fi.width # approximate number for gastroflex, 250 ym field of view, 576 px - self.mpp_y = self.fovy/self.fi.height - - - # generate circular mask for this file - self.circMask = circularMask(self.fi.width,self.fi.height, self.fi.width-2).mask - #print('Circular mask shape:',self.circMask.shape) - self.properties = { openslide.PROPERTY_NAME_BACKGROUND_COLOR: '000000', - openslide.PROPERTY_NAME_MPP_X: self.mpp_x, - openslide.PROPERTY_NAME_MPP_Y: self.mpp_y, - openslide.PROPERTY_NAME_OBJECTIVE_POWER:20, - openslide.PROPERTY_NAME_VENDOR: 'MKT'} - - - def getMetaInfo(self): - # the meta information at the MKT file always starts with "00616C6C 6F776564 5F656761 696E5F65 6F666673 65745F70 61697273 3D" or allowed_egain_eoffset_pairs= in the hex code - # we need to read the whole file to find the meta information - with open(self.fileName, 'rb') as f: - fileContent = f.read() - # find start position with fileContent.find(b'allowed_egain_eoffset_pairs=') that searches for the byte sequence - metaStart = fileContent.rfind(b'allowed_egain_eoffset_pairs=') - metaEnd = fileContent.find(b'f', fFPSByte)[0] - - self.fi = fileinfo() - self.fileHandle.seek(10) # we find the image size at position 10 - fSizeByte = self.fileHandle.read(4) - self.fi.size = int.from_bytes(fSizeByte, byteorder='big', signed=True) - self.fi.nImages=1000 - - self.fi.width = 576 - if ((self.fi.size / (2 * self.fi.width)) % 2 != 0): - self.fi.width=512 - self.fi.height=int(self.fi.size/(2*self.fi.width)) - else: - self.fi.height=int(self.fi.size/(2*self.fi.width)) - - self.filestats = stat(self.fileName) - self.fi.nImages = int((self.filestats.st_size-self.fi.offset) / (self.fi.size+self.fi.gapBetweenImages)) - - - self.numberOfFrames = self.fi.nImages - - self.geometry_imsize = [self.fi.height, self.fi.width] - self.imsize = [self.fi.height, self.fi.width] - self.geometry_tilesize = [(self.fi.height, self.fi.width)] - self.geometry_rows = [1] - self.geometry_columns = [1] - self.levels = [1] - self.channels = 1 - self.mpp_x = 250/576 # approximate number for gastroflex, 250 ym field of view, 576 px - self.mpp_y = 250/576 - - self.properties = { openslide.PROPERTY_NAME_BACKGROUND_COLOR: '000000', - openslide.PROPERTY_NAME_MPP_X: self.mpp_x, - openslide.PROPERTY_NAME_MPP_Y: self.mpp_y, - openslide.PROPERTY_NAME_OBJECTIVE_POWER:20, - openslide.PROPERTY_NAME_VENDOR: 'MKT'} - - self.circMask = circularMask(self.fi.width,self.fi.height, self.fi.width-2).mask + self.mpp_x = self.fovx / self.fi.width + self.mpp_y = self.fovy / self.fi.height + + self.circMask = circularMask(self.fi.width, self.fi.height, self.fi.width - 2).mask + + self.properties = { + openslide.PROPERTY_NAME_BACKGROUND_COLOR: '000000', + openslide.PROPERTY_NAME_MPP_X: self.mpp_x, + openslide.PROPERTY_NAME_MPP_Y: self.mpp_y, + openslide.PROPERTY_NAME_OBJECTIVE_POWER: 20, + openslide.PROPERTY_NAME_VENDOR: 'MKT', + } + + def getMetaInfo(self) -> dict[str, str]: + """Return all key=value metadata from the type-5 chunk.""" + return self._kv + + def getMostRelevantMetaInfo(self) -> dict[str, str]: + """Return metadata values for all labeled fields (keys defined in meta_data_dict).""" + result = {k: self._kv[k] for k in self.meta_data_dict if k in self._kv} + if self.fps > 0: + result['duration_seconds'] = str(self.numberOfFrames / self.fps) + return result @property def meta_data(self) -> dict: - # Cache metadata to avoid repeated file I/O and parsing. - if not hasattr(self, "_meta_data_cache"): - self._meta_data_cache = self.getMostRelevantMetaInfo() - return self._meta_data_cache - + return self.getMostRelevantMetaInfo() + @property def meta_data_dict(self) -> dict: - meta_data_dict = { - 'width': 'Image width (pixels)', - 'height': 'Image height (pixels)', - 'framerate': 'Frame rate (fps)', - 'duration_seconds': 'Duration (seconds)', - 'patient_id': 'Patient ID', - 'fovx': 'Field of view x (micrometers)', - 'fovy': 'Field of view y (micrometers)' - } - return meta_data_dict + labels = { + # Patient / recording + 'patient_id': 'Patient ID', + 'biopsy': 'Biopsy taken', + 'duration_seconds': 'Duration (s)', + 'framerate': 'Frame rate (fps)', + 'mosaicing_enabled': 'Mosaicing enabled', + 'utc_offset_seconds': 'UTC offset (s)', + # Image geometry + 'width': 'Image width (px)', + 'height': 'Image height (px)', + 'fovx': 'FOV x (µm)', + 'fovy': 'FOV y (µm)', + 'pfovx': 'Probe FOV x (µm)', + 'pfovy': 'Probe FOV y (µm)', + 'bbox_min_x': 'Mosaic bbox min x (px)', + 'bbox_min_y': 'Mosaic bbox min y (px)', + 'bbox_max_x': 'Mosaic bbox max x (px)', + 'bbox_max_y': 'Mosaic bbox max y (px)', + # Display + 'pal_cropMin': 'Display window min', + 'pal_cropMax': 'Display window max', + 'mask_level': 'Mask level', + 'compression_type': 'Compression', + # Laser / acquisition + 'eocu_mode': 'Imaging mode', + 'eocu_laser_wavelength': 'Laser wavelength (nm)', + 'eocu_lpwr': 'Laser power (%)', + 'eocu_serial_number': 'EOCU serial number', + 'eocu_uid': 'EOCU UID', + # Probe + 'proflex_type': 'Probe type', + 'proflex_uid': 'Probe UID', + 'proflex_diameter': 'Probe diameter (mm)', + 'proflex_length': 'Probe length (m)', + 'proflex_lateral_res': 'Lateral resolution (µm)', + 'proflex_axial_res': 'Axial resolution (µm)', + 'proflex_working_dist': 'Working distance (µm)', + 'proflex_sensitivity': 'Probe sensitivity', + # Software + 'mzversion_str': 'Software version', + } + # Only return labels for fields present in this file + return {k: v for k, v in labels.items() if k in self._kv or k == 'duration_seconds'} @property def seriesInstanceUID(self) -> str: @@ -208,132 +213,109 @@ def seriesInstanceUID(self) -> str: def level_downsamples(self): return [1] - @property + @property def level_dimensions(self): return [self.geometry_imsize] - @property + @property def nFrames(self): return self.numberOfFrames - - @property - def frame_descriptors(self) -> list[str]: - """ returns a list of strings, used as descriptor for each frame - """ - return ['%.2f s (%d)' % (float(frame_id)/float(self.fps), frame_id) for frame_id in range(self.nFrames)] @property def frame_descriptors(self) -> list[str]: - return 0 - + return ['%.2f s (%d)' % (float(frame_id) / float(self.fps), frame_id) + for frame_id in range(self.nFrames)] - @property + @property def nLayers(self): return 1 @property def layer_descriptors(self) -> list[str]: - """ returns a list of strings, used as descriptor for each layer - """ return [''] @property def frame_type(self): return FrameType.TIMESERIES - def get_thumbnail(self, size): - return self.read_region((0,0),0, self.dimensions).resize(size) + return self.read_region((0, 0), 0, self.dimensions).resize(size) def readImage(self, position=0): + seekpos = self._frame_offsets[position] + image = np.fromfile(self.fileName, offset=seekpos, dtype=np.int16, + count=int(self.fi.size / 2)) - seekpos=self.fi.offset + self.fi.size*position + self.fi.gapBetweenImages*position - image = np.fromfile(self.fileName, offset=seekpos, dtype=np.int16, count=int(self.fi.size/2)) - - if (image.size > 0): + if image.size > 0: image = np.clip(image, 0, np.max(image)) - if (image.shape[0]!=self.fi.height* self.fi.width): + if image.shape[0] != self.fi.height * self.fi.width: image = np.zeros((self.fi.height, self.fi.width)) - image=np.reshape(image, newshape=(self.fi.height, self.fi.width)) - - return image + return np.reshape(image, (self.fi.height, self.fi.width)) - def scaleImageUINT8(self, image, mask = None): - # read image and scale to uint8 [0;255] format + def scaleImageUINT8(self, image, mask=None): + if mask is None: + mask = self.circMask - if (mask is None): - mask = self.circMask + maskedImage = image[mask] + cmin, cmax = np.percentile(maskedImage, 0.5), np.percentile(maskedImage, 99.5) + if cmax > 5000: + cmax = 5000 + dyn = cmax - cmin - maskedImage = image[mask] + compr = 255 / dyn + image = (image - cmin) * compr + return np.uint8(np.clip(np.round(image), 0, 255)) - cmin,cmax = np.percentile(maskedImage,0.5), np.percentile(maskedImage,99.5) - if (cmax>5000): - cmax=5000 - dyn=cmax-cmin - - # compress - compr=255/dyn - image = image-cmin - image = image*compr - - # limit to 0 - image = np.clip(np.round(image),0,255) - image=np.uint8(image) - - return image def readImageUINT8(self, position=0): - # read image and scale to uint8 [0;255] format - image=self.readImage(position) - - image = self.scaleImageUINT8(image) - return image + return self.scaleImageUINT8(self.readImage(position)) - def read_region(self, location: tuple, level:int, size:tuple, frame:int=0): - img = np.zeros((size[1],size[0],4), np.uint8) - img[:,:,3]=255 - offset=[0,0] - if (location[1]<0): + def read_region(self, location: tuple, level: int, size: tuple, frame: int = 0): + img = np.zeros((size[1], size[0], 4), np.uint8) + img[:, :, 3] = 255 + offset = [0, 0] + if location[1] < 0: offset[0] = -location[1] - location = (location[0],0) - if (location[0]<0): + location = (location[0], 0) + if location[0] < 0: offset[1] = -location[0] - location = (0,location[1]) + location = (0, location[1]) pixel_array = self.readImageUINT8(position=frame) - imgcut = pixel_array[location[1]:location[1]+size[1]-offset[0],location[0]:location[0]+size[0]-offset[1]] - imgcut = np.uint8(np.clip(np.float32(imgcut),0,255)) + imgcut = pixel_array[ + location[1] : location[1] + size[1] - offset[0], + location[0] : location[0] + size[0] - offset[1], + ] + imgcut = np.uint8(np.clip(np.float32(imgcut), 0, 255)) for k in range(3): - img[offset[0]:imgcut.shape[0]+offset[0],offset[1]:offset[1]+imgcut.shape[1],k] = imgcut + img[offset[0] : imgcut.shape[0] + offset[0], + offset[1] : offset[1] + imgcut.shape[1], k] = imgcut return Image.fromarray(img) - @property def dimensions(self): return self.level_dimensions[0] - def get_best_level_for_downsample(self,downsample): - return np.argmin(np.abs(np.asarray(self.level_downsamples)-downsample)) + def get_best_level_for_downsample(self, downsample): + return np.argmin(np.abs(np.asarray(self.level_downsamples) - downsample)) @property def level_count(self): return len(self.levels) - def imagePos_to_id(self, imagePos:tuple, level:int): + def imagePos_to_id(self, imagePos: tuple, level: int): id_x, id_y = imagePos - if (id_y>=self.geometry_rows[level]): - id_x=self.geometry_columns[level] # out of range - - if (id_x>=self.geometry_columns[level]): - id_y=self.geometry_rows[level] # out of range - return (id_x+(id_y*self.geometry_columns[level])) - - - def get_id(self, pixelX:int, pixelY:int, level:int) -> (int, int, int): - - id_x = round(-0.5+(pixelX/self.geometry_tilesize[level][1])) - id_y = round(-0.5+(pixelY/self.geometry_tilesize[level][0])) - - return (id_x,id_y), pixelX-(id_x*self.geometry_tilesize[level][0]), pixelY-(id_y*self.geometry_tilesize[level][1]), \ No newline at end of file + if id_y >= self.geometry_rows[level]: + id_x = self.geometry_columns[level] + if id_x >= self.geometry_columns[level]: + id_y = self.geometry_rows[level] + return id_x + (id_y * self.geometry_columns[level]) + + def get_id(self, pixelX: int, pixelY: int, level: int): + id_x = round(-0.5 + (pixelX / self.geometry_tilesize[level][1])) + id_y = round(-0.5 + (pixelY / self.geometry_tilesize[level][0])) + return ((id_x, id_y), + pixelX - (id_x * self.geometry_tilesize[level][0]), + pixelY - (id_y * self.geometry_tilesize[level][1])) diff --git a/exact/util/slideio.py b/exact/util/slideio.py index 604b2188..fd7d09f8 100644 --- a/exact/util/slideio.py +++ b/exact/util/slideio.py @@ -45,9 +45,19 @@ def __init__(self,filename): # Zeiss is funny and thinks BGR is the proper way to store images ... if (os.path.splitext(filename.upper())[-1] == '.CZI'): self.mode = 'BGRA' - - - #print('Circular mask shape:',self.circMask.shape) + + # SlideIO sometimes returns 0 for resolution on DICOM files that lack + # standard PixelSpacing tags. Fall back to pydicom for those tags. + if self.mpp_x == 0.0: + _dcm = self._read_dcm_meta() + if _dcm is not None: + ps = (getattr(_dcm, 'PixelSpacing', None) + or getattr(_dcm, 'ImagerPixelSpacing', None) + or getattr(_dcm, 'NominalScannedPixelSpacing', None)) + if ps and float(ps[0]) > 0: + self.mpp_x = float(ps[1]) * 1000 # mm → µm + self.mpp_y = float(ps[0]) * 1000 + self.properties = { openslide.PROPERTY_NAME_BACKGROUND_COLOR: '000000', openslide.PROPERTY_NAME_MPP_X: self.mpp_x, openslide.PROPERTY_NAME_MPP_Y: self.mpp_y, @@ -71,11 +81,35 @@ def level_dimensions(self): def nFrames(self): return self.scene.num_z_slices + def _read_dcm_meta(self): + """Read DICOM metadata via pydicom; returns dataset or None.""" + if not hasattr(self, '_dcm_meta_cache'): + ext = os.path.splitext(self.fileName.lower())[1] + if ext == '.dcm': + try: + import pydicom + self._dcm_meta_cache = pydicom.dcmread(self.fileName, stop_before_pixels=True) + except Exception: + self._dcm_meta_cache = None + else: + self._dcm_meta_cache = None + return self._dcm_meta_cache + @property def frame_descriptors(self) -> list[str]: """ returns a list of strings, used as descriptor for each frame """ - return ['%.2f µ' % (self.scene.z_resolution*1E6*x) for x in range(self.scene.num_z_slices)] + z_res = self.scene.z_resolution * 1e6 + if z_res > 0: + return ['%.2f µ' % (z_res * x) for x in range(self.scene.num_z_slices)] + # For DICOM temporal sequences, fall back to FrameTime tag + _dcm = self._read_dcm_meta() + if _dcm is not None: + ft = getattr(_dcm, 'FrameTime', None) + if ft is not None: + frame_time_s = float(ft) / 1000.0 + return ['%.2f s (%d)' % (frame_time_s * x, x) for x in range(self.scene.num_z_slices)] + return ['Frame %d' % x for x in range(self.scene.num_z_slices)] @property def default_frame(self): @@ -93,6 +127,11 @@ def layer_descriptors(self) -> list[str]: @property def frame_type(self): + _dcm = self._read_dcm_meta() + if _dcm is not None: + modality = str(getattr(_dcm, 'Modality', '')) + if modality in ('CT', 'MR', 'PT', 'NM'): + return FrameType.ZSTACK return FrameType.TIMESERIES @@ -105,7 +144,15 @@ def read_region(self, location: tuple, level:int, size:tuple, frame:int=0): img = self.scene.read_block(rect=[*location, int(size[0]/ds), int(size[1]/ds)], size=size, slices=(frame, frame+1)) img_4ch = np.zeros([size[1], size[0],4], dtype=np.uint8) img_4ch[:,:,3] = 255 - img_4ch[:,:,0:3] = img if self.mode=='RGBA' else img[:,:,::-1] + if img.ndim == 2: + # Grayscale (single-channel) — broadcast to RGB + img_4ch[:,:,0] = img + img_4ch[:,:,1] = img + img_4ch[:,:,2] = img + elif self.mode == 'RGBA': + img_4ch[:,:,0:3] = img + else: + img_4ch[:,:,0:3] = img[:,:,::-1] return Image.fromarray(img_4ch) @property @@ -129,7 +176,7 @@ def imagePos_to_id(self, imagePos:tuple, level:int): return (id_x+(id_y*self.geometry_columns[level])) - def get_id(self, pixelX:int, pixelY:int, level:int) -> (int, int, int): + def get_id(self, pixelX:int, pixelY:int, level:int) -> tuple: id_x = round(-0.5+(pixelX/self.geometry_tilesize[level][1])) id_y = round(-0.5+(pixelY/self.geometry_tilesize[level][0]))