Examples

If you don’t have PWS experimental data to run these examples with you can find a dataset that is used for automated testing on Zenodo here: https://zenodo.org/record/5976039#.Yf6aPt_MJPY

Running basic PWS analysis on a single image.

This example runs the PWS analysis on a single measurement. The results of the analysis will be saved alongside the raw data under the analyses folder. Running this example may be required in order for other examples which require that analysis results already be available.

 1# -*- coding: utf-8 -*-
 2# Copyright 2018-2021 Nick Anthony, Backman Biophotonics Lab, Northwestern University
 3#
 4# This file is part of PWSpy.
 5#
 6# PWSpy is free software: you can redistribute it and/or modify
 7# it under the terms of the GNU General Public License as published by
 8# the Free Software Foundation, either version 3 of the License, or
 9# (at your option) any later version.
10#
11# PWSpy is distributed in the hope that it will be useful,
12# but WITHOUT ANY WARRANTY; without even the implied warranty of
13# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14# GNU General Public License for more details.
15#
16# You should have received a copy of the GNU General Public License
17# along with PWSpy.  If not, see <https://www.gnu.org/licenses/>.
18
19"""
20This script saves pws analysis results to the test data. This must be run before many of the other examples will work.
21"""
22
23from pwspy import analysis
24from pwspy import dataTypes as pwsdt
25import pathlib as pl
26
27### User Variables ###
28PWSExperimentPath: pl.Path = ...  # Set this to a folder containing multiple "Cell{x}" acquisition folders. If you have downloaded the test dataset you can use `pl.Path(__file__).parent.parent / 'tests' / 'resources' / 'test_data' / 'sequencer'`
29######################
30
31settings = analysis.pws.PWSAnalysisSettings.loadDefaultSettings("Recommended")
32
33# Load our blank reference image which will be used for normalization.
34refAcq = pwsdt.Acquisition(PWSExperimentPath / 'Cell3')
35ref = refAcq.pws.toDataClass()
36
37anls = analysis.pws.PWSAnalysis(settings=settings, extraReflectance=None, ref=ref)  # Create a new analysis for the given reference image and analysis settings. The "ExtraReflection" calibration is ignored in this case.
38
39acq = pwsdt.Acquisition(PWSExperimentPath / "Cell1")  # Create an "Acquisition" object to handle operations for the data associated with a single acquisition
40cube = acq.pws.toDataClass()  # Request that the PWS metadata object load the full data.
41results, warnings = anls.run(cube)  # Run the pre-setup analysis on our data. Get the analysis results and potentially a list of warnings from the analysis.
42
43acq.pws.saveAnalysis(results, 'myAnalysis', overwrite=True)  # Save our analysis results to file in the default location alongside the raw data under the `analyses` folder.
44
45loadedResults = acq.pws.loadAnalysis('myAnalysis')  # This is just going to be a copy of `results`, loaded from file.
46
47print(f"Saved anlysis `myAnalysis` to {acq.filePath}")

Performing FFT to visualize the estimated depth of cell features

 1# -*- coding: utf-8 -*-
 2# Copyright 2018-2021 Nick Anthony, Backman Biophotonics Lab, Northwestern University
 3#
 4# This file is part of PWSpy.
 5#
 6# PWSpy is free software: you can redistribute it and/or modify
 7# it under the terms of the GNU General Public License as published by
 8# the Free Software Foundation, either version 3 of the License, or
 9# (at your option) any later version.
10#
11# PWSpy is distributed in the hope that it will be useful,
12# but WITHOUT ANY WARRANTY; without even the implied warranty of
13# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14# GNU General Public License for more details.
15#
16# You should have received a copy of the GNU General Public License
17# along with PWSpy.  If not, see <https://www.gnu.org/licenses/>.
18
19"""
20Load the OPD from a previously saved analysis result and plot it using a special multi-dimensional plotting widget.
21
22@author: Nick Anthony
23"""
24
25import pwspy.dataTypes as pwsdt
26from mpl_qt_viz.visualizers import PlotNd  # The mpl_qt_viz library is not a dependency of PWSpy and will need to be installed separately. It can be found on Conda-Forge (conda install -c conda-forge mpl_qt_viz) or on PyPi (pip install mpl_qt_viz)
27import matplotlib.pyplot as plt
28import pathlib as pl
29from PyQt5.QtWidgets import QApplication
30
31### User Variables ###
32PWSImagePath: pl.Path = ... # Set this to the path of a "Cell{X}" folder of a PWS acquisition.
33######################
34
35plt.ion()  # Without this we will get a crash when trying to open the PlotNd window because a Qt application loop must be running.
36plt.figure()
37
38acquisition = pwsdt.Acquisition(PWSImagePath)  # Get a reference to the top-level folder for the measurement.
39
40roiSpecs = acquisition.getRois()  # Get a list of the (name, number, file format) of the available ROIs.
41print("ROIs:\n", roiSpecs)
42
43analysis = acquisition.pws.loadAnalysis(acquisition.pws.getAnalyses()[0])  # Load a reference to an analysis file.
44kCube = analysis.reflectance  # Load the processed `reflectance` array from the analysis file.
45
46opd, opdValues = kCube.getOpd(useHannWindow=False, indexOpdStop=50)  # Use FFT to transform the reflectance array to OPD
47
48# Scale the opdValues to give estimated depth instead of raw OPD. Factor of 2 because light is making a round-trip.
49ri = 1.37  # Estimated RI of livecell chromatin
50opdValues = opdValues / (2 * ri)
51
52app = QApplication([])  # Start a QApplication instance to run the PlotNd qwidget
53plotWindow = PlotNd(opd, names=('y', 'x', 'depth'),
54                    indices=(None, None, opdValues), title="Estimated Depth")
55app.exec()
_images/opdExample.gif

Compile analysis results to a table of average values within each ROI

 1# Copyright 2018-2021 Nick Anthony, Backman Biophotonics Lab, Northwestern University
 2#
 3# This file is part of PWSpy.
 4#
 5# PWSpy is free software: you can redistribute it and/or modify
 6# it under the terms of the GNU General Public License as published by
 7# the Free Software Foundation, either version 3 of the License, or
 8# (at your option) any later version.
 9#
10# PWSpy is distributed in the hope that it will be useful,
11# but WITHOUT ANY WARRANTY; without even the implied warranty of
12# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13# GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License
16# along with PWSpy.  If not, see <https://www.gnu.org/licenses/>.
17
18"""
19This script loads analysis results from a list of PWS acquisitions and uses the "PWSRoiCompiler" class to get average output
20values within the ROIs for each acquisition. The compiled results are then placed into a Pandas dataframe.
21"""
22import pandas
23from pwspy.analysis.compilation import PWSRoiCompiler, PWSCompilerSettings
24import pwspy.dataTypes as pwsdt
25import pathlib as pl
26
27### User Variables ###
28PWSExperimentPath: pl.Path = ...  # Set this to a folder containing multiple "Cell{x}" acquisition folders. If you have downloaded the test dataset you can use `pl.Path(__file__).parent.parent / 'tests' / 'resources' / 'test_data' / 'sequencer'`
29######################
30
31
32tmpList = []
33compiler = PWSRoiCompiler(PWSCompilerSettings(reflectance=True, rms=True))
34
35listOfAcquisitions = [pwsdt.Acquisition(i) for i in PWSExperimentPath.glob("Cell[0-9]")]
36for acquisition in listOfAcquisitions:
37    for analysisName in acquisition.pws.getAnalyses():
38        analysisResults = acquisition.pws.loadAnalysis(analysisName)
39        for roiSpec in acquisition.getRois():
40            roiFile = acquisition.loadRoi(*roiSpec)
41            results, warnings = compiler.run(analysisResults, roiFile.getRoi())
42
43            if len(warnings) > 0:
44                print(warnings)
45
46            tmpList.append(dict(
47                acquisition=acquisition,
48                cellNumber=acquisition.getNumber(),
49                analysisResults=analysisResults,
50                rms=results.rms,
51                reflectance=results.reflectance,
52                roiNum=roiFile.number,
53                roiName=roiFile.name
54            ))
55
56dataFrame = pandas.DataFrame(tmpList)
57print(dataFrame)

Basic loading of ROI’s to extract data from specific regions

 1# -*- coding: utf-8 -*-
 2# Copyright 2018-2021 Nick Anthony, Backman Biophotonics Lab, Northwestern University
 3#
 4# This file is part of PWSpy.
 5#
 6# PWSpy is free software: you can redistribute it and/or modify
 7# it under the terms of the GNU General Public License as published by
 8# the Free Software Foundation, either version 3 of the License, or
 9# (at your option) any later version.
10#
11# PWSpy is distributed in the hope that it will be useful,
12# but WITHOUT ANY WARRANTY; without even the implied warranty of
13# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14# GNU General Public License for more details.
15#
16# You should have received a copy of the GNU General Public License
17# along with PWSpy.  If not, see <https://www.gnu.org/licenses/>.
18
19"""
20Loop through all ROIs for all acquisitions in a directory and plot a histogram of the RMS values within the ROI.
21"""
22
23import pwspy.dataTypes as pwsdt
24import matplotlib.pyplot as plt
25import pathlib
26import numpy as np
27
28
29### User Variables ###
30PWSExperimentPath = ...  # Set this to a folder containing multiple "Cell{x}" acquisition folders. If you have downloaded the test dataset you can use `pl.Path(__file__).parent.parent / 'tests' / 'resources' / 'test_data' / 'sequencer'`
31analysisName = 'script'  # This will often be "p0"
32######################
33
34plt.ion()
35
36
37def plotHist(roi, rms):
38    """
39    This function takes an ROI
40    and a 2D RMS image and plots a histogram of the RMS values inside the ROI
41    """
42    # Check input values just to be safe.
43    assert isinstance(roi, pwsdt.Roi)  # Make sure roiFile variable is actually an ROI
44    assert isinstance(rms, np.ndarray)  # Make sure the RMS image is an numpy array
45    assert roi.mask.shape == rms.shape  # Make sure the ROI and RMS arrays have the same dimensions.
46
47    vals = rms[roi.mask]  # A 1D array of the values inside the ROI
48    plt.hist(vals)  # Plot a histogram
49
50
51cellFolderIterator = pathlib.Path(PWSExperimentPath).glob("Cell[0-9]")  # An iterator for all folders that are below workingDirectory and match the "regex" pattern "Cell[0-9]"
52for folder in cellFolderIterator:
53    acq = pwsdt.Acquisition(folder)  # An object handling the contents of a single "Cell{X}" folder
54
55    try:
56        anls = acq.pws.loadAnalysis(analysisName)  # Load the analysis results from file.
57    except:
58        print(f"Analysis loading failed for {acq.filePath}")
59        continue  # Skip to the next loop iteration
60
61    roiSpecs = acq.getRois()  # A list of the names, numbers, and fileFormats of the ROIs in this acquisition
62
63    for name, number, fformat in roiSpecs:  # Loop through every ROI.
64        roiFile = acq.loadRoi(name, number, fformat)  # Load the ROI from file.
65        plotHist(roiFile.getRoi(), anls.rms)  # Use the function defined above to plot a histogram

Using a hand-drawn ROI to generate a reference pseudo-measurement

 1
 2# -*- coding: utf-8 -*-
 3# Copyright 2018-2021 Nick Anthony, Backman Biophotonics Lab, Northwestern University
 4#
 5# This file is part of PWSpy.
 6#
 7# PWSpy is free software: you can redistribute it and/or modify
 8# it under the terms of the GNU General Public License as published by
 9# the Free Software Foundation, either version 3 of the License, or
10# (at your option) any later version.
11#
12# PWSpy is distributed in the hope that it will be useful,
13# but WITHOUT ANY WARRANTY; without even the implied warranty of
14# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15# GNU General Public License for more details.
16#
17# You should have received a copy of the GNU General Public License
18# along with PWSpy.  If not, see <https://www.gnu.org/licenses/>.
19
20"""
21This script allows the user to select a region of an PwsCube. the spectra of this
22region is then averaged over the X and Y dimensions. This spectra is then saved
23as a reference dataTypes with the same initial dimensions.
24Can help to make a reference when you don't actually have one for some reason
25"""
26
27import pwspy.dataTypes as pwsdt
28import matplotlib.pyplot as plt
29import numpy as np
30import pathlib as pl
31
32### User Variables ###
33PWSImagePath: pl.Path = ... # Set this to the path of a "Cell{X}" folder of a PWS acquisition.
34######################
35
36plt.ion()
37a = pwsdt.Acquisition(PWSImagePath).pws.toDataClass()  # Load a measurement from file.
38
39roi = a.selectLassoROI()  # Prompt the user for a hand-drawn ROI
40spec, std = a.getMeanSpectra(mask=roi)  # Get the average spectra within the ROI
41newData = np.zeros(a.data.shape)
42newData[:, :, :] = spec[np.newaxis, np.newaxis, :]  # Extend the averaged spectrum along the full dimensions of the original measurement.
43ref = pwsdt.PwsCube(newData, a.metadata)  # Create a new synthetic measurement using the averaged spectrum
44plt.plot(a.wavelengths, spec)
45

Blurring data laterally to smooth a reference image.

 1# -*- coding: utf-8 -*-
 2# Copyright 2018-2021 Nick Anthony, Backman Biophotonics Lab, Northwestern University
 3#
 4# This file is part of PWSpy.
 5#
 6# PWSpy is free software: you can redistribute it and/or modify
 7# it under the terms of the GNU General Public License as published by
 8# the Free Software Foundation, either version 3 of the License, or
 9# (at your option) any later version.
10#
11# PWSpy is distributed in the hope that it will be useful,
12# but WITHOUT ANY WARRANTY; without even the implied warranty of
13# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14# GNU General Public License for more details.
15#
16# You should have received a copy of the GNU General Public License
17# along with PWSpy.  If not, see <https://www.gnu.org/licenses/>.
18
19"""
20This script blurs an image cube in the xy direction. Allows you to turn an
21image of cells into something that can be used as a reference image, assuming
22most of the the FOV is glass. In reality you should just have a good reference image to use and not resort to something
23like this.
24"""
25
26import copy
27import matplotlib.pyplot as plt
28import pwspy.dataTypes as pwsdt
29
30### User Variables ###
31PWSImagePath = ...  # Set this to the path of a "Cell{X}" folder of a PWS acquisition.
32######################
33
34plt.ion()
35
36acq = pwsdt.Acquisition(PWSImagePath)
37a = acq.pws.toDataClass()
38
39a.correctCameraEffects()  # Correct for dark counts and potentially for camera nonlinearity using metadata stored with the original measurement.
40a.normalizeByExposure()  # Divide by exposure time to get data in units of `counts/ms`. This isn't strictly necessary in this case.
41
42mirror = copy.deepcopy(a)
43mirror.filterDust(10)  # Apply a gaussian blurring with sigma=10 microns along the XY plane.
44
45a.plotMean()  # Plot the mean reflectance of the original
46mirror.plotMean()  # Plot the mean reflectance after filtering.
47a.normalizeByReference(mirror)  # Normalize raw by reference
48a.plotMean()  # Plot the measurement after normalization.
49plt.figure()
50plt.imshow(a.data.std(axis=2))  # Plot RMS after normalization.

Measuring Sigma using only a limited range of the OPD signal.

  1
  2# -*- coding: utf-8 -*-
  3# Copyright 2018-2021 Nick Anthony, Backman Biophotonics Lab, Northwestern University
  4#
  5# This file is part of PWSpy.
  6#
  7# PWSpy is free software: you can redistribute it and/or modify
  8# it under the terms of the GNU General Public License as published by
  9# the Free Software Foundation, either version 3 of the License, or
 10# (at your option) any later version.
 11#
 12# PWSpy is distributed in the hope that it will be useful,
 13# but WITHOUT ANY WARRANTY; without even the implied warranty of
 14# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 15# GNU General Public License for more details.
 16#
 17# You should have received a copy of the GNU General Public License
 18# along with PWSpy.  If not, see <https://www.gnu.org/licenses/>.
 19
 20"""
 21This script is based on a matlab script written by Lusik Cherkezyan for NC.
 22Nano uses this method to extract rms from phantom make from ChromEM cells embedded in resin.
 23The phantom has a strong thin-film spectrum. This script is meant to filter out the thin film components
 24of the fourier transfrom and extract RMS from what is left.
 25"""
 26from pwspy.dataTypes import CameraCorrection, Acquisition, PwsCube, KCube
 27import matplotlib.pyplot as plt
 28import scipy.signal as sps
 29import os
 30import numpy as np
 31import pathlib as pl
 32
 33'''User Input'''
 34path: pl.Path = ...
 35refName = 'Cell3'  # This is an PwsCube of glass, used for normalization.
 36cellNames = ['Cell1', 'Cell2']  # , 'Cell3', 'Cell4','Cell5']
 37maskSuffix = 'resin'
 38
 39# identify the depth in um to which the OPD spectra need to be integrated
 40integrationDepth = 2.0  # in um
 41isHannWindow = True  # Should Hann windowing be applied to eliminate edge artifacts?
 42subtractResinOpd = True
 43resetResinMasks = False
 44wvStart = 510  # start wavelength for poly subtraction
 45wvEnd = 690  # end wavelength for poly subtraction
 46sampleRI = 1.545  # The refractive index of the resin. This is taken from matlab code, I don't know if it's correct.
 47orderPolyFit = 0
 48wv_step = 2
 49correction = CameraCorrection(2000, (0.977241216, 1.73E-06, 1.70E-11))
 50
 51'''************'''
 52
 53b, a = sps.butter(6, 0.1 * wv_step)
 54opdIntegralEnd = integrationDepth * 2 * sampleRI  # We need to convert from our desired depth into an opd value. There are some questions about having a 2 here but that's how it is in the matlab code so I'm keeping it.
 55
 56### load and save mirror or glass image cube
 57ref = PwsCube.fromMetadata(Acquisition(os.path.join(path, refName)).pws)
 58ref.correctCameraEffects(correction)
 59ref.filterDust(6, pixelSize=1)
 60ref.normalizeByExposure()
 61
 62if subtractResinOpd:
 63    ### load and save reference empty resin image cube
 64    fig, ax = plt.subplots()
 65    resinOpds = {}
 66    for cellName in cellNames:
 67        resin = PwsCube.fromMetadata(Acquisition(os.path.join(path, cellName)).pws)
 68        resin.correctCameraEffects(correction)
 69        resin.normalizeByExposure()
 70        resin /= ref
 71        resin = KCube.fromPwsCube(resin)
 72        if resetResinMasks:
 73            [resin.metadata.acquisitionDirectory.deleteRoi(name, num) for name, num, fformat in resin.metadata.acquisitionDirectory.getRois() if name == maskSuffix]
 74        if maskSuffix in [name for name, number, fformat in resin.metadata.acquisitionDirectory.getRois()]:
 75            resinRoi = resin.metadata.acquisitionDirectory.loadRoi(maskSuffix, 1).getRoi()
 76        else:
 77            print('Select a region containing only resin.')
 78            resinRoi = resin.selectLassoROI()
 79            resin.metadata.acquisitionDirectory.saveRoi(maskSuffix, 1, resinRoi)
 80        resin.data -= resin.data.mean(axis=2)[:, :, np.newaxis]
 81        opdResin, xvals = resin.getOpd(isHannWindow, indexOpdStop=None, mask=resinRoi.mask)
 82        resinOpds[cellName] = opdResin
 83        ax.plot(xvals, opdResin, label=cellName)
 84        ax.vlines([opdIntegralEnd], ymin=opdResin.min(), ymax=opdResin.max())
 85        ax.set_xlabel('OPD')
 86        ax.set_ylabel("Amplitude")
 87    ax.legend()
 88    plt.pause(0.2)
 89
 90print("Beginning processing.")
 91rmses = {}  # Store the rms maps for later saving
 92for cellName in cellNames:
 93    cube = PwsCube.fromMetadata(Acquisition(os.path.join(path, cellName)).pws)
 94    cube.correctCameraEffects(correction)
 95    cube.normalizeByExposure()
 96    cube /= ref
 97    cube.data = sps.filtfilt(b, a, cube.data, axis=2)
 98    cube = KCube.fromPwsCube(cube)
 99
100    ## -- Polynomial Fit
101    print("Subtracting Polynomial")
102    polydata = cube.data.reshape((cube.data.shape[0] * cube.data.shape[1], cube.data.shape[2]))
103    polydata = np.rollaxis(polydata, 1)  # Flatten the array to 2d and put the wavenumber axis first.
104    cubePoly = np.zeros(polydata.shape)  # make an empty array to hold the fit values.
105    polydata = np.polyfit(cube.wavenumbers, polydata,
106                          orderPolyFit)  # At this point polydata goes from holding the cube data to holding the polynomial values for each pixel. still 2d.
107    for i in range(orderPolyFit + 1):
108        cubePoly += (np.array(cube.wavenumbers)[:, np.newaxis] ** i) * polydata[i,
109                                                                       :]  # Populate cubePoly with the fit values.
110    cubePoly = np.moveaxis(cubePoly, 0, 1)
111    cubePoly = cubePoly.reshape(cube.data.shape)  # reshape back to a cube.
112    # Remove the polynomial fit from filtered cubeCell.
113    cube.data = cube.data - cubePoly
114
115    rmsData = np.sqrt(np.mean(cube.data ** 2, axis=2)) #This can be compared to rmsOPDIntData, when the integralStopIdx is high and we don't do spectral subtraction they should be equivalent.
116
117    # Find the fft for each signal in the desired wavelength range
118    opdData, xvals = cube.getOpd(isHannWindow, None)
119
120    if subtractResinOpd:
121        opdData = opdData - resinOpds[cellName]
122
123    try:
124        integralStopIdx = np.where(xvals >= opdIntegralEnd)[0][0]
125    except IndexError:  # If we get an index error here then our opdIntegralEnd is probably bigger than we can achieve. Just use the biggest value we have.
126        integralStopIdx = None
127        opdIntegralEnd = max(xvals)
128        print(f'Integrating to OPD {opdIntegralEnd}')
129
130    opdSquared = np.sum(opdData[:, :, :integralStopIdx] ** 2,axis=2)  # Parseval's theorem tells us that this is equivalent to the sum of the squares of our original signal
131    opdSquared *= len(cube.wavenumbers) / opdData.shape[2]  # If the original data and opd were of the same length then the above line would be correct. Since the fft has been upsampled. we need to normalize.
132    rmsOpdIntData = np.sqrt(opdSquared)  # this should be equivalent to normal RMS if our stop index is high and resin subtraction is disabled.
133
134    cmap = plt.get_cmap('jet')
135    fig, axs = plt.subplots(1, 2, sharex=True, sharey=True)
136    im = axs[0].imshow(rmsData, cmap=cmap, clim=[np.percentile(rmsData, 0.5), np.percentile(rmsData, 99.5)])
137    fig.colorbar(im, ax=axs[0])
138    axs[0].set_title('RMS')
139    im = axs[1].imshow(rmsOpdIntData, cmap=cmap,
140                       clim=[np.percentile(rmsOpdIntData, 0.5), np.percentile(rmsOpdIntData, 99.5)])
141    fig.colorbar(im, ax=axs[1])
142    axs[1].set_title(f'RMS from OPD below {opdIntegralEnd} after resin OPD subtraction')
143    fig.suptitle(cellName)
144    rmses[cellName] = rmsOpdIntData
145    plt.pause(0.2)
146plt.pause(0.5)

Generating new position lists to enable colocalized measurements on multiple systems.

 1# Copyright 2018-2020 Nick Anthony, Backman Biophotonics Lab, Northwestern University
 2#
 3# This file is part of PWSpy.
 4#
 5# PWSpy is free software: you can redistribute it and/or modify
 6# it under the terms of the GNU General Public License as published by
 7# the Free Software Foundation, either version 3 of the License, or
 8# (at your option) any later version.
 9#
10# PWSpy is distributed in the hope that it will be useful,
11# but WITHOUT ANY WARRANTY; without even the implied warranty of
12# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13# GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License
16# along with PWSpy.  If not, see <https://www.gnu.org/licenses/>.
17from pwspy.utility.micromanager import PropertyMap
18
19from pwspy.utility.micromanager.positions import PositionList
20
21"""This example demonstrates how to use generate new cell positions from a set of positions after the sample has been picked up and likely shifted or rotated.
22This method relies on measuring a set (at least 3) of reference positions before and after moving the dish. You can then use these positions to generate an 
23affine transform. This affine transform can then be applied to your original cell positions in order to generate a new set of positions for the same cells.
24In the case of a standard cell culture dish it is best to use the corners of the glass coverslip as your reference locations.
25
26@author: Nick Anthony
27"""
28
29import pathlib as pl
30import matplotlib.pyplot as plt
31from skimage.transform import AffineTransform
32import numpy as np
33
34plt.ion()
35NCPath = pl.Path(r'Z:\Nick\NU-NC_EtOH fixation comparison\NC\Buccal cell')
36NUPath = pl.Path(r'Z:\Nick\NU-NC_EtOH fixation comparison\NU\Buccal')
37
38# Load the position list of the coverslip corners taken at the beginning of the experiment.
39preTreatRefPositions = PositionList.fromNanoMatFile(NCPath / 'corners_list2.mat', 'TIXYDrive')
40# Load the position list of the coverslip corners after placing the dish back on the microscope after treatment.
41postTreatRefPositions = PositionList.fromPropertyMap(PropertyMap.loadFromFile(NUPath / 'corners.pos'))
42# Generate an affine transform describing the difference between the two position lists.
43transformMatrix = preTreatRefPositions.getAffineTransform(postTreatRefPositions)
44# Load the positions of the cells we are measuring before the dish was removed.
45preTreatCellPositions = PositionList.fromNanoMatFile(NCPath / 'cell_list.mat', 'TIXYDrive')
46# Transform the cell positions to the new expected locations.
47postTreatCellPositions = preTreatCellPositions.applyAffineTransform(transformMatrix)
48# Save the new positions to a file that can be loaded by Micro-Manager.
49postTreatCellPositions.toPropertyMap().saveToFile(NUPath / 'transformedPositions.pos')
50
51# Plot the reference and cell positions before treatment
52fig, ax = plt.subplots()
53preTreatRefPositions.plot(fig, ax)
54preTreatCellPositions.plot(fig, ax)
55
56# Plot the reference and cell positions after treatment
57fig2, ax2 = plt.subplots()
58postTreatRefPositions.plot(fig2, ax2)
59postTreatCellPositions.plot(fig2, ax2)
60
61af = AffineTransform(np.vstack([transformMatrix, [0, 0, 1]]))
62print(f"Scale: {af.scale}, Shear: {af.shear}, Rotation: {af.rotation / (2*np.pi) * 360}, Translation: {af.translation}")