Merge pull request #1 from cdtools-developers/signal_estimate

Add the tools to quickly estimate the signal level from raw data
This commit is contained in:
allevitan
2024-09-20 13:44:51 +02:00
committed by GitHub
4 changed files with 535 additions and 43 deletions
+38 -5
View File
@@ -7,6 +7,8 @@ from cdtools.datasets import CDataset
from cdtools.datasets.random_selection import random_selection
from cdtools.tools import data as cdtdata
from cdtools.tools import plotting
from matplotlib import pyplot as plt
from cdtools.tools import analysis
from copy import deepcopy
__all__ = ['Ptycho2DDataset']
@@ -207,7 +209,13 @@ class Ptycho2DDataset(CDataset):
cdtdata.add_shot_to_shot_info(cxi_file, self.intensities, 'intensities')
def inspect(self, logarithmic=True, units='um', log_offset=1):
def inspect(
self,
logarithmic=True,
units='um',
log_offset=1,
plot_mean_pattern=True
):
"""Launches an interactive plot for perusing the data
This launches an interactive plotting tool in matplotlib that
@@ -248,14 +256,39 @@ class Ptycho2DDataset(CDataset):
# nanomap_values = (self.mask * self.patterns).sum(dim=(1,2)).detach().cpu().numpy()
if logarithmic:
cbar_title = ('Log Base 10 of Diffraction Intensity + %0.2f'
% log_offset)
cbar_title = f'Log Base 10 of Intensity + {log_offset}'
else:
cbar_title = 'Diffraction Intensity'
cbar_title = 'Intensity'
if plot_mean_pattern:
self.plot_mean_pattern(log_offset=log_offset)
return plotting.plot_nanomap_with_images(self.translations.detach().cpu(), get_images, values=nanomap_values, nanomap_units=units, image_title='Diffraction Pattern', image_colorbar_title=cbar_title)
def plot_mean_pattern(self, log_offset=1):
"""Plots the mean diffraction pattern across the dataset
The output is normalized so that the summed intensity on the
detector is equal to the total intensity of light that passed
through the sample within each detector conjugate field of view.
The plot is plotted as log base 10 of the output plus log_offset.
By default, log_offset is set equal to 1, which is a good level for
shot-noise limited data captured in units of photons. More
generally, log_offset should be set roughly at the background noise
level.
"""
mean_pattern, bins, ssnr = analysis.calc_spectral_info(self)
cmap_label = f'Log Base 10 of Intensity + {log_offset}'
title = 'Scaled mean diffraction pattern'
return plotting.plot_real(
t.log10(t.as_tensor(mean_pattern + log_offset)),
cmap_label=cmap_label,
title=title,
)
def split(self):
"""Splits a dataset into two pseudorandomly selected sub-datasets
"""
+29 -6
View File
@@ -25,6 +25,7 @@ class FancyPtycho(CDIModel):
background=None,
probe_basis=None,
translation_offsets=None,
probe_fourier_shifts=None,
mask=None,
weights=None,
translation_scale=1,
@@ -134,6 +135,13 @@ class FancyPtycho(CDIModel):
t_o = t.as_tensor(translation_offsets, dtype=t.float32)
t_o = t_o / translation_scale
self.translation_offsets = t.nn.Parameter(t_o)
if probe_fourier_shifts is None:
self.probe_fourier_shifts = None
else:
self.probe_fourier_shifts = t.nn.Parameter(
t.as_tensor(translation_offsets, dtype=t.float32)
)
self.register_buffer('translation_scale',
t.as_tensor(translation_scale, dtype=dtype))
@@ -152,7 +160,7 @@ class FancyPtycho(CDIModel):
t.as_tensor(simulate_probe_translation, dtype=bool)
)
if simulate_probe_translation:
if simulate_probe_translation or (self.probe_fourier_shifts is not None):
Is = t.arange(self.probe.shape[-2], dtype=dtype)
Js = t.arange(self.probe.shape[-1], dtype=dtype)
Is, Js = t.meshgrid(Is/t.max(Is), Js/t.max(Js))
@@ -195,6 +203,7 @@ class FancyPtycho(CDIModel):
fourier_probe=False,
loss='amplitude mse',
units='um',
allow_probe_fourier_shifts=False,
simulate_probe_translation=False,
simulate_finite_pixels=False,
exponentiate_obj=False,
@@ -327,6 +336,11 @@ class FancyPtycho(CDIModel):
translation_offsets = 0 * (t.rand((len(dataset), 2)) - 0.5)
if allow_probe_fourier_shifts:
probe_fourier_shifts = t.zeros((len(dataset), 2), dtype=t.float32)
else:
probe_fourier_shifts = None
if dm_rank is not None and dm_rank != 0:
if dm_rank > n_modes:
raise KeyError('Density matrix rank cannot be greater than the number of modes. Use dm_rank = -1 to use a full rank matrix.')
@@ -347,8 +361,9 @@ class FancyPtycho(CDIModel):
Ws = t.ones(len(dataset))
if hasattr(dataset, 'intensities') and dataset.intensities is not None:
Ws *= (dataset.intensities.to(dtype=Ws.dtype)[:,...]
/ t.mean(dataset.intensities))
intensities = dataset.intensities.to(dtype=Ws.dtype)[:,...]
weights = t.sqrt(intensities)
Ws *= (weights / t.mean(weights))
if hasattr(dataset, 'mask') and dataset.mask is not None:
mask = dataset.mask.to(t.bool)
@@ -380,6 +395,7 @@ class FancyPtycho(CDIModel):
fourier_probe=fourier_probe,
oversampling=oversampling,
loss=loss, units=units,
probe_fourier_shifts=probe_fourier_shifts,
simulate_probe_translation=simulate_probe_translation,
simulate_finite_pixels=simulate_finite_pixels,
phase_only=phase_only,
@@ -431,12 +447,19 @@ class FancyPtycho(CDIModel):
# Maybe this can be done with a matmul now?
prs = t.sum(Ws[..., None, None] * basis_prs, axis=-3)
if self.simulate_probe_translation:
det_pix_trans = tools.interactions.translations_to_pixel(
if self.simulate_probe_translation or (self.probe_fourier_shifts is not None):
if self.probe_fourier_shifts is not None:
det_pix_trans = self.probe_fourier_shifts[index]
else:
det_pix_trans = t.zeros_like(translations)
if self.simulate_probe_translation:
det_pix_trans = det_pix_trans + tools.interactions.translations_to_pixel(
self.det_basis,
translations,
surface_normal=self.surface_normal)
probe_masks = t.exp(1j* (det_pix_trans[:,0,None,None] *
self.I_phase[None,...] +
det_pix_trans[:,1,None,None] *
+98 -9
View File
@@ -9,9 +9,11 @@ data has been stored in numpy arrays.
import torch as t
import numpy as np
from cdtools.tools import image_processing as ip
import cdtools
from scipy import linalg as sla
from scipy import special
from scipy import optimize as opt
from scipy import spatial
__all__ = [
'product_svd',
@@ -30,6 +32,7 @@ __all__ = [
'remove_amplitude_exponent',
'standardize_reconstruction_set',
'standardize_reconstruction_pair',
'calc_spectral_info',
]
@@ -1338,16 +1341,15 @@ def standardize_reconstruction_pair(
obj_2, probe_2 = remove_phase_ramp(
half_2['obj'], window, probe=half_2['probe'])
# TODO weights are not included
if correct_amplitude_exponent:
obj_1, probe_1, weights_1 = remove_amplitude_exponent(
obj_1, probe_1 = remove_amplitude_exponent(
obj_1, window, probe=probe_1,
weights=half_1['weights'],
basis=half_1['basis'],
basis=half_1['obj_basis'],
translations=half_1['translations'])
obj_2, probe_2, weights_2 = remove_amplitude_exponent(
obj_2, probe_2 = remove_amplitude_exponent(
obj_2, window, probe=probe_2,
weights=half_2['weights'],
basis=half_2['basis'],
basis=half_2['obj_basis'],
translations=half_2['translations'])
@@ -1355,7 +1357,6 @@ def standardize_reconstruction_pair(
obj_1 = np.exp(-1j* np.angle(np.sum(obj_1[window]))) * obj_1
obj_2 = np.exp(-1j* np.angle(np.sum(obj_2[window]))) * obj_2
# Todo update the translations to account for the determined shift
shift = ip.find_shift(
t.as_tensor(ip.hann_window(obj_1[window])),
@@ -1367,6 +1368,7 @@ def standardize_reconstruction_pair(
t.as_tensor(probe_1[0]),
t.as_tensor(probe_2[0]),
)
for idx in range(probe_2.shape[0]):
probe_2[idx] = ip.sinc_subpixel_shift(
t.as_tensor(probe_2[idx]), probe_shift).numpy()
@@ -1388,9 +1390,16 @@ def standardize_reconstruction_pair(
limit=frc_limit,
)
probe_1_intensity = np.sum(np.abs(probe_1)**2)
probe_2_intensity = np.sum(np.abs(probe_2)**2)
probe_nmse = 1 - (calc_fidelity(probe_1, probe_2)
/ (probe_1_intensity * probe_2_intensity))
probe_nrms_error = calc_generalized_rms_error(
probe_1,
probe_2,
probe_1[0:],
probe_2[0:],
normalize=True
)
@@ -1420,7 +1429,87 @@ def standardize_reconstruction_pair(
'probe_frc': probe_frc,
'probe_frc_threshold': probe_frc_threshold,
'probe_nrms_error': probe_nrms_error,
'probe_nmse': probe_nmse,
}
return results
def calc_spectral_info(dataset, nbins=50):
"""Makes a properly normalized sum diffraction pattern
This returns a scaled version of sum of all the diffraction patterns
within the dataset. The scaling is defined so that the total intensity
in the final image is equal to the intensity arising from a region of
the scan pattern whose area matches one detector conjugate field of
view.
Parameters
----------
dataset : Ptycho2DDataset
A ptychography dataset to use
nbins : int
The number of bins to use for the SNR curve
Returns
-------
spectrum : t.tensor
An image of the spectral signal rate
freqs : t.tensor
The frequencies at which the SSNR is estimated
SSNR : t.tensor
The estimated SSNR
"""
scan_hull = spatial.ConvexHull(dataset.translations[:,:2].cpu().numpy())
scan_area = scan_hull.volume
ewg = cdtools.tools.initializers.exit_wave_geometry
obj_basis = ewg(
dataset.detector_geometry['basis'],
dataset[0][1].shape,
dataset.wavelength,
dataset.detector_geometry['distance'],
)
det_conj_fov_area = np.linalg.norm(
np.cross(obj_basis[:,0]*dataset.patterns.shape[-2],
obj_basis[:,1]*dataset.patterns.shape[-1])
)
scale_factor = det_conj_fov_area / scan_area
mask = dataset.mask.cpu().numpy().astype(int)
sum_pattern = dataset.mask * t.sum(dataset.patterns, dim=0) * scale_factor
sum_pattern = sum_pattern.cpu().numpy()
# TODO this assumes orthogonal axes
pix_sizes = np.linalg.norm(obj_basis, axis=0)
i_freqs = np.fft.fftshift(np.fft.fftfreq(
sum_pattern.shape[0],d=pix_sizes[0]))
j_freqs = np.fft.fftshift(np.fft.fftfreq(
sum_pattern.shape[1],d=pix_sizes[1]))
Js,Is = np.meshgrid(j_freqs,i_freqs)
Rs = np.sqrt(Is**2+Js**2)
max_i = np.max(i_freqs)
max_j = np.max(j_freqs)
frc_range = [0, max(max_i,max_j)]
sum_spectrum, frc_bins = np.histogram(Rs, bins=nbins, range=frc_range,
weights=sum_pattern)
sum_spectrum_sq, frc_bins = np.histogram(Rs, bins=nbins, range=frc_range,
weights=sum_pattern**2)
n_pix, frc_bins = np.histogram(Rs, bins=nbins, range=frc_range,
weights=mask)
mean_spectrum = sum_spectrum / n_pix
pattern_snr = sum_spectrum_sq / sum_spectrum
return sum_pattern, frc_bins[:-1], mean_spectrum
+370 -23
View File
@@ -13,12 +13,22 @@ from matplotlib.colors import hsv_to_rgb
from matplotlib.widgets import Slider
from matplotlib import ticker, patheffects
from matplotlib import transforms as mtransforms
from matplotlib import colors
__all__ = ['colorize', 'plot_amplitude', 'plot_phase',
'plot_colorized', 'plot_translations', 'get_units_factor',
'plot_nanomap', 'plot_real', 'plot_imag',
'plot_nanomap_with_images']
__all__ = [
'colorize',
'plot_amplitude',
'plot_phase',
'plot_colorized',
'plot_translations',
'get_units_factor',
'plot_nanomap',
'plot_real',
'plot_imag',
'plot_nanomap_with_images',
'cmocean_phase'
]
def colorize(z):
@@ -83,7 +93,21 @@ def get_units_factor(units):
return factor
def plot_image(im, plot_func=lambda x: x, fig=None, basis=None, view_basis='ortho', units='$\\mu$m', cmap='viridis', cmap_label=None, interpolation=None, **kwargs):
def plot_image(
im,
plot_func=lambda x: x,
fig=None,
basis=None,
view_basis='ortho',
units='$\\mu$m',
cmap='viridis',
cmap_label=None,
show_cbar=True,
vmin=None,
vmax=None,
interpolation=None,
**kwargs
):
"""Plots an image with a colorbar and on an appropriate spatial grid
If a figure is given explicitly, it will clear that existing figure and
@@ -113,7 +137,13 @@ def plot_image(im, plot_func=lambda x: x, fig=None, basis=None, view_basis='orth
cmap : str
Default is 'viridis', the colormap to plot with
cmap_label : str
What to label the colorbar when plotting
What to label the colorbar when plotting.
show_cbar : bool
Default is True, whether or not to show the colorbar
vmin : int
Default is min(plot_func(im)), the minimum value for the colormap
vmax : int
Default is max(plot_func(im)), the maximum value for the colormap
interpolation : str
What interpolation to use for imshow
\\**kwargs
@@ -161,7 +191,9 @@ def plot_image(im, plot_func=lambda x: x, fig=None, basis=None, view_basis='orth
mpl_im = plt.imshow(
to_plot,
cmap = cmap,
interpolation = interpolation
interpolation = interpolation,
vmin=vmin,
vmax=vmax,
)
plt.gca().set_facecolor('k')
@@ -224,10 +256,11 @@ def plot_image(im, plot_func=lambda x: x, fig=None, basis=None, view_basis='orth
plt.gca().set_xlim([mins[0], maxes[0]])
plt.gca().set_ylim([mins[1], maxes[1]])
plt.gca().invert_yaxis()
cbar = plt.colorbar()
if cmap_label is not None:
cbar.set_label(cmap_label)
if show_cbar:
cbar = plt.colorbar()
if cmap_label is not None:
cbar.set_label(cmap_label)
if basis is not None:
plt.xlabel('X (' + units + ')')
@@ -388,8 +421,18 @@ def plot_amplitude(im, fig = None, basis=None, units='$\\mu$m', cmap='viridis',
**kwargs)
def plot_phase(im, fig=None, basis=None, units='$\\mu$m', cmap='auto', cmap_label='Phase (rad)', **kwargs):
""" Plots the phase of a complex array with dimensions NxMx2
def plot_phase(
im,
fig=None,
basis=None,
units='$\\mu$m',
cmap='cividis',
cmap_label='Phase (rad)',
vmin=None,
vmax=None,
**kwargs
):
""" Plots the phase of a complex array with dimensions NxM
If a figure is given explicitly, it will clear that existing figure and
plot over it. Otherwise, it will generate a new figure.
@@ -397,6 +440,9 @@ def plot_phase(im, fig=None, basis=None, units='$\\mu$m', cmap='auto', cmap_labe
If a basis is explicitly passed, the image will be plotted in real-space
coordinates
If the cmap is entered as 'phase', it will plot the cmocean phase colormap,
and by default set the limits to [-pi,pi].
Parameters
----------
im : array
@@ -408,9 +454,14 @@ def plot_phase(im, fig=None, basis=None, units='$\\mu$m', cmap='auto', cmap_labe
units : str
The length units to mark on the plot, default is um
cmap : str
Default is 'viridis', the colormap to plot with
Default is 'cividis', the colormap to plot with.
cmap_label : str
What to label the colorbar when plotting
vmin : int
Default is min(angle(im)), the minimum value for the colormap
vmax : int
Default is max(angle(im)), the maximum value for the colormap
\\**kwargs
All other args are passed to fig.add_subplot(111, \\**kwargs)
@@ -419,17 +470,17 @@ def plot_phase(im, fig=None, basis=None, units='$\\mu$m', cmap='auto', cmap_labe
used_fig : matplotlib.figure.Figure
The figure object that was actually plotted to.
"""
if cmap == 'auto':
if 'twilight' in plt.colormaps():
cmap = 'twilight'
elif 'hsv' in plt.colormaps():
cmap = 'hsv'
else:
raise AttributeError('Neither twilight or hsv colormap exists in this screwed up matplotlib install')
plot_func = lambda x: np.angle(x)
if cmap == 'cyclic' or cmap == 'phase' or cmap == 'cmocean_phase':
cmap = cmocean_phase
vmin = (-np.pi if (vmin is None) else vmin)
vmax = (np.pi if (vmax is None) else vmax)
return plot_image(im, plot_func=plot_func, fig=fig, basis=basis,
units=units, cmap=cmap, cmap_label=cmap_label,
vmin=vmin,vmax=vmax,
**kwargs)
@@ -468,7 +519,7 @@ def plot_colorized(im, fig=None, basis=None, units='$\\mu$m', **kwargs):
"""
plot_func = lambda x: colorize(x)
return plot_image(im, plot_func=plot_func, fig=fig, basis=basis,
units=units, **kwargs)
units=units, show_cbar=False, **kwargs)
def plot_translations(translations, fig=None, units='$\\mu$m', lines=True, invert_xaxis=True, **kwargs):
@@ -837,3 +888,299 @@ def plot_nanomap_with_images(translations, get_image_func, values=None, mask=Non
update(0)
return fig
#
# Some code to include the "phase" colormap from cmocean, which is
# beautiful, without having to add a dependency on the whole cmocean
# package
#
# License and authorship info for the cmocean package, which this code
# is adapted from:
#
# The MIT License (MIT)
#
# Copyright (c) 2015 Kristen M. Thyng
# 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.
#
#
cm_data = [[ 0.65830839, 0.46993917, 0.04941288],
[ 0.66433742, 0.4662019 , 0.05766473],
[ 0.67020869, 0.46248014, 0.0653456 ],
[ 0.67604299, 0.45869838, 0.07273174],
[ 0.68175228, 0.45491407, 0.07979262],
[ 0.6874028 , 0.45108417, 0.08667103],
[ 0.6929505 , 0.44723893, 0.09335869],
[ 0.69842619, 0.44335768, 0.09992839],
[ 0.7038123 , 0.43945328, 0.1063871 ],
[ 0.70912069, 0.43551765, 0.11277174],
[ 0.71434524, 0.43155576, 0.11909348],
[ 0.71949289, 0.42756272, 0.12537606],
[ 0.72455619, 0.4235447 , 0.13162325],
[ 0.72954895, 0.41949098, 0.13786305],
[ 0.73445172, 0.41541774, 0.14408039],
[ 0.73929496, 0.41129973, 0.15032217],
[ 0.74403834, 0.40717158, 0.15654335],
[ 0.74873695, 0.40298519, 0.16282282],
[ 0.75332319, 0.39880107, 0.16907566],
[ 0.75788083, 0.39454245, 0.17542179],
[ 0.7623326 , 0.39028096, 0.18175915],
[ 0.76673205, 0.38596549, 0.18816819],
[ 0.77105247, 0.38162141, 0.19461532],
[ 0.77529528, 0.37724732, 0.20110652],
[ 0.77948666, 0.37281509, 0.2076873 ],
[ 0.78358534, 0.36836772, 0.21429736],
[ 0.78763763, 0.363854 , 0.22101648],
[ 0.79161134, 0.35930804, 0.2277974 ],
[ 0.79550606, 0.3547299 , 0.23464353],
[ 0.79935398, 0.35007959, 0.24161832],
[ 0.80311671, 0.34540152, 0.24865892],
[ 0.80681033, 0.34067452, 0.25580075],
[ 0.8104452 , 0.33588248, 0.26307222],
[ 0.8139968 , 0.33105538, 0.27043183],
[ 0.81747689, 0.32617526, 0.27791096],
[ 0.82089415, 0.32122629, 0.28553846],
[ 0.82422713, 0.3162362 , 0.29327617],
[ 0.82747661, 0.31120154, 0.30113388],
[ 0.83066399, 0.30608459, 0.30917579],
[ 0.83376307, 0.30092244, 0.31734921],
[ 0.83677286, 0.29571346, 0.32566199],
[ 0.83969693, 0.29044723, 0.33413665],
[ 0.84253873, 0.28511151, 0.34279962],
[ 0.84528297, 0.27972917, 0.35162078],
[ 0.84792704, 0.27430045, 0.36060681],
[ 0.85046793, 0.26882624, 0.36976395],
[ 0.85291056, 0.26328859, 0.37913116],
[ 0.855242 , 0.25770888, 0.38868217],
[ 0.85745673, 0.25209367, 0.39841601],
[ 0.85955023, 0.24644737, 0.40833625],
[ 0.86151767, 0.24077563, 0.41844557],
[ 0.86335392, 0.23508521, 0.42874606],
[ 0.86505685, 0.22937288, 0.43926008],
[ 0.86661606, 0.22366308, 0.44996127],
[ 0.86802578, 0.21796785, 0.46084758],
[ 0.86928003, 0.21230132, 0.47191554],
[ 0.87037274, 0.20667988, 0.48316015],
[ 0.87129781, 0.2011224 , 0.49457479],
[ 0.87204914, 0.19565041, 0.50615118],
[ 0.87262076, 0.19028829, 0.51787932],
[ 0.87300686, 0.18506334, 0.5297475 ],
[ 0.8732019 , 0.18000588, 0.54174232],
[ 0.87320066, 0.1751492 , 0.55384874],
[ 0.87299833, 0.17052942, 0.56605016],
[ 0.87259058, 0.16618514, 0.57832856],
[ 0.87197361, 0.16215698, 0.59066466],
[ 0.87114414, 0.15848667, 0.60303881],
[ 0.87009966, 0.15521687, 0.61542844],
[ 0.86883823, 0.15238892, 0.62781175],
[ 0.86735858, 0.15004199, 0.64016651],
[ 0.8656601 , 0.14821149, 0.65247022],
[ 0.86374282, 0.14692762, 0.66470043],
[ 0.86160744, 0.14621386, 0.67683495],
[ 0.85925523, 0.14608582, 0.68885204],
[ 0.85668805, 0.14655046, 0.70073065],
[ 0.85390829, 0.14760576, 0.71245054],
[ 0.85091881, 0.14924094, 0.7239925 ],
[ 0.84772287, 0.15143717, 0.73533849],
[ 0.84432409, 0.15416865, 0.74647174],
[ 0.84072639, 0.15740403, 0.75737678],
[ 0.83693394, 0.16110786, 0.76803952],
[ 0.83295108, 0.16524205, 0.77844723],
[ 0.82878232, 0.16976729, 0.78858858],
[ 0.82443225, 0.17464414, 0.7984536 ],
[ 0.81990551, 0.179834 , 0.80803365],
[ 0.81520674, 0.18529984, 0.8173214 ],
[ 0.81034059, 0.19100664, 0.82631073],
[ 0.80531176, 0.1969216 , 0.83499645],
[ 0.80012467, 0.20301465, 0.84337486],
[ 0.79478367, 0.20925826, 0.8514432 ],
[ 0.78929302, 0.21562737, 0.85919957],
[ 0.78365681, 0.22209936, 0.86664294],
[ 0.77787898, 0.22865386, 0.87377308],
[ 0.7719633 , 0.23527265, 0.88059043],
[ 0.76591335, 0.24193947, 0.88709606],
[ 0.7597325 , 0.24863985, 0.89329158],
[ 0.75342394, 0.25536094, 0.89917908],
[ 0.74699063, 0.26209137, 0.90476105],
[ 0.74043533, 0.2688211 , 0.91004033],
[ 0.73376055, 0.27554128, 0.91502 ],
[ 0.72696862, 0.28224415, 0.91970339],
[ 0.7200616 , 0.2889229 , 0.92409395],
[ 0.71304134, 0.29557159, 0.92819525],
[ 0.70590945, 0.30218508, 0.9320109 ],
[ 0.69866732, 0.30875887, 0.93554451],
[ 0.69131609, 0.31528914, 0.93879964],
[ 0.68385669, 0.32177259, 0.94177976],
[ 0.6762898 , 0.32820641, 0.94448822],
[ 0.6686159 , 0.33458824, 0.94692818],
[ 0.66083524, 0.3409161 , 0.94910264],
[ 0.65294785, 0.34718834, 0.95101432],
[ 0.64495358, 0.35340362, 0.95266571],
[ 0.63685208, 0.35956083, 0.954059 ],
[ 0.62864284, 0.3656591 , 0.95519608],
[ 0.62032517, 0.3716977 , 0.95607853],
[ 0.61189825, 0.37767607, 0.95670757],
[ 0.60336117, 0.38359374, 0.95708408],
[ 0.59471291, 0.3894503 , 0.95720861],
[ 0.58595242, 0.39524541, 0.95708134],
[ 0.5770786 , 0.40097871, 0.95670212],
[ 0.56809041, 0.40664983, 0.95607045],
[ 0.55898686, 0.41225834, 0.95518556],
[ 0.54976709, 0.41780374, 0.95404636],
[ 0.5404304 , 0.42328541, 0.95265153],
[ 0.53097635, 0.42870263, 0.95099953],
[ 0.52140479, 0.43405447, 0.94908866],
[ 0.51171597, 0.43933988, 0.94691713],
[ 0.50191056, 0.44455757, 0.94448311],
[ 0.49198981, 0.44970607, 0.94178481],
[ 0.48195555, 0.45478367, 0.93882055],
[ 0.47181035, 0.45978843, 0.93558888],
[ 0.46155756, 0.46471821, 0.93208866],
[ 0.45119801, 0.46957218, 0.92831786],
[ 0.44073852, 0.47434688, 0.92427669],
[ 0.43018722, 0.47903864, 0.9199662 ],
[ 0.41955166, 0.4836444 , 0.91538759],
[ 0.40884063, 0.48816094, 0.91054293],
[ 0.39806421, 0.49258494, 0.90543523],
[ 0.38723377, 0.49691301, 0.90006852],
[ 0.37636206, 0.50114173, 0.89444794],
[ 0.36546127, 0.5052684 , 0.88857877],
[ 0.35454654, 0.5092898 , 0.88246819],
[ 0.34363779, 0.51320158, 0.87612664],
[ 0.33275309, 0.51700082, 0.86956409],
[ 0.32191166, 0.52068487, 0.86279166],
[ 0.31113372, 0.52425144, 0.85582152],
[ 0.3004404 , 0.52769862, 0.84866679],
[ 0.28985326, 0.53102505, 0.84134123],
[ 0.27939616, 0.53422931, 0.83386051],
[ 0.26909181, 0.53731099, 0.82623984],
[ 0.258963 , 0.5402702 , 0.81849475],
[ 0.24903239, 0.54310763, 0.8106409 ],
[ 0.23932229, 0.54582448, 0.80269392],
[ 0.22985664, 0.54842189, 0.79467122],
[ 0.2206551 , 0.55090241, 0.78658706],
[ 0.21173641, 0.55326901, 0.77845533],
[ 0.20311843, 0.55552489, 0.77028973],
[ 0.1948172 , 0.55767365, 0.76210318],
[ 0.1868466 , 0.55971922, 0.75390763],
[ 0.17921799, 0.56166586, 0.74571407],
[ 0.1719422 , 0.56351747, 0.73753498],
[ 0.16502295, 0.56527915, 0.72937754],
[ 0.15846116, 0.566956 , 0.72124819],
[ 0.15225499, 0.56855297, 0.71315321],
[ 0.14639876, 0.57007506, 0.70509769],
[ 0.14088284, 0.57152729, 0.69708554],
[ 0.13569366, 0.57291467, 0.68911948],
[ 0.13081385, 0.57424211, 0.68120108],
[ 0.12622247, 0.57551447, 0.67333078],
[ 0.12189539, 0.57673644, 0.66550792],
[ 0.11780654, 0.57791235, 0.65773233],
[ 0.11392613, 0.5790468 , 0.64999984],
[ 0.11022348, 0.58014398, 0.64230637],
[ 0.10666732, 0.58120782, 0.63464733],
[ 0.10322631, 0.58224198, 0.62701729],
[ 0.0998697 , 0.58324982, 0.61941001],
[ 0.09656813, 0.58423445, 0.61181853],
[ 0.09329429, 0.58519864, 0.60423523],
[ 0.09002364, 0.58614483, 0.5966519 ],
[ 0.08673514, 0.58707512, 0.58905979],
[ 0.08341199, 0.58799127, 0.58144971],
[ 0.08004245, 0.58889466, 0.57381211],
[ 0.07662083, 0.58978633, 0.56613714],
[ 0.07314852, 0.59066692, 0.55841474],
[ 0.06963541, 0.5915367 , 0.55063471],
[ 0.06610144, 0.59239556, 0.54278681],
[ 0.06257861, 0.59324304, 0.53486082],
[ 0.05911304, 0.59407833, 0.52684614],
[ 0.05576765, 0.5949003 , 0.5187322 ],
[ 0.05262511, 0.59570732, 0.51050978],
[ 0.04978881, 0.5964975 , 0.50216936],
[ 0.04738319, 0.59726862, 0.49370174],
[ 0.04555067, 0.59801813, 0.48509809],
[ 0.04444396, 0.59874316, 0.47635 ],
[ 0.04421323, 0.59944056, 0.46744951],
[ 0.04498918, 0.60010687, 0.45838913],
[ 0.04686604, 0.60073837, 0.44916187],
[ 0.04988979, 0.60133103, 0.43976125],
[ 0.05405573, 0.60188055, 0.4301812 ],
[ 0.05932209, 0.60238289, 0.42040543],
[ 0.06560774, 0.60283258, 0.41043772],
[ 0.07281962, 0.60322442, 0.40027363],
[ 0.08086177, 0.60355283, 0.38990941],
[ 0.08964366, 0.60381194, 0.37934208],
[ 0.09908952, 0.60399554, 0.36856412],
[ 0.10914617, 0.60409695, 0.35755799],
[ 0.11974119, 0.60410858, 0.34634096],
[ 0.13082746, 0.6040228 , 0.33491416],
[ 0.14238003, 0.60383119, 0.323267 ],
[ 0.1543847 , 0.60352425, 0.31138823],
[ 0.16679093, 0.60309301, 0.29931029],
[ 0.17959757, 0.60252668, 0.2870237 ],
[ 0.19279966, 0.60181364, 0.27452964],
[ 0.20634465, 0.60094466, 0.2618794 ],
[ 0.22027287, 0.5999043 , 0.24904251],
[ 0.23449833, 0.59868591, 0.23611022],
[ 0.24904416, 0.5972746 , 0.2230778 ],
[ 0.26382006, 0.59566656, 0.21004673],
[ 0.2788104 , 0.5938521 , 0.19705484],
[ 0.29391494, 0.59183348, 0.18421621],
[ 0.3090634 , 0.58961302, 0.17161942],
[ 0.32415577, 0.58720132, 0.15937753],
[ 0.3391059 , 0.58461164, 0.14759012],
[ 0.35379624, 0.58186793, 0.13637734],
[ 0.36817905, 0.5789861 , 0.12580054],
[ 0.38215966, 0.57599512, 0.1159504 ],
[ 0.39572824, 0.57290928, 0.10685038],
[ 0.40881926, 0.56975727, 0.09855521],
[ 0.42148106, 0.56654159, 0.09104002],
[ 0.43364953, 0.56329296, 0.08434116],
[ 0.44538908, 0.56000859, 0.07841305],
[ 0.45672421, 0.5566943 , 0.07322913],
[ 0.46765017, 0.55336373, 0.06876762],
[ 0.47819138, 0.5500213 , 0.06498436],
[ 0.48839686, 0.54666195, 0.06182163],
[ 0.49828924, 0.5432874 , 0.05922726],
[ 0.50789114, 0.53989827, 0.05714466],
[ 0.51722475, 0.53649429, 0.05551476],
[ 0.5263115 , 0.53307443, 0.05427793],
[ 0.53517186, 0.52963707, 0.05337567],
[ 0.54382515, 0.52618009, 0.05275208],
[ 0.55228947, 0.52270103, 0.05235479],
[ 0.56058163, 0.51919713, 0.0521356 ],
[ 0.56871719, 0.51566545, 0.05205062],
[ 0.57671045, 0.51210292, 0.0520602 ],
[ 0.5845745 , 0.50850636, 0.05212851],
[ 0.59232129, 0.50487256, 0.05222299],
[ 0.5999617 , 0.50119827, 0.05231367],
[ 0.60750568, 0.49748022, 0.05237234],
[ 0.61496232, 0.49371512, 0.05237168],
[ 0.62233999, 0.48989963, 0.05228423],
[ 0.62964652, 0.48603032, 0.05208127],
[ 0.63688935, 0.48210362, 0.05173155],
[ 0.64407572, 0.4781157 , 0.0511996 ],
[ 0.65121289, 0.47406244, 0.05044367],
[ 0.65830839, 0.46993917, 0.04941288]]
rgb = np.array(cm_data)
rgb_with_alpha = np.zeros((rgb.shape[0],4))
rgb_with_alpha[:,:3] = rgb
rgb_with_alpha[:,3] = 1. #set alpha channel to 1
cmocean_phase = colors.ListedColormap(rgb_with_alpha, N=rgb.shape[0])