diff --git a/.github/workflows/test_and_deploy.yml b/.github/workflows/test_and_deploy.yml index 3eb44c5..a44155a 100644 --- a/.github/workflows/test_and_deploy.yml +++ b/.github/workflows/test_and_deploy.yml @@ -23,7 +23,7 @@ jobs: strategy: matrix: platform: [ubuntu-latest, windows-latest, macos-latest] - python-version: ['3.8', '3.9', '3.10'] + python-version: ['3.9', '3.10'] steps: - uses: actions/checkout@v3 diff --git a/MANIFEST.in b/MANIFEST.in index 5760a99..c7ce461 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -4,4 +4,3 @@ include README.md recursive-exclude * __pycache__ recursive-exclude * *.py[co] include src/napari_flim_phasor_plotter/data/* -include src/napari_flim_phasor_plotter/data/hazelnut_FLIM_z_stack/* diff --git a/setup.cfg b/setup.cfg index b31ba28..b17d966 100644 --- a/setup.cfg +++ b/setup.cfg @@ -18,7 +18,6 @@ classifiers = Programming Language :: Python Programming Language :: Python :: 3 Programming Language :: Python :: 3 :: Only - Programming Language :: Python :: 3.8 Programming Language :: Python :: 3.9 Programming Language :: Python :: 3.10 Topic :: Scientific/Engineering :: Image Processing @@ -34,7 +33,8 @@ install_requires = numpy magicgui qtpy - napari-clusters-plotter + napari-clusters-plotter>=0.8.0 + ptufile sdtfile natsort rocket-fft @@ -43,7 +43,7 @@ install_requires = napari-segment-blobs-and-things-with-membranes napari-skimage-regionprops -python_requires = >=3.8 +python_requires = >=3.9 include_package_data = True package_dir = =src diff --git a/src/napari_flim_phasor_plotter/__init__.py b/src/napari_flim_phasor_plotter/__init__.py index b741884..4d02687 100644 --- a/src/napari_flim_phasor_plotter/__init__.py +++ b/src/napari_flim_phasor_plotter/__init__.py @@ -2,14 +2,14 @@ from ._reader import napari_get_reader from ._sample_data import load_seminal_receptacle_image, load_hazelnut_image, load_hazelnut_z_stack, load_lifetime_cat_synthtetic_single_image -from ._io.readPTU_FLIM import PTUreader +from ._io import convert_to_zarr from . import phasor, filters, _plotting, _widget __all__ = ( "napari_get_reader", "_sample_data", - "PTUreader", + "convert_to_zarr", "phasor", "filters", "_plotting", diff --git a/src/napari_flim_phasor_plotter/_io/convert_to_zarr.py b/src/napari_flim_phasor_plotter/_io/convert_to_zarr.py index f3ddbfb..c6605a4 100644 --- a/src/napari_flim_phasor_plotter/_io/convert_to_zarr.py +++ b/src/napari_flim_phasor_plotter/_io/convert_to_zarr.py @@ -33,6 +33,7 @@ def convert_folder_to_zarr(folder_path: pathlib.Path, cancel_button: bool = Fals Path to the folder containing the FLIM images. """ import zarr + from numcodecs import Blosc import dask.array as da import numpy as np from pathlib import Path @@ -76,9 +77,11 @@ def convert_folder_to_zarr(folder_path: pathlib.Path, cancel_button: bool = Fals # zarr file will be saved in the same folder as the input folder output_path = folder_path / (folder_path.stem + '.zarr') # Using zarr to automatically guess chunk sizes + # Use Blosc compressor for better compression and speed + compressor = Blosc(cname='zstd', clevel=3, shuffle=Blosc.BITSHUFFLE) # Create an empty zarr array of a specified shape and dtype filled with zeros zarr_array = zarr.open(output_path, mode='w', - shape=stack_shape, dtype=image_dtype) + shape=stack_shape, dtype=image_dtype, compressor=compressor) # Using dask to rechunk micro-time axis in single chunk (for fft calculation afterwards) dask_array = da.from_zarr(output_path) # Rechunk axis 1 (micro-time axis) to a single chunk diff --git a/src/napari_flim_phasor_plotter/_io/readPTU_FLIM.py b/src/napari_flim_phasor_plotter/_io/readPTU_FLIM.py deleted file mode 100644 index 97c876d..0000000 --- a/src/napari_flim_phasor_plotter/_io/readPTU_FLIM.py +++ /dev/null @@ -1,601 +0,0 @@ -'''' -This package was copied from readPTU_FLIM library -(https://github.com/SumeetRohilla/readPTU_FLIM) - -MIT License - -Copyright (c) 2019 SumeetRohilla - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. - -''' - - -# Read PTU Library and FLIM data -# Author: Sumeet Rohilla -# sumeetrohilla@gmail.com - -# PTU Reader Library - -""" -Created on Tue, 14 May 2019 -Updated on Sun, 20 Dec, 2020 -@author: SumeetRohilla, sumeetrohilla@gmail.com - -"Good artists copy; great artists steal." - -Largely inspired from examples: -- PicoQuant demo codes - https://github.com/PicoQuant/PicoQuant-Time-Tagged-File-Format-Demos -- from a jupyter notebook by tritemio on GitHub: - https://gist.github.com/tritemio/734347586bc999f39f9ffe0ac5ba0e66 - -Aim : Open and convert Picoquant .ptu image files for FLIM analysis - * Use : Simply select your .ptu file and library will provide: - * 1) Lifetime image stack for each channel (1-8). - * flim_data_stack = [numPixX numPixY numDetectors numTCSPCbins] - * Fluorescence decays in each pixel/num_detection_channel are available for the whole acquisition - * Frame-wise info is not implemented here, but in principle is pretty straightforwad to implement - * 2) Intensity is just acquisition flim_data_stack(in photons) is accessed by binning across axis = numTCSPCbins and numDetectors - * 3) get_flim_data_stack class method is numba accelarated (using @jit decorator) to gain speed in building flim_data_stack from raw_tttr_data - -""" - -import time -import sys -import struct -import numpy as np -from numba import njit - -@njit -def get_flim_data_stack_static(sync, tcspc, channel, special, header_variables): - - ImgHdr_Ident = header_variables[0] - MeasDesc_Resolution = header_variables[1] - MeasDesc_GlobalResolution = header_variables[2] - ImgHdr_PixX = header_variables[3] - ImgHdr_PixY = header_variables[4] - ImgHdr_LineStart = header_variables[5] - ImgHdr_LineStop = header_variables[6] - ImgHdr_Frame = header_variables[7] - - if (ImgHdr_Ident == 9) or (ImgHdr_Ident == 3): - - tcspc_bin_resolution = 1e9*MeasDesc_Resolution # in Nanoseconds - sync_rate = np.ceil(MeasDesc_GlobalResolution*1e9) # in Nanoseconds - -# num_of_detectors = np.max(channel)+1 -# num_tcspc_channel = np.max(tcspc)+1 -# num_tcspc_channel = floor(sync_rate/tcspc_bin_resolution)+1 - - num_of_detectors = np.unique(channel).size - num_tcspc_channel = np.unique(tcspc).size - num_pixel_X = ImgHdr_PixX - num_pixel_Y = ImgHdr_PixY - - - flim_data_stack = np.zeros((num_pixel_Y, num_pixel_X, num_of_detectors,num_tcspc_channel), dtype = np.uint16) - - # Markers necessary to make FLIM image stack - LineStartMarker = 2**(ImgHdr_LineStart-1) - LineStopMarker = 2**(ImgHdr_LineStop-1) - FrameMarker = 2**(ImgHdr_Frame-1) - - # Get Number of Frames - FrameSyncVal = sync[np.where(special == FrameMarker)] - num_of_Frames = FrameSyncVal.size - read_data_range = np.where(sync == FrameSyncVal[num_of_Frames-1])[0][0] - - L1 = sync[np.where(special == LineStartMarker)] # Get Line start marker sync values - L2 = sync[np.where(special == LineStopMarker)] # Get Line start marker sync values - - syncPulsesPerLine = np.floor(np.mean(L2[10:] - L1[10:])) - -# Get pixel dwell time values from header for PicoQuant_FLIMBee or Zeiss_LSM scanner - -# if 'StartedByRemoteInterface' in head.keys(): - -# #syncPulsesPerLine = round((head.TimePerPixel/head.MeasDesc_GlobalResolution)*num_pixel_X); -# syncPulsesPerLine = np.floor(np.mean(L2[10:] - L1[10:])) -# else: -# #syncPulsesPerLine = floor(((head.ImgHdr_TimePerPixel*1e-3)/head.MeasDesc_GlobalResolution)*num_pixel_X); -# syncPulsesPerLine = np.floor(np.mean(L2[10:] - L1[10:])) - - - # Initialize Variable - currentLine = 0 - currentSync = 0 - syncStart = 0 - currentPixel = 0 - unNoticed_events = 0 - countFrame = 0 - insideLine = False - insideFrame = False - isPhoton = False - - for event in range(read_data_range+1): - - if num_of_Frames == 1: - # when only zero/one frame marker is present in TTTR file - insideFrame = True - - currentSync = sync[event] - special_event = special[event] - - # is the record a valid photon event or a special marker type event - if special[event] == 0: - isPhoton = True - else: - isPhoton = False - - if not(isPhoton): - #This is not needed once inside the first Frame marker - if (special_event == FrameMarker): - insideFrame = True - countFrame += 1 - currentLine = 0 - - if (special_event == LineStartMarker): - - insideLine = True - syncStart = currentSync - - elif (special_event == LineStopMarker): - - insideLine = False - currentLine += 1 - syncStart = 0 - - if (currentLine >= num_pixel_Y): - insideFrame = False - currentLine = 0 - - # Build FLIM image data stack here for N-spectral/tcspc-input channels - - if (isPhoton and insideLine and insideFrame): - - currentPixel = int(np.floor((((currentSync - syncStart)/syncPulsesPerLine)*num_pixel_X))) - tmpchan = channel[event] - tmptcspc = tcspc[event] - - if (currentPixel < num_pixel_X) and (tmptcspc 32 bit record - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # | | | | | | | | | | | | | | | | | |x|x|x|x|x|x|x|x|x|x|x|x|x|x|x|x| --> Sync - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # | | | | |x|x|x|x|x|x|x|x|x|x|x|x| | | | | | | | | | | | | | | | | | --> TCSPC bin - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # |x|x|x|x| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | --> Spectral/TCSPC input Channel - # +-------------------------------+ +-------------------------------+ - - WRAPAROUND = 65536 # After this sync counter will overflow - sync = np.bitwise_and(t3records, 65535) # Lowest 16 bits - tcspc = np.bitwise_and(np.right_shift(t3records, 16), 4095) # Next 12 bits, dtime can be obtained from header - chan = np.bitwise_and(np.right_shift(t3records, 28),15) # Next 4 bits - special = ((chan==15)*1)*(np.bitwise_and(tcspc,15)*1) # Last bit for special markers - - index = ((chan==15)*1)*((np.bitwise_and(tcspc,15)==0)*1) # Find overflow locations - chan[chan==15] = special[chan==15] - chan[chan==1] = 0 - chan[chan==2] = 1 - chan[chan==4] = 0 - - elif record_type == 'rtPicoHarpT2': - - print('TCSPC Hardware: {}'.format(record_type[2:])) - - # +----------------------+ T2 32 bit record +---------------------+ - # |x|x|x|x|x|x|x|x|x|x|x|x|x|x|x|x| |x|x|x|x|x|x|x|x|x|x|x|x|x|x|x|x| --> 32 bit record - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # | | | | |x|x|x|x|x|x|x|x|x|x|x|x| |x|x|x|x|x|x|x|x|x|x|x|x|x|x|x|x| --> Sync - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # |x|x|x|x| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | --> Spectral/TCSPC input Channel - # +-------------------------------+ +-------------------------------+ - - WRAPAROUND = 210698240 # After this sync counter will overflow - sync = np.bitwise_and(t3records, 268435455) # Lowest 28 bits - tcspc = np.bitwise_and(t3records, 15) # Next 4 bits, dtime can be obtained from header - chan = np.bitwise_and(np.right_shift(t3records, 28),15) # Next 4 bits - special = ((chan==15)*1)*np.bitwise_and(tcspc,15) # Last bit for special markers - index = ((chan==15)*1)*((np.bitwise_and(tcspc,15)==0)*1) # Find overflow locations - - elif record_type in ['rtHydraHarpT3', 'rtHydraHarp2T3', 'rtTimeHarp260NT3', 'rtTimeHarp260PT3','rtMultiHarpNT3']: - - print('TCSPC Hardware: {}'.format(record_type[2:])) - - # +----------------------+ T3 32 bit record +---------------------+ - # |x|x|x|x|x|x|x|x|x|x|x|x|x|x|x|x| |x|x|x|x|x|x|x|x|x|x|x|x|x|x|x|x| --> 32 bit record - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # | | | | | | | | | | | | | | | | | | | | | | | |x|x|x|x|x|x|x|x|x|x| --> Sync - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # | | | | | | | |x|x|x|x|x|x|x|x|x| |x|x|x|x|x|x| | | | | | | | | | | --> TCSPC bin - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # | |x|x|x|x|x|x| | | | | | | | | | | | | | | | | | | | | | | | | | | --> Spectral/TCSPC input Channel - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # |x| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | --> Special markers - # +-------------------------------+ +-------------------------------+ - WRAPAROUND = 1024 # After this sync counter will overflow - sync = np.bitwise_and(t3records, 1023) # Lowest 10 bits - tcspc = np.bitwise_and(np.right_shift(t3records, 10), 32767) # Next 15 bits, dtime can be obtained from header - chan = np.bitwise_and(np.right_shift(t3records, 25),63) # Next 8 bits - special = np.bitwise_and(t3records,2147483648)>0 # Last bit for special markers - index = (special*1)*((chan==63)*1) # Find overflow locations - special = (special*1)*chan - - elif record_type == 'rtHydraHarpT2': - - print('TCSPC Hardware: {}'.format(record_type[2:])) - - # +----------------------+ T3 32 bit record +---------------------+ - # |x|x|x|x|x|x|x|x|x|x|x|x|x|x|x|x| |x|x|x|x|x|x|x|x|x|x|x|x|x|x|x|x| --> 32 bit record - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # | | | | | | | |x|x|x|x|x|x|x|x|x| |x|x|x|x|x|x|x|x|x|x|x|x|x|x|x|x| --> Sync - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # | |x|x|x|x|x|x| | | | | | | | | | | | | | | | | | | | | | | | | | | --> Spectral/TCSPC input Channel - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # |x| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | --> Special markers - # +-------------------------------+ +-------------------------------+ - WRAPAROUND = 33552000 # After this sync counter will overflow - sync = np.bitwise_and(t3records, 33554431) # Lowest 25 bits - chan = np.bitwise_and(np.right_shift(t3records, 25),63) # Next 6 bits - tcspc = np.bitwise_and(chan, 15) - special = np.bitwise_and(np.right_shift(t3records, 31),1) # Last bit for special markers - index = (special*1) * ((chan==63)*1) # Find overflow locations - special = (special*1)*chan - - elif record_type in ['rtHydraHarp2T2', 'rtTimeHarp260NT2', 'rtTimeHarp260PT2','rtMultiHarpNT2']: - - print('TCSPC Hardware: {}'.format(record_type[2:])) - - # +----------------------+ T3 32 bit record +---------------------+ - # |x|x|x|x|x|x|x|x|x|x|x|x|x|x|x|x| |x|x|x|x|x|x|x|x|x|x|x|x|x|x|x|x| --> 32 bit record - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # | | | | | | | |x|x|x|x|x|x|x|x|x| |x|x|x|x|x|x|x|x|x|x|x|x|x|x|x|x| --> Sync - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # | |x|x|x|x|x|x| | | | | | | | | | | | | | | | | | | | | | | | | | | --> Spectral/TCSPC input Channel - # +-------------------------------+ +-------------------------------+ - # +-------------------------------+ +-------------------------------+ - # |x| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | --> Special markers - # +-------------------------------+ +-------------------------------+ - - WRAPAROUND = 33554432 # After this sync counter will overflow - sync = np.bitwise_and(t3records, 33554431) # Lowest 25 bits - chan = np.bitwise_and(np.right_shift(t3records, 25),63) # Next 6 bits - tcspc = np.bitwise_and(chan, 15) - special = np.bitwise_and(np.right_shift(t3records, 31),1) # Last bit for special markers - index = (special*1) * ((chan==63)*1) # Find overflow locations - special = (special*1)*chan - - else: - print('Illegal RecordType!') - exit(0) - - - - # Fill in the correct sync values for overflow location - #sync[np.where(index == 1)] = 1 # assert values of sync = 1 just right after overflow to avoid any overflow-correction instability in next step - if record_type in ['rtHydraHarp2T3','rtTimeHarp260PT3','rtMultiHarpNT3']: - sync = sync + (WRAPAROUND*np.cumsum(index*sync)) # For HydraHarp V1 and TimeHarp260 V1 overflow corrections - else: - sync = sync + (WRAPAROUND*np.cumsum(index)) # correction for overflow to sync varibale - - sync = np.delete(sync, np.where(index == 1), axis = 0) - tcspc = np.delete(tcspc, np.where(index == 1), axis = 0) - chan = np.delete(chan, np.where(index == 1), axis = 0) - special = np.delete(special, np.where(index == 1), axis = 0) - - del index - - # Convert to appropriate data type to save memory - - self.sync = sync.astype(np.uint64, copy=False) - self.tcspc = tcspc.astype(np.uint16, copy=False) - self.channel = chan.astype(np.uint8, copy=False) - self.special = special.astype(np.uint8, copy=False) - self.x_size = self.head['ImgHdr_PixX'] - self.y_size = self.head['ImgHdr_PixY'] - print("Raw Data has been Read!\n") - - return None - - def get_flim_data_stack(self): - - # Check if it's FLIM image - if self.head["Measurement_SubMode"] == 0: - raise IOError("This is not a FLIM PTU file.!!! \n") - sys.exit() - elif (self.head["ImgHdr_Ident"] == 1) or (self.head["ImgHdr_Ident"] == 5): - raise IOError("Piezo Scanner Data Reader Not Implemented Yet!!! \n") - sys.exit() - else: - # Create numpy array of important variables to be passed into numba accelaratd get_flim_data_stack_static function - # as numba doesn't recognizes python dict type files - # Check - - header_variables = np.array([self.head["ImgHdr_Ident"], self.head["MeasDesc_Resolution"],self.head["MeasDesc_GlobalResolution"],self.head["ImgHdr_PixX"], self.head["ImgHdr_PixY"], self.head["ImgHdr_LineStart"],self.head["ImgHdr_LineStop"], self.head["ImgHdr_Frame"]],dtype = np.uint64) - - sync = self.sync - tcspc = self.tcspc - channel = self.channel - special = self.special - - del self.sync, self.special# self.tcspc self.channel - - flim_data_stack = get_flim_data_stack_static(sync, tcspc, channel, special, header_variables) - - if flim_data_stack.ndim == 4: - - tmp_intensity_image = np.sum(flim_data_stack, axis = 3) # sum across tcspc bin - intensity_image = np.sum(tmp_intensity_image, axis = 2)# sum across spectral channels - - elif flim_data_stack.ndim == 3: - - intensity_image = np.sum(flim_data_stack, axis = 3) # sum across tcspc bin, only 1 detection channel - - return flim_data_stack, intensity_image - -@njit -def get_lifetime_image(flim_data_stack,channel_number,timegating_start1,timegating_stop1,meas_resolution,estimated_irf): - - work_data = flim_data_stack[:,:,channel_number,timegating_start1:timegating_stop1] - bin_range = np.reshape(np.linspace(0,timegating_stop1,timegating_stop1),(1,1,timegating_stop1)) - - fast_flim = (np.sum(work_data*bin_range,axis = 2)*meas_resolution)/np.sum(work_data,axis = 2) - - return fast_flim diff --git a/src/napari_flim_phasor_plotter/_reader.py b/src/napari_flim_phasor_plotter/_reader.py index 45f229b..7f21370 100644 --- a/src/napari_flim_phasor_plotter/_reader.py +++ b/src/napari_flim_phasor_plotter/_reader.py @@ -30,18 +30,98 @@ def napari_get_reader(path): return None +def read_single_ptu_file_2d_timelapse(path, *args, **kwargs): + """Read a 2D timelapse PTU file. + + Parameters + ---------- + path : str + Path to the PTU file. + + Returns + ------- + data : np.ndarray + The data array with dimensions (ch, ut, t, z, y, x), where z has length 1 (following napari convention). + metadata_per_channel : list + List of metadata dictionaries for each channel. + """ + from ptufile import PtuFile + import numpy as np + + ptu = PtuFile(path) + ptu_dims = list(ptu.dims) + data = ptu[:] + + # Re-order dimensions to standard order + standard_dims_order = ['C', 'H', 'T', 'Y', 'X'] # (ch, ut, t, y, x) + argsorted = [ptu_dims.index(dim) for dim in standard_dims_order] + data = np.moveaxis(data, argsorted, range(len(argsorted))) + # Add unitary dimension for z + data = np.expand_dims(data, axis=3) # (ch, ut, t, z, y, x) + + # Get metadata + # metadata per channel/detector + metadata_per_channel = [] + metadata = ptu.tags + metadata['file_type'] = 'ptu' + metadata['frequency'] = ptu.frequency + metadata['tcspc_resolution'] = ptu.tcspc_resolution + metadata['x_pixel_size'] = ptu.coords['X'][1] + metadata['y_pixel_size'] = ptu.coords['Y'][1] + metadata['frame_time'] = ptu.frame_time + # Add same metadata to each channel + for channel in range(data.shape[0]): + metadata_per_channel.append(metadata) + return data, metadata_per_channel + def read_single_ptu_file(path, *args, **kwargs): - """Read a single ptu file.""" - from napari_flim_phasor_plotter._io.readPTU_FLIM import PTUreader + """Read a single ptu file. - ptu_file = PTUreader(path, print_header_data=False) - data, _ = ptu_file.get_flim_data_stack() - # from (x, y, ch, ut) to (ch, ut, y, x) - data = np.moveaxis(data, [0, 1], [-2, -1]) + Only accepts single frame PTU files. If PTU file is a time-lapse, split into single frame PTU files first, otherwise just first frame gets read. + + Parameters + ---------- + path : str + Path to the PTU file. + + Returns + ------- + data : np.ndarray + The data array with dimensions (ch, ut, y, x). + metadata_per_channel : list + List of metadata dictionaries for each channel. + """ + from ptufile import PtuFile + import numpy as np + from warnings import warn + + ptu = PtuFile(path) + ptu_dims = list(ptu.dims) + if 'T' in ptu_dims: + t_pos = ptu_dims.index('T') + n_frames = ptu.shape[t_pos] + if n_frames > 1: + warn("Timelapse found in single ptu file. This function is for single frame PTU files only. Returning first frame.\nUse read_single_ptu_file_2d_timelapse() for timelapse PTU files.") + # read single frame from axis where time is present + data = np.take(ptu[:], 0, axis=t_pos) + ptu_dims.pop(t_pos) + else: + data = ptu[:] + + # Re-order dimensions to standard order + standard_dims_order = ['C', 'H', 'Y', 'X'] # (ch, ut, y, x) + argsorted = [ptu_dims.index(dim) for dim in standard_dims_order] + data = np.moveaxis(data, argsorted, range(len(argsorted))) + + # Get metadata # metadata per channel/detector metadata_per_channel = [] - metadata = ptu_file.head + metadata = ptu.tags metadata['file_type'] = 'ptu' + metadata['frequency'] = ptu.frequency + metadata['tcspc_resolution'] = ptu.tcspc_resolution + metadata['x_pixel_size'] = ptu.coords['X'][1] + metadata['y_pixel_size'] = ptu.coords['Y'][1] # Add same metadata to each channel for channel in range(data.shape[0]): metadata_per_channel.append(metadata) @@ -224,6 +304,7 @@ def flim_file_reader(path): """ from pathlib import Path import tifffile + from ptufile import PtuFile import numpy as np from napari.utils. notifications import show_warning # handle both a string and a list of strings @@ -231,6 +312,7 @@ def flim_file_reader(path): # Use Path from pathlib paths = [Path(path) for path in paths] + imread = None layer_data = [] for path in paths: # Assume stack if paths are folders (which includes .zarr here) @@ -256,8 +338,18 @@ def flim_file_reader(path): channel_axis = None if len(shape) < 4: # single 2D image (ut, y, x) channel_axis = None - imread = get_read_function_from_extension[file_extension] - # (ch, ut, y, x) or (ch, ut, t, z, y, x) in case of single tif stack + # if .ptu, check if it is a 2D timelapse + if file_extension == '.ptu': + ptu = PtuFile(path) + ptu_dims = list(ptu.dims) + if 'T' in ptu_dims: + t_pos = ptu_dims.index('T') + n_frames = ptu.shape[t_pos] + if n_frames > 1: + imread = read_single_ptu_file_2d_timelapse + if imread is None: + imread = get_read_function_from_extension[file_extension] + # (ch, ut, y, x) or (ch, ut, t, z, y, x) in case of single tif stack data, metadata_list = imread(file_path, channel_axis=channel_axis, viewer_exists=True) if data.ndim == 4: # expand dims if not a stack already data = np.expand_dims(data, axis=(2, 3)) # (ch, ut, t, z, y, x) @@ -322,9 +414,10 @@ def read_stack(folder_path): 'Stack is larger than 4GB, please convert to .zarr') print('Stack is larger than 4GB, please convert to .zarr') return None, None - # TO DO: remove print - print('stack = True\n', 'data type: ', file_extension, - '\ndata_shape = ', data.shape, '\n') + print('stack = True', + '\ndata type: ', file_extension, + '\ndata_shape = ', data.shape, + '\n') return data, metadata_list @@ -347,13 +440,17 @@ def get_max_slice_shape_and_dtype(file_paths, file_extension): Max shape and data type. """ # TO DO: offer fast reading option by calculating max shape from metadata (array may become bigger) + from tqdm import tqdm + progress_bar = tqdm(total=len(file_paths), desc='Calculating stack shape from individual files', unit='files') shapes_list = [] for file_path in file_paths: if file_path.suffix == file_extension: imread = get_read_function_from_extension[file_extension] image_slice, _ = imread(file_path) shapes_list.append(image_slice.shape) # (ch, ut, y, x) - # Get slice max shape (ch, mt, y, x) + progress_bar.update(1) + progress_bar.close() + # Get slice max shape (ch, ut, y, x) return max(shapes_list), image_slice.dtype @@ -372,6 +469,7 @@ def make_full_numpy_stack(file_paths, file_extension): numpy_stack, metadata_per_channel : Tuple(numpy array, List(dict)) A numpy array of shape (ch, ut, t, z, y, x) and a metadata list (one metadata per channel). """ + from tqdm import tqdm # Read all images to get max slice shape image_slice_shape, image_dtype = get_max_slice_shape_and_dtype( file_paths, file_extension) @@ -379,6 +477,9 @@ def make_full_numpy_stack(file_paths, file_extension): list_of_time_point_paths = get_structured_list_of_paths( file_paths, file_extension) + progress_bar = tqdm(total=len(list_of_time_point_paths) * len(list_of_time_point_paths[0]), + desc='Reading stack', unit='slices') + z_list, t_list = [], [] for list_of_zslice_paths in list_of_time_point_paths: for zslice_path in list_of_zslice_paths: @@ -387,10 +488,12 @@ def make_full_numpy_stack(file_paths, file_extension): z_slice[:data.shape[0], :data.shape[1], :data.shape[2], :data.shape[3]] = data z_list.append(z_slice) + progress_bar.update(1) z_stack = np.stack(z_list) t_list.append(z_stack) z_list = [] stack = np.stack(t_list) + progress_bar.close() # from (t, z, ch, ut, y, x) to (ch, ut, t, z, y, x) stack = np.moveaxis(stack, [-4, -3], [0, 1]) return stack, metadata_per_channel @@ -485,6 +588,8 @@ def get_stack_estimated_size(file_paths, file_extension, from_file_size=False): stack_size : float Stack size in MB. """ + from tqdm import tqdm + progress_bar = tqdm(total=len(file_paths), desc='Calculating decompressed stack size', unit='files') stack_size = 0 for file_path in file_paths: if file_path.suffix == file_extension: @@ -495,6 +600,9 @@ def get_stack_estimated_size(file_paths, file_extension, from_file_size=False): data, metadata_list = imread(file_path) # (ch, ut, y, x) file_size = data.nbytes / 1e6 # in MB stack_size += file_size + progress_bar.update(1) + progress_bar.close() + print(f"Estimated decompressed stack size: {stack_size} MB") return stack_size diff --git a/src/napari_flim_phasor_plotter/_sample_data.py b/src/napari_flim_phasor_plotter/_sample_data.py index 66b08a3..da40ee1 100644 --- a/src/napari_flim_phasor_plotter/_sample_data.py +++ b/src/napari_flim_phasor_plotter/_sample_data.py @@ -85,16 +85,43 @@ def load_hazelnut_z_stack(): """ import numpy as np import zipfile + import requests + import shutil from pathlib import Path + from tqdm import tqdm from napari_flim_phasor_plotter._reader import read_stack - - zip_file_path = Path(DATA_ROOT / "hazelnut_FLIM_z_stack.zip") extracted_folder_path = Path(DATA_ROOT / "unzipped_hazelnut_FLIM_z_stack") - # Create the target directory if it doesn't exist - extracted_folder_path.mkdir(parents=True, exist_ok=True) - # Extract the zip file - with zipfile.ZipFile(zip_file_path, 'r') as zip_ref: - zip_ref.extractall(extracted_folder_path) + # If extracted folder does not exist or is empty, download and extract the zip file + if not extracted_folder_path.exists() or (extracted_folder_path.exists() and not any(extracted_folder_path.iterdir())): + + zip_url = 'https://github.com/zoccoler/hazelnut_FLIM_z_stack_data/raw/main/hazelnut_FLIM_z_stack.zip' + zip_file_path = Path(DATA_ROOT / "hazelnut_FLIM_z_stack.zip") + # Download the zip file + response = requests.get(zip_url) + + # Total size in bytes. + total_size = int(response.headers.get('content-length', 0)) + print(f"Total download size: {total_size/1e6} MBytes") + print(f"Downloading to {zip_file_path}" ) + block_size = 1024 + progress_bar = tqdm(total=total_size, unit='iB', unit_scale=True, desc="Downloading zip file") + with open(zip_file_path, 'wb') as zip_file: + for block in response.iter_content(block_size): + progress_bar.update(len(block)) + zip_file.write(block) + progress_bar.close() + if total_size != 0 and progress_bar.n != total_size: + print("ERROR: Something went wrong with the download") + + + # Create the target directory + extracted_folder_path.mkdir(parents=True, exist_ok=True) + # Extract the zip file + with zipfile.ZipFile(zip_file_path, 'r') as zip_ref: + zip_ref.extractall(extracted_folder_path) + + # Delete the zip file after extraction + Path(zip_file_path).unlink() folder_path = extracted_folder_path / "hazelnut_FLIM_z_stack" image, metadata = read_stack(folder_path) diff --git a/src/napari_flim_phasor_plotter/_tests/test_sample_data.py b/src/napari_flim_phasor_plotter/_tests/test_sample_data.py index 637a223..791d307 100644 --- a/src/napari_flim_phasor_plotter/_tests/test_sample_data.py +++ b/src/napari_flim_phasor_plotter/_tests/test_sample_data.py @@ -22,10 +22,10 @@ def test_loading_sample_data(): assert receptacle_raw_flim_image.shape == (256, 1, 1, 512, 512) assert receptacle_intensity_image.shape == (1, 1, 512, 512) - assert hazelnut_raw_flim_image.shape == (137, 1, 1, 256, 256) + assert hazelnut_raw_flim_image.shape == (139, 1, 1, 256, 256) assert hazelnut_intensity_image.shape == (1, 1, 256, 256) - assert hazelnut_z_stack_raw_flim_image.shape == (138, 1, 11, 256, 256) + assert hazelnut_z_stack_raw_flim_image.shape == (139, 1, 11, 256, 256) assert hazelnut_z_stack_intensity_image.shape == (1, 11, 256, 256) assert cat_lifetime_raw_flim_image.shape == (256, 1, 1, 256, 256) diff --git a/src/napari_flim_phasor_plotter/_widget.py b/src/napari_flim_phasor_plotter/_widget.py index b437d8b..1efef94 100644 --- a/src/napari_flim_phasor_plotter/_widget.py +++ b/src/napari_flim_phasor_plotter/_widget.py @@ -117,7 +117,8 @@ def make_flim_phasor_plot(image_layer: "napari.layers.Image", name='Labelled_pixels_from_' + image_layer.name, features=table, scale=image_layer.scale[1:], - visible=False) + visible=True, + opacity=0.2) # Check if plotter was alrerady added to dock_widgets # TODO: avoid using private method access to napari_viewer.window._dock_widgets (will be deprecated) @@ -134,15 +135,10 @@ def make_flim_phasor_plot(image_layer: "napari.layers.Image", plotter_widget = widgets.findChild(PhasorPlotterWidget) # Get labels layer with labelled pixels (labels) - # Commented code below will work with napari-clusters-plotter 0.7.4 (not released yet) or 0.8.0 depending on the next version number - # plotter_widget.layer_select.value = [ - # choice for choice in plotter_widget.layer_select.choices if choice.name.startswith("Labelled_pixels")][0] - for choice in plotter_widget.labels_select.choices: + for choice in plotter_widget.layer_select.choices: if choice.name == 'Labelled_pixels_from_' + image_layer.name: - plotter_widget.labels_select.value = choice + plotter_widget.layer_select.value = choice break - # plotter_widget.labels_select.value = [ - # choice for choice in plotter_widget.labels_select.choices if choice.name.startswith("Labelled_pixels")][0] # Set G and S as features to plot (update_axes_list method clears Comboboxes) plotter_widget.plot_x_axis.setCurrentIndex(1) plotter_widget.plot_y_axis.setCurrentIndex(2) diff --git a/src/napari_flim_phasor_plotter/data/hazelnut_FLIM_z_stack.zip b/src/napari_flim_phasor_plotter/data/hazelnut_FLIM_z_stack.zip deleted file mode 100644 index 0b7474b..0000000 Binary files a/src/napari_flim_phasor_plotter/data/hazelnut_FLIM_z_stack.zip and /dev/null differ diff --git a/tox.ini b/tox.ini index edde9b0..cbc978d 100644 --- a/tox.ini +++ b/tox.ini @@ -1,11 +1,10 @@ # For more information about tox, see https://tox.readthedocs.io/en/latest/ [tox] -envlist = py{38,39,310}-{linux,macos,windows} +envlist = py{39,310}-{linux,macos,windows} isolated_build=true [gh-actions] python = - 3.8: py38 3.9: py39 3.10: py310 @@ -44,4 +43,13 @@ deps = napari-tools-menu napari-skimage-regionprops>=0.2.0 hdbscan + napari-clusters-plotter>=0.8.0 + ptufile + sdtfile + natsort + rocket-fft + dask + zarr + napari-segment-blobs-and-things-with-membranes + napari-skimage-regionprops commands = pytest -v --color=yes --cov=napari_flim_phasor_plotter --cov-report=xml