diff --git a/LICENSE b/LICENSE index 1007679..1d74164 100644 --- a/LICENSE +++ b/LICENSE @@ -1,5 +1,5 @@ -Copyright (c) 2022, Napari Developer +Copyright (c) 2022, Marcelo L. Zoccoler All rights reserved. Redistribution and use in source and binary forms, with or without diff --git a/README.md b/README.md index 20c8aae..aafec12 100644 --- a/README.md +++ b/README.md @@ -7,33 +7,32 @@ [![codecov](https://codecov.io/gh/zoccoler/napari-flim-phasor-calculator/branch/main/graph/badge.svg)](https://codecov.io/gh/zoccoler/napari-flim-phasor-calculator) [![napari hub](https://img.shields.io/endpoint?url=https://api.napari-hub.org/shields/napari-flim-phasor-calculator)](https://napari-hub.org/plugins/napari-flim-phasor-calculator) -A simple plugin to use FooBar segmentation within napari +A napari plugin to generate a phasor plot for TCSPC FLIM data. ---------------------------------- -This [napari] plugin was generated with [Cookiecutter] using [@napari]'s [cookiecutter-napari-plugin] template. +## Usage - + -## Installation +Manually draw curves on the plot to get the corresponding pixels highlighted in a new labels layer. Hold 'SHIFT' while drawing to add more than two colors. -You can install `napari-flim-phasor-calculator` via [pip]: +This plugin integrates with [napari-clusters-plotter plugin](https://github.com/BiAPoL/napari-clusters-plotter). - pip install napari-flim-phasor-calculator +This plugin can read the following FLIM file types: + - ".ptu" +This plugin works with the following data shapes: + - 2D timelapse, where time is the first dimension. +## Installation -To install latest development version : +You can install `napari-flim-phasor-calculator` via [pip]. Currently, only the development version is available: pip install git+https://github.com/zoccoler/napari-flim-phasor-calculator.git - ## Contributing Contributions are very welcome. Tests can be run with [tox], please ensure diff --git a/images/napari_FLIM_phasor_calculator_Demo.gif b/images/napari_FLIM_phasor_calculator_Demo.gif new file mode 100644 index 0000000..fd39111 Binary files /dev/null and b/images/napari_FLIM_phasor_calculator_Demo.gif differ diff --git a/setup.cfg b/setup.cfg index 9aaaa2c..4ab7f1b 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,12 +1,12 @@ [metadata] name = napari-flim-phasor-calculator version = 0.0.1 -description = A simple plugin to use FooBar segmentation within napari +description = A plugin that performs phasor plot from TCSPC FLIM data. long_description = file: README.md long_description_content_type = text/markdown url = https://github.com/zoccoler/napari-flim-phasor-calculator -author = Napari Developer -author_email = yourname@example.com +author = Marcelo L. Zoccoler and Cornelia Wetzker +author_email = marzoccoler@gmail.com license = BSD-3-Clause license_files = LICENSE classifiers = diff --git a/src/napari_flim_phasor_calculator/__init__.py b/src/napari_flim_phasor_calculator/__init__.py index f888bf7..c8f7206 100644 --- a/src/napari_flim_phasor_calculator/__init__.py +++ b/src/napari_flim_phasor_calculator/__init__.py @@ -2,9 +2,15 @@ from ._reader import napari_get_reader from ._sample_data import make_sample_data +from ._io.readPTU_FLIM import PTUreader +from . import phasor, filters, _plotting __all__ = ( "napari_get_reader", "make_sample_data", + "PTUreader", + "phasor", + "filters", + "_plotting" ) diff --git a/src/napari_flim_phasor_calculator/_io/__init__.py b/src/napari_flim_phasor_calculator/_io/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/napari_flim_phasor_calculator/_io/readPTU_FLIM.py b/src/napari_flim_phasor_calculator/_io/readPTU_FLIM.py new file mode 100644 index 0000000..9118050 --- /dev/null +++ b/src/napari_flim_phasor_calculator/_io/readPTU_FLIM.py @@ -0,0 +1,599 @@ +'''' +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) + 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.tcspc, self.channel , self.special + + 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_calculator/_plotting.py b/src/napari_flim_phasor_calculator/_plotting.py new file mode 100644 index 0000000..88712d1 --- /dev/null +++ b/src/napari_flim_phasor_calculator/_plotting.py @@ -0,0 +1,62 @@ +from napari_clusters_plotter._plotter import PlotterWidget +from qtpy.QtCore import QSize + +def add_phasor_circle(ax): + ''' + Generate FLIM universal semi-circle plot + ''' + import numpy as np + angles = np.linspace(0, np.pi, 180) + x =(np.cos(angles) + 1) / 2 + y = np.sin(angles) / 2 + ax.plot(x,y, 'yellow', alpha=0.3) + return ax + +def add_tau_lines(ax, tau_list, frequency): + import numpy as np + if not isinstance(tau_list, list): + tau_list = [tau_list] + frequency = frequency * 1E6 # MHz to Hz + w = 2*np.pi*frequency # Hz to radians/s + for tau in tau_list: + tau = tau * 1E-9 # nanoseconds to seconds + g = 1 / ( 1 + ( (w * tau)**2) ) + s = (w * tau) / ( 1 + ( (w * tau)**2) ) + dot, = ax.plot(g, s, marker='o', mfc='none') + array = np.linspace(0, g, 50) + y = (array*s/g) + ax.plot(array, y, color = dot.get_color()) + +def add_2d_histogram(ax, x, y): + import matplotlib.pyplot as plt + output = ax.hist2d( + x=x, + y=y, + bins=10, + cmap='jet', + norm='log', + alpha=0.5 + ) + return ax + +class PhasorPlotterWidget(PlotterWidget): + def __init__(self, napari_viewer): + super().__init__(napari_viewer) + self.setMinimumSize(QSize(100,300)) + def run(self, + features, + plot_x_axis_name, + plot_y_axis_name, + plot_cluster_name=None, + redraw_cluster_image=True,): + super().run(features=features, + plot_x_axis_name=plot_x_axis_name, + plot_y_axis_name=plot_y_axis_name, + plot_cluster_name=plot_cluster_name, + redraw_cluster_image=redraw_cluster_image,) + add_phasor_circle(self.graphics_widget.axes) + self.graphics_widget.draw() + # To Do: Add 2d histogram + # add_2d_histogram(self.graphics_widget.axes, + # features[plot_x_axis_name], + # features[plot_y_axis_name]) diff --git a/src/napari_flim_phasor_calculator/_reader.py b/src/napari_flim_phasor_calculator/_reader.py index 8d93ff9..aa2d8ef 100644 --- a/src/napari_flim_phasor_calculator/_reader.py +++ b/src/napari_flim_phasor_calculator/_reader.py @@ -6,7 +6,7 @@ https://napari.org/stable/plugins/guides.html?#readers """ import numpy as np - +from napari_flim_phasor_calculator._io.readPTU_FLIM import PTUreader def napari_get_reader(path): """A basic implementation of a Reader contribution. @@ -28,10 +28,6 @@ def napari_get_reader(path): # so we are only going to look at the first file. path = path[0] - # if we know we cannot read the file, we immediately return None. - if not path.endswith(".npy"): - return None - # otherwise we return the *function* that can read ``path``. return reader_function @@ -60,13 +56,15 @@ def reader_function(path): """ # handle both a string and a list of strings paths = [path] if isinstance(path, str) else path - # load all files into array - arrays = [np.load(_path) for _path in paths] - # stack arrays into single array - data = np.squeeze(np.stack(arrays)) - - # optional kwargs for the corresponding viewer.add_* method - add_kwargs = {} - - layer_type = "image" # optional, default is "image" - return [(data, add_kwargs, layer_type)] + layer_data = [] + for path in paths: + ptu_file = PTUreader(path, print_header_data = False) + data, _ = ptu_file.get_flim_data_stack() + # Move xy dimensions to the end + # TO DO: handle 3D images + data = np.moveaxis(data, [0, 1], [-2, -1]) + # optional kwargs for the corresponding viewer.add_* method + add_kwargs = {'channel_axis': 0, 'metadata': ptu_file.head} + layer_type = "image" + layer_data.append((data, add_kwargs, layer_type)) + return layer_data diff --git a/src/napari_flim_phasor_calculator/_synthetic.py b/src/napari_flim_phasor_calculator/_synthetic.py new file mode 100644 index 0000000..b480079 --- /dev/null +++ b/src/napari_flim_phasor_calculator/_synthetic.py @@ -0,0 +1,49 @@ +def monoexp(x, A, tau): + import numpy as np + return A * np.exp(-(1/tau)*x) + +def create_time_array(frequency, n_points=100): + ''' + Create time array from laser frequency + + Parameters + ---------- + frequency: float + Frequency of the pulsed laser (in MHz) + n_points: int, optional + The number of samples collected between each aser shooting + Returns + ------- + time_array : array + Time array (in nanoseconds) + ''' + import numpy as np + time_window = 1 / (frequency * 10**6) + time_window_ns = time_window * 10**9 # in nanoseconds + time_step = time_window_ns / n_points # ns + array = np.arange(0, n_points) + time_array = array * time_step + return time_array + +def make_synthetic_flim_data(time_array, amplitude_list, tau_list): + """ + Create a synthetic FLIM image from amplitudes and tau + + Each different tau in the list adds a new pixel to the image + """ + import numpy as np + # Handle input types + if not isinstance(tau_list, list): + tau_list = [tau_list] + if not isinstance(amplitude_list, list): + amplitude_list = [amplitude_list] + if len(amplitude_list)==1 and len(amplitude_list) 0 - layer_data_tuple = layer_data_list[0] - assert isinstance(layer_data_tuple, tuple) and len(layer_data_tuple) > 0 +# # make sure we're delivering the right format +# layer_data_list = reader(my_test_file) +# assert isinstance(layer_data_list, list) and len(layer_data_list) > 0 +# layer_data_tuple = layer_data_list[0] +# assert isinstance(layer_data_tuple, tuple) and len(layer_data_tuple) > 0 - # make sure it's the same as it started - np.testing.assert_allclose(original_data, layer_data_tuple[0]) +# # make sure it's the same as it started +# np.testing.assert_allclose(original_data, layer_data_tuple[0]) -def test_get_reader_pass(): - reader = napari_get_reader("fake.file") - assert reader is None +# def test_get_reader_pass(): +# reader = napari_get_reader("fake.file") +# assert reader is None diff --git a/src/napari_flim_phasor_calculator/_tests/test_sample_data.py b/src/napari_flim_phasor_calculator/_tests/test_sample_data.py index 9a477a3..363b3e2 100644 --- a/src/napari_flim_phasor_calculator/_tests/test_sample_data.py +++ b/src/napari_flim_phasor_calculator/_tests/test_sample_data.py @@ -1,7 +1,2 @@ -# from napari_flim_phasor_calculator import make_sample_data - -# add your tests here... - - def test_something(): pass diff --git a/src/napari_flim_phasor_calculator/_tests/test_widget.py b/src/napari_flim_phasor_calculator/_tests/test_widget.py index 226797a..51345e0 100644 --- a/src/napari_flim_phasor_calculator/_tests/test_widget.py +++ b/src/napari_flim_phasor_calculator/_tests/test_widget.py @@ -1,36 +1,29 @@ -import numpy as np - -from napari_flim_phasor_calculator import ExampleQWidget, example_magic_widget - +from napari_flim_phasor_calculator._widget import make_flim_phasor_plot +from napari_flim_phasor_calculator._synthetic import make_synthetic_flim_data +from napari_flim_phasor_calculator._synthetic import create_time_array # make_napari_viewer is a pytest fixture that returns a napari viewer object # capsys is a pytest fixture that captures stdout and stderr output streams -def test_example_q_widget(make_napari_viewer, capsys): - # make viewer and add an image layer using our fixture +def test_make_flim_phasor_plot(make_napari_viewer, capsys): viewer = make_napari_viewer() - viewer.add_image(np.random.random((100, 100))) - - # create our widget, passing in the viewer - my_widget = ExampleQWidget(viewer) - - # call our widget method - my_widget._on_click() - - # read captured output and check that it's as we expected - captured = capsys.readouterr() - assert captured.out == "napari has 1 layers\n" - - -def test_example_magic_widget(make_napari_viewer, capsys): - viewer = make_napari_viewer() - layer = viewer.add_image(np.random.random((100, 100))) + laser_frequency = 40 # MHz + amplitude = 1 + tau_list = [0.1, 0.2, 0.5, 1, 2, 5, 10, 25, 40] # ns + number_of_time_points = 1000 + + time_array = create_time_array(laser_frequency, number_of_time_points) + flim_data = make_synthetic_flim_data(time_array, amplitude, tau_list) + flim_data = flim_data.reshape(number_of_time_points, 3, 3) + + layer = viewer.add_image(flim_data, rgb = False) # this time, our widget will be a MagicFactory or FunctionGui instance - my_widget = example_magic_widget() + my_widget = make_flim_phasor_plot() # if we "call" this object, it'll execute our function - my_widget(viewer.layers[0]) - - # read captured output and check that it's as we expected - captured = capsys.readouterr() - assert captured.out == f"you have selected {layer}\n" + my_widget() + + labels_layer = viewer.layers[-1] + + assert len(viewer.layers) == 2 + assert labels_layer.features.shape == (9, 3) diff --git a/src/napari_flim_phasor_calculator/_widget.py b/src/napari_flim_phasor_calculator/_widget.py index 6199a07..6595608 100644 --- a/src/napari_flim_phasor_calculator/_widget.py +++ b/src/napari_flim_phasor_calculator/_widget.py @@ -9,78 +9,113 @@ from typing import TYPE_CHECKING from magicgui import magic_factory -from qtpy.QtWidgets import QHBoxLayout, QPushButton, QWidget -from napari.types import ImageData, LayerDataTuple +from napari.types import LayerDataTuple +from napari.layers import Image, Labels +from napari import Viewer import numpy as np import pandas as pd +from qtpy.QtWidgets import QWidget +# import napari_clusters_plotter + + +from .phasor import get_phasor_components +from .filters import make_time_mask, make_space_mask_from_manual_threshold +from .filters import apply_median_filter +from ._plotting import PhasorPlotterWidget if TYPE_CHECKING: import napari -# Uses the `autogenerate: true` flag in the plugin manifest -# to indicate it should be wrapped as a magicgui to autogenerate -# a widget. - -def get_phasor_components(flim_data, harmonic=1): +def connect_events(widget): ''' - Calculate phasor components G and S from the fourier transform + Connect widget events to make some visible/invisible depending on others ''' - import numpy as np - flim_data_fft = np.fft.fft(flim_data, axis=0) - dc = flim_data_fft[0].real - # change the zeros to the img average - # dc = np.where(dc != 0, dc, int(np.mean(dc))) - dc = np.where(dc != 0, dc, 1) - g = flim_data_fft[harmonic].real - g = g / dc - s = abs(flim_data_fft[harmonic].imag) - s /= dc - return g, s, dc + def toggle_median_n_widget(event): + widget.median_n.visible = event + # Connect events + widget.apply_median.changed.connect(toggle_median_n_widget) + # Intial visibility states + widget.median_n.visible = False + widget.laser_frequency.label = 'Laser Frequency (MHz)' -def create_time_array(frequency, n_points=100): - ''' - Create time array from laser frequency +@magic_factory(widget_init=connect_events) +def make_flim_phasor_plot(image_layer : Image, + laser_frequency : float = 40, + harmonic : int = 1, + threshold : int = 0, + apply_median : bool = False, + median_n : int = 1, + napari_viewer : Viewer = None) -> None: + from skimage.segmentation import relabel_sequential + image = image_layer.data + if 'TTResult_SyncRate' in image_layer.metadata: + laser_frequency = image_layer.metadata['TTResult_SyncRate'] *1E-6 #MHz - Parameters - ---------- - frequency: float - Frquency of the pulsed laser (in MHz) - n_points: int, optional - The number of samples collected between each aser shooting - Returns - ------- - time_array : array - Time array (in nanoseconds) - ''' - import numpy as np - time_window = 1 / (frequency * 10**6) - time_window_ns = time_window * 10**9 # in nanoseconds - time_step = time_window_ns / n_points # ns - array = np.arange(0, n_points) - time_array = array * time_step - return time_array - - -def make_flim_label_layer(image : ImageData, laser_frequency : float, harmonic : int = 1, threshold : int = 0, apply_median : bool = False, viewer=None) -> LayerDataTuple: + time_mask = make_time_mask(image, laser_frequency) + + space_mask = make_space_mask_from_manual_threshold(image, threshold) - # create time array based on laser frequency - time_array = create_time_array(laser_frequency, n_points = image.shape[0]) - time_step = time_array[1] - # choose starting index based on maximum value of image histogram - heights, bin_edges = np.histogram(np.ravel(np.argmax(image, axis=0) * time_step), bins=time_array) - start_index = np.argmax(heights[1:]) + 1 + image = image[time_mask] - time_mask = time_array >= time_array[start_index] + if apply_median: + image = apply_median_filter(image, median_n) - g, s, dc = get_phasor_components(image[time_mask], harmonic = harmonic) + g, s, dc = get_phasor_components(image, harmonic = harmonic) label_image = np.arange(dc.shape[0]*dc.shape[1]).reshape(dc.shape) + 1 + label_image[~space_mask] = 0 + label_image = relabel_sequential(label_image)[0] - phasor_components = {'label': np.ravel(label_image), 'G': np.ravel(g), 'S': np.ravel(s)} + phasor_components = {'label': np.ravel(label_image[space_mask]), + 'G': np.ravel(g[space_mask]), + 'S': np.ravel(s[space_mask])} table = pd.DataFrame(phasor_components) + # The layer has to be created here so the plotter can be filled properly + # below. Overwrite layer if it already exists. + for layer in napari_viewer.layers: + if (isinstance(layer, Labels)) & (layer.name=='Label_' + image_layer.name): + labels_layer = layer + labels_layer.data = label_image + labels_layer.features = table + break + else: + labels_layer = napari_viewer.add_labels(label_image, + name = 'Label_' + image_layer.name, + features = table) + + # Check if plotter was alrerady added to dock_widgets + dock_widgets_names = [key for key, value in napari_viewer.window._dock_widgets.items()] + if 'Plotter Widget' not in dock_widgets_names: + plotter_widget = PhasorPlotterWidget(napari_viewer) + napari_viewer.window.add_dock_widget(plotter_widget, name = 'Plotter Widget') + else: + widgets = napari_viewer.window._dock_widgets['Plotter Widget'] + plotter_widget = widgets.findChild(PhasorPlotterWidget) + # If we were to use the original plotter, we could add it as below + # _, plotter_widget = napari_viewer.window.add_plugin_dock_widget( + # 'napari-clusters-plotter', + # widget_name='Plotter Widget') + # Set G and S as features to plot (update_axes_list method clears Comboboxes) + plotter_widget.update_axes_list() + plotter_widget.plot_x_axis.setCurrentIndex(1) + plotter_widget.plot_y_axis.setCurrentIndex(2) + # Show parent (PlotterWidget) so that run function can run properly + plotter_widget.parent().show() + # Disconnect selector to reset collection of points in plotter + # (it gets reconnected when 'run' method is run) + plotter_widget.graphics_widget.selector.disconnect() + plotter_widget.run(labels_layer.features, + plotter_widget.plot_x_axis.currentText(), + plotter_widget.plot_y_axis.currentText()) - return (label_image, {'features' : table}, 'labels') + # Update laser frequency spinbox + # To Do: access and update widget in a better way + if 'Make FLIM Phasor Plot (napari-flim-phasor-calculator)' in dock_widgets_names: + widgets = napari_viewer.window._dock_widgets['Make FLIM Phasor Plot (napari-flim-phasor-calculator)'] + laser_frequency_spinbox = widgets.children()[4].children()[2].children()[-1] + laser_frequency_spinbox.setValue(laser_frequency) + + return + -def example_function_widget(img_layer: "napari.layers.Image"): - print(f"you have selected {img_layer}") diff --git a/src/napari_flim_phasor_calculator/filters.py b/src/napari_flim_phasor_calculator/filters.py new file mode 100644 index 0000000..52cc236 --- /dev/null +++ b/src/napari_flim_phasor_calculator/filters.py @@ -0,0 +1,60 @@ +import numpy as np + +def make_time_mask(image, laser_frequency): + ''' + Create a time mask from the image histogram maximum onwards + + Parameters + ---------- + image: array + The flim timelapse image + laser_frequency: float + Frequency of the pulsed laser (in MHz) + Returns + ------- + time_mask : boolean array + Time mask + ''' + from napari_flim_phasor_calculator._synthetic import create_time_array + # create time array based on laser frequency + time_array = create_time_array(laser_frequency, n_points = image.shape[0]) + time_step = time_array[1] + # choose starting index based on maximum value of image histogram + heights, bin_edges = np.histogram( + np.ravel(np.argmax(image, axis=0) * time_step), bins=time_array + ) + + start_index = np.argmax(heights[1:]) + 1 + time_mask = time_array >= time_array[start_index] + + return time_mask + +def make_space_mask_from_manual_threshold(image, threshold): + ''' + Create a space mask from the summed intensity image over time, keeping + pixels whose value is above threshold. + + Parameters + ---------- + image: array + The flim timelapse image + threshold: float + An integer threshold value. + Returns + ------- + space_mask : boolean array + A boolean mask representing pixels to keep. + ''' + intensity_image = np.sum(image, axis=0) + space_mask = intensity_image >= threshold + + return space_mask + +def apply_median_filter(image, n = 1): + from skimage.filters import median + from skimage.morphology import cube + footprint = cube(3) + image_filt = np.copy(image) + for i in range(n): + image_filt = median(image_filt, footprint) + return image_filt \ No newline at end of file diff --git a/src/napari_flim_phasor_calculator/napari.yaml b/src/napari_flim_phasor_calculator/napari.yaml index 9645d1a..29cc73d 100644 --- a/src/napari_flim_phasor_calculator/napari.yaml +++ b/src/napari_flim_phasor_calculator/napari.yaml @@ -9,17 +9,16 @@ contributions: python_name: napari_flim_phasor_calculator._sample_data:make_sample_data title: Load sample data from FLIM phasor calculator - id: napari-flim-phasor-calculator.make_func_widget - python_name: napari_flim_phasor_calculator._widget:make_flim_label_layer + python_name: napari_flim_phasor_calculator._widget:make_flim_phasor_plot title: Make FLIM Label Layer readers: - command: napari-flim-phasor-calculator.get_reader accepts_directories: false - filename_patterns: ['*.npy'] + filename_patterns: ['*.ptu'] sample_data: - command: napari-flim-phasor-calculator.make_sample_data display_name: FLIM phasor calculator key: unique_id.1 widgets: - command: napari-flim-phasor-calculator.make_func_widget - autogenerate: true - display_name: Make FLIM Label Layer + display_name: Make FLIM Phasor Plot diff --git a/src/napari_flim_phasor_calculator/phasor.py b/src/napari_flim_phasor_calculator/phasor.py new file mode 100644 index 0000000..7d7d6d8 --- /dev/null +++ b/src/napari_flim_phasor_calculator/phasor.py @@ -0,0 +1,43 @@ +''' +This function was adapted from PhasorPy v1.0.2 +(https://pypi.org/project/PhasorPy) + +MIT License + +Copyright (c) 2022 Bruno Schuty + +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. +''' + +def get_phasor_components(flim_data, harmonic=1): + ''' + Calculate phasor components G and S from the Fourier transform. + ''' + import numpy as np + flim_data_fft = np.fft.fft(flim_data, axis=0) + dc = flim_data_fft[0].real + # change the zeros to the img average + dc = np.where(dc != 0, dc, int(np.mean(dc))) + + g = flim_data_fft[harmonic].real + g /= dc + s = flim_data_fft[harmonic].imag + s /= -dc + + return g, s, dc \ No newline at end of file diff --git a/src/napari_flim_phasor_calculator/testing.py b/src/napari_flim_phasor_calculator/testing.py new file mode 100644 index 0000000..9fbd94d --- /dev/null +++ b/src/napari_flim_phasor_calculator/testing.py @@ -0,0 +1,36 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Dec 19 09:11:43 2022 + +@author: mazo260d +""" + +import napari +from napari_flim_phasor_calculator._widget import make_flim_phasor_plot +from napari_flim_phasor_calculator._synthetic import create_time_array, make_synthetic_flim_data +viewer = napari.Viewer() + +laser_frequency = 60 # MHz +amplitude = 1 +tau_list = [0.1, 0.2, 0.5, 1, 2, 5, 10, 25, 40] # ns +number_of_harmonics = 5 +number_of_time_points = 1000 + +time_array = create_time_array(laser_frequency, number_of_time_points) +flim_data = make_synthetic_flim_data(time_array, amplitude, tau_list) +flim_data = flim_data.reshape(number_of_time_points, 3, 3) + + +layer = viewer.add_image(flim_data, rgb = False)#, metadata={'TTResult_SyncRate': 39100000}) + +# this time, our widget will be a MagicFactory or FunctionGui instance +my_widget = make_flim_phasor_plot() + +viewer.window.add_plugin_dock_widget(plugin_name = 'napari-flim-phasor-calculator', + widget_name='Make FLIM Phasor Plot') + +# if we "call" this object, it'll execute our function +my_widget(layer, laser_frequency = laser_frequency, threshold = 30) + + +