diff --git a/.gitignore b/.gitignore index f40e32a..de24304 100644 --- a/.gitignore +++ b/.gitignore @@ -7,4 +7,4 @@ *.tif build/ dist/ -*.egg_info/ +*.egg-info/ diff --git a/README.md b/README.md index 6e5f50f..7858899 100644 --- a/README.md +++ b/README.md @@ -1,34 +1,109 @@ -Bsisb (Xi-cam PluginMaker) -=============================== +BSISB (A Xi-cam plugin for FTIR data analysis) +============================================== -version number: 0.1 -author: Liang Chen +Version number: 0.1 -Overview --------- +Author: Liang Chen +## Installation +--------------------- +### Install C++ build tools (Windows OS only): +Download and install "Build Tools for Visual Studio" from [here](https://visualstudio.microsoft.com/downloads/). +See the screenshot below for reference: -Installation / Usage --------------------- -Install Xi-CAM: +![build](images/buildtools.png) -Follow the guide [here](https://xi-cam2.readthedocs.io/en/latest/install.html) +### Install git and python3.7.7: +You will need to ensure that you have both git and python3 installed on your system for Xi-cam installation. -Clone the repo: +**MacOS** + +Open the Terminal application (in Applications/Utilities). In the terminal, check to see if git is installed by typing `git --version`. Either a version number will be printed, indicating git is already installed, or a dialog will open asking The “git” command requires the command line developer tools. Would you like to install the tools now? Click the **Install** button to install the developer tools, which will install git for you. + +You will also need to install python3.7.7. You can download python3.7.7 at [python.org](https://www.python.org/downloads/mac-osx/). You will want to get the macOS 64-bit installer if you are running any macOS version since Mavericks. + +**Windows OS** + +Download git [here](https://git-scm.com/download/win) and follow the installer’s instructions. This will install **Git for Windows**, which provides a **Git Bash** command line interface as well as a **Git GUI** graphical interface for git. + +You will want to go to the python3.7.7 download page at [python.org](https://www.python.org/downloads/release/python-377/). For modern systems, install the **Windows x86-64 executable installer** at the bottom. + +When you run the installer, make sure sure to check the box that says **Add Python 3.x to PATH** to ensure that the interpreter will be placed in your execution path. + +### Create and Activate a Virtual Environment: + +The latest python3 version comes with the venv module, which can be used to create a virtual environment. A virtual environment is a sequestered space where you can install and uninstall packages without modifying your system’s installed packages. + +Create a virtual environment for installing the Xi-cam components and dependencies. You will then want to activate the virtual environment you created so that any packages you install with python’s package manager, **pip**, will be installed into that active virtual environment. In the commands below, create a virtual environment called **venv_xicam** and activate it: + +**Linux/macOS** + + $ python3 -m venv venv_xicam + $ source venv_xicam/bin/actviate + +**Windows OS** + + $ python -m venv venv_xicam + $ venv_xicam\Scripts\activate + +### Install Xi-CAM: +Before installing Xi-cam, run the following commands: + + $ python -m pip install --upgrade pip + $ pip install wheel + $ pip install --upgrade setuptools + $ pip install numcodecs + +Then run the following commands to install Xi-cam: + + $ pip install xicam + +### Install lbl_ir and Xi-cam.BSISB plugin: +Run the following commands to install lbl_ir package and Xi-cam.BSISB plugin: $ git clone https://github.com/lchen23/Xi-cam.BSISB.git - $ cd Xi-CAM.BSISB + $ cd Xi-cam.BSISB/lbl_ir $ pip install -e . - $ cd .. + $ cd .. + $ pip install -e . + +### Run the Xi-cam.BSISB program: + +Run the following command: + $ xicam + +### Run the Xi-cam.BSISB GUI from a new terminal: +First, activate the virtual environment, then run xicam: + +**Linux/macOS** + + $ source venv_xicam/bin/actviate + $ xicam + +**Windows OS** + +In Git bash: + + $ source venv_xicam/Scripts/activate + $ xicam + +In Windows CMD: + + $ venv_xicam\Scripts\activate + $ xicam + +After the Xi-cam GUI program launches, open the +`Xi-cam.BSISB/tests/PC12-NGF-3h.h5` file by double clicking the h5 file +in the file browser as shown below: + +![example](images/example.png) + Contributing ------------ TBD -Example -------- -TBD diff --git a/images/buildtools.png b/images/buildtools.png new file mode 100644 index 0000000..6f81aed Binary files /dev/null and b/images/buildtools.png differ diff --git a/images/example.png b/images/example.png new file mode 100644 index 0000000..516b3d8 Binary files /dev/null and b/images/example.png differ diff --git a/lbl_ir/.gitignore b/lbl_ir/.gitignore new file mode 100644 index 0000000..cc715cc --- /dev/null +++ b/lbl_ir/.gitignore @@ -0,0 +1,2 @@ +__pycache__/ +.ipynb_checkpoints/ \ No newline at end of file diff --git a/lbl_ir/README.md b/lbl_ir/README.md new file mode 100644 index 0000000..234c43d --- /dev/null +++ b/lbl_ir/README.md @@ -0,0 +1,30 @@ +# README # + + +This README would normally document whatever steps are necessary to get your application up and running. + +### What is this repository for? ### + +* Quick summary +* Version +* [Learn Markdown](https://bitbucket.org/tutorials/markdowndemo) + +### How do I get set up? ### + +* Summary of set up +* Configuration +* Dependencies +* Database configuration +* How to run tests +* Deployment instructions + +### Contribution guidelines ### + +* Writing tests +* Code review +* Other guidelines + +### Who do I talk to? ### + +* Repo owner or admin +* Other community or team contact diff --git a/lbl_ir/lbl_ir/GPR/GPR_engine.py b/lbl_ir/lbl_ir/GPR/GPR_engine.py new file mode 100644 index 0000000..e1c7f75 --- /dev/null +++ b/lbl_ir/lbl_ir/GPR/GPR_engine.py @@ -0,0 +1,143 @@ +import numpy as np +import matplotlib.pyplot as plt +import numba +from scipy.optimize import minimize +from scipy.optimize import basinhopping + +@numba.jit(nopython=True) #,cache=True) +def dmat(X1,X2=None): + if X2 is None: + X2 = X1 + N1,M1 = X1.shape + N2,M2 = X2.shape + + dd = np.zeros( (N2,N1) ) + for ii in range(N2): + for jj in range(N1): + tmp = 0 + for kk in range(M1): + delta = X2[ii,kk]-X1[jj,kk] + tmp += delta*delta + tmp = np.sqrt( tmp ) + dd[ii,jj]=tmp + return dd + +@numba.jit(nopython=True) #,cache=True) +def rbf_kern(dmat, variance, length): + result = np.exp( -dmat*dmat/(1e-8+length*length*2.0))*(1e-8+variance) + return result + +@numba.jit(nopython=True) #,cache=True) +def d_rbf_kern_FD(X2, X, variance, length, h=1e-4): + dmat_plus = dmat(X,X2+h/2) + dmat_minus= dmat(X,X2-h/2) + Kplus = np.exp( -dmat_plus*dmat_plus/(length*length*2.0) )*variance + Kminus = np.exp( -dmat_minus*dmat_minus/(length*length*2.0) )*variance + result = (Kplus - Kminus)/h + return result + +def d_rbf_kern(X2, X, dxx2, Kxx2, variance, length): + # first derivative + result1 = [] + for xx in X2.flatten(): + result1.append( X.flatten()-xx ) + result1 = np.vstack(result1) + result1 = result1*Kxx2 + + # second derivative + result2 = [] + N2,N1 = result1.shape + for ii in range(N2): + xx = X2[ii,:] + dK = result1[ii,:] + K = Kxx2[ii,:] + tmp = -K/(length*length) + tmp2 = (X-xx).flatten()*dK.flatten()/(length*length) + result2.append( tmp.flatten()+tmp2.flatten() ) + result2 = np.vstack(result2) + return result1, result2 + + + +class GPR_exp_fitter(object): + def __init__(self,sigma,length, mu=0,niter=3): + """ + :param sigma: The sigma + :param length: The length scale + :param mu: The mean value + :param niter: the number of iterations in basin hopping to fit the GPR model + """ + + self.sigma = sigma + self.length = length + self.mu = mu + + # to compute later + self.dX1X1 = None + + + def marginal_likelihood(self,hparams): + this_mu = hparams[0] + this_sigma = hparams[1] + this_length= hparams[2] + KX1X1 = rbf_kern( self.dX1X1, this_sigma*this_sigma, this_length ) + + II = np.diag(self.sY*self.sY) + KX1X1_inv = np.linalg.pinv(KX1X1+II) + tY = self.Y - this_mu + tY = tY.reshape(-1,1) + KiY = KX1X1_inv.dot(tY) + term1 = 0.5*tY.transpose().dot(KiY) + term2 = 0.5*np.linalg.slogdet(KX1X1+II)[1] + + result = term1+term2 + result = result.flatten()[0] + return result + + + def fit(self,X,Y,sY): + # first we need to build a distance matrix + self.sY = sY + if type(self.sY) is float: + self.sY = self.Y*0+sY + self.X = X + self.Y = Y + self.dX1X1 = dmat(X.reshape(-1,1)) + fitter = basinhopping( func = self.marginal_likelihood, + x0 = np.array([self.mu, self.sigma, self.length]), + niter = 13) + self.mu = fitter.x[0] + self.sigma = abs(fitter.x[1]) + self.length = abs(fitter.x[2]) + + + self.tY = (Y - self.mu).reshape(-1,1) + self.KX1X1 = rbf_kern( self.dX1X1, self.sigma**2.0, self.length ) + np.eye( self.dX1X1.shape[0])*self.sY + self.KX1X1_inv = np.linalg.pinv(self.KX1X1) + self.KiY = self.KX1X1_inv.dot(self.tY) + + + + def predict(self,X): + dX1X2 = dmat(self.X.reshape(-1,1), X.reshape(-1,1)) + kX1X2 = rbf_kern( dX1X2, self.sigma**2, self.length ) + tmp = kX1X2.dot( self.KiY ) + return tmp+self.mu + + + + +def tst(): + x = np.linspace(-10,10,20, True) + xx = np.linspace(-10,10,200,True) + y = np.sin(x) #10 + x*x + x + + plt.plot(x,y,'.');plt.show() + obj = GPR_exp_fitter(1,1,1) + obj.fit(x,y,0.01) + yy = obj.predict( xx ).flatten() + plt.plot(x,y,'.'); plt.plot(xx,yy);plt.show() + + +if __name__ == "__main__": + tst() diff --git a/lbl_ir/lbl_ir/GPR/GPR_peaks.py b/lbl_ir/lbl_ir/GPR/GPR_peaks.py new file mode 100644 index 0000000..32c7069 --- /dev/null +++ b/lbl_ir/lbl_ir/GPR/GPR_peaks.py @@ -0,0 +1,227 @@ +import numpy as np +import matplotlib.pyplot as plt +import numba +from numpy.polynomial.polynomial import Polynomial +from scipy.optimize import minimize +from lbl_ir.GPR import GPR_engine + +""" +Here I try to use Gaussian process regression for identifying peaks +and estimate their standard deviation. +I build a simple GPR fitter myself because I needed more control over +kernel and hyper parameters. Things need to be cleaned up. + +""" + + +#@numba.jit(nopython=True) +def dmat(X1,X2=None): + if X2 is None: + X2 = X1 + N1,M1 = X1.shape + N2,M2 = X2.shape + + dd = np.zeros( (N2,N1) ) + for ii in range(N2): + for jj in range(N1): + tmp = 0 + for kk in range(M1): + delta = X2[ii,kk]-X1[jj,kk] + tmp += delta*delta + tmp = np.sqrt( tmp ) + dd[ii,jj]=tmp + return dd + +@numba.jit(nopython=True) +def rbf_kern(dmat, variance, length): + result = np.exp( -dmat*dmat/(length*length*2.0))*variance + return result + +def d_rbf_kerni_FD(X2, X, variance, length, h=1e-4): + dmat_plus = dmat(X,X2+h/2) + dmat_minus= dmat(X,X2-h/2) + Kplus = np.exp( -dmat_plus*dmat_plus/(length*length*2.0) )*variance + Kminus = np.exp( -dmat_minus*dmat_minus/(length*length*2.0) )*variance + result = (Kplus - Kminus)/h + return result + +def d_rbf_kern(X2, X, dxx2, Kxx2, variance, length): + # first derivative + result1 = [] + for xx in X2.flatten(): + result1.append( X.flatten()-xx ) + result1 = np.vstack(result1) + result1 = result1*Kxx2 + + # second derivative + result2 = [] + N2,N1 = result1.shape + for ii in range(N2): + xx = X2[ii,:] + dK = result1[ii,:] + K = Kxx2[ii,:] + tmp = -K/(length*length) + tmp2 = (X-xx).flatten()*dK.flatten()/(length*length) + result2.append( tmp.flatten()+tmp2.flatten() ) + result2 = np.vstack(result2) + return result1, result2 + + + +class peak_picker(object): + def __init__(self,X,Y,sY): + self.X = X + self.Y = Y + self.N = X.shape[0] + self.sY = sY + init_scale = np.std( self.X.flatten()) + obj = GPR_engine.GPR_exp_fitter( sigma=np.std(Y.flatten()), length=init_scale) + obj.fit( self.X, self.Y, self.sY ) + + self.lengthscale = obj.length + self.variance = obj.sigma**2.0 + self.mu = obj.mu + + # set things up for peak picking purposes + self.dXX = dmat(self.X) + self.Kxx = rbf_kern(self.dXX,self.variance,self.lengthscale) + self.Keff = np.eye( self.N )*self.sY*self.sY + self.Kxx + self.Keff_inv = np.linalg.pinv(self.Keff,rcond=1e-10) + self.Keff_invY= self.Keff_inv.dot(self.Y-self.mu) + + self.dx = np.mean( np.sort( self.dXX.flatten() )[self.N:self.N*2] ) + + + + + def predict(self,x, grads=True): + # setup stuff + dxx2 = dmat(self.X, x ) + dx2x2 = dmat(x) + Kxx2 = rbf_kern( dxx2, self.variance,self.lengthscale ) + Kx2x2 = rbf_kern( dx2x2, self.variance, self.lengthscale ) + + # mean and variance + cond_mean = Kxx2.dot( self.Keff_invY ) + self.mu + cond_variance = Kx2x2 - Kxx2.dot(self.Keff_inv).dot(Kxx2.transpose()) + + if grads: + # derivatives if requested + dK,ddK = d_rbf_kern(x, self.X, dxx2, Kxx2, self.variance, self.lengthscale) + dmu = dK.dot( self.Keff_invY ) + ddmu = ddK.dot( self.Keff_invY ) + return cond_mean, cond_variance, dmu, ddmu + + else: + return cond_mean, cond_variance + + def f(self,x): + xx = np.array([x]).reshape(-1,1) + result = self.predict(xx)[0][0,0]+self.mu + return -result + + + def find_peak(self, x, eps=1e-2, max_iter=1000): + init_simplex = np.array( [x, x+eps] ).reshape(-1,1) + x_opt = minimize(fun=self.f, method='Nelder-Mead',x0=x, options={'initial_simplex':init_simplex} ) + return x_opt.x + + def find_peak_oof(self, x, eps=1e-8,max_iter=1000, damp=1.0, peak_range=None): + # do a simple root finding starting from x + converged=False + x_in = np.array(x).reshape(-1,1) + iter = 0 + restart = 0 + max_delta = 1.0 + if peak_range is not None: + max_delta = np.max_peak(range) - np.min(peak_range) + max_delta = max_delta / 10.0 + while not converged: + m,v,dm,ddm = self.predict(x_in,True) + tmp = np.random.uniform(0.95,1.05,1)[0] + #print("HERE",m,v,dm,ddm) + delta = tmp*dm/ddm + #if np.abs(delta) > max_delta: + # delta = max_delta*(np.sign(delta)) + #print(x_in, delta, dm, ddm,"DELTA", x_in - delta*damp ) + x_in = x_in - delta*damp + iter+=1 + if np.abs(dm) < eps: + converged=True + if iter > max_iter: + converged = False + print("Early termination in peak finding. Apply damping") + restart += 1 + iter = 0 + x_in = np.array(x).reshape(-1,1) + damp = damp * 0.5 + if restart > 5: + print("Early termination in peak finding. Damping doesn't work") + converged = True + return None + return x_in[0,0] + + + def fst_der_prop(self,x): + x0 = np.array([x]).reshape(-1,1) + # we need to get the vector + M = self.X.flatten()-x + M = np.diag(M/self.lengthscale*self.lengthscale) + + dxx2 = dmat(self.X, x0 ) + dx2x2 = dmat(x0) + Kxx2 = rbf_kern( dxx2, self.variance,self.lengthscale ) + Kx2x2 = rbf_kern( dx2x2, self.variance, self.lengthscale ) + nmu = M.dot(Kxx2.transpose()).transpose().dot(self.Keff_invY) + tmp = M.dot(Kxx2.transpose()).transpose() + nmu = tmp.dot(self.Keff_invY) + nvar = Kx2x2 - tmp.dot(self.Keff_inv).dot(tmp.transpose()) + return nmu[0][0], nvar[0][0] + + def peak_and_std_via_resample(self,x_start,N=10, factor = 0.75): + x_peak = self.find_peak( x_start ) + h = self.dx*factor + x_around = np.array( [x_peak-h, x_peak, x_peak+h] ).reshape(-1,1) + # collect the new posterior mean and variance + mu, var = self.predict(x_around, grads=False) + fs = np.random.multivariate_normal( mu.flatten(), var, N ) + ds = [] + for ii in range(N): + c = Polynomial.fit( x_around.flatten()-x_peak, fs[ii,:], deg = 2 ) + c = c.coef + dx = -0.5*c[1]/c[2] + ds.append(dx) + return x_peak, np.std( ds ), mu[1], np.sqrt(np.abs(var[1][1])) + + +def tst(P=8,S=0.1, show=False): + X = np.random.uniform(-3,3,P).reshape(-1,1) # np.linspace(-3,3,P).reshape((P,1)) + X = np.linspace(-3,3,P).reshape((P,1)) + Y = 10.0*np.exp(-X*X) + np.random.normal(0,S, (P,1) ) + + this_one = np.argmax( Y.flatten() ) + this_x = X.flatten()[this_one] + S = np.zeros( P )+S + obj = peak_picker(X,Y,S) + peak,sigma,val,sig = obj.peak_and_std_via_resample(this_x, N=100) + assert abs(peak/sigma) < 4 + print( "OK" ) + if show: + Xstar = np.linspace(-3,3,1024).reshape(-1,1) + y,v,dm,ddm = obj.predict(Xstar, True) + plt.plot( X.flatten() , Y.flatten(), '.' ) + plt.plot( Xstar.flatten(), y.flatten(), '-' ) + #plt.plot( Xstar.flatten(), dm.flatten(), '--', lw=3) + #plt.plot( Xstar.flatten(), ddm.flatten(), '--', lw=2) + print(peak) + peak_val, peak_std = obj.predict(np.array([[peak[0]]]), False) + print(peak_val, peak_std) + plt.errorbar( peak, peak_val, peak_std, sigma )# '.', markersize=10) + plt.show() + + + + + +if __name__ =="__main__": + tst(P=35,S=0.25,show=True) diff --git a/tests/__init__.py b/lbl_ir/lbl_ir/GPR/__init__.py similarity index 100% rename from tests/__init__.py rename to lbl_ir/lbl_ir/GPR/__init__.py diff --git a/lbl_ir/lbl_ir/GPR/spectral_peak_picker.py b/lbl_ir/lbl_ir/GPR/spectral_peak_picker.py new file mode 100644 index 0000000..45d8fc3 --- /dev/null +++ b/lbl_ir/lbl_ir/GPR/spectral_peak_picker.py @@ -0,0 +1,116 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.signal import find_peaks +import sys + +from lbl_ir.io_tools.read_map import read_all_formats +from lbl_ir.tasks.preprocessing.transform import to_absorbance +from lbl_ir.GPR.GPR_peaks import peak_picker + + +class spectral_peak_picker(object): + def __init__(self, + wavenumbers, + spectrum, + sigma=None, + window=5, + peak_height_threshold=0.1, + peak_separation_threshold=3, + peak_prominence=0.001, + ): + """ + + :param wavenumbers: The wavenumbers + :param spectrum: The IR spectrum (Absorbance) + :param sigma: The associated standard deviation. If None, a suitable default will be chosen. + :param window: A parameter that determines the window size when fitting peak using a GPR approach. + The default should be fine. + :param peak_height_threshold: minimum peak height + :param peak_separation_threshold: minimum distance between peaks + :param peak_prominence: peak prominance + """ + + self.wavenumbers = wavenumbers + self.Nwavs = len(wavenumbers) + self.spectrum = spectrum + self.sigma = sigma + self.window = window + self.peak_height_threshold = peak_height_threshold + self.peak_separation_threshold = peak_separation_threshold + self.peak_prominence = peak_prominence + + # first pass: find the peak using a simple approach + self.raw_peaks_indx, self.peak_props = find_peaks(self.spectrum, + height=self.peak_height_threshold, + distance=self.peak_separation_threshold, + prominence=self.peak_prominence ) + self.prominances = self.peak_props['prominences'] + + def refine_peaks(self, level=0.05, sigma_multi=0.0005): + these_peaks = self.prominances > level + these_peak_indx = self.raw_peaks_indx[ these_peaks ] + # now loop over these peaks and do the GPR stuff + peak_locations = [] + peak_sigmas = [] + peak_vals = [] + val_sigmas = [] + ok_flags = [] + for ii in these_peak_indx: + min_indx = max(ii-self.window,0) + max_indx = min(ii+self.window+1,self.Nwavs) + these_waves = self.wavenumbers[min_indx:max_indx] + these_specs = self.spectrum[min_indx:max_indx] + obs_sigma = np.sqrt(np.abs( these_specs ))*sigma_multi + + + + obj = peak_picker(X=these_waves.reshape(-1,1),Y=these_specs.reshape(-1,1),sY=obs_sigma) + x_start = np.mean(these_waves.flatten()) + + peak, sigma, val, val_sig = obj.peak_and_std_via_resample(x_start = x_start ) + if abs(peak - x_start) > 3*(these_waves[1]- these_waves[0]): + peak = x_start + sigma = 2*(these_waves[1]- these_waves[0]) + val = self.spectrum[ii] + val_sig = -1 + ok_flags.append(False) + else: + ok_flags.append(True) + + peak_locations.append(peak) + peak_sigmas.append(sigma) + peak_vals.append(val) + val_sigmas.append(val_sig) + return np.array(peak_locations).flatten(), \ + np.array(peak_sigmas).flatten(), \ + np.array(peak_vals).flatten(), \ + np.array(val_sigmas).flatten(), \ + np.array(ok_flags).flatten() + + + + + + + + + + + +def tst(filename): + map, fmt = read_all_formats(filename) + data,bg = to_absorbance( map.data,map.wavenumbers ) + waves = map.wavenumbers + spec = np.mean(data, axis=0) + #plt.plot(spec);plt.show() + obj = spectral_peak_picker(waves, spec, sigma=0.05, peak_height_threshold=0.1) + peaks, sigma, val, vs, ok = obj.refine_peaks( 0.005) + plt.plot(waves, spec,'.-', markersize=4) + plt.plot(peaks,val,'x',markersize=8);plt.show() + plt.plot( peaks, sigma, '.', markersize=4); plt.show() + + + + +if __name__ =="__main__": + tst(sys.argv[1]) diff --git a/lbl_ir/lbl_ir/__init__.py b/lbl_ir/lbl_ir/__init__.py new file mode 100644 index 0000000..4f1fbed --- /dev/null +++ b/lbl_ir/lbl_ir/__init__.py @@ -0,0 +1 @@ +from . import data_objects, GPR, gui_tools, io_tools, math_tools, simulations, tasks diff --git a/xicam.BSISB.egg-info/dependency_links.txt b/lbl_ir/lbl_ir/data_objects/__init__.py similarity index 50% rename from xicam.BSISB.egg-info/dependency_links.txt rename to lbl_ir/lbl_ir/data_objects/__init__.py index 8b13789..139597f 100644 --- a/xicam.BSISB.egg-info/dependency_links.txt +++ b/lbl_ir/lbl_ir/data_objects/__init__.py @@ -1 +1,2 @@ + diff --git a/lbl_ir/lbl_ir/data_objects/ir_map.py b/lbl_ir/lbl_ir/data_objects/ir_map.py new file mode 100644 index 0000000..d5cc26e --- /dev/null +++ b/lbl_ir/lbl_ir/data_objects/ir_map.py @@ -0,0 +1,504 @@ +import numpy as np +import h5py +import datetime +import sys +import os +import matplotlib.pyplot as plt + +def val2ind(val, an_array): + return np.argmin(abs(an_array-val), axis=0) + +class sample_info(object): + """ + A simple class that contains sample info. + + Arguments: + ---------- + + sample_id : A string that identifies the sample, spaces will be + substituted for underscores. + Specifying this is mandatory. + + sample_meta_data : More verbose description of the data. Having this + structured isn't a bad idea, but not enforced at this + level. Not requiered but highly encouraged. + + sample_date : The date at which the data was taken. Default is + today. + + Attributes: + ----------- + + show(out) : prints the contents of id, meta data and date to out. + If out is None, it defaults to sys.stdout + + Examples: + --------- + + si = sample_info(sample_id = 'C_elegans', + sample_date='2004_04_18', + sample_meta_data='Some details.') + si.show() + + """ + def __init__(self, sample_id="Unknown", sample_meta_data="None", sample_date=None): + assert ' ' not in sample_id + self.sample_id = sample_id + + self.sample_meta_data = sample_meta_data + if sample_date is None: + sample_date = str(datetime.date.today().year)+ \ + '_'+str(datetime.date.today().month)+ \ + '_'+str(datetime.date.today().day) + + self.sample_date = sample_date + + def show(self, out=None): + if out is None: + out = sys.stdout + print("Sample_id : %s "%self.sample_id, file=out) + print("Sample date : %s "%self.sample_date, file=out ) + print("Sample description:\n ", self.sample_meta_data, file=out ) + + +class ir_map(object): + """ A simple data object that contains IR data. + + Arguments: + ---------- + wavenumbers : An array of wavenumbers + + sample_info : A sample_info object + + filename : The hdf5 filename where data will be read from. Required for + the 'hdf5' mode. + + data_type : Either 'transmission', 'reflection', or 'absorbance'(default) + + _mode : Choice between 'memory' or 'hdf5' + If 'memory', the dataset currenly resides in memory + If 'hdf5', the dataset currenly resides in an hdf5 file + + _N_obs : The number of rows in the 2D spectral matrix. + + Attributes: + ----------- + As above plus + + add_data(self, spectrum, xy) : add data, provide a numpy array of spectra + and xy positions + + _allocate_space(self) : this is a function that allocates space for + an hdf5 file. no need to call it yourself. + + write_as_hdf5(self, filename) : Writes in-memory data as an hdf5 file. + + """ + + def __init__(self, + wavenumbers = [], + sample_info = sample_info(), + filename = None, + data_type = 'absorbance', + with_image_cube = False, + with_factorization = False): + self.wavenumbers = wavenumbers + self.N_w = len(wavenumbers) + self.sample_info = sample_info + self._h5_filename = filename + self.xy = np.empty( (0,2) ) + self.data = np.empty( (0,self.N_w) ) + self._N_obs = 0 + self._with_image_cube = with_image_cube + self._with_factorization = with_factorization + + assert data_type in ['transmission', 'reflection', 'absorbance'] + self.data_type = data_type + + if self._h5_filename is not None: # hdf5 mode, load data from hdf5 file + self._mode = 'hdf5' + + self._h5= h5py.File(self._h5_filename,'r') + with self._h5: + self._root = list(self._h5.keys())[0] #get sample root group name + self.sample_info.sample_id = self._h5[self._root+'/info/sample_id'][()] + self.sample_info.sample_meta_data = self._h5[self._root+'/info/sample_meta_data'][()] + self.sample_info.sample_date = self._h5[self._root+'/info/sample_date'][()] + else: # memory mode, data is in memory + self._mode = 'memory' + self._h5 = None # this stays None until we allocate space + +# assert self.sample_info.sample_id != "Unknown", "Missing 'sample_info' keyword parameter" + self._root = str(self.sample_info.sample_id) + '_' + str(self.sample_info.sample_date) + + def add_data(self, spectrum=None, xy=None, ind=[]): + """When working in memory mode, append the data in memory into the ir_map object. + + When working in hdf5 mode, load the data from the hdf5 file. + + Arguments: + + ---------- + spectrum : The 2D spectral data matrix + + xy : An array of xy positions + + ind : An array of indices for selecting specific rows in the 2D spectra matrix. + If it is not given(default), full spectra matrix with all data points are loaded + into the ir_map object, 1D int array + + """ + if self._mode == 'memory': + assert spectrum is not None, "please provide a spectrum matrix" + assert xy is not None, "please provide a xy position array" + + if len(ind) == 0: + self.xy = np.append(self.xy, xy, axis = 0) + self.data = np.append(self.data, spectrum, axis = 0) + else: + self.xy = np.append(self.xy, xy[ind,:], axis = 0) + self.data = np.append(self.data, spectrum[ind,:], axis = 0) + + if self._mode == 'hdf5': + self._h5= h5py.File(self._h5_filename,'r') + with self._h5: + self.wavenumbers = self._h5[self._root+'/data/wavenumbers'][:] + if len(ind) == 0: + self.data = self._h5[self._root+'/data/spectra'][:,:] + self.xy = self._h5[self._root+'/data/xy'][:,:] + else: + self.data = self._h5[self._root+'/data/spectra'][ind,:] + self.xy = self._h5[self._root+'/data/xy'][ind,:] + + def add_image_cube(self, imageCube=None, imageMask=None, image_grid_param=None, ind=[]): + """When working in memory mode, load the image cube data in memory into the ir_map object + and flatten the image cube to 2d spectrum matrix. + + When working in hdf5 mode, load the data from the hdf5 file and flatten the image cube + to 2d spectrum matrix. + + Arguments: + + ---------- + imageCube : The spectral image cube, 3D float array + + imageMask : An image mask where non-blank pixels = True, 2D bool array + + image_grid_param : [x0, y0, dx, dy], 1D float list or array + + ind : An array of indices for selecting wavenumber range. If it is not given(default), + full spectrum are loaded into the ir_map object, 1D int array + """ + self._with_image_cube = True + + if self._mode == 'memory': + assert imageCube is not None, "please provide an image cube" + assert imageMask is not None, "please provide an image mask matrix" + assert image_grid_param is not None, "please provide image grid parameters : [x0, y0, dx, dy]" + + self.imageMask = imageMask + self.image_grid_param = image_grid_param + self.N_y, self.N_x = imageMask.shape[0], imageMask.shape[1] + + if len(ind) == 0:# read in full spectra + self.imageCube = imageCube + else:# read in partial spectra + self.imageCube = imageCube[:,:,ind] + assert len(ind) <= len(self.wavenumbers), "The selected wavenumber indices is longer than the full wavenumber range" + self.wavenumbers = self.wavenumbers[ind] + self.N_w = len(self.wavenumbers) + # convert image cube to 2d data matrix and load the data into self.data, self.xy + self.flatten_image_cube(imageCube, imageMask, image_grid_param) + + if self._mode == 'hdf5': + self._h5= h5py.File(self._h5_filename,'r') + with self._h5: + self.imageMask = self._h5[self._root+'/data/image/image_mask'][:,:] + self.image_grid_param = self._h5[self._root+'/data/image/image_grid_param'][:] + self.wavenumbers = self._h5[self._root+'/data/wavenumbers'][:] + if len(ind) == 0:# read in full spectra + self.imageCube = self._h5[self._root+'/data/image/image_cube'][:,:,:] + else:# read in partial spectrum + self.imageCube = self._h5[self._root+'/data/image/image_cube'][:,:,ind] + assert len(ind) <= len(self.wavenumbers), "The selected wavenumber indices is longer than the full wavenumber range" + self.wavenumbers = self.wavenumbers[ind] + self.N_w = len(self.wavenumbers) + # convert image cube to 2d data matrix and load the data into self.data, self.xy + self.flatten_image_cube(self.imageCube, self.imageMask, self.image_grid_param) + + def add_factorization(self, component=None, component_coef=None, prefix='PCA', ind=[]): + """When working in memory mode, load the factorized components data in memory into the ir_map object. + + When working in hdf5 mode, load the data from the hdf5 file. + + Arguments: + + ---------- + component : The spectral image cube, 3D float array + + component_coef : An image mask where non-blank pixels = True, 2D bool array + + prefix : Name of the factorized components, e.g. PCA, MCR + + ind : An array of indices for selecting specific rows in the 2D spectra matrix. + If it is not given(default), full spectra matrix with all data points are loaded + into the ir_map object, 1D int array + """ + self._with_factorization = True + self._factor_prefix = prefix + '_' + + if self._mode == 'memory': + assert component is not None, "please provide a component matrix" + assert component_coef is not None, "please provide a component_coef matrix" + + self.component = component + self.N_component = component.shape[0] + if len(ind) == 0:# read in all data + self.component_coef = component_coef + else:# read in partial data points + self.component_coef = component_coef[ind,:] + + if self._mode == 'hdf5': + self._h5= h5py.File(self._h5_filename,'r') + with self._h5: + self.component = self._h5[self._root + '/data/factorization/' + self._factor_prefix + 'component'][:,:] + self.N_component = self.component.shape[0] + if len(ind) == 0: # read in all data + self.component_coef = self._h5[self._root + '/data/factorization/' + self._factor_prefix + 'component_coef'][:,:] + else: # read in partial data points + self.component_coef = self._h5[self._root + '/data/factorization/' + self._factor_prefix + 'component_coef'][ind,:] + # check the component dimensions match self.data dimensions + assert self.component.shape[1] == self.data.shape[1], "number of wavenumbers in component does not match that of spectra matrix" + assert self.component_coef.shape[0] == self.data.shape[0], "number of rows in component_coef does not match that of spectra matrix" + + def _allocate_space(self): + + self._h5= h5py.File(self._h5_filename,'w') + data_group = self._h5.create_group( self._root ) + data_group.create_dataset( 'data/xy', + (self._N_obs, 2), + dtype='float32') # we just allocate space + data_group.create_dataset('data/wavenumbers', + data = self.wavenumbers, + dtype='float32') # this we can keep in memory + data_group.create_dataset('data/spectra', + (self._N_obs,self.N_w), + dtype='float32') # we just allocate space + + dt = h5py.special_dtype(vlen=str) + data_group.create_dataset('info/sample_id', + data = self.sample_info.sample_id, + dtype= dt ) + + data_group.create_dataset('info/sample_meta_data', + data = self.sample_info.sample_meta_data, + dtype= dt ) + + data_group.create_dataset('info/sample_date', + data = self.sample_info.sample_date, + dtype= dt ) + if self._with_image_cube: + data_group.create_dataset('data/image/image_cube', + (self.N_y, self.N_x, self.N_w), + dtype='float32') # we just allocate space + data_group.create_dataset('data/image/image_mask', + (self.N_y, self.N_x), + dtype='bool') # we just allocate space + data_group.create_dataset('data/image/ind_rc_map', + (self._N_obs, 3), + dtype='int') # we just allocate space + data_group.create_dataset('data/image/image_grid_param', + data = self.image_grid_param, + dtype='float32') # this we keep in memory + + if self._with_factorization: + data_group.create_dataset('data/factorization/' + self._factor_prefix + 'component', + (self.N_component, self.N_w), + dtype='float32') # we just allocate space + data_group.create_dataset('data/factorization/' + self._factor_prefix + 'component_coef', + (self._N_obs, self.N_component), + dtype='float32') # we just allocate space + + def write_as_hdf5(self, filename): + """Save the object out as an hdf5 file. + + Arguments: + + ---------- + filename : The hdf5 filename where data will be written to. + + """ + # prevent overwriting the existing hdf5 files + assert self._h5_filename != filename, \ + "The given hdf5 filename already exists. Please provide a different filename" + + self._h5_filename = filename + self._N_obs = self.data.shape[0] + self._allocate_space( ) + + with self._h5: + self._h5[self._root + '/data/xy'][:,:] = self.xy + self._h5[self._root + '/data/spectra'][:,:] = self.data + + if self._with_image_cube: #save image cube + self._h5[self._root + '/data/image/image_cube'][:,:,:] = self.imageCube + self._h5[self._root + '/data/image/image_mask'][:,:] = self.imageMask + self._h5[self._root + '/data/image/ind_rc_map'][:,:] = self.ind_rc_map + self._h5[self._root + '/data/image/image_grid_param'][:] = self.image_grid_param + + if self._with_factorization: #save factorization + self._h5[self._root + '/data/factorization/' + self._factor_prefix +'component'][:,:] = self.component + self._h5[self._root + '/data/factorization/' + self._factor_prefix +'component_coef'][:,:] = self.component_coef + + print(f'Data is saved as an HDF5 file. Filename : {filename}') + + def to_image_cube(self, N_x=64, N_y=64, x0=0, y0=0, dx=1, dy=1): + """Transform the spectra matrix to 3D image cube. + + The third dimension of the image cube is the spectra data. + + Arguments: + ---------- + (x0, y0) : The starting location of the ir image. + + (N_x, N_y) : The size of the image. + + (dx, dy) : The step size of the image. + + Returns: + -------- + imageCube : The spectral image cube, 3D float array + + imageMask : An image mask where non-blank pixels = True, 2D bool array + + pointCounts : A mask that shows measurement counts in each pixel, 2D int array + """ + + self._with_image_cube = True + + self.N_x = N_x + self.N_y = N_y + self.image_grid_param = np.array([x0, y0, dx, dy], dtype='float32') + + self.imageCube = np.zeros((self.N_y, self.N_x, self.N_w), dtype='float32') + self.pointCounts = np.zeros((self.N_y, self.N_x), dtype='int') + ind_rc_map = np.zeros((self.xy.shape[0], 3), dtype='int') # ind to row-col mapping [i, row, col] + + x = np.arange(self.N_x)*dx + x0 + y = np.arange(self.N_y)*dy + y0 + + for i in range(self.xy.shape[0]): + ind_rc_map[i, 0] = i + ind_rc_map[i, 1] = val2ind(self.xy[i, 1], y) # align y coordinate to get row + ind_rc_map[i, 2] = val2ind(self.xy[i, 0], x) # align x coordinate to get col + self.pointCounts[ind_rc_map[i, 1], ind_rc_map[i, 2]] += 1 # count how many measurements fall in a pixel + self.imageCube[ind_rc_map[i, 1], ind_rc_map[i, 2], :] += self.data[i,:]# add all spectra that fall in a pixel + + self.imageCube /= np.where(self.pointCounts != 0, self.pointCounts,1)[:,:,np.newaxis]# get average spectra per pixel + self.imageMask = self.pointCounts.astype('bool') # convert to boolean matrix + self.ind_rc_map = ind_rc_map + + return self.imageCube, self.imageMask, self.pointCounts + + def flatten_image_cube(self, imageCube, imageMask, image_grid_param): + """Transform a 3D image cube into a spectra matrix using imageMask to filter out blank data points, + and load the matrix into self.data, load the xy positions of the data points into self.xy + + Arguments: + ---------- + imageCube : The spectral image cube, 3D float array + + imageMask : An image mask where non-blank pixels = True, 2D bool array + + image_grid_param : [x0, y0, dx, dy], 1D float list or array + """ + self.N_x = N_x = imageCube.shape[1] + self.N_y = N_y = imageCube.shape[0] + x0 = image_grid_param[0] + y0 = image_grid_param[1] + dx = image_grid_param[2] + dy = image_grid_param[3] + # set up xy grid and use imageMask to pull out non-blank pixel xy-coordinate + x = np.linspace(x0, x0+dx*(N_x-1), N_x) + y = np.linspace(y0, y0+dy*(N_y-1), N_y) + xv, yv = np.meshgrid(x, y) + xy_grid = np.zeros((N_y, N_x, 2)) + xy_grid[:,:,0] = xv + xy_grid[:,:,1] = yv + + # set up image grid and use imageMask to pull out non-blank pixel row, col position + x = np.arange(N_x) + y = np.arange(N_y) + X, Y = np.meshgrid(x,y) + ind_rc_map = np.zeros((len(X[imageMask]), 3), dtype='int') # ind to row-col mapping [i, row, col] + for i, (r,c) in enumerate(zip(Y[imageMask], X[imageMask])): + ind_rc_map[i,:] = [i, r, c] + + self.xy = xy_grid[imageMask,:] + self.data = imageCube[imageMask,:] + self.ind_rc_map = ind_rc_map + +if __name__ == "__main__": + si = sample_info(sample_id = 'C_elegans') + si.show() + #prespare sample data + np.random.seed(3) + N_wav = 100 + N_obs = 100 + waves = np.linspace(500,4000,N_wav) + data = np.random.uniform(0,1, (N_obs, N_wav) ) + xy = np.random.uniform(-5,5, (N_obs, 2) ) + imageCube = data.reshape(-1, 10, N_wav) + imageMask = np.random.random((N_obs//10, 10)) > 0.5 + x0, y0, dx, dy = 0, 0, 1, 1 + image_grid_param = [x0, y0, dx, dy] + component = np.random.random((3, N_wav)) + component_coef = np.random.random((imageMask.sum(), 3)) + # test loading 2d spectra matrix and writing into hdf5 file + ir_data = ir_map( waves, si) + ir_data.add_data( data, xy ) + ir_data.write_as_hdf5('tst_file.h5') + + # test loading 2d spectra matrix from a hdf5 file + ir_data2 = ir_map(filename ='tst_file.h5' ) + ir_data2.add_data() #load full data matrix + print(ir_data2.data.shape) + ir_data2.add_data(ind=np.arange(20)) #load partial data matrix + print(ir_data2.data.shape) + os.remove('tst_file.h5') + + # test loading image cube data + ir_data3 = ir_map( waves, si) + ir_data3.add_image_cube(imageCube, imageMask, image_grid_param) + assert ir_data3.data.shape[0] == imageMask.sum(), "number of rows in ir_data3.data doesn't match non-blank pixels in imageMask" + r , c= np.where(imageMask) + assert np.all(ir_data3.xy[:,0] == c), "x coordinates in ir_data3.xy doesn't match that of non-blank pixels in imageMask" + assert np.all(ir_data3.xy[:,1] == r), "y coordinates in ir_data3.xy doesn't match that of non-blank pixels in imageMask" + + # loading factorization component + ir_data3.add_factorization(component, component_coef) + + # writing into hdf5 file + ir_data3.write_as_hdf5('tst_file2.h5') + # show hdf5 file dataset structure + lst=[] + with h5py.File('tst_file2.h5','r') as f: + root_name = list(f.keys())[0] + h = f[root_name] + f.visit(lst.append) + ind_rc_map = f[root_name + '/data/image/ind_rc_map'][:,:] + print(*lst, sep='\n') + print(ind_rc_map[:20, :]) + plt.imshow(imageMask) + + # test loading image cube and factorization components from a hdf5 file + ir_data4 = ir_map(filename='tst_file2.h5') + ir_data4.add_image_cube() + ir_data4.add_factorization() + print(ir_data4.data.shape) + print(ir_data4.component.shape) + print(ir_data4.component_coef.shape) + + os.remove('tst_file2.h5') + + print('OK') diff --git a/lbl_ir/lbl_ir/gui_tools/__init__.py b/lbl_ir/lbl_ir/gui_tools/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/lbl_ir/lbl_ir/gui_tools/ispectrum.py b/lbl_ir/lbl_ir/gui_tools/ispectrum.py new file mode 100644 index 0000000..a0f0bd2 --- /dev/null +++ b/lbl_ir/lbl_ir/gui_tools/ispectrum.py @@ -0,0 +1,85 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 22 17:34:27 2019 + +@author: lchen43 +""" +import matplotlib.pyplot as plt +from matplotlib.patches import Circle +from lbl_ir.data_objects.ir_map import val2ind + +class ispectrum: + + def __init__(self, wavenumbers, y): + self.f, self.ax = plt.subplots(figsize=(8,4)) + self.wavenumbers, self.y = wavenumbers, y + self.ax.plot(wavenumbers, y) + plt.xlim([4000,500]) + plt.title('Click on the spectrum to select peaks; right click to cancel selections'); + plt.show() + + self.pos = [] + self.connect() + + def connect(self): + '''connect to all the events + ''' + self.cidpress = self.f.canvas.mpl_connect('button_press_event', self.onclick) + self.cidmotion = self.f.canvas.mpl_connect('motion_notify_event', self.hover) + self.cidenter_axes = self.f.canvas.mpl_connect('axes_enter_event', self.in_axes) + self.cidleave_axes = self.f.canvas.mpl_connect('axes_leave_event', self.leave_axes) + + def disconnect(self): + '''disconnect all the stored connection ids + ''' + self.f.canvas.mpl_disconnect(self.cidpress) + self.f.canvas.mpl_disconnect(self.cidmotion) + self.f.canvas.mpl_disconnect(self.cidenter_axes) + self.f.canvas.mpl_disconnect(self.cidleave_axes) + return self.pos + + def onclick(self, event): + ind = val2ind(event.xdata, self.wavenumbers) + if event.button == 1: + self.pos.append([self.wavenumbers[ind], self.y[ind]]) + self.ax.plot(self.wavenumbers[ind], self.y[ind],'ro') + self.ax.text(self.wavenumbers[ind], self.y[ind]+0.03,str(len(self.pos))) + self.ax.texts[-1], self.ax.texts[-2] = self.ax.texts[-2], self.ax.texts[-1] + else: + if len(self.pos) > 0: + self.pos.pop() + self.ax.lines[-1].remove() + self.ax.texts[-2].remove() + + def hover(self, event): + ind = val2ind(event.xdata, self.wavenumbers) + self.ax.patches[-1].set_center((self.wavenumbers[ind],self.y[ind])) + self.ax.texts[-1].set_position((self.wavenumbers[ind],self.y[ind]+0.08)) + self.ax.texts[-1].set_text(str(self.wavenumbers[ind]) + ', ' +str(self.y[ind])) + + + def in_axes(self, event): + self.ax.texts = [] + self.ax.patches = [] + if event.inaxes: + ind = val2ind(event.xdata, self.wavenumbers) + self.ax.add_patch(Circle((self.wavenumbers[ind], self.y[ind]), radius = 0.05, color = 'r')) + if len(self.pos) > 0: + for i in range(len(self.pos)): + self.ax.text(self.pos[i][0], self.pos[i][1]+0.03, str(i+1)) + self.ax.text(self.wavenumbers[ind], self.y[ind], str(len(self.ax.texts))+ '-' + str(self.wavenumbers[ind]) + ', ' +str(self.y[ind])) + + def leave_axes(self, event): + self.ax.texts = [] + self.ax.patches = [] + +if __name__ == "__main__": + + import os + from ..io_tools.map_IO import read_spa + + test_data_home = '../test_irdata/' + + spa_file = os.path.join(test_data_home, 'test_data0001.spa') + wavenumbers, y, _, _ = read_spa(spa_file) + spec = ispectrum(wavenumbers, y) \ No newline at end of file diff --git a/lbl_ir/lbl_ir/gui_tools/rc2ind_mapping.py b/lbl_ir/lbl_ir/gui_tools/rc2ind_mapping.py new file mode 100644 index 0000000..a283ef4 --- /dev/null +++ b/lbl_ir/lbl_ir/gui_tools/rc2ind_mapping.py @@ -0,0 +1,23 @@ +import numpy as np + + +def getRC2Ind(imgShape, imageMask=None): + """ + In a 2D image, get a dictionary mapping from (row, col) to linear index i of flattened image, and vice versa, + :param imgShape: the shape of the image + :param imageMask: if the image is sparse, use imageMask to generate the mapping + :return: + """ + if imageMask == None: + imageMask = np.ones(imgShape) > 0 + N_x, N_y = imgShape[1], imgShape[0] + x = np.arange(N_x) + y = np.arange(N_y) + X, Y = np.meshgrid(x, y) + ind_rc_map = np.zeros((imgShape[0] * imgShape[1], 3), dtype='int') + for i, (r, c) in enumerate(zip(Y[imageMask], X[imageMask])): + ind_rc_map[i, :] = [i, r, c] + # make a dictionary + ind2rc = {x[0]: tuple(x[1:]) for x in ind_rc_map} + rc2ind = {tuple(x[1:]): x[0] for x in ind_rc_map} + return ind2rc, rc2ind \ No newline at end of file diff --git a/lbl_ir/lbl_ir/io_tools/Omnic_PyMca5/DataObject.py b/lbl_ir/lbl_ir/io_tools/Omnic_PyMca5/DataObject.py new file mode 100644 index 0000000..d330440 --- /dev/null +++ b/lbl_ir/lbl_ir/io_tools/Omnic_PyMca5/DataObject.py @@ -0,0 +1,166 @@ +#/*########################################################################## +# +# The PyMca X-Ray Fluorescence Toolkit +# +# Copyright (c) 2004-2014 European Synchrotron Radiation Facility +# +# This file is part of the PyMca X-ray Fluorescence Toolkit developed at +# the ESRF by the Software group. +# +# 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. +# +#############################################################################*/ +__author__ = "V. Armando Sole - ESRF Data Analysis" +__contact__ = "sole@esrf.fr" +__license__ = "MIT" +__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" +import numpy + +class DataObject(object): + ''' + Simple container of an array and associated information. + Basically it has the members: + info: A dictionnary + data: An array, usually 2D, 3D, ... + + In the past also incorporated selection methods. + Now each different data source implements its selection methods. + + Plotting routines may add additional members + + x: A list containing arrays to be considered axes + y: A list of data to be considered as signals + m: A list containing the monitor data + ''' + GETINFO_DEPRECATION_WARNING = True + GETDATA_DEPRECATION_WARNING = True + SELECT_DEPRECATION_WARNING = True + + def __init__(self): + ''' + Defaut Constructor + ''' + self.info = {} + self.data = numpy.array([]) + + # all the following methods are here for compatibility purposes + # they are obsolete and bound to disappear. + + def getInfo(self): + """ + Deprecated method + """ + if DataObject.GETINFO_DEPRECATION_WARNING: + print("DEPRECATION WARNING: DataObject.getInfo()") + DataObject.GETINFO_DEPRECATION_WARNING = False + return self.info + + def getData(self): + """ + Deprecated method + """ + if DataObject.GETDATA_DEPRECATION_WARNING: + print("DEPRECATION WARNING: DataObject.getData()") + DataObject.GETDATA_DEPRECATION_WARNING = False + return self.data + + def select(self, selection=None): + """ + Deprecated method + """ + if DataObject.SELECT_DEPRECATION_WARNING: + print("DEPRECATION WARNING: DataObject.select(selection=None)") + DataObject.SELECT_DEPRECATION_WARNING = False + dataObject = DataObject() + dataObject.info = self.info + dataObject.info['selection'] = selection + if selection is None: + dataObject.data = self.data + return dataObject + if type(selection) == dict: + #dataObject.data = self.data #should I set it to none??? + dataObject.data = None + if 'rows' in selection: + dataObject.x = None + dataObject.y = None + dataObject.m = None + if 'x' in selection['rows']: + for rownumber in selection['rows']['x']: + if rownumber is None: + continue + if dataObject.x is None: + dataObject.x = [] + dataObject.x.append(self.data[rownumber, :]) + + if 'y' in selection['rows']: + for rownumber in selection['rows']['y']: + if rownumber is None: + continue + if dataObject.y is None: + dataObject.y = [] + dataObject.y.append(self.data[rownumber, :]) + + if 'm' in selection['rows']: + for rownumber in selection['rows']['m']: + if rownumber is None: + continue + if dataObject.m is None: + dataObject.m = [] + dataObject.m.append(self.data[rownumber, :]) + elif ('cols' in selection) or ('columns' in selection): + if 'cols' in selection: + key = 'cols' + else: + key = 'columns' + dataObject.x = None + dataObject.y = None + dataObject.m = None + if 'x' in selection[key]: + for rownumber in selection[key]['x']: + if rownumber is None: + continue + if dataObject.x is None: + dataObject.x = [] + dataObject.x.append(self.data[:, rownumber]) + + if 'y' in selection[key]: + for rownumber in selection[key]['y']: + if rownumber is None: + continue + if dataObject.y is None: + dataObject.y = [] + dataObject.y.append(self.data[:, rownumber]) + + if 'm' in selection[key]: + for rownumber in selection[key]['m']: + if rownumber is None: + continue + if dataObject.m is None: + dataObject.m = [] + dataObject.m.append(self.data[:, rownumber]) + if dataObject.x is None: + if 'Channel0' in dataObject.info: + ch0 = int(dataObject.info['Channel0']) + else: + ch0 = 0 + dataObject.x = [numpy.arange(ch0, + ch0 + len(dataObject.y[0])).astype(numpy.float)] + if not ("selectiontype" in dataObject.info): + dataObject.info["selectiontype"] = "%dD" % len(dataObject.y) + return dataObject diff --git a/lbl_ir/lbl_ir/io_tools/Omnic_PyMca5/OmnicMap.py b/lbl_ir/lbl_ir/io_tools/Omnic_PyMca5/OmnicMap.py new file mode 100644 index 0000000..a5bf3b3 --- /dev/null +++ b/lbl_ir/lbl_ir/io_tools/Omnic_PyMca5/OmnicMap.py @@ -0,0 +1,363 @@ +#/*########################################################################## +# +# The PyMca X-Ray Fluorescence Toolkit +# +# Copyright (c) 2004-2014 European Synchrotron Radiation Facility +# +# This file is part of the PyMca X-ray Fluorescence Toolkit developed at +# the ESRF by the Software group. +# +# 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. +# +# Note: Last modified by Liang Chen 5/14/2019 +#############################################################################*/ +__author__ = "V.A. Sole - ESRF Data Analysis" +__contact__ = "sole@esrf.fr" +__license__ = "MIT" +__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" +import os +import sys +import re +import struct +import numpy +import copy + +from lbl_ir.io_tools.Omnic_PyMca5 import DataObject #modified + +DEBUG = 0 +SOURCE_TYPE = "EdfFileStack" + + +class OmnicMap(DataObject.DataObject): + ''' + Class to read OMNIC .map files + + It reads the spectra into a DataObject instance. + This class info member contains all the parsed information. + This class data member contains the map itself as a 3D array. + ''' + def __init__(self, filename): + ''' + Parameters: + ----------- + filename : str + Name of the .map file. + It is expected to work with OMNIC versions 7.x and 8.x + ''' + DataObject.DataObject.__init__(self) + if sys.platform == 'win32' or 1: #modified, "added or 1" + fid = open(filename, 'rb') + else: + fid = open(filename, 'rb') + data = fid.read() + fid.close() + + try: + omnicInfo = self._getOmnicInfo(data) + except: + omnicInfo = None + self.sourceName = [filename] + if sys.version < '3.0': + searchedChain = "Spectrum " + else: + searchedChain = bytes("Spectrum", 'utf-8') + firstByte = data.index(searchedChain) + s = data[firstByte:(firstByte + 100 - 16)] + if sys.version >= '3.0': + s = str(s) + if DEBUG: + print("firstByte = %d" % firstByte) + print("s1 = %s " % s) + exp = re.compile('(-?[0-9]+\.?[0-9]*)') + tmpValues = exp.findall(s) + spectrumIndex = int(tmpValues[0]) + self.nSpectra = int(tmpValues[1]) + if "X = " in s: + xPosition = float(tmpValues[2]) + yPosition = float(tmpValues[3]) + else: + # I have to calculate them from the scan + xPosition, yPosition = self.getPositionFromIndexAndInfo(0, omnicInfo) + if DEBUG: + print("spectrumIndex, nSpectra, xPosition, yPosition = %d %d %f %f" %\ + (spectrumIndex, self.nSpectra, xPosition, yPosition)) + if sys.version < '3.0': + chain = "Spectrum" + else: + chain = bytes("Spectrum", 'utf-8') + secondByte = data[(firstByte + 1):].index(chain) + secondByte += firstByte + 1 + if DEBUG: + print("secondByte = ", secondByte) + self.nChannels = int((secondByte - firstByte - 100) / 4) + if DEBUG: + print("nChannels = %d" % self.nChannels) + self.firstSpectrumOffset = firstByte - 16 + + #fill the header + self.header = [] + oldXPosition = xPosition + oldYPosition = yPosition + self.nRows = 0 + for i in range(self.nSpectra): + offset = int(firstByte + i * (100 + self.nChannels * 4)) + if sys.version < '3.0': + s = data[offset:(offset + 100 - 16)] + else: + s = str(data[offset:(offset + 100 - 16)]) + tmpValues = exp.findall(s) + spectrumIndex = int(tmpValues[0]) + if "X = " in s: + xPosition = float(tmpValues[2]) + yPosition = float(tmpValues[3]) + else: + #I have to calculate them from the scan + xPosition, yPosition = self.getPositionFromIndexAndInfo(i, omnicInfo) + if (abs(yPosition - oldYPosition) > 1.0e-6) and\ + (abs(xPosition - oldXPosition) < 1.0e-6): + break + self.nRows = self.nRows + 1 + if DEBUG: + print("DIMENSIONS X = %f Y=%d" %\ + ((self.nSpectra * 1.0) / self.nRows, self.nRows)) + + #arrange as an EDF Stack + self.info = {} + self.__nFiles = int(self.nSpectra / self.nRows) + self.data = numpy.zeros((self.__nFiles, self.nRows, self.nChannels), + dtype=numpy.float32) + + self.__nImagesPerFile = 1 + offset = firstByte - 16 + 100 # starting position of the data + delta = 100 + self.nChannels * 4 + fmt = "%df" % self.nChannels + for i in range(self.__nFiles): + for j in range(self.nRows): + # this approach is inneficient when compared to a direct + # data readout, but it allows to deal with nan at the source + tmpData = numpy.zeros((self.nChannels,), dtype=numpy.float32) + tmpData[:] = struct.unpack(fmt,\ + data[offset:(offset + delta - 100)]) + finiteData = numpy.isfinite(tmpData) + self.data[i, j, finiteData] = tmpData[finiteData] + offset = int(offset + delta) + shape = self.data.shape + for i in range(len(shape)): + key = 'Dim_%d' % (i + 1,) + self.info[key] = shape[i] + + self.info["SourceType"] = SOURCE_TYPE + self.info["SourceName"] = self.sourceName + self.info["Size"] = self.__nFiles * self.__nImagesPerFile + self.info["NumberOfFiles"] = self.__nFiles * 1 + self.info["FileIndex"] = 0 + self.info["Channel0"] = 0.0 + if omnicInfo is not None: + self.info['McaCalib'] = [omnicInfo['First X value'] * 1.0, + omnicInfo['Data spacing'] * 1.0, + 0.0] + else: + self.info["McaCalib"] = [0.0, 1.0, 0.0] + self.info['OmnicInfo'] = omnicInfo + + def _getOmnicInfo(self, data): + ''' + Parameters: + ----------- + data : The contents of the .map file + + Returns: + -------- + A dictionnary with acquisition information + ''' + #additional information + fmt = "I" # unsigned long in 32-bit + offset = 372 # 93*4 unsigned integers + infoBlockIndex = (struct.unpack(fmt, data[offset:(offset + 4)])[0] - 204) / 4. + infoBlockIndex = int(infoBlockIndex) + #infoblock is the position of the information block + offset = infoBlockIndex * 4 + #read 13 unsigned integers + nValues = 13 + fmt = "%dI" % nValues + values = struct.unpack(fmt, data[offset:(offset + 4 * nValues)]) + ddict = {} + ddict['Number of points'] = values[0] + ddict['Number of scan points'] = values[6] + ddict['Interferogram peak position'] = values[7] + ddict['Number of sample scans'] = values[8] + ddict['Number of FFT points'] = values[10] + ddict['Number of background scans'] = values[12] + offset = (infoBlockIndex + 3) * 4 + nFloats = 47 + fmt = "%df" % nFloats + vFloats = struct.unpack(fmt, data[offset:(offset + 4 * nFloats)]) + lastX = vFloats[0] + firstX = vFloats[1] + ddict['First X value'] = firstX + ddict['Last X value'] = lastX + ddict['Identifier for start indices of spectra'] = vFloats[14] + ddict['Laser frequency'] = vFloats[16] + ddict['Data spacing'] = (lastX - firstX) / (ddict['Number of points'] - 1.0) + ddict['Background gain'] = vFloats[10] + if DEBUG: + for key in ddict.keys(): + print(key, ddict[key]) + ddict.update(self.getMapInformation(data)) + return ddict + + def getMapInformation(self, data): + ''' + Internal method to help finding spectra coordinates + Parameters: + ----------- + data : Contents of the .map file + + Returns: + -------- + Dictionnary with map gemoetrical acquisition parameters + ''' + #look for the chain 'Position' (this is the old version, modified by Liang Chen 3/11/2019) +# if sys.version < '3.0': +# chain = 'Position' +# else: +# chain = bytes('Position', 'utf-8') +# offset = data.index(chain) +# positions = [offset] +# while True: +# try: +# a = data[(offset + 1):].index(chain) +# offset = a + offset + 1 +# positions.append(offset) +# except ValueError: +# break + + #look for the chain 'micrometers' (this is the new version, modified by Liang Chen 3/11/2019) + if sys.version < '3.0': + chain = 'micrometers' + else: + chain = bytes('micrometers', 'utf-8') + offset = data.index(chain) + x0_trial = struct.unpack('f', data[offset - 88 : offset - 84])[0] + + #look for the chain 'Spectrum ' + if sys.version < '3.0': + searchedChain = "Spectrum " + else: + searchedChain = bytes("Spectrum", 'utf-8') + firstSpectrum = data.index(searchedChain) + s = str(data[firstSpectrum : (firstSpectrum + 100 - 16)]) + exp = re.compile('(-?[0-9]+\.?[0-9]*)') + tmpValues = exp.findall(s) + x0 = float(tmpValues[2]) + + #look for the chain 'GMT' + if sys.version < '3.0': + chain = "GMT" + else: + chain = bytes("GMT", 'utf-8') + offset_year = data.index(chain) + year_collected = int(data[offset_year - 6 :offset_year - 2].decode('utf-8')) + ddict = {} + #map description position (this is the old version, modified by Liang Chen 3/11/2019) +# if (positions[1] - positions[0]) == 66: # reverse engineered magic number :-) +# mapDescriptionOffset = positions[0] - 90 +# mapDescription = struct.unpack('6f', data[mapDescriptionOffset:mapDescriptionOffset + 24]) +# y0, y1, deltaY, x0, x1, deltaX = mapDescription +# ddict['First map location'] = [x0, y0] +# ddict['Last map location'] = [x1, y1] +# ddict['Mapping stage X step size'] = deltaX +# ddict['Mapping stage Y step size'] = deltaY +# ddict['Number of spectra'] = abs((1 + ((y1 - y0) / deltaY)) * (1 + ((x1 - x0) / deltaX))) + + #map description position (this is the new version, modified by Liang Chen 5/14/2019) + if (round(x0_trial) == round(x0)) or (year_collected >= 2018): # if map was collected in 2018 or later, use default: offset - 100 + mapDescriptionOffset = offset - 100 + else: + mapDescriptionOffset = offset - 96 + + mapDescription = struct.unpack('6f', data[mapDescriptionOffset:mapDescriptionOffset + 24]) + y0, y1, deltaY, x0, x1, deltaX = mapDescription + ddict['First map location'] = [x0, y0] + ddict['Last map location'] = [x1, y1] + ddict['Mapping stage X step size'] = deltaX + ddict['Mapping stage Y step size'] = deltaY + ddict['Mapping stage parameters'] = [x0, y0, deltaX, deltaY] + ddict['Number of spectra'] = round(abs((1 + ((y1 - y0) / deltaY)) * (1 + ((x1 - x0) / deltaX)))) + + if DEBUG: + for key in ddict.keys(): + print(key, ddict[key]) + return ddict + + def getOmnicInfo(self): + """ + Returns a dictionnary with the parsed OMNIC information + """ + return copy.deepcopy(self.info['OmnicInfo']) + + def getPositionFromIndexAndInfo(self, index, info=None): + ''' + Internal method to obtain the position at which a spectrum + was acquired + Parameters: + ----------- + index : int + Index of spectrum + info : Dictionnary + Information recovered with _getOmnicInfo + Returns: + -------- + x, y : floats + Position at which the spectrum was acquired. + ''' + if info is None: + return 0.0, 0.0 + ddict = info + #first variation on X and then on Y + try: + x0, y0 = ddict['First map location'] + except KeyError: + return 0.0, 0.0 + x1, y1 = ddict['Last map location'] + deltaX = ddict['Mapping stage X step size'] + deltaY = ddict['Mapping stage Y step size'] + nX = round(1 + ((x1 - x0) / deltaX)) + x = x0 + (index % nX) * deltaX + y = y0 + int(index / nX) * deltaY + return x, y + +if __name__ == "__main__": + filename = None + if len(sys.argv) > 2: + DEBUG = int(sys.argv[2]) + if len(sys.argv) > 1: + filename = sys.argv[1] + elif os.path.exists("SambaPhg_IR.map"): + filename = "SambaPhg_IR.map" + if filename is not None: + w = OmnicMap(filename) + print(type(w)) + print(type(w.data[0:10])) + print(w.data[0:10]) + print("shape = ", w.data.shape) + print(type(w.info)) + print("INFO = ", w.info['OmnicInfo']) + else: + print("Please supply input filename") diff --git a/lbl_ir/lbl_ir/io_tools/Omnic_PyMca5/__init__.py b/lbl_ir/lbl_ir/io_tools/Omnic_PyMca5/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/lbl_ir/lbl_ir/io_tools/__init__.py b/lbl_ir/lbl_ir/io_tools/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/lbl_ir/lbl_ir/io_tools/basic_parser.py b/lbl_ir/lbl_ir/io_tools/basic_parser.py new file mode 100644 index 0000000..e61a347 --- /dev/null +++ b/lbl_ir/lbl_ir/io_tools/basic_parser.py @@ -0,0 +1,142 @@ +import configparser +import sys +import ast +import io + +""" +I put the output from configparser through some helper classes and create objects +with names and key / value pairs as laid out in the input file. +I find this a lot easier to work with then the raw config parser objects. + +""" + +class section_object(object): + def __init__(self, items): + for item in items: + key = item[0] + val = item[1].strip() + try: + val = ast.literal_eval(val.strip()) + except: pass + setattr(self, key, val) + + def as_txt(self): + txt = """""" + keys = [] + for key in self.__dict__.keys(): + if '__' not in key: + keys.append(key) + keys.sort() + for key in keys: + txt += key+'='+str(self.__dict__[key])+' \n' + return txt + + +class config_object(object): + def __init__(self, config): + for section in config.sections(): + items = config.items(section) + this_section_object = section_object(items) + setattr(self, section, this_section_object ) + + def as_txt(self): + txt = """""" + keys = [] + for key in self.__dict__.keys(): + if '__' not in key: + keys.append(key) + for key in keys: + txt+='[%s]\n'%key + txt+=self.__dict__[key].as_txt()+' \n' + return txt + + def show(self,f=None): + if f is None: + f = sys.stdout + print(self.as_txt(), file=f) + + +def read_and_parse(inputs,defaults=None): + # make config object + config = configparser.ConfigParser() + + # not sure how the defaults work + if defaults is not None: + default_config = None + if type(defaults) is str: + default_config = configparser.ConfigParser() + default_config.read_string( defaults ) + #default_config.readfp(io.BytesIO(defaults)) + else: + default_config = configparser.ConfigParser() + default_config.read_file(defaults) + # if this aint a config object, we're f-ed anyway + config._sections = default_config._sections + + if inputs is None: + co = config_object( default_config ) + co.show() + raise SystemExit('##--- No inputs provide, use template shown above ---##') + + + if type(inputs) is type(""""""): + config.read_string(inputs) #inputs = inputs + else: + # update the config object accoring to the specified input + config.read_file(inputs) + + # build the object. + co = config_object( config ) + return co + +def tst(): + default_instructions = """ +[data] +experiment = None +run = None +index_start = 0 +index_stop = 500 +index_stride = 1 + +[output] +filename = output.h5 +comments = "No comments" +""" + + instructions = """ +[data] +experiment = amox26916 +run = 56 +index_start = 0 +index_stop = 5000 + +[output] +filename = output_59.h5 +comments = "No soup for you!" +ooops = 9 +""" + + obj = read_and_parse(instructions, default_instructions) + + instructions=""" +[stuff] +coordinates = (1,4) +""" + obj = read_and_parse(instructions) + assert obj.stuff.coordinates[0] == 1 + assert obj.stuff.coordinates[1] == 4 + obj.show() + + can_we_do_comments = """ +[breakfast] +eggs = True +bacon = True +spam = False # i don't like spam +""" + obj = read_and_parse(can_we_do_comments) + obj.show() + print('OK') + + +if __name__ == "__main__": + tst() diff --git a/lbl_ir/lbl_ir/io_tools/map_IO.py b/lbl_ir/lbl_ir/io_tools/map_IO.py new file mode 100644 index 0000000..4b588b8 --- /dev/null +++ b/lbl_ir/lbl_ir/io_tools/map_IO.py @@ -0,0 +1,249 @@ +""" +Created on Tue Dec 18 16:14:43 2018 + +@author: Liang Chen +""" + +import spectral.io.envi as envi +import numpy as np +import re + + +def read_envi(hdr_file): + """Load an ENVI map from the .hdr header file + + Parameters: + ----------- + hdr_file: string + file path for a .hdr file + + Returns: + -------- + out : object + a spectral.io.bipfile.BipFile object + + Examples: + --------- + >>> hdr_file = os.path.join(test_data_home, 'test_envi.hdr') + >>> img = envi.open(hdr_file) + (17, 32, 1738) + """ + + # read in parameter list in the header file + params = [] + + with open(hdr_file, 'r', encoding='utf-8', errors='ignore') as header: + for line in header: + params.append(line.split('=')[0].strip()) + + # if parameter 'byte order' is not found, append 'byte order = 0' to the header file + if 'byte order' not in params: + with open(hdr_file, 'a', encoding='utf-8', errors='ignore') as header: + header.write('byte order = 0\n') + + img = envi.open(hdr_file) + + return img + + +def read_binary(fileObj, byteType='uint8', size=1): + """A helper function to readin values from a binary file + + Parameters: + ----------- + fileObj : object + a binary file object + + byteType : string, optional, default 'uint8' + the type of readin values + + size : int, optional + the number of bytes to readin. Default is 1. + + Returns: + -------- + out : a value or a tuple of values or a string + the readout value from the bytes + """ + import struct + + typeNames = { + 'int8': ('b', struct.calcsize('b')), + 'uint8': ('B', struct.calcsize('B')), + 'int16': ('h', struct.calcsize('h')), + 'uint16': ('H', struct.calcsize('H')), + 'int32': ('i', struct.calcsize('i')), + 'uint32': ('I', struct.calcsize('I')), + 'int64': ('q', struct.calcsize('q')), + 'uint64': ('Q', struct.calcsize('Q')), + 'float': ('f', struct.calcsize('f')), + 'double': ('d', struct.calcsize('d')), + 'char': ('s', struct.calcsize('s'))} + + if size == 1: + return struct.unpack(typeNames[byteType][0], fileObj.read(typeNames[byteType][1]))[0] + elif size > 1: + return struct.unpack(typeNames[byteType][0] * size, fileObj.read(typeNames[byteType][1] * size)) + else: + return None + + +def read_spa(spa_file): + """Load a spectrum from a .spa file + + Parameters: + ----------- + spa_file: string + file path for a .spa file + + Returns: + -------- + wavenumbers : float array + the wavenumber values of the spectrum (x-axis) + + spectrum : float array + the transmission/reflection/absorption coefficient + at wavenumber of the spectrum (y-axis) + + title : string + the title of the spectrum + + comment : string + the comment section of the spectrum, if it exists + + Examples: + --------- + >>> spa_file = os.path.join(test_data_home, 'test_data0001.spa') + >>> wavenumbers, spectrum, title, comment = read_spa(spa_file) + >>> print(title) + C:\\Users\\lchen43\\Documents\\CDIPS_2017\\lbl-ir\\test_data\\test_data.map - Spectrum #1 + Position (X,Y): 499.51, 5715.69 Thu Apr 22 13:57:44 2010 (GMT-07:00) + >>> print(comment) + Split map from: C:\\Users\\lchen43\\Documents\\CDIPS_2017\\lbl-ir\\test_data\\test_data.map + X Range: 499.51, 654.51 + Y Range: 5715.69, 5795.69 + Position (X,Y): 499.51, 5715.69 + Thu Apr 22 13:57:44 2010 (GMT-07:00) + """ + + with open(spa_file, 'rb') as f: + f.seek(30) + readChar = read_binary(f, 'uint8', 255) + title = ''.join([chr(i) for i in filter(lambda x: x > 0, readChar)]) + + f.seek(564) + spectrumPts = read_binary(f, 'int32') + + f.seek(576) + maxWavenum = read_binary(f, 'float') + minWavenum = read_binary(f, 'float') + wavenumbers = np.linspace(maxWavenum, minWavenum, spectrumPts) + + # The starting byte location of the spectrum data is stored in the + # header. It immediately follows a flag value of 3. + # If there is a comment section, the flag value is 27 and it should be in front of the flag value of 3. + f.seek(338) + Flag = 0 + while (Flag != 3 and Flag != 27): + Flag = read_binary(f, 'uint16') + + if Flag == 3: # no comment section, look for data starting position + dataPosition = read_binary(f, 'uint16') + f.seek(dataPosition) + spectrum = read_binary(f, 'float', spectrumPts) + comment = '' + elif Flag == 27: # there is a comment section, look for comment starting position + commentPosition = read_binary(f, 'uint16') + # move forward 14 bytes, that's the data starting position + f.seek(14, 1) + dataPosition = read_binary(f, 'uint16') + commentLength = dataPosition - commentPosition + + f.seek(commentPosition) + readChar = read_binary(f, 'uint8', commentLength) + comment = ''.join([chr(i) + for i in filter(lambda x: x > 0, readChar)]) + f.seek(dataPosition) + spectrum = read_binary(f, 'float', spectrumPts) + else: + raise Exception('The flag value of 3 or 27 cannot be found') + + spectrum = np.array(spectrum) + return wavenumbers, spectrum, title, comment + + +def read_series(file_name, wavLen=1738): + """ + read Ominc series map + :param file_name: path of minc series map + :param wavLen: the length of the wavenumbers vector + :return: + wav: the wavenumbers vector + spectra: all spectra in the series map + xy: all xy coordinate in the series map + """ + + wav = spectra = xy = None + + with open(file_name, 'rb') as fid: + data = fid.read() + # read out firstWav, lastWav, construct wav + fid.seek(0) + s = [] + for i in range(1000): + s.append((read_binary(fid, 'uint32'))) + offset = s.index(wavLen) * 4 + fid.seek(offset + 12) + val = read_binary(fid, 'float', 2) + firstWav, lastWav = val[1], val[0] + wav = np.linspace(firstWav, lastWav, wavLen) + + # read out num of spectra + Chain = bytes('Spectrum', 'utf-8') + firstByte = data.index(Chain) + s1 = str(data[firstByte:(firstByte + 100 - 16)]) + exp = re.compile('(-?[0-9]+\.?[0-9]*)') + tmpValues = exp.findall(s1) + nSpectra = int(tmpValues[1]) + + # read out all spectra + spectra = np.zeros((nSpectra, wavLen)) + delta = wavLen * 4 + 96 + for i in range(nSpectra): + fid.seek(firstByte + delta * i + 80) + spectra[i, :] = read_binary(fid, 'float', wavLen) + + # find xy positions + chain = bytes('Position', 'utf-8') + firstPos = data.index(chain) + secondPos = data[(firstPos + 1):].index(chain) + secondPos += firstPos + 1 + fid.seek(secondPos + 48) + val = read_binary(fid, 'float', 2 * nSpectra) + xy = np.zeros((nSpectra, 2)) + xy[:, 0], xy[:, 1] = val[::2], val[1::2] + + return wav, spectra, xy + + +if __name__ == "__main__": + import os + + test_data_home = '../../test_irdata/' + + hdr_file = os.path.join(test_data_home, 'test_envi.hdr') + img = read_envi(hdr_file) + print('====Envi file====') + print(img.shape) + + spa_file = os.path.join(test_data_home, 'test_data0001.spa') + wavenumbers, spectrum, title, comment = read_spa(spa_file) + print('====' + title + '====') + print(comment) + + map_file = os.path.join(test_data_home, 'typeII-010_12x9.map') + wavenumbers, spectra, xy = read_series(map_file) + print('====Series map====') + print(spectra.shape) + print(xy[:3,:]) + print(spectra[:4, -3:]) diff --git a/lbl_ir/lbl_ir/io_tools/read_XAS.py b/lbl_ir/lbl_ir/io_tools/read_XAS.py new file mode 100644 index 0000000..71e7aec --- /dev/null +++ b/lbl_ir/lbl_ir/io_tools/read_XAS.py @@ -0,0 +1,54 @@ +import os +import numpy as np +import h5py +from lbl_ir.data_objects import ir_map + + +def get_grid_info(coords): + xsorted = sorted(set(coords[:,0])) + ysorted = sorted(set(coords[:,1])) + x0, xmax, y0, ymax = xsorted[0], xsorted[-1], ysorted[0], ysorted[-1] + dx = np.mean(np.diff(xsorted)) + dy = np.mean(np.diff(ysorted)) + if dx > dy: + step = dy / 2 + else: + step = dx / 2 + Nx = int(round((xmax - x0)/step) + 1) + Ny = int(round((ymax - y0)/step) + 1) + return x0, y0, step, Nx, Ny + +def read_xasH5(filePath): + xasTypes = [] + energy = {} + dataSets = {} + coords = {} + xas_maps = {} + + with h5py.File(filePath, 'r') as f: + xasSpectra = f['xas/'] + for k in xasSpectra: + xasTypes.append(k) + for k1 in xasSpectra[k]: + k1 = '/' + k1 + spectra = xasSpectra[k + k1 + '/raw'][:, :, :] + energy[k] = spectra[0, 0, :] + dataSets[k] = spectra[:, 1, :] + dataSets[k] = np.where(dataSets[k] != np.inf, dataSets[k], 0) #filter np.inf + + samples = f['maps/samples/'] + for i, k in enumerate(samples): + coords[xasTypes[i]] = samples[k + '/xas_coords'][:, :] + + for _type in xasTypes: + assert coords[_type].shape[0] == dataSets[_type].shape[0], 'xas and coords sample sequences were mis-aligned.' + + fileName = os.path.basename(filePath) + sample_info = ir_map.sample_info(fileName[:-3]) + xas_maps[_type] = ir_map.ir_map(wavenumbers=energy[_type], sample_info=sample_info) + xas_maps[_type].add_data(spectrum=dataSets[_type], xy=coords[_type]) + x0, y0, step, Nx, Ny = get_grid_info(coords[_type]) + xas_maps[_type].to_image_cube(Nx, Ny, x0, y0, step, step) + + return xas_maps + diff --git a/lbl_ir/lbl_ir/io_tools/read_map.py b/lbl_ir/lbl_ir/io_tools/read_map.py new file mode 100644 index 0000000..0c74394 --- /dev/null +++ b/lbl_ir/lbl_ir/io_tools/read_map.py @@ -0,0 +1,47 @@ +import h5py +import sys + +from lbl_ir.io_tools import read_omnic +from lbl_ir.data_objects import ir_map + +def read_all_formats(filename, sample_info=None): + + ok = False + format = None + if not ok: + try: + data = read_omnic.read_and_convert(filename, sample_info=None) + ok= True + format = "Omnic" + except: pass + + if not ok: + try: + data = ir_map.ir_map(filename=filename) + ok = True + format = "hdf5" + + with h5py.File(filename,'r') as f: + root_name = list(f.keys())[0] + # if there is an image group, load imagecube and data, otherwise load data group + if 'image' in f[root_name + '/data']: + data.add_image_cube() + else: + data.add_data() + # if there is an factorization group, load factorization data + if 'factorization' in f[root_name + '/data']: + data.add_factorization() + except: pass + + + + if not ok: + print("Could not read file; Check input or file formats and header integrity.") + else: + return data, format + + + +if __name__ == "__main__": + data,fmt = read_all_formats(sys.argv[1]) + print("Data read in with format", fmt) diff --git a/lbl_ir/lbl_ir/io_tools/read_numpy.py b/lbl_ir/lbl_ir/io_tools/read_numpy.py new file mode 100644 index 0000000..d6d8c34 --- /dev/null +++ b/lbl_ir/lbl_ir/io_tools/read_numpy.py @@ -0,0 +1,42 @@ +import numpy as np +from lbl_ir.data_objects import ir_map + + +def read_npy(filename, wavenumbers=None, data_type="absorbance", sample_info=None): + if sample_info is None: + sample_info = ir_map.sample_info() + + data = np.load(filename) + + if wavenumbers is None: + wavenumbers = np.arange(data.shape[2]) + + image_grid_param = [0, 0, 1, 1] + image_mask = np.ones(data.shape[0:2]) > 0.5 + + this_ir_map = ir_map.ir_map(wavenumbers=wavenumbers, + sample_info=sample_info, + data_type=data_type) + this_ir_map.add_image_cube(data, image_mask, image_grid_param) + return this_ir_map + +def read_npz(filename, data_type="absorbance", sample_info=None): + if sample_info is None: + sample_info = ir_map.sample_info() + + npzFile = np.load(filename) + + for k in npzFile.files: + if 'energy' in k: + wavenumbers = npzFile[k] + else: + data = npzFile[k] + + image_grid_param = [0, 0, 1, 1] + image_mask = np.ones(data.shape[0:2]) > 0.5 + + this_ir_map = ir_map.ir_map(wavenumbers=wavenumbers, + sample_info=sample_info, + data_type=data_type) + this_ir_map.add_image_cube(data, image_mask, image_grid_param) + return this_ir_map \ No newline at end of file diff --git a/lbl_ir/lbl_ir/io_tools/read_omnic.py b/lbl_ir/lbl_ir/io_tools/read_omnic.py new file mode 100644 index 0000000..e64f649 --- /dev/null +++ b/lbl_ir/lbl_ir/io_tools/read_omnic.py @@ -0,0 +1,51 @@ +import numpy as np +import sys + +from lbl_ir.io_tools.Omnic_PyMca5 import OmnicMap +from lbl_ir.data_objects import ir_map + + +def read_and_convert(filename, start_wav=None, stop_wav=None, data_type="absorbance", sample_info=None): + if sample_info is None: + sample_info = ir_map.sample_info() + + omnic_object = OmnicMap.OmnicMap( filename ) + wavenumbers = None + if wavenumbers is None: + n_wav = omnic_object.data.shape[2] + if (start_wav is None) or ( stop_wav is None) : + if omnic_object.info['OmnicInfo'] is not None: + start_wav = omnic_object.info['OmnicInfo']['First X value'] + stop_wav = omnic_object.info['OmnicInfo']['Last X value'] + n_wav = omnic_object.info['OmnicInfo']['Number of points'] + else: + raise ValueError('There is an issue with the Ominc file or its parser. Make sure that the meta info from the header is present and parsed correctly.') + + wavenumbers = np.linspace( start_wav, stop_wav, n_wav ) + image_grid_param = omnic_object.info['OmnicInfo']['Mapping stage parameters'] + image_mask = np.ones( omnic_object.data.shape[0:2]) > 0.5 + + # build the basis object + this_ir_map = ir_map.ir_map( wavenumbers = wavenumbers, + sample_info = sample_info, + data_type = data_type + ) + this_ir_map.add_image_cube(omnic_object.data, image_mask, image_grid_param) + return this_ir_map + +if __name__ == "__main__": + + map_name = '190519_N2_L2w1_mp2' + '.map' + file_name = '../../ir_data/' + map_name + sample_id = file_name + sample_id = sample_id.replace("../","") + sample_id = sample_id.replace(".map","") + sample_id = sample_id.replace("/","_") + sample_id = sample_id.replace(" ","_") + print(sample_id) + + sample_info = ir_map.sample_info( sample_id = sample_id, sample_meta_data='Hello World') + data = read_and_convert(file_name, sample_info=sample_info) + data.write_as_hdf5('../../ir_data/' + map_name + '.h5') + + diff --git a/lbl_ir/lbl_ir/math_tools/__init__.py b/lbl_ir/lbl_ir/math_tools/__init__.py new file mode 100644 index 0000000..b28b04f --- /dev/null +++ b/lbl_ir/lbl_ir/math_tools/__init__.py @@ -0,0 +1,3 @@ + + + diff --git a/lbl_ir/lbl_ir/math_tools/batched_SVD.py b/lbl_ir/lbl_ir/math_tools/batched_SVD.py new file mode 100644 index 0000000..70ca0f4 --- /dev/null +++ b/lbl_ir/lbl_ir/math_tools/batched_SVD.py @@ -0,0 +1,467 @@ +from scipy.sparse.linalg import svds, eigsh +from scipy.linalg.interpolative import svd as isvd + + +import os +import numpy as np +import matplotlib.pyplot as plt +import time + +from mpi4py import MPI +import h5py + +from tqdm import tqdm + + +def test_function( N, M , noise=0.001): + """ + Generate some test data we can use. The resulting matrix should have a rank of 4. + """ + + x = np.linspace(-1,1,M) + p0 = x*0+1.0 + p1 = x + p2 = 0.5*(3*x*x-1) + p3 = 0.5*(5*x*x*x-3*x) + + result = [] + for ii in range(N): + tmp = np.random.uniform(-1,1,4) + tmp = tmp[0]*p0 + tmp[1]*p1 + tmp[2]*p2 + tmp[3]*p3 + tmp = tmp + np.random.normal(0,1.0,M)*noise + result.append(tmp) + result = np.vstack(result) + + return result + +def invert_permutation(p): + """Given a permutation, provide an array that undoes the permutation. + """ + s = np.empty(p.size, p.dtype) + s[p] = np.arange(p.size) + return s + + +class batched_SVD(object): + """Compute an SVD of all data, but in small batches to overcome memory issues. + + Parameters: + ---------- + + data: A pointer to a data object, say a N,M matrix. + N experimental observations of dimension M + + N_max: The maximum number of entries in sub block of data + + k_singular: The number of singular values to consider + + randomize : A flag which determines if the data will be split in a random fashion. True by default + + + Attributes: + ----------- + self.order : The order in which the data will be examined + + self.inv_order : The inverse of the above array + + self.N_split : The number of batches of data + + self.parts : A list of selection arrays + + self.partial_svd_u : A list of (truncated) svd matrices (U) from individual batches of data + + self.partial_svd_s : A list of (truncated) svd matrices (S) from individual batches of data + + self.partial_svd_vt : A list of (truncated) svd matrices (V^T) from individual batches of data + + self.partial_bases : A list of (truncated) svd matrices (SV^T) from individual batches of data + + + Examples: + --------- + + data = test_function(10000,3000,1.1) + bSVD = batched_SVD(data, 2000, 5, randomize=True) + u,s,vt = bSVD.go_svd() + + """ + + def __init__(self, data, N_max, k_singular, randomize=True): + self.data = data + self.N_max = N_max + self.k_singular = k_singular + + self.randomize = randomize + + self.order = np.arange( self.data.shape[0] ) + if self.randomize: + np.random.shuffle( self.order ) + + self.N_split = int( np.floor(self.data.shape[0] / self.N_max) ) +1 + self.parts = np.array_split( self.order, self.N_split) + + # if we have these numbers in order, we can use them as slices in a hdf5 setting + # this is only requiered when we randomize the lot + if self.randomize: + tmp = [] + for part in self.parts: + part = np.sort(part) + tmp.append(part) + self.parts = tmp + self.order = np.concatenate(self.parts) + self.inv_order = invert_permutation(self.order) + + # here the partial svds are stored + self.partial_svd_u = [] + self.partial_svd_s = [] + self.partial_svd_vt = [] + self.partial_bases = [] + + + def SVD_on_chunk(self, this_chunk): + """Perform an SVD on a subset of the data + + Parameters: + ----------- + + this_chunk: a list of indices that maps back to the self.data array + """ + + tmp_data = self.data[ this_chunk, : ] + + u,s,v = isvd(tmp_data.astype('float64'), self.k_singular) + vt = v.transpose() + self.partial_svd_u.append( u ) + self.partial_svd_s.append( s ) + self.partial_svd_vt.append( vt ) + self.partial_bases.append( np.diag(s).dot(vt) ) + + + def get_u_given_basis(self, sigma, v_transpose, this_chunk): + inv_bases = v_transpose.transpose().dot( np.diag( 1.0 / sigma ) ) + tmp_data = self.data[ this_chunk, : ] + new_u = tmp_data.dot( inv_bases ) + return new_u + + def go_svd(self): + """ Seperate SVD's for individual data chunks are computed and combined into a single, best estimate + of the matrix S and V^T. The data is subsequently revisited to get a new estimate of the matrix U. + + The estimate of U on the basis of the chunked SVD approach is possible, but not as accurate. + """ + + for this_chunk in self.parts: + self.SVD_on_chunk( this_chunk ) + + # do an SVD on the SVD results of the individual chunks. + tmp_bases = np.vstack( self.partial_bases ) + ub,sb,vb = isvd(tmp_bases.astype('float64'), self.k_singular) + vbt = vb.transpose() + ubs = np.array_split(ub,self.N_split) + new_s = sb + new_vt = vbt + + # revisit the data to get the U matrix again + new_us = [] + for part in self.parts: + this_new_u = self.get_u_given_basis(new_s, new_vt, part) + new_us.append(this_new_u) + + # stack it up and unmix the data + new_us = np.vstack(new_us) + if self.randomize: + new_us = new_us[ self.inv_order ] + return new_us, sb, vbt + + +class parallel_SVD_MPI(object): + """ + THSI NEEDS WORK!!! + An MPI based version of an SVD method, splitting the data among a number of cores. + The outline is a bit different than the batched_SVD version. + + The idea is to first split all data equally across all cores. In each core, we compute + the an SVD using the procedure coded up for the batched_SVD approach, and end up with a + SVD estimate on each core. We subsequently pool all these SVD's and reestimate it on + the root node, and scatter these functions back to all subsequent nodes. Subsequently, + we need to compute U on on each node for all data associated with this node. The final + results are stored in a pointer, which is either a hdf5 file or a numpy array. + + + Parameters: + ---------- + + data: A pointer to a data object, say a N,M matrix. + N experimental observations of dimension M + + N_max: The maximum number of entries in sub block of data + + k_singular: The number of singular values to consider + + MPI_COMM_WORLD: an instance of MPI.COMM_WORLD + + randomize : A flag which determines if the data will be split in a random fashion. True by default + + Attributes: + ----------- + + go_svd() : this function performs the svd and returns the SVD results for each single core, including + an array that allows one to map the data back to the original order in which it was presented. + + Example: + -------- + + + + + + """ + + + def __init__(self, data, N_max, k_singular, MPI_COMM_WORLD, randomize=True, selection = None ): + # randomize the order in which we analyze the data + self.randomize = randomize + + # sort out MPI stuff + self.mpi_comm = MPI_COMM_WORLD + self.mpi_rank = self.mpi_comm.Get_rank() + self.mpi_size = self.mpi_comm.Get_size() + + # data etc + self.data = data # the data + self.k_singular = k_singular # the number of singular vectors + self.N_max = N_max # the maximum number of data points per chunk per core + + # split the data amond cores + self.N_obs, self.N_dim = self.data.shape + if selection is not None: + self.Nobs = len(selection) + self.order = np.arange(self.N_obs) + self.selection = selection + if self.selection is not None: + self.order = self.order[selection] + self.inv_order = None + + if self.randomize: + np.random.shuffle( self.order ) + self.rank_splits = np.array_split( self.order, self.mpi_size ) + + if self.randomize: + tmp = [] + for part in self.rank_splits: + part = np.sort(part) # we want this ordered from low to high to enable slicing an hdf5 array + tmp.append(part) + self.rank_splits = tmp + self.order = np.concatenate(self.rank_splits) + self.inv_order = invert_permutation(self.order) + self.inv_order_split = np.array_split( self.inv_order, self.mpi_size ) + self.mpi_comm.Barrier() + + def go_svd(self): + """Do the SVD ihn each core, on several chunks. + """ + + partial_u = [] + partial_s = [] + partial_vt = [] + partial_bases = [] + + rank_selection = self.rank_splits[ self.mpi_rank ] + N_chunks = int( len(rank_selection) / self.N_max ) + 1 + chunks = np.array_split( rank_selection, N_chunks ) + u = None + s = None + vt = None + this_rank_bases = None + for chunk in tqdm(chunks, position=self.mpi_rank): + partial_data = self.data[ chunk, : ] + u,s,vt = isvd(partial_data.astype('float64'), self.k_singular) + partial_u.append( u ) + partial_s.append( s ) + partial_vt.append( vt ) + + partial_bases.append( np.diag(s).dot(vt) ) + # now that we have the partial svd results, we bnring stuff together + if N_chunks > 1: + all_bases = np.vstack( partial_bases ) + uc,sc,vct = isvd(all_bases.astype('float64'), self.k_singular) + # now we need to pass this guy to the main rank + this_rank_bases = np.diag(sc).dot(vct) + else: + this_rank_bases = np.diag(s).dot(vt) + self.mpi_comm.Barrier() + + + gathered_bases = None + if self.mpi_rank == 0: + gathered_bases = np.empty( [self.mpi_size, self.k_singular, self.N_dim], dtype='d' ) + gathered_bases = self.mpi_comm.gather( this_rank_bases, root=0 ) + + final_sg = np.zeros( [ self.k_singular ], dtype='float32' ) + final_vgt = np.zeros( [ self.k_singular, self.N_dim], dtype='float32' ) + + if self.mpi_rank == 0: + gathered_bases = np.vstack( gathered_bases ) + ug,final_sg,final_vgt = isvd(gathered_bases.astype('float64'), self.k_singular) + # we now need to scatter back the sg and vgt matrices + self.mpi_comm.Barrier() + self.mpi_comm.Bcast( final_sg, root=0 ) + self.mpi_comm.Barrier() + self.mpi_comm.Bcast( final_vgt , root=0 ) + self.mpi_comm.Barrier() + # now that we have the final and best sigma and Vt, we need to go back to the data and + # reestimate the the U matrices + inv_multi = final_vgt.transpose().dot( np.diag(final_sg) ) + chunks_of_u = [] + for chunk in chunks: + partial_data = self.data[ chunk, : ] + this_u = partial_data.dot( inv_multi ) + chunks_of_u.append( this_u) + + chunks_of_u = np.vstack( chunks_of_u ) + # we need to return this, including a placement array + return chunks_of_u, final_sg, final_vgt, self.inv_order_split + + + + + + +def tst_batched(): + N = 10000 + M = 2000 + K = 200 + P = 4 + + data = test_function(N,M,0.0001) + print("Data constructed") + e0 = time.time() + u,s,vt = isvd(data.astype('float64'), P) + vt = vt.transpose() + e1 = time.time() + + bSVD = batched_SVD(data, K, P, randomize=False) + e4 = time.time() + us, ss, vst = bSVD.go_svd() + e5 = time.time() + + assert np.std( (s-ss)/ss ) < 1e-3 + print("Singular values match between batch and full approach") + + # checking the reconstructions + da = us.dot(np.diag(ss).dot( vst )) + dc = u.dot(np.diag(s).dot( vt )) + delta_both = np.std( (da-dc) ) + assert np.abs( delta_both ) < 1e-2 + + + print('Reconstruction Error is similar in batched versus full approach') + print ('Time for full: %4.2f batched: %4.2f'%(e1-e0, e5-e4)) + + # Now we want to do this using a random order or data + + bSVD = batched_SVD(data, K, P, randomize=True) + e4 = time.time() + us, ss, vst = bSVD.go_svd() + e5 = time.time() + + assert np.std( (s-ss)/s ) < 1e-3 + print("Singular values match between randomized batch and full approach") + delta_both = np.std( (dc - da) ) + assert np.abs( delta_both ) < 1e-2 + print('Reconstruction Error is decent') + print ('Time for full: %4.2f batched: %4.2f'%(e1-e0, e5-e4)) + + +def tst_MPI(): + N = 10000 # number of observations + M = 2000 # dimension of an observations + P = 4 # number of singular values + + # make the data + mpi_comm = MPI.COMM_WORLD + if mpi_comm.Get_rank() == 0: + data = test_function(N,M,0.0001) + + e0 = time.time() + u,s,vt = isvd(data.astype('float64'),P) + e1 = time.time() + print('Standard svd takes %12.3f seconds'%(e1-e0)) + print('\n') + + e0 = time.time() + # lets quickly see if this works in the batched setup + bSVD = batched_SVD(data, 1000, P, randomize=True) + ub, sb, vbt = bSVD.go_svd() + e1 = time.time() + print('Single core batched version with data in memory takes %12.3f seconds'%(e1-e0) ) + + + # write this to an hdf5 file + print('Creating h5 data file') + f = h5py.File('test_data.h5','w') + dset = f.create_dataset("data", data=data, dtype='float32') + f.close() + del data + + # now we read the data + f = h5py.File('test_data.h5','r') + data = f['/data'] + print(data.shape) + e0 = time.time() + # lets quickly see if this works in the batched setup + bSVD = batched_SVD(data, 1000, P, randomize=False) + ub, sb, vbt = bSVD.go_svd() + e1 = time.time() + print('Single core batched version while reading data from HDF5 takes %12.3f seconds'%( e1-e0 ) ) + print(sb) + f.close() + + mpi_comm.Barrier() + f = h5py.File('test_data.h5','r' ) # lets see if we can get away with this + data = f['data'] + e2 = time.time() + print('build it') + mSVD = parallel_SVD_MPI( data, 1000, P, mpi_comm, False, None) + print('go') + u,s,v,sel = mSVD.go_svd() + e3 = time.time() + print('Core %i with batches version takes %12.3f seconds, while reading data from hdf5'%( mpi_comm.rank , e3-e2) ) + print(s, mpi_comm.rank ) + + #lets read the data into memory + print('Reading data into memory') + data2 = f['data'].value #[:,:] + print(data2.shape) + print( 'done') + + mpi_comm.Barrier() + e2 = time.time() + mSVD = parallel_SVD_MPI( data2, 500, P, mpi_comm, True) + u,s,v,sel = mSVD.go_svd() + e3 = time.time() + print('Core %i with batches version takes %12.3f seconds, with data in memory'%(mpi_comm.rank , e3-e2) ) + + mpi_comm.Barrier() + selection = np.arange(1000) + print('Testing setup with an included selection array') + data2[1000,:]=0. + mSVD = parallel_SVD_MPI( data2, 50000, P, mpi_comm, False, selection = selection) + SVDs = batched_SVD.batched_SVD( ) + up,sp,vp, order_split = mSVD.go_svd() + uf,sf,vf = isvd(data2[0:1000,:].astype('float64'), P) + assert np.mean( ( np.abs(sp-sf)/sf ) ) < 2e-2 + + mpi_comm.Barrier() + if mpi_comm.Get_rank() == 0: + # now remove this file again + print('Removing h5 data file') + os.remove('test_data.h5') + # done + + + + + +if __name__ == "__main__": + tst_batched() diff --git a/lbl_ir/lbl_ir/simulations/__init__.py b/lbl_ir/lbl_ir/simulations/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/lbl_ir/lbl_ir/simulations/spectra_map_simulator.py b/lbl_ir/lbl_ir/simulations/spectra_map_simulator.py new file mode 100644 index 0000000..091bda6 --- /dev/null +++ b/lbl_ir/lbl_ir/simulations/spectra_map_simulator.py @@ -0,0 +1,287 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Jan 2 16:14:47 2019 + +@author: Liang Chen +""" +import numpy as np +import matplotlib.pyplot as plt +from lbl_ir.data_objects.ir_map import sample_info, ir_map + +def lorentzian(x,x0,gamma=10): + + return 1/(np.power((x-x0)/gamma,2)+1) + +def gaussian(x,x0=0,sigma=10): + + return 1/np.exp(np.power((x-x0)/sigma,2)) + +class spectra_map_simulator: + + """Generate a simulated spectral map + + Parameters: + ----------- + NbaseSpectra: int, optional + The number of basis spectrum. Default is 3. + + Nclusters: int, optional + The number of data point clusters in the spectral map. Default is 4. + + ptsPerCluster: int, optional + The number of data points per cluster. Default is 50. + + Nx: int, optional + The pixel numbers of x-axis. Default is 64. + + Ny: int, optional + The pixel numbers of y-axis. Default is 64. + + cov: int or float matrix, optional + The covariance matrix of 2D gaussian distribution used to generate random data point clusters. + Default is [[20, 10], [10, 25]]. + + sigma: int or float, optional + The standard deviation of 1D gaussian distribution used to generate spectral weights. Default is 10. + + startWavenumber: int or float, optional + The beginning wavenumber of spectral range. Default is 400. + + endWavenumber: int or float, optional + The ending wavenumber of spectral range. Default is 4000. + + Nwavenumber: int, optional + The number of wavenumber values in the spectrum. Default is 1600. + + random_state: int, optional + The seed of the pseudo random number generator to use when generating random numbers. Default to 17. + + Attributes: + -------- + data : float matrix + The final spectral data matrix. Number of rows equals number of non-zero data points; + Number of columns equals number of wavenumbers. + + densityMatCondense : float matrix + The distribution weight matrix of each basis spectrum (only stores non-zero data points, see nonZeroInd). + Number of rows equals number of non-zero data points; + Number of columns equals number of basis spectrum. + + spectraMat : float matrix + The basis spectra matrix. Number of rows equals number of basis spectrum; + Number of columns equals number of wavenumbers. + + wavenumber : float array + The wavenumber values of the spectrum (x-axis). + + nonZeroInd : int array + The linear index of all non-zero data points in the flatterned 2D map. + + Examples: + --------- + >>> s = spectra_map_simulator(random_state=3) + >>> s.spectra_map_gen() + >>> print(s.data.shape) + (446, 1600) + """ + + def __init__(self, NbaseSpectra=3, Nclusters=4, ptsPerCluster=50, Nx=64, cov=[[20, 10], [10, 25]], sigma=10, + startWavenumber=400, endWavenumber=4000, Nwavenumber=1600, random_state=17): + + self.NbaseSpectra = NbaseSpectra + self.Nclusters = Nclusters + self.ptsPerCluster = ptsPerCluster + self.Nx = Nx + self.Ny = Nx + self.cov = cov + self.sigma = sigma + self.startWavenumber = startWavenumber + self.endWavenumber = endWavenumber + self.Nwavenumber = Nwavenumber + self.random_state = random_state + + def cluster_placement(self, random_state=17): + """Generate simulated data point locations and weights in a 2D map + + Returns: + -------- + points : float matrix + The data point location matrix. The first two columns are x-y coordiantes of data points. + The third column are the weight coefficients of data points. + + densityVec : float array + The 1D array generated from flattened 2D weight coefficients matrix (densityMat) + """ + + np.random.seed(seed=random_state) + centroids = np.random.randint(int(self.Nx*0.8), size=(self.Nclusters,2)) + self.points = np.zeros((self.ptsPerCluster*self.Nclusters, 3)) + densityMat = np.zeros((self.Ny, self.Nx)) + + for i in range(self.Nclusters): + self.points[i*self.ptsPerCluster:(i+1)*self.ptsPerCluster, :2] = np.random.multivariate_normal(centroids[i,:], self.cov, self.ptsPerCluster) + self.points[i*self.ptsPerCluster:(i+1)*self.ptsPerCluster, 2] = gaussian(np.sqrt(np.power(self.points[i*self.ptsPerCluster:(i+1)*self.ptsPerCluster, :2] + -centroids[i,:],2).sum(axis = 1)),sigma=self.sigma) + + for i in self.points: + idx = tuple(np.clip(i[1::-1].astype(int), 0, self.Nx-1)) + densityMat[idx] += i[2] + + densityVec = densityMat.flatten() + + return densityVec + + def spectrum_gen(self, Npeaks=3, firstPeakPosition=1200, peakWidth=[30,60,150], random_state=17): + """Generate a simulated spectrum + + Parameters: + ----------- + Npeaks: int, optional + The number of peaks. Default is 3. + + firstPeakPosition: int, optional + The position of the first peak. Default is 1200. + + peakWidth: int or float list, optional + The list of peak widths. Default is [30,60,150]. + + Returns: + -------- + wavenumber : float array + The wavenumber values of the spectrum (x-axis). + + spectrum : float array + The transmission/reflection/absorption coefficients of the spectrum (y-axis) + """ + + if len(peakWidth) != Npeaks: + raise Exception("The number of peak width values doesn't match the number of peaks") + + np.random.seed(seed=random_state) + self.wavenumber = np.linspace(self.startWavenumber,self.endWavenumber,self.Nwavenumber) + peakPositions = np.linspace(firstPeakPosition,self.endWavenumber,Npeaks,endpoint=False) + spectrum = np.zeros(len(self.wavenumber)) + weights = np.random.rand(Npeaks) + + for i in range(Npeaks): + singlePeak = lorentzian(self.wavenumber,peakPositions[i],peakWidth[i]) + spectrum += singlePeak * weights[i] + spectrum /= weights.sum() + + return spectrum + + def spectra_map_gen(self, Npeaks = [3,4,5], positions = [800, 1000, 1200]): + """Generate a spectral map data cube + + Parameters: + ----------- + Npeaks: int list, optional + The list of number of peaks. Default is [3,4,5]. + + positions: int or float list, optional + The list of first peak positions. Default is [800, 1000, 1200]. + + Returns + ------- + self : object + """ + + if len(positions) != len(Npeaks): + raise Exception("The number of first peak positions doesn't match the number of basis spectrum") + + np.random.seed(seed=self.random_state) + random_seeds = np.random.randint(50, size=len(Npeaks)) + + densityVecs = np.zeros((self.Nx*self.Ny,self.NbaseSpectra)) # component coefficients + self.spectraMat = np.zeros((self.NbaseSpectra, self.Nwavenumber)) # components spectra matrix + + for i in range(self.NbaseSpectra): + + densityVecs[:,i] = self.cluster_placement(random_state = random_seeds[i]) + + np.random.seed(seed=random_seeds[i]) + peakWidth = np.random.randint(4,20, Npeaks[i])*10 + self.spectraMat[i,:] = self.spectrum_gen(Npeaks=Npeaks[i], firstPeakPosition = positions[i], + peakWidth=peakWidth, random_state = random_seeds[i]) + + self.nonZeroInd = ~np.all(densityVecs==0, axis=1) + self.densityVecsCondense = densityVecs[self.nonZeroInd] # remove all zero rows + self.data = self.densityVecsCondense.dot(self.spectraMat) # matrix mulplication + + mask = np.zeros((self.Ny, self.Nx), dtype='bool').flatten() + mask[self.nonZeroInd] = True + self.mask = mask.reshape(self.Ny, self.Nx) + + def save(self, sample_id='simulated_dataset'): + """Save the simulated dataset as an hdf5 file. + + Arguments: + + ---------- + sample_id : A string that identifies the sample, spaces will be substituted for underscores. + + """ + si = sample_info(sample_id=sample_id, sample_meta_data=f"This dataset has {self.NbaseSpectra} components and map size is {self.Ny} * {self.Nx}") + si.show() + + self.spectra_map_gen() # generate simulated dataset + y, x = np.where(self.mask) + self.xy = np.c_[x,y] + + ir_data = ir_map(self.wavenumber, si, with_factorization=True) + ir_data.add_data(spectrum=self.data, xy=self.xy) + ir_data.to_image_cube() + ir_data.add_factorization(component=self.spectraMat, component_coef=self.densityVecsCondense, prefix='MCR') + ir_data.write_as_hdf5(f'{sample_id}.h5') + + def load(self, filename='simulated_dataset.h5'): + """load the simulated dataset from an hdf5 file. + + Arguments: + + ---------- + filename : The hdf5 filename where data will be read from. + """ + ir_data = ir_map(filename=filename) + ir_data.add_image_cube() + ir_data.add_factorization(prefix='MCR') + + return ir_data + + def plot_spectra_map(self): + """Plot the spectral distribution map and basis spectra + + """ + + self.spectra_map_gen() # generate simulated dataset + + plt.figure(figsize=(12, 6)) + for i in range(self.NbaseSpectra): + + n_row = 2 + + plt.subplot(n_row,3,i+1) + densityMat = np.zeros((self.Ny, self.Nx)).flatten() + densityMat[self.nonZeroInd] = self.densityVecsCondense[:,i] + densityMat = densityMat.reshape(self.Ny, self.Nx) + plt.imshow(np.flipud(densityMat)) + plt.title(f'Distribution of Component {i+1}') + plt.colorbar() + plt.clim([0, 2]) + + if n_row == 3: + plt.subplot(n_row,3,i+7) + plt.imshow(np.flipud(densityMat>0)) + + plt.subplot(n_row,3,i+4) + plt.subplots_adjust(hspace=0.3) + plt.plot(self.wavenumber, self.spectraMat[i,:]) + plt.title(f'Spectrum of Component {i+1}') + plt.xlim([4000,400]) + + return None + +if __name__ == "__main__": + + s = spectra_map_simulator(random_state=3) + s.plot_spectra_map() \ No newline at end of file diff --git a/lbl_ir/lbl_ir/tasks/NMF/__init__.py b/lbl_ir/lbl_ir/tasks/NMF/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/lbl_ir/lbl_ir/tasks/NMF/basic_NMF.py b/lbl_ir/lbl_ir/tasks/NMF/basic_NMF.py new file mode 100644 index 0000000..fa9d2cf --- /dev/null +++ b/lbl_ir/lbl_ir/tasks/NMF/basic_NMF.py @@ -0,0 +1,40 @@ +from sklearn.decomposition import NMF +import matplotlib.pyplot as plt +import pyqtgraph as pq + +import sys +import numpy as np +from lbl_ir.io_tools import read_map +from lbl_ir.tasks.preprocessing import data_prep + + +def simple_NMF(ir_map, wmask, smask,components=40): + """ + A tools to do NMF on a selected region in a map + + :param ir_map: Image data cube + :param wmask: wavelength mask + :param smask: spatial mask + :param components: number of components + :return: spatial maps of components and associated spectra + """ + ori_shape = ir_map.imageCube.shape + data = ir_map.imageCube[:,:,wmask] + data = data.reshape( ori_shape[0]*ori_shape[1], data.shape[2] ) + position_selection = np.where(smask.flatten()>0.5)[0] + data = data[ position_selection, : ] + + + NMF_obj = NMF(n_components=components) + W = NMF_obj.fit_transform( data ) + H = NMF_obj.components_ + + maps = [] + for nn in range(components): + this_map = np.zeros(ori_shape[0:2]) + this_map = this_map.flatten() + this_map[position_selection] = W[:,nn] + this_map = this_map.reshape( ori_shape[0:2]) + maps.append(this_map) + + return maps, H diff --git a/lbl_ir/lbl_ir/tasks/NMF/multi_set_analyses.py b/lbl_ir/lbl_ir/tasks/NMF/multi_set_analyses.py new file mode 100644 index 0000000..edbebcf --- /dev/null +++ b/lbl_ir/lbl_ir/tasks/NMF/multi_set_analyses.py @@ -0,0 +1,124 @@ +from sklearn.decomposition import NMF +import matplotlib.pyplot as plt +#import pyqtgraph as pq + + +import sys +import numpy as np +from lbl_ir.io_tools import read_map +from lbl_ir.tasks.preprocessing import data_prep, transform +#from umap import UMAP + + +class aggregate_data(object): + def __init__(self, names, data, wmask, components=40): + self.names = names + self.Nset = len(names) + self.wmask = wmask + self.master_wmask = self.get_wavenumber_intersection() + self.wavenumbers = data[0].wavenumbers[self.master_wmask] + self.data, self.dims, self.labels = self.get_data(data) + del data + self.components = components + print(self.dims) + + + def get_wavenumber_intersection(self): + result = self.wmask[0] + for ii in range(1,self.Nset): + result = np.intersect1d(result,self.wmask[ii]) + return result + + def get_data(self,data): + result = [] + dims = [] + labels = [] + for kk,cube in enumerate(data): + sel_cube = cube.imageCube[:,:, self.master_wmask ] + dim = sel_cube.shape + sel_cube = sel_cube.reshape( dim[0]*dim[1],dim[2] ) + result.append( sel_cube ) + labs = np.zeros( dim[0]*dim[1] ) + kk + labels.append(labs) + dims.append( dim[:2] ) + result = np.vstack( result ) + labels = np.concatenate(labels).flatten() + return transform.to_absorbance(result,False)[0], dims, labels + + def splitter(self,X,to_map=False): + sets = [] + for kk in range( self.Nset ): + sel = self.labels == kk + sel = np.where(sel)[0] + tmp = X[sel,:] + if to_map: + tmp = tmp.reshape( self.dims[kk] ) + sets.append( tmp ) + return sets + + + + +def multi_set_umap(agg_data, fraction=0.25): + umap_object = UMAP(n_components=2, n_neighbors=5) + Nobs = agg_data.data.shape[0] + these_ones = np.arange(Nobs) + np.random.shuffle(these_ones) + these_ones = these_ones[:int(Nobs*fraction)] + these_ones = np.sort(these_ones) + + low_dim = umap_object.fit_transform( agg_data.data[these_ones,:] ) + low_dim_all = umap_object.transform( agg_data.data) + low_dim_all_split = agg_data.splitter(low_dim_all,False) + + for set in low_dim_all_split: + plt.plot(set[:,0],set[:,1],'.' , markersize=1.5 ) + plt.savefig('UMAP.png') + +def multi_set_NMF(agg_data,components=40): + NMF_obj = NMF(n_components=components) + W = NMF_obj.fit_transform( agg_data.data ) + H = NMF_obj.components_ + spectra = [] + maps = [] + for ii in range(components): + Ws = W[:,ii].reshape(-1,1) + Ws_split = agg_data.splitter( Ws,True ) + spectra.append(H[ii,:]) + maps.append( Ws_split ) + return maps, spectra, agg_data.wavenumbers + + +def single_set_NMF(data, wav_mask, components): + A = transform.fix_data( data.data[:,wav_mask] ) + NMF_obj = NMF(n_components=components) + W = NMF_obj.fit_transform( A ) + H = NMF_obj + return data.wavenumbers[wav_mask],W,H + + + +if __name__ == "__main__": + data_files = [] + wav_masks = [] + names = ['../../../ir_data/test_data.map'] + if len(names)>1: + for fname in names: + data,fmt = read_map.read_all_formats( fname ) + data_files.append( data ) + ds = data_prep.data_prepper(data) + wav_masks.append( ds.decent_bands ) + + ad = aggregate_data(names, data_files, wav_masks) + #multi_set_umap( ad ) + multi_set_NMF( ad ) + else: + data,fmt = read_map.read_all_formats( names[0] ) + print(fmt) + wmask = data_prep.data_prepper(data).decent_bands + wavs,W,H = single_set_NMF(data, wmask, 4 ) + for i in range(4): + plt.plot(wavs, H.components_[i], label='NMF'+str(i)) + plt.xlim([4000, 500]) + plt.legend() + diff --git a/lbl_ir/lbl_ir/tasks/Worms/Shai_Hulut.py b/lbl_ir/lbl_ir/tasks/Worms/Shai_Hulut.py new file mode 100644 index 0000000..55d7b1a --- /dev/null +++ b/lbl_ir/lbl_ir/tasks/Worms/Shai_Hulut.py @@ -0,0 +1,227 @@ +""" +Here are some specific worm processing tools. + +""" + +import numpy as np +import matplotlib.pyplot as plt + +# reading an omnic map +from lbl_ir.io_tools.Omnic_PyMca5.OmnicMap import OmnicMap + +# doing SVD +from lbl_ir.math_tools.batched_SVD import batched_SVD + +# for segmentating of the image +from scipy.ndimage.morphology import binary_closing, binary_dilation, binary_fill_holes +from scipy.ndimage import gaussian_filter + +from skimage.morphology import skeletonize +from skimage import draw +import skimage +import pandas +from skan import csr +from scipy.ndimage.morphology import distance_transform_edt + + + +class little_maker( object ): + def __init__(self, filename, k_singular=10, z_threshold =14): + self.filename = filename + self.k_singular = k_singular + self.z_threshold = z_threshold + + self.omnic_object = OmnicMap(filename) + if self.omnic_object.info['OmnicInfo'] is not None: + start_wav = self.omnic_object.info['OmnicInfo']['First X value'] + stop_wav = self.omnic_object.info['OmnicInfo']['Last X value'] + n_wav = self.omnic_object.info['OmnicInfo']['Number of points'] + else: + start_wav = 699 + stop_wav = 3999 + n_wav = self.omnic_object.data.shape[-1] + self.waves = np.linspace(start_wav, stop_wav, n_wav ) # I'm not sure where to find this in the data object! + + self.data_shape = self.omnic_object.data.shape + self.flattened_data = self.omnic_object.data.reshape(-1, self.data_shape[2] ) + + self.full_frame_svd() + self.z_scores = self.get_background_map() + + def full_frame_svd(self,N_max=1000): + svd_obj = batched_SVD( self.flattened_data, + N_max=N_max, + k_singular=self.k_singular, + randomize=True) + self.U,self.S,self.VT = svd_obj.go_svd() + + + def skeletonize(self, threshold): + fin_mask = (self.z_scores > threshold).astype(int) + skel = skeletonize( fin_mask ).astype(int) + sk_obj = csr.Skeleton(skel, spacing=1, keep_images=True ) + + path_lengths = sk_obj.path_lengths() + this_one = np.argmax(path_lengths) + path_coordinates = sk_obj.path_coordinates(this_one).astype(int) + + pruned_skel = np.zeros( fin_mask.shape ) + for pair in path_coordinates: + pruned_skel[ pair[0], pair[1] ] = 1 + return pruned_skel, path_coordinates[::-1,:] + + def blobber(self): + None + + def extract_extrema(self, N_pixels=20, radius=20, threshold=14): + skel, coords = self.skeletonize(threshold) + side_0_coords = coords[ 0:N_pixels ] + dxy = side_0_coords[0,:]- side_0_coords[-1,:] + angle_0 = np.arctan2(dxy[1], dxy[0]) + print(angle_0*180.0/3.14) + print( side_0_coords ) + nimage = skimage.transform.rotate( self.z_scores , -angle_0*180.0/np.pi, clip=False, resize=True ) + plt.imshow(nimage); plt.show() + + side_1_coords = coords[ -(N_pixels): ] + dxy = side_1_coords[0,:]- side_1_coords[-1,:] + angle_1 = np.arctan2(dxy[1], dxy[0]) + print(angle_1*180.0/3.14) + print( side_1_coords ) + nimage = skimage.transform.rotate( self.z_scores , -angle_1*180.0/np.pi, clip=False, resize=True ) + plt.imshow(nimage); plt.show() + + + + + + + + def get_background_map(self, percentile=75.0, safety1=10, safety2=10, window_MM=6): + """ + We first use a local roughness check to determine edges. We do this on the mean images, + but could easly use the full hypercube if need be. + """ + + + # get a mean image + mean_img = self.U[:,0].reshape( self.data_shape[0:2] ) + mean_img = mean_img - np.mean(mean_img) + + # lets padd the image to avoid FFT periodicity issues + X,Y = mean_img.shape + # we need to fill this with some decent values + fill_sigma = np.std(mean_img) + fill_mean = np.median( mean_img ) + new_img = np.random.normal( fill_mean, fill_sigma, (X+2*safety1,Y+2*safety1)) + # blur it a bit + new_img = gaussian_filter( new_img, sigma=10) + + # paste the real image in here + new_img[ safety1:safety1+X,safety1:safety1+Y ] = mean_img + mean_img = new_img+0 + + + # define a windowing function of MM by MM pixels + MM = window_MM + kernel = np.zeros( mean_img.shape ) + kernel[0:MM,0:MM] = 1.0 + + # compute a local mean + FT_img = np.fft.fft2(mean_img) + FT_kernel = np.fft.fft2(kernel) + local_mean= np.fft.ifft2( FT_img* FT_kernel.conjugate() ).real / ( MM*MM) + + # compute a local variance + mean_img_sq = mean_img * mean_img + FT_img_sq = np.fft.fft2(mean_img_sq) + local_var = np.fft.ifft2( FT_img_sq* FT_kernel.conjugate() ).real / (MM*MM) + local_var = local_var - local_mean*local_mean + local_sigma = np.sqrt( local_var ) + + ############################################################## + ## Now we build a rough mask + ############################################################## + threshold = np.percentile( local_sigma.flatten(), percentile ) + sel = local_sigma > threshold + mask = np.zeros( mean_img.shape ) + mask[sel] = 1.0 + + # The morphological operators need some room on the sides + safety2 = 0 + V,W = mask.shape + nV = V+safety2*2 + nW = W+safety2*2 + new_mask = np.zeros( (nV,nW) ) + # place the mask + new_mask[safety2:safety2+V, safety2:safety2+W ] = mask + + BM=5 + structure = np.ones((BM,BM)) + closed = binary_closing(new_mask,structure) + + done=False + BM=5 + while not done: + BM = BM + 1 + structure = np.ones((BM,BM)) + new_closed = binary_closing(closed,structure).astype(int) + delta = np.sum( np.abs(new_closed.astype(int) - closed.astype(int)) ) + if delta == 0: + done = True + if BM >20: + done = True + closed = new_closed+0 + + new_mask = closed+0 + + # for some reason, the whole mask is shifted along one axis. This has likely to do with some + # origin definition of the structuring elements, but lets solve this using an FFT based shift calculation + # The origin option in the scipy morphology toolbox kill my kernel + + mask = new_mask[safety2:safety2+V, safety2:safety2+W] + ft_mask = np.fft.fft2(mask) + ft_mean = np.fft.fft2(np.abs(mean_img) ) + TF = np.fft.ifft2( ft_mean*ft_mask.conjugate() ).real + dX,dY = np.meshgrid( np.arange(mask.shape[1]), np.arange(mask.shape[0]) ) + here = np.argmax( TF ) + dX = dX.flatten()[here] + dY = dY.flatten()[here] + + # no wild shifts please + if dX > 10: + dX = 0 + if dY > 10: + dY = 0 + + mask = np.roll( mask, dX, axis=0) + mask = np.roll( mask, dY, axis=1) + # lift out the section of the mask that we are interested in + mask = mask[safety1:safety1+X, safety1:safety1+Y] + + + ###################################################### + # Here we build a more fine tuned mask / z_score map + ###################################################### + + + # now we use this mask to define the background + bg_sel = mask.flatten() < 0.5 + background = self.U[bg_sel,:] + mean_bg = np.mean( background, axis=0) + var_covar = np.cov( (background-mean_bg).transpose() ) + inv_vcv = np.linalg.pinv( var_covar ) + t = self.U - mean_bg + z_scores = [] + for tt in t: + z = tt.reshape(1,-1) + z_scores.append( np.sqrt( z.dot(inv_vcv).dot(z.transpose()) ) ) + z_scores = np.array(z_scores).reshape( self.data_shape[0:2] ) + return z_scores + + + + + + + diff --git a/lbl_ir/lbl_ir/tasks/Worms/__init__.py b/lbl_ir/lbl_ir/tasks/Worms/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/lbl_ir/lbl_ir/tasks/__init__.py b/lbl_ir/lbl_ir/tasks/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/lbl_ir/lbl_ir/tasks/baseline/__init__.py b/lbl_ir/lbl_ir/tasks/baseline/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/lbl_ir/lbl_ir/tasks/baseline/rubberband.py b/lbl_ir/lbl_ir/tasks/baseline/rubberband.py new file mode 100644 index 0000000..5d39c33 --- /dev/null +++ b/lbl_ir/lbl_ir/tasks/baseline/rubberband.py @@ -0,0 +1,29 @@ +import numpy as np +from scipy.spatial import ConvexHull + +def rubberband(wavenumbers, spectrum): + """ + A rubberband baseline is essentially the lower half of tyhe convex hull + of the data. See this discussion for more details and code. + + https://dsp.stackexchange.com/questions/2725/how-to-perform-a-rubberband-correction-on-spectroscopic-data + """ + + # First find the convex hull using scipy.spatial tools + tmp = np.vstack([ wavenumbers.flatten(),spectrum.flatten()]).transpose() + v = ConvexHull(tmp).vertices + + # rool it untill the first point is the first point measured + v = np.roll(v, -v.argmin()) + # We don't care aboput the top half of the convex hull + v = v[:v.argmax()] + + # Use interpolation to get the baseline across the whole + # spectrum + return np.interp(wavenumbers, wavenumbers[v], spectrum[v]) + + + + + + diff --git a/lbl_ir/lbl_ir/tasks/preprocessing/EMSC.py b/lbl_ir/lbl_ir/tasks/preprocessing/EMSC.py new file mode 100644 index 0000000..bc518f8 --- /dev/null +++ b/lbl_ir/lbl_ir/tasks/preprocessing/EMSC.py @@ -0,0 +1,510 @@ +import numpy as np +import scipy.optimize +import sklearn.decomposition as skl_decomposition +from scipy.signal import hilbert + + +def konevskikh_parameters(a, n0, f): + """ + Compute parameters for Konevskikh algorithm + :param a: cell radius + :param n0: refractive index + :param f: scaling factor + :return: parameters alpha0 and gamma + """ + alpha0 = 4.0 * np.pi * a * (n0 - 1.0) + gamma = np.divide(f, n0 - 1.0) + return alpha0, gamma + + +def GramSchmidt(V): + """ + Perform Gram-Schmidt normalization for the matrix V + :param V: matrix + :return: nGram-Schmidt normalized matrix + """ + V = np.array(V) + U = np.zeros(np.shape(V)) + + for k in range(len(V)): + sum1 = 0 + for j in range(k): + sum1 += np.dot(V[k], U[j]) / np.dot(U[j], U[j]) * U[j] + U[k] = V[k] - sum1 + return U + + +def check_orthogonality(U): + """ + Check orthogonality of a matrix + :param U: matrix + """ + for i in range(len(U)): + for j in range(i, len(U)): + if i != j: + print(np.dot(U[i], U[j])) + + +def find_nearest_number_index(array, value): + """ + Find the nearest number in an array and return its index + :param array: + :param value: value to be found inside the array + :return: position of the number closest to value in array + """ + array = np.array(array) # Convert to numpy array + if np.shape(np.array(value)) is (): # If only one value wants to be found: + index = (np.abs(array - value)).argmin() # Get the index of item closest to the value + else: # If value is a list: + value = np.array(value) + index = np.zeros(np.shape(value)) + k = 0 + # Find the indexes for all values in value + for val in value: + index[k] = (np.abs(array - val)).argmin() + k += 1 + index = index.astype(int) # Convert the indexes to integers + return index + + +def Q_ext_kohler(wn, alpha): + """ + Compute the scattering extinction values for a given alpha and a range of wavenumbers + :param wn: array of wavenumbers + :param alpha: scalar alpha + :return: array of scattering extinctions calculated for alpha in the given wavenumbers + """ + rho = alpha * wn + Q = 2.0 - (4.0 / rho) * np.sin(rho) + (2.0 / rho) ** 2.0 * (1.0 - np.cos(rho)) + return Q + + +def apparent_spectrum_fit_function(wn, Z_ref, p, b, c, g): + """ + Function used to fit the apparent spectrum + :param wn: wavenumbers + :param Z_ref: reference spectrum + :param p: principal components of the extinction matrix + :param b: Reference's linear factor + :param c: Offset + :param g: Extinction matrix's PCA scores (to be fitted) + :return: fitting of the apparent specrum + """ + A = b * Z_ref + c + np.dot(g, p) # Extended multiplicative scattering correction formula + return A + + +def reference_spectrum_fit_function(wn, p, c, g): + """ + Function used to fit a reference spectrum (without using another spectrum as reference). + :param wn: wavenumbers + :param p: principal components of the extinction matrix + :param c: offset + :param g: PCA scores (to be fitted) + :return: fitting of the reference spectrum + """ + A = c + np.dot(g, p) + return A + + +def apparent_spectrum_fit_function_Bassan(wn, Z_ref, p, c, m, h, g): + """ + Function used to fit the apparent spectrum in Bassan's algorithm + :param wn: wavenumbers + :param Z_ref: reference spectrum + :param p: principal componetns of the extinction matrix + :param c: offset + :param m: linear baseline + :param h: reference's linear factor + :param g: PCA scores to be fitted + :return: fitting of the apparent spectrum + """ + A = c + m * wn + h * Z_ref + np.dot(g, p) + return A + + +def correct_reference(m, wn, a, d, w_regions): + """ + Correct reference spectrum as in Kohler's method + :param m: reference spectrum + :param wn: wavenumbers + :param a: Average refractive index range + :param d: Cell diameter range + :param w_regions: Weighted regions + :return: corrected reference spectrum + """ + n_components = 6 # Set the number of principal components + + # Copy the input variables + m = np.copy(m) + wn = np.copy(wn) + + # Compute the alpha range: + alpha = 4.0 * np.pi * 0.5 * np.linspace(np.min(d) * (np.min(a) - 1.0), np.max(d) * (np.max(a) - 1.0), 150) + + p0 = np.ones(1 + n_components) # Initial guess for the fitting + + # Compute extinction matrix + Q_ext = np.zeros((np.size(alpha), np.size(wn))) + for i in range(np.size(alpha)): + Q_ext[i][:] = Q_ext_kohler(wn, alpha=alpha[i]) + + # Perform PCA to Q_ext + pca = skl_decomposition.IncrementalPCA(n_components=n_components) + pca.fit(Q_ext) + p_i = pca.components_ # Get the principal components of the extinction matrix + + # Get the weighted regions of the wavenumbers, the reference spectrum and the principal components + w_indexes = [] + for pair in w_regions: + min_pair = min(pair) + max_pair = max(pair) + ii1 = find_nearest_number_index(wn, min_pair) + ii2 = find_nearest_number_index(wn, max_pair) + w_indexes.extend(np.arange(ii1, ii2)) + wn_w = np.copy(wn[w_indexes]) + m_w = np.copy(m[w_indexes]) + p_i_w = np.copy(p_i[:, w_indexes]) + + def min_fun(x): + """ + Function to be minimized for the fitting + :param x: offset and PCA scores + :return: difference between the spectrum and its fitting + """ + cc, g = x[0], x[1:] + # Return the squared norm of the difference between the reference spectrum and its fitting: + return np.linalg.norm(m_w - reference_spectrum_fit_function(wn_w, p_i_w, cc, g)) ** 2.0 + + # Perform the minimization using Powell method + res = scipy.optimize.minimize(min_fun, p0, bounds=None, method='Powell') + + c, g_i = res.x[0], res.x[1:] # Obtain the fitted parameters + + # Apply the correction: + m_corr = np.zeros(np.shape(m)) + for i in range(len(wn)): + sum1 = 0 + for j in range(len(g_i)): + sum1 += g_i[j] * p_i[j][i] + m_corr[i] = (m[i] - c - sum1) + + return m_corr # Return the corrected spectrum + + +def Kohler(wavenumbers, App, m0, n_components=8): + """ + Correct scattered spectra using Kohler's algorithm + :param wavenumbers: array of wavenumbers + :param App: apparent spectrum + :param m0: reference spectrum + :param n_components: number of principal components to be calculated + :return: corrected data + """ + # Make copies of all input data: + wn = np.copy(wavenumbers) + A_app = np.copy(App) + m_0 = np.copy(m0) + ii = np.argsort(wn) # Sort the wavenumbers from smallest to largest + # Sort all the input variables accordingly + wn = wn[ii] + A_app = A_app[ii] + m_0 = m_0[ii] + + # Initialize the alpha parameter: + alpha = np.linspace(3.14, 49.95, 150) * 1.0e-4 # alpha = 2 * pi * d * (n - 1) * wavenumber + p0 = np.ones(2 + n_components) # Initialize the initial guess for the fitting + + # # Initialize the extinction matrix: + Q_ext = np.zeros((np.size(alpha), np.size(wn))) + for i in range(np.size(alpha)): + Q_ext[i][:] = Q_ext_kohler(wn, alpha=alpha[i]) + + # Perform PCA of Q_ext: + pca = skl_decomposition.IncrementalPCA(n_components=n_components) + pca.fit(Q_ext) + p_i = pca.components_ # Extract the principal components + + # print(np.sum(pca.explained_variance_ratio_)*100) # Print th explained variance ratio in percentage + + def min_fun(x): + """ + Function to be minimized by the fitting + :param x: array containing the reference linear factor, the offset, and the PCA scores + :return: function to be minimized + """ + bb, cc, g = x[0], x[1], x[2:] + # Return the squared norm of the difference between the apparent spectrum and the fit + return np.linalg.norm(A_app - apparent_spectrum_fit_function(wn, m_0, p_i, bb, cc, g)) ** 2.0 + + # Minimize the function using Powell method + res = scipy.optimize.minimize(min_fun, p0, bounds=None, method='Powell') + # print(res) # Print the minimization result + # assert(res.success) # Raise AssertionError if res.success == False + + b, c, g_i = res.x[0], res.x[1], res.x[2:] # Obtain the fitted parameters + + # Apply the correction to the apparent spectrum + Z_corr = np.zeros(np.shape(m_0)) + for i in range(len(wavenumbers)): + sum1 = 0 + for j in range(len(g_i)): + sum1 += g_i[j] * p_i[j][i] + Z_corr[i] = (A_app[i] - c - sum1) / b + + return Z_corr[::-1] # Return the correction in reverse order for compatibility + +def Kohler_zero(wavenumbers, App, w_regions, n_components=8): + """ + Correct scattered spectra using Kohler's algorithm + :param wavenumbers: array of wavenumbers + :param App: apparent spectrum + :param m0: reference spectrum + :param n_components: number of principal components to be calculated + :return: corrected data + """ + # Make copies of all input data: + wn = np.copy(wavenumbers) + A_app = np.copy(App) + m_0 = np.zeros(len(wn)) + ii = np.argsort(wn) # Sort the wavenumbers from smallest to largest + # Sort all the input variables accordingly + wn = wn[ii] + A_app = A_app[ii] + m_0 = m_0[ii] + + # Initialize the alpha parameter: + alpha = np.linspace(1.25, 49.95, 150) * 1.0e-4 # alpha = 2 * pi * d * (n - 1) * wavenumber + p0 = np.ones(2 + n_components) # Initialize the initial guess for the fitting + + # # Initialize the extinction matrix: + Q_ext = np.zeros((np.size(alpha), np.size(wn))) + for i in range(np.size(alpha)): + Q_ext[i][:] = Q_ext_kohler(wn, alpha=alpha[i]) + + # Perform PCA of Q_ext: + pca = skl_decomposition.IncrementalPCA(n_components=n_components) + pca.fit(Q_ext) + p_i = pca.components_ # Extract the principal components + + # print(np.sum(pca.explained_variance_ratio_)*100) # Print th explained variance ratio in percentage + w_indexes = [] + for pair in w_regions: + min_pair = min(pair) + max_pair = max(pair) + ii1 = find_nearest_number_index(wn, min_pair) + ii2 = find_nearest_number_index(wn, max_pair) + w_indexes.extend(np.arange(ii1, ii2)) + wn_w = np.copy(wn[w_indexes]) + A_app_w = np.copy(A_app[w_indexes]) + m_w = np.copy(m_0[w_indexes]) + p_i_w = np.copy(p_i[:, w_indexes]) + + def min_fun(x): + """ + Function to be minimized by the fitting + :param x: array containing the reference linear factor, the offset, and the PCA scores + :return: function to be minimized + """ + bb, cc, g = x[0], x[1], x[2:] + # Return the squared norm of the difference between the apparent spectrum and the fit + return np.linalg.norm(A_app_w - apparent_spectrum_fit_function(wn_w, m_w, p_i_w, bb, cc, g)) ** 2.0 + + # Minimize the function using Powell method + res = scipy.optimize.minimize(min_fun, p0, bounds=None, method='Powell') + # print(res) # Print the minimization result + # assert(res.success) # Raise AssertionError if res.success == False + + b, c, g_i = res.x[0], res.x[1], res.x[2:] # Obtain the fitted parameters + + # Apply the correction to the apparent spectrum + Z_corr = (A_app - c - np.dot(g_i, p_i)) # Apply the correction + base = np.dot(g_i, p_i) + + return Z_corr, base + +def Bassan(wavenumbers, App, m0, n_components=8, iterations=1, w_regions=None): + """ + Correct scattered spectra using Bassan's algorithm. + :param wavenumbers: array of wavenumbers + :param App: apparent spectrum + :param m0: reference spectrum + :param n_components: number of principal components to be calculated for the extinction matrix + :param iterations: number of iterations of the algorithm + :param w_regions: the regions to be taken into account for the fitting + :return: corrected apparent spectrum + """ + # Copy the input data + wn = np.copy(wavenumbers) + A_app = np.copy(App) + m_0 = np.copy(m0) + ii = np.argsort(wn) # Sort the wavenumbers + # Apply the sorting to the input variables + wn = wn[ii] + A_app = A_app[ii] + m_0 = m_0[ii] + + # Define the weighted regions: + if w_regions is not None: + m_0 = correct_reference(np.copy(m_0), wn, a, d, w_regions) # Correct the reference spectrum as in Kohler method + w_indexes = [] + # Get the indexes of the regions to be taken into account + for pair in w_regions: + min_pair = min(pair) + max_pair = max(pair) + ii1 = find_nearest_number_index(wn, min_pair) + ii2 = find_nearest_number_index(wn, max_pair) + w_indexes.extend(np.arange(ii1, ii2)) + # Take the weighted regions of wavenumbers, apparent and reference spectrum + wn_w = np.copy(wn[w_indexes]) + A_app_w = np.copy(A_app[w_indexes]) + m_0_w = np.copy(m_0[w_indexes]) + + n_loadings = 10 # Number of values to be computed for each parameter (a, b, d) + a = np.linspace(1.1, 1.5, n_loadings) # Average refractive index + d = np.linspace(2.0, 8.0, n_loadings) * 1.0e-4 # Cell diameter + Q = np.zeros((n_loadings ** 3, len(wn))) # Initialize the extinction matrix + m_n = np.copy(m_0) # Initialize the reference spectrum, that will be updated after each iteration + for iteration in range(iterations): + # Compute the scaled real part of the refractive index by Kramers-Kronig transform: + nkk = -1.0 * np.imag(hilbert(m_n)) + # Build the extinction matrix + n_row = 0 + for i in range(n_loadings): + b = np.linspace(0.0, a[i] - 1.0, 10) # Range of amplification factors of nkk + for j in range(n_loadings): + for k in range(n_loadings): + n = a[i] + b[j] * nkk # Compute the real refractive index + alpha = 2.0 * np.pi * d[k] * (n - 1.0) + rho = alpha * wn + # Compute the extinction coefficients for each combination of a, b and d: + Q[n_row] = 2.0 - np.divide(4.0, rho) * np.sin(rho) + \ + np.divide(4.0, rho ** 2.0) * (1.0 - np.cos(rho)) + n_row += 1 + + # Orthogonalization of th extinction matrix with respect to the reference spectrum: + for i in range(n_loadings ** 3): + Q[i] -= np.dot(Q[i], m_0) / np.linalg.norm(m_0) ** 2.0 * m_0 + + # Perform PCA of the extinction matrix + pca = skl_decomposition.IncrementalPCA(n_components=n_components) + pca.fit(Q) + p_i = pca.components_ # Get the principal components + + if w_regions is None: # If all regions have to be taken into account: + def min_fun(x): + """ + Function to be minimized for the fitting + :param x: fitting parameters (offset, baseline, reference's linear factor, PCA scores) + :return: squared norm of the difference between the apparent spectrum and its fitting + """ + cc, mm, hh, g = x[0], x[1], x[2], x[3:] + return np.linalg.norm(A_app - apparent_spectrum_fit_function_Bassan(wn, m_0, p_i, cc, mm, hh, g)) ** 2.0 + else: # If only the specified regions have to be taken into account: + # Take the indexes of the specified regions + w_indexes = [] + for pair in w_regions: + min_pair = min(pair) + max_pair = max(pair) + ii1 = find_nearest_number_index(wn, min_pair) + ii2 = find_nearest_number_index(wn, max_pair) + w_indexes.extend(np.arange(ii1, ii2)) + p_i_w = np.copy(p_i[:, w_indexes]) # Get the principal components of the extinction matrix at the + + # specified regions + + def min_fun(x): + """ + Function to be minimized for the fitting + :param x: fitting parameters (offset, baseline, reference's linear factor, PCA scores) + :return: squared norm of the difference between the apparent spectrum and its fitting + """ + cc, mm, hh, g = x[0], x[1], x[2], x[3:] + return np.linalg.norm(A_app_w - + apparent_spectrum_fit_function_Bassan(wn_w, m_0_w, p_i_w, cc, mm, hh, g)) ** 2.0 + + p0 = np.append([1.0, 0.0005, 0.9], np.ones(n_components)) # Initial guess for the fitting + res = scipy.optimize.minimize(min_fun, p0, method='Powell') # Perform the fitting + + # print(res) # Print the result of the minimization + # assert(res.success) # Raise AssertionError if res.success == False + + c, m, h, g_i = res.x[0], res.x[1], res.x[2], res.x[3:] # Take the fitted parameters + + Z_corr = (A_app - c - m * wn - np.dot(g_i, p_i)) / h # Apply the correction + + m_n = np.copy(Z_corr) # Take the corrected spectrum as the reference for the next iteration + + return np.copy(Z_corr[::-1]) # Return the corrected spectrum in inverted order for compatibility + + +def Konevskikh(wavenumbers, App, m0, n_components=8, iterations=1): + """ + Correct scattered spectra using Konevskikh algorithm + :param wavenumbers: array of wavenumbers + :param App: apparent spectrum + :param m0: reference spectrum + :param n_components: number of components + :param iterations: number of iterations + :return: corrected spectrum + """ + # Copy the input variables + wn = np.copy(wavenumbers) + A_app = np.copy(App) + m_0 = np.copy(m0) + ii = np.argsort(wn) # Sort the wavenumbers + wn = wn[ii] + A_app = A_app[ii] + m_0 = m_0[ii] + + # Initialize parameters range: + alpha_0, gamma = np.array([np.logspace(np.log10(0.1), np.log10(2.2), num=10) * 4.0e-4 * np.pi, + np.logspace(np.log10(0.05e4), np.log10(0.05e5), num=10) * 1.0e-2]) + p0 = np.ones(2 + n_components) + Q_ext = np.zeros((len(alpha_0) * len(gamma), len(wn))) # Initialize extinction matrix + + m_n = np.copy(m_0) # Copy the reference spectrum + for n_iteration in range(iterations): + ns_im = np.divide(m_n, wn) # Compute the imaginary part of the refractive index + # Compute the real part of the refractive index by Kramers-Kronig transform + ns_re = -1.0 * np.imag(hilbert(ns_im)) + + # Compute the extinction matrix + n_index = 0 + for i in range(len(alpha_0)): + for j in range(len(gamma)): + for k in range(len(A_app)): + rho = alpha_0[i] * (1.0 + gamma[j] * ns_re[k]) * wn[k] + beta = np.arctan(ns_im[k] / (1.0 / gamma[j] + ns_re[k])) + Q_ext[n_index][k] = 2.0 - 4.0 * np.exp(-1.0 * rho * np.tan(beta)) * (np.cos(beta) / rho) * \ + np.sin(rho - beta) - 4.0 * np.exp(-1.0 * rho * np.tan(beta)) * (np.cos(beta) / rho) ** 2.0 * \ + np.cos(rho - 2.0 * beta) + 4.0 * (np.cos(beta) / rho) ** 2.0 * np.cos(2.0 * beta) + # TODO: rewrite this in a simpler way + + n_index += 1 + + # Orthogonalize the extinction matrix with respect to the reference: + for i in range(n_index): + Q_ext[i][:] -= np.dot(Q_ext[i][:], m_0) / np.linalg.norm(m_0) ** 2.0 * m_0 + # Q_ext = GramSchmidt(np.copy(Q_ext)) # Apply Gram-Schmidt othogonalization to Q_ext (don't uncomment this) + + # Compute PCA of the extinction matrix + pca = skl_decomposition.IncrementalPCA(n_components=n_components) + pca.fit(Q_ext) + p_i = pca.components_ # Get the principal components + + def min_fun(x): + bb, cc, g = x[0], x[1], x[2:] + return np.linalg.norm(A_app - apparent_spectrum_fit_function(wn, m_0, p_i, bb, cc, g)) ** 2.0 + + res = scipy.optimize.minimize(min_fun, p0, method='Powell') + # print(res) # Print the minimization results + # assert(res.success) # Raise AssertionError if res.success == False + + b, c, g_i = res.x[0], res.x[1], res.x[2:] # Get the fitted parameters + + Z_corr = (A_app - c - np.dot(g_i, p_i)) / b # Apply the correction + + m_n = np.copy(Z_corr) # Update te reference with the correction + + return Z_corr[::-1] # Return the corrected spectrum diff --git a/lbl_ir/lbl_ir/tasks/preprocessing/__init__.py b/lbl_ir/lbl_ir/tasks/preprocessing/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/lbl_ir/lbl_ir/tasks/preprocessing/data_prep.py b/lbl_ir/lbl_ir/tasks/preprocessing/data_prep.py new file mode 100644 index 0000000..afc9901 --- /dev/null +++ b/lbl_ir/lbl_ir/tasks/preprocessing/data_prep.py @@ -0,0 +1,96 @@ +import matplotlib.pyplot as plt +#import pyqtgraph as pq + + +import sys +import numpy as np +from lbl_ir.io_tools import read_map + +from scipy.stats import iqr +import scipy.signal +import scipy.ndimage + +from numba import jit + + +@jit#(nopython=True) +def band_score_numba(map,subsample=1,band=10): + Nx,Ny,Nwav = map.shape + result = np.zeros((Nwav,Nwav)) + norma = np.zeros((Nwav,Nwav)) + for ii in range(Nwav): + for jj in range( ii+1 ,min(ii+band,Nwav)): + slab_0 = map[::subsample,::subsample,ii].flatten() + slab_1 = map[::subsample,::subsample,jj].flatten() + cc = np.corrcoef(slab_0, slab_1)[0][1] + norma[ii,jj]=1.0 + norma[jj,ii]=1.0 + result[ii,jj]=cc + result[jj,ii]=cc + result = np.sum( result, axis=1) + norma = np.sum( norma, axis=1) + return result/norma + + +class data_prepper(object): + """ + Prep the data for further data analyses. + We try to detect systematic issues, such as bad bands and synchortron noise spikes + + + """ + def __init__(self, data_map, band_limit=0.99, threshold=6.0,band=5, additional_selection = None): + self.data_map = data_map + self.waves = self.data_map.wavenumbers + self.threshold = threshold + self.additional_selection = additional_selection + self.band_limit = band_limit + self.band_scores = None + self.bad_bands = None + self.decent_bands = None + self.score_bands(band) + + def score_bands(self,band=10): + # we need to loop over every single frame and detect spikes + self.band_scores = band_score_numba( self.data_map.imageCube,subsample=1,band=band ) + self.bad_bands = np.where( self.band_scores < self.band_limit )[0] + self.decent_bands = np.where( self.band_scores >= self.band_limit)[0] + + def wavenumber_mask(self): + # we have to merge these auto score bands with the manual selection + selections = self.band_scores >= self.band_limit + tmp = np.zeros( len(self.waves) ) > 1.0 + if self.additional_selection is not None: + for this_sel in self.additional_selection: + pair = np.sort( this_sel ) + sel1 = (self.waves > pair[0]) & (self.waves or < + :param delta: selects band size around provided peak values. + :return: + """ + sel1 = (ir_map.wavenumbers > wave1-delta) & (ir_map.wavenumbers < wave1+delta) + sel2 = (ir_map.wavenumbers > wave2-delta) & (ir_map.wavenumbers < wave2+delta) + map1 = np.mean( ir_map.imageCube[:,:,sel1], axis=2 ) + map2 = np.mean( ir_map.imageCube[:,:,sel2], axis=2 ) + + mask1 = map1 > bg_threshold + mask2 = map2 > bg_threshold + bg_mask = mask1*mask2 + + + fin = map1 / map2 + if operator == '<': + fin = fin < peak_threshold + if operator == '>': + fin = fin > peak_threshold + fin = fin & bg_mask + mask = np.zeros( ir_map.imageCube.shape[0:2] ) + mask[fin] = 1.0 + return mask + + + + + + diff --git a/lbl_ir/lbl_ir/tasks/preprocessing/svd_data.py b/lbl_ir/lbl_ir/tasks/preprocessing/svd_data.py new file mode 100644 index 0000000..49ef021 --- /dev/null +++ b/lbl_ir/lbl_ir/tasks/preprocessing/svd_data.py @@ -0,0 +1,16 @@ +import numpy as np + +from lbl_ir.math_tools.batched_SVD import batched_SVD + + + +def svd(data, k_singular, N_max=10000): + svd_obj = batched_SVD( data.data, + N_max=N_max, + k_singular=k_singular, + randomize=True) + U,S,VT = svd_obj.go_svd() + return U,S,VT + + + diff --git a/lbl_ir/lbl_ir/tasks/preprocessing/transform.py b/lbl_ir/lbl_ir/tasks/preprocessing/transform.py new file mode 100644 index 0000000..772e87c --- /dev/null +++ b/lbl_ir/lbl_ir/tasks/preprocessing/transform.py @@ -0,0 +1,22 @@ +import numpy as np +from lbl_ir.tasks.baseline import rubberband + +def to_absorbance(data, wavenumbers, eps=1e-5, rubber_band = False): + tmp = data /100.0 + sel_low = tmp < eps + sel_high = tmp > 1-eps + tmp[sel_low] = eps + tmp[sel_high] = 1.0-eps + result = -np.log(tmp) + bg = result*0 + if rubber_band: + shape = result.shape + for ii in range(shape[0]): + for jj in range(shape[1]): + tmp = result[ii,jj,:] + this_bg = rubberband.rubberband( wavenumbers, tmp ) + bg[ii,jj,:]=this_bg + return result, bg + else: + return result,None + diff --git a/lbl_ir/setup.py b/lbl_ir/setup.py new file mode 100644 index 0000000..322a1b1 --- /dev/null +++ b/lbl_ir/setup.py @@ -0,0 +1,100 @@ +""" +Usage: pip install -e . + python setup.py install + python setup.py bdist_wheel + python setup.py sdist bdist_egg + twine upload dist/* +""" + +from setuptools import setup, find_packages + +setup( + name='lbl_ir', + + # Versions should comply with PEP440. For a discussion on single-sourcing + # the version across setup.py and the project code, see + # https://packaging.python.org/en/latest/single_source_version.html + version='0.1.0', + + description='A toolset to analyze OMNIC map files.', + + # The project's main homepage. + url='https://bitbucket.org/berkeleylab/lbl_ir/src/master/', + + # Author details + author='Liang Chen; Petrus Zwart', + author_email='lchen2@lbl.gov', + + # Choose your license + license='BSD', + + # See https://pypi.python.org/pypi?%3Aaction=list_classifiers + classifiers=[ + # How mature is this project? Common values are + # 3 - Alpha + # 4 - Beta + # 5 - Production/Stable + 'Development Status :: 4 - Beta', + + # Indicate who your project is intended for + 'Intended Audience :: Science/Research', + 'Topic :: Scientific/Engineering :: Physics', + + # Pick your license as you wish (should match "license" above) + 'License :: OSI Approved :: BSD License', + + # Specify the Python versions you support here. In particular, ensure + # that you indicate whether you support Python 2, Python 3 or both. + 'Programming Language :: Python :: 3.6' + ], + + # What does your project relate to? + keywords='Infrared Microscopy Analysis', + + # You can just specify the packages manually here if your project is + # simple. Or you can use find_packages(). + packages=find_packages(), + + package_dir={}, + + # Alternatively, if you want to distribute just a my_module.py, uncomment + # this: + # py_modules=["__init__"], + + # List run-time dependencies here. These will be installed by pip when + # your project is installed. For an analysis of "install_requires" vs pip's + # requirements files see: + # https://packaging.python.org/en/latest/requirements.html + install_requires=['numpy', 'scipy', 'seaborn', 'pandas', 'matplotlib', 'spectral', 'h5py', + 'scikit-learn', 'umap-learn'], + + setup_requires=[], + + # List additional groups of dependencies here (e.g. development + # dependencies). You can install these using the following syntax, + # for example: + # $ pip install -e .[dev,tests] + extras_require={ + # 'dev': ['check-manifest'], + }, + + # If there are data files included in your packages that need to be + # installed, specify them here. If using Python 2.6 or less, then these + # have to be included in MANIFEST.in as well. + package_data={}, + + # Although 'package_data' is the preferred approach, in some case you may + # need to place data files outside of your packages. See: + # http://docs.python.org/3.4/distutils/setupscript.html#installing-additional-files # noqa + # In this case, 'data_file' will be installed into '/my_data' + # data_files=[#('lib/python2.7/site-packages/gui', glob.glob('gui/*')), + # ('lib/python2.7/site-packages/yaml/tomography',glob.glob('yaml/tomography/*'))], + + # To provide executable scripts, use entry points in preference to the + # "scripts" keyword. Entry points provide cross-platform support and allow + # pip to create the appropriate form of executable for the target platform. + entry_points={}, + + ext_modules=[], + include_package_data=True +) diff --git a/requirements.txt b/requirements.txt index 2429bca..63a273e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -10,4 +10,7 @@ matplotlib seaborn scikit-learn umap-learn - +pyMCR +xlrd +xlwt +openpyxl \ No newline at end of file diff --git a/setup.py b/setup.py index d44db80..1ef8e3f 100644 --- a/setup.py +++ b/setup.py @@ -35,5 +35,10 @@ author='Liang Chen', install_requires=install_requires, dependency_links=dependency_links, - author_email='lchen2@lbl.gov' + author_email='lchen2@lbl.gov', + entry_points={"xicam.plugins.GUIPlugin": + ["BSISB = xicam.BSISB:BSISB"], + "xicam.plugins.DataHandlerPlugin": + ["mapfile = xicam.BSISB.formats.mapfile:MapFilePlugin"] + } ) diff --git a/tests/PC12-NGF-3h.h5 b/tests/PC12-NGF-3h.h5 new file mode 100644 index 0000000..b03090a Binary files /dev/null and b/tests/PC12-NGF-3h.h5 differ diff --git a/xicam.BSISB.egg-info/PKG-INFO b/xicam.BSISB.egg-info/PKG-INFO deleted file mode 100644 index eed8d82..0000000 --- a/xicam.BSISB.egg-info/PKG-INFO +++ /dev/null @@ -1,46 +0,0 @@ -Metadata-Version: 1.1 -Name: xicam.BSISB -Version: 0.1 -Summary: UNKNOWN -Home-page: UNKNOWN -Author: Liang Chen -Author-email: lchen2@lbl.gov -License: BSD -Description: Bsisb (Xi-cam PluginMaker) - =============================== - - version number: 0.1 - author: Liang Chen - - Overview - -------- - - - - Installation / Usage - -------------------- - - To install use pip: - - $ pip install Xi-cam.plugins.BSISB - - - Or clone the repo: - - $ git clone - $ python setup.py install - - Contributing - ------------ - - TBD - - Example - ------- - - TBD - -Platform: UNKNOWN -Classifier: Development Status :: 3 - Alpha -Classifier: Intended Audience :: Developers -Classifier: Programming Language :: Python :: 3 diff --git a/xicam.BSISB.egg-info/SOURCES.txt b/xicam.BSISB.egg-info/SOURCES.txt deleted file mode 100644 index bf7fb28..0000000 --- a/xicam.BSISB.egg-info/SOURCES.txt +++ /dev/null @@ -1,16 +0,0 @@ -MANIFEST.in -README.md -requirements.txt -setup.cfg -setup.py -xicam.BSISB.egg-info/PKG-INFO -xicam.BSISB.egg-info/SOURCES.txt -xicam.BSISB.egg-info/dependency_links.txt -xicam.BSISB.egg-info/requires.txt -xicam.BSISB.egg-info/top_level.txt -xicam/BSISB/__init__.py -xicam/BSISB/formats/__init__.py -xicam/BSISB/formats/mapfile.py -xicam/BSISB/widgets/__init__.py -xicam/BSISB/widgets/mapviewwidget.py -xicam/BSISB/widgets/spectraplotwidget.py \ No newline at end of file diff --git a/xicam.BSISB.egg-info/requires.txt b/xicam.BSISB.egg-info/requires.txt deleted file mode 100644 index 4d9416b..0000000 --- a/xicam.BSISB.egg-info/requires.txt +++ /dev/null @@ -1,3 +0,0 @@ -nose -coverage -pypi-publisher diff --git a/xicam.BSISB.egg-info/top_level.txt b/xicam.BSISB.egg-info/top_level.txt deleted file mode 100644 index db1e266..0000000 --- a/xicam.BSISB.egg-info/top_level.txt +++ /dev/null @@ -1 +0,0 @@ -xicam diff --git a/xicam/BSISB/BSISB.yapsy-plugin b/xicam/BSISB/BSISB.yapsy-plugin deleted file mode 100644 index 4d507f8..0000000 --- a/xicam/BSISB/BSISB.yapsy-plugin +++ /dev/null @@ -1,9 +0,0 @@ -[Core] -Name = xicam.BSISB -Module = __init__.py - -[Documentation] -Author = Liang Chen -Version = 0.1 -Website = -Description = diff --git a/xicam/BSISB/__init__.py b/xicam/BSISB/__init__.py index 9f054a6..fa70d18 100644 --- a/xicam/BSISB/__init__.py +++ b/xicam/BSISB/__init__.py @@ -1,266 +1,33 @@ -import os from functools import partial from qtpy.QtCore import * from qtpy.QtGui import * -from qtpy.QtWidgets import * -import pickle -import pyqtgraph as pg -import numpy as np +from xicam.core import msg from xicam.core.data import NonDBHeader -from xicam.BSISB.widgets.uiwidget import MsgBox, uiSaveFile, uiGetFile from xicam.BSISB.widgets.mapconvertwidget import mapToH5 from xicam.BSISB.widgets.mapviewwidget import MapViewWidget +from xicam.BSISB.widgets.spectramaproiwidget import MapView from xicam.BSISB.widgets.spectraplotwidget import SpectraPlotWidget from xicam.BSISB.widgets.factorizationwidget import FactorizationWidget - +from xicam.BSISB.widgets.preprocesswidget import PreprocessWidget +from xicam.BSISB.widgets.clusteringwidget import ClusteringWidget from xicam.plugins import GUIPlugin, GUILayout from xicam.gui.widgets.tabview import TabView -from pyqtgraph.parametertree import ParameterTree, Parameter - - - -class MapView(QSplitter): - sigRoiPixels = Signal(object) - sigRoiState = Signal(object) - sigAutoMaskState = Signal(object) - sigSelectMaskState = Signal(object) - - def __init__(self, header: NonDBHeader = None, field: str = 'primary', ): - super(MapView, self).__init__() - # layout set up - self.setOrientation(Qt.Vertical) - self.imageview = MapViewWidget() - self.spectra = SpectraPlotWidget() - - self.imageview_and_toolbar = QSplitter() - self.imageview_and_toolbar.setOrientation(Qt.Horizontal) - self.toolbar_and_param = QSplitter() - self.toolbar_and_param.setOrientation(Qt.Vertical) - #define tool bar - self.toolBar = QWidget() - self.gridlayout = QGridLayout() - self.toolBar.setLayout(self.gridlayout) - #add tool bar buttons - self.roiBtn = QToolButton() - self.roiBtn.setText('Manual ROI') - self.roiBtn.setCheckable(True) - self.roiMeanBtn = QToolButton() - self.roiMeanBtn.setText('ROI Mean') - self.autoMaskBtn = QToolButton() - self.autoMaskBtn.setText('Auto ROI') - self.autoMaskBtn.setCheckable(True) - self.selectMaskBtn = QToolButton() - self.selectMaskBtn.setText('Mark Select') - self.selectMaskBtn.setCheckable(True) - self.saveRoiBtn = QToolButton() - self.saveRoiBtn.setText('Save ROI') - self.saveRoiBtn.setCheckable(False) - self.loadRoiBtn = QToolButton() - self.loadRoiBtn.setText('Load ROI') - self.loadRoiBtn.setCheckable(False) - self.gridlayout.addWidget(self.roiBtn, 0, 0, 1, 1) - self.gridlayout.addWidget(self.autoMaskBtn, 0, 1, 1, 1) - self.gridlayout.addWidget(self.selectMaskBtn, 1, 0, 1, 1) - self.gridlayout.addWidget(self.roiMeanBtn, 1, 1, 1, 1) - self.gridlayout.addWidget(self.saveRoiBtn, 2, 0, 1, 1) - self.gridlayout.addWidget(self.loadRoiBtn, 2, 1, 1, 1) - - self.parameterTree = ParameterTree() - self.parameter = Parameter(name='Threshhold', type='group', - children=[{'name': 'Amide II', - 'value': 0, - 'type': 'float'}, - {'name': "ROI type", - 'values': ['+', '-'], - 'value': '+', - 'type': 'list'}, - ]) - self.parameter.child('Amide II').setOpts(step=0.1) - self.parameterTree.setParameters(self.parameter, showTop=False) - self.parameterTree.setHeaderLabels(['Params','Value']) - self.parameterTree.setIndentation(0) - - # Assemble widgets - self.toolbar_and_param.addWidget(self.toolBar) - self.toolbar_and_param.addWidget(self.parameterTree) - self.toolbar_and_param.setSizes([1000, 1]) #adjust initial splitter size - self.imageview_and_toolbar.addWidget(self.toolbar_and_param) - self.imageview_and_toolbar.addWidget(self.imageview) - self.imageview_and_toolbar.setSizes([1, 1000])#adjust initial splitter size - self.addWidget(self.imageview_and_toolbar) - self.addWidget(self.spectra) - self.setSizes([1000, 1000]) # adjust initial splitter size - # readin header - self.imageview.setHeader(header, field='image') - self.spectra.setHeader(header, field='spectra') - self.header = header - # init pixel selection dict - self.pixSelection = {'ROI': None, 'Mask': None} - - #setup ROI item - sideLen = 10 - self.roi = pg.PolyLineROI(positions=[[0, 0], [sideLen, 0], [sideLen, sideLen], [0, sideLen]], closed=True) - self.imageview.view.addItem(self.roi) - self.roiInitState = self.roi.getState() - self.roi.hide() - - #constants - self.path = os.path.expanduser('~/') - - # Connect signals - self.imageview.sigShowSpectra.connect(self.spectra.showSpectra) - self.spectra.sigEnergyChanged.connect(self.imageview.setEnergy) - self.roiBtn.clicked.connect(self.roiBtnClicked) - self.roi.sigRegionChangeFinished.connect(self.roiSelectPixel) - self.roi.sigRegionChangeFinished.connect(self.showSelectMask) - self.sigRoiPixels.connect(self.spectra.getSelectedPixels) - self.roiMeanBtn.clicked.connect(self.spectra.showMeanSpectra) - self.autoMaskBtn.clicked.connect(self.showAutoMask) - self.selectMaskBtn.clicked.connect(self.showSelectMask) - self.saveRoiBtn.clicked.connect(self.saveRoi) - self.loadRoiBtn.clicked.connect(self.loadRoi) - self.parameter.child('Amide II').sigValueChanged.connect(self.showAutoMask) - self.parameter.child('Amide II').sigValueChanged.connect(self.intersectSelection) - self.parameter.child('ROI type').sigValueChanged.connect(self.intersectSelection) - - def roiBtnClicked(self): - self.roiSelectPixel() - if self.roiBtn.isChecked(): - self.imageview.cross.hide() - self.roi.show() - self.sigRoiState.emit((True, self.roi.getState())) - else: - self.roi.hide() - self.roi.setState(self.roiInitState) - self.sigRoiState.emit((False, self.roi.getState())) - - def saveRoi(self): - parameterDict = {name: self.parameter[name] for name in self.parameter.names.keys()} - roiStates = {'roiBtn': self.roiBtn.isChecked(), 'maskBtn': self.autoMaskBtn.isChecked(), - 'roiState': self.roi.getState(), 'parameter': parameterDict} - filePath, fileName, canceled = uiSaveFile('Save ROI state', self.path, "Pickle Files (*.pkl)") - if not canceled: - with open(filePath + fileName, 'wb') as f: - pickle.dump(roiStates, f) - MsgBox(f'ROI state file was saved! \nFile Location: {filePath + fileName}') - - def loadRoi(self): - filePath, fileName, canceled = uiGetFile('Open ROI state file', self.path, "Pickle Files (*.pkl)") - if not canceled: - with open(filePath + fileName, 'rb') as f: - roiStates = pickle.load(f) - self.roiBtn.setChecked(roiStates['roiBtn']) - self.roi.setState(roiStates['roiState']) - if roiStates['roiBtn']: - self.roi.show() - self.autoMaskBtn.setChecked(roiStates['maskBtn']) - self.selectMaskBtn.setChecked(True) - for k, v in roiStates['parameter'].items(): - self.parameter[k] = v - MsgBox(f'ROI states were loaded from: \n{filePath + fileName}') - else: - return +class BSISBTabview(TabView): - def roiMove(self, roi): - roiState = roi.getState() - self.roi.setState(roiState) - - def getImgShape(self, imgShape): - self.row, self.col = imgShape[0], imgShape[1] - #set up X,Y grid - x = np.linspace(0, self.col - 1, self.col) - y = np.linspace(self.row - 1, 0, self.row) - self.X, self.Y = np.meshgrid(x, y) - self.fullMap = list(zip(self.Y.ravel(), self.X.ravel())) - # setup automask item - self.autoMask = np.ones((self.row, self.col)) - self.autoMaskItem = pg.ImageItem(self.autoMask, axisOrder="row-major", autoLevels=True, opacity=0.3) - self.imageview.view.addItem(self.autoMaskItem) - self.autoMaskItem.hide() - # setup selctmask item to mark selected pixels - self.selectMask = np.ones((self.row, self.col)) - self.selectMaskItem = pg.ImageItem(self.selectMask, axisOrder="row-major", autoLevels=True, opacity=0.3, - lut = np.array([[0, 0, 0], [255, 0, 0]])) - self.imageview.view.addItem(self.selectMaskItem) - self.selectMaskItem.hide() - - def roiSelectPixel(self): - if self.roiBtn.isChecked(): - #get x,y positions list - xPos = self.roi.getArrayRegion(self.X, self.imageview.imageItem) - xPos = np.round(xPos[xPos > 0]) - yPos = self.roi.getArrayRegion(self.Y, self.imageview.imageItem) - yPos = np.round(yPos[yPos > 0]) - - # extract x,y coordinate from selected region - selectedPixels = list(zip(yPos, xPos)) - self.intersectSelection('ROI', selectedPixels) - self.sigRoiState.emit((True, self.roi.getState())) - else: - self.intersectSelection('ROI', None) # no ROI, select all pixels - self.sigRoiState.emit((False, self.roi.getState())) - - def showSelectMask(self): - if self.selectMaskBtn.isChecked(): - # update and show mask - self.selectMaskItem.setImage(self.selectMask) - self.selectMaskItem.show() - self.sigSelectMaskState.emit((True, self.selectMask)) - else: - self.selectMaskItem.hide() - self.sigSelectMaskState.emit((False, self.selectMask)) - - def showAutoMask(self): - if self.autoMaskBtn.isChecked(): - # update and show mask - self.autoMask = self.imageview.makeMask([self.parameter['Amide II']]) - self.autoMaskItem.setImage(self.autoMask) - self.autoMaskItem.show() - # select pixels - mask = self.autoMask.astype(np.bool) - selectedPixels = list(zip(self.Y[mask], self.X[mask])) - self.intersectSelection('Mask', selectedPixels) - self.sigAutoMaskState.emit((True, self.autoMask)) - else: - self.autoMaskItem.hide() - self.autoMask[:, :] = 1 - self.intersectSelection('Mask', None) # no mask, select all pixels - self.sigAutoMaskState.emit((False, self.autoMask)) - - def intersectSelection(self, selector, selectedPixels): - # update pixel selection dict - if (selector == 'ROI') or (selector == 'Mask'): - self.pixSelection[selector] = selectedPixels - # reverse ROI selection - if (self.parameter['ROI type'] == '-') and (self.pixSelection['ROI'] is not None): - roi_copy = self.pixSelection['ROI'] - reverseROI = set(self.fullMap) - set(self.pixSelection['ROI']) - self.pixSelection['ROI'] = list(reverseROI) - - if (self.pixSelection['ROI'] is None) and (self.pixSelection['Mask'] is None): - self.sigRoiPixels.emit(None) # no ROI, select all pixels - self.selectMask = np.ones((self.row, self.col)) - return - elif self.pixSelection['ROI'] is None: - allSelected = set(self.pixSelection['Mask']) #de-duplication of pixels - elif self.pixSelection['Mask'] is None: - allSelected = set(self.pixSelection['ROI']) #de-duplication of pixels - else: - allSelected = set(self.pixSelection['ROI']) & set(self.pixSelection['Mask']) + def __init__(self, *args, **kwargs): + super(BSISBTabview, self).__init__(*args, **kwargs) - allSelected = np.array(list(allSelected), dtype='int') # convert to array - self.selectMask = np.zeros((self.row, self.col)) - if len(allSelected) > 0: - self.selectMask[allSelected[:, 0], allSelected[:, 1]] = 1 - self.selectMask = np.flipud(self.selectMask) - self.sigRoiPixels.emit(allSelected) - # show SelectMask - self.showSelectMask() - #recover ROI selection - if (self.parameter['ROI type'] == '-') and (self.pixSelection['ROI'] is not None): - self.pixSelection['ROI'] = roi_copy + def closeTab(self, i): + newindex = self.currentIndex() + if (i <= self.currentIndex()) and (newindex > 0): + newindex -= 1 + self.removeTab(i) + self.catalogmodel.removeRow(i) + self.setCurrentIndex(newindex) + self.selectionmodel.setCurrentIndex(self.catalogmodel.index(newindex, 0), QItemSelectionModel.Rows + | QItemSelectionModel.ClearAndSelect) class BSISB(GUIPlugin): name = 'BSISB' @@ -273,26 +40,30 @@ def __init__(self, *args, **kwargs): # Selection model self.selectionmodel = QItemSelectionModel(self.headermodel) - - self.PCA_widget = FactorizationWidget(self.headermodel, self.selectionmodel) - self.NMF_widget = FactorizationWidget(self.headermodel, self.selectionmodel) + self.preprocess = PreprocessWidget(self.headermodel, self.selectionmodel) + self.FA_widget = FactorizationWidget(self.headermodel, self.selectionmodel) + self.clusterwidget = ClusteringWidget(self.headermodel, self.selectionmodel) # update headers list when a tab window is closed - self.headermodel.rowsRemoved.connect(partial(self.PCA_widget.setHeader, 'spectra')) - self.headermodel.rowsRemoved.connect(partial(self.NMF_widget.setHeader, 'volume')) + self.headermodel.rowsRemoved.connect(partial(self.FA_widget.setHeader, 'spectra')) # Setup tabviews and update map selection - self.imageview = TabView(self.headermodel, self.selectionmodel, MapView, 'image') + self.imageview = BSISBTabview(self.headermodel, self.selectionmodel, MapView, 'image') + self.imageview.currentChanged.connect(self.updateTab) self.stages = {"MapToH5": GUILayout(self.mapToH5), "Image View": GUILayout(self.imageview), - "PCA": GUILayout(self.PCA_widget), - "NMF": GUILayout(self.NMF_widget)} + "Preprocess": GUILayout(self.preprocess), + "Decomposition": GUILayout(self.FA_widget), + "Clustering": GUILayout(self.clusterwidget)} super(BSISB, self).__init__(*args, **kwargs) def appendHeader(self, header: NonDBHeader, **kwargs): + # get fileName and update status bar + fileName = header.startdoc.get('sample_name', '????') + msg.showMessage(f'Opening {fileName}.h5') # init item - item = QStandardItem(header.startdoc.get('sample_name', '????') + '_' + str(self.headermodel.rowCount())) + item = QStandardItem(fileName + '_' + str(self.headermodel.rowCount())) item.header = header item.selectedPixels = None @@ -302,21 +73,24 @@ def appendHeader(self, header: NonDBHeader, **kwargs): # read out image shape imageEvent = next(header.events(fields=['image'])) imgShape = imageEvent['imgShape'] + rc2ind = imageEvent['rc_index'] # get current MapView widget currentMapView = self.imageview.currentWidget() # transmit imgshape to currentMapView - currentMapView.getImgShape(imgShape) + currentMapView.getImgShape(imgShape, rc2ind) # get xy coordinates of ROI selected pixels currentMapView.sigRoiPixels.connect(partial(self.appendSelection, 'pixel')) currentMapView.sigRoiState.connect(partial(self.appendSelection, 'ROI')) currentMapView.sigAutoMaskState.connect(partial(self.appendSelection, 'autoMask')) currentMapView.sigSelectMaskState.connect(partial(self.appendSelection, 'select')) - self.PCA_widget.setHeader(field='spectra') - self.NMF_widget.setHeader(field='volume') + self.preprocess.setHeader(field='spectra') + self.FA_widget.setHeader(field='spectra') + self.clusterwidget.setHeader(field='spectra') for i in range(4): - self.PCA_widget.roiList[i].sigRegionChangeFinished.connect(self.updateROI) + self.FA_widget.roiList[i].sigRegionChangeFinished.connect(self.updateROI) + self.clusterwidget.roi.sigRegionChangeFinished.connect(self.updateROI) def appendSelection(self, sigCase, sigContent): # get current widget and append selectedPixels to item @@ -325,17 +99,16 @@ def appendSelection(self, sigCase, sigContent): self.headermodel.item(currentItemIdx).selectedPixels = sigContent elif sigCase == 'ROI': self.headermodel.item(currentItemIdx).roiState = sigContent - self.PCA_widget.updateRoiMask() - self.NMF_widget.updateRoiMask() + self.FA_widget.updateRoiMask() + self.clusterwidget.updateRoiMask() elif sigCase == 'autoMask': self.headermodel.item(currentItemIdx).maskState = sigContent - self.PCA_widget.updateRoiMask() - self.NMF_widget.updateRoiMask() + self.FA_widget.updateRoiMask() + self.clusterwidget.updateRoiMask() elif sigCase == 'select': self.headermodel.item(currentItemIdx).selectState = sigContent - self.PCA_widget.updateRoiMask() - self.NMF_widget.updateRoiMask() - + self.FA_widget.updateRoiMask() + self.clusterwidget.updateRoiMask() def updateROI(self, roi): if self.selectionmodel.hasSelection(): @@ -344,7 +117,6 @@ def updateROI(self, roi): selectMapIdx = 0 self.imageview.widget(selectMapIdx).roiMove(roi) - def updateTab(self, tabIdx): - if tabIdx >= 0: - self.selectionmodel.select(self.headermodel.index(tabIdx, 0), QItemSelectionModel.ClearAndSelect) - + def updateTab(self): + # clean up all widgets + self.preprocess.cleanUp() diff --git a/xicam/BSISB/__pycache__/__init__.cpython-37.pyc b/xicam/BSISB/__pycache__/__init__.cpython-37.pyc deleted file mode 100644 index 87b763d..0000000 Binary files a/xicam/BSISB/__pycache__/__init__.cpython-37.pyc and /dev/null differ diff --git a/xicam/BSISB/formats/__pycache__/mapfile.cpython-37.pyc b/xicam/BSISB/formats/__pycache__/mapfile.cpython-37.pyc deleted file mode 100644 index ff5c4d8..0000000 Binary files a/xicam/BSISB/formats/__pycache__/mapfile.cpython-37.pyc and /dev/null differ diff --git a/xicam/BSISB/formats/mapfile.py b/xicam/BSISB/formats/mapfile.py index f41731b..b47073f 100644 --- a/xicam/BSISB/formats/mapfile.py +++ b/xicam/BSISB/formats/mapfile.py @@ -86,7 +86,7 @@ def getImageEvents(cls, path, descriptor_uid): for i in range(n): yield embedded_local_event_doc(descriptor_uid, 'image', cls, (path,), resource_kwargs={'E': i}, - metadata={'wavenumbers': wavenumbers, 'rc_index': rc2ind, 'index_rc': ind2rc, 'imgShape': imgShape}) + metadata={'path': path, 'wavenumbers': wavenumbers, 'rc_index': rc2ind, 'index_rc': ind2rc, 'imgShape': imgShape}) @classmethod def getSpectraDescriptor(cls, path, start_uid): @@ -112,7 +112,7 @@ def getSpectraEvents(cls, path, descriptor_uid): for i in range(n): yield embedded_local_event_doc(descriptor_uid, 'spectra', cls, (path,), resource_kwargs={'i': i}, - metadata={'wavenumbers':wavenumbers, 'rc_index': rc2ind, 'index_rc': ind2rc, 'imgShape':imgShape}) + metadata={'path': path, 'wavenumbers':wavenumbers, 'rc_index': rc2ind, 'index_rc': ind2rc, 'imgShape':imgShape}) @classmethod def ingest(cls, paths): diff --git a/xicam/BSISB/formats/mapfile.yapsy-plugin b/xicam/BSISB/formats/mapfile.yapsy-plugin deleted file mode 100644 index 39da3f9..0000000 --- a/xicam/BSISB/formats/mapfile.yapsy-plugin +++ /dev/null @@ -1,9 +0,0 @@ -[Core] -Name = BSISB Map File Format -Module = mapfile.py - -[Documentation] -Author = Liang Chen -Version = 0.1 -Website = -Description = diff --git a/xicam/BSISB/widgets/__pycache__/__init__.cpython-37.pyc b/xicam/BSISB/widgets/__pycache__/__init__.cpython-37.pyc deleted file mode 100644 index 5c9ad30..0000000 Binary files a/xicam/BSISB/widgets/__pycache__/__init__.cpython-37.pyc and /dev/null differ diff --git a/xicam/BSISB/widgets/__pycache__/mapviewwidget.cpython-37.pyc b/xicam/BSISB/widgets/__pycache__/mapviewwidget.cpython-37.pyc deleted file mode 100644 index ef8c62c..0000000 Binary files a/xicam/BSISB/widgets/__pycache__/mapviewwidget.cpython-37.pyc and /dev/null differ diff --git a/xicam/BSISB/widgets/__pycache__/spectraplotwidget.cpython-37.pyc b/xicam/BSISB/widgets/__pycache__/spectraplotwidget.cpython-37.pyc deleted file mode 100644 index 0a5d9fa..0000000 Binary files a/xicam/BSISB/widgets/__pycache__/spectraplotwidget.cpython-37.pyc and /dev/null differ diff --git a/xicam/BSISB/widgets/clusteringwidget.py b/xicam/BSISB/widgets/clusteringwidget.py new file mode 100644 index 0000000..fd992aa --- /dev/null +++ b/xicam/BSISB/widgets/clusteringwidget.py @@ -0,0 +1,664 @@ +from functools import partial +import os +import numpy as np +import pandas as pd +from lbl_ir.data_objects.ir_map import val2ind +from matplotlib import cm +from pyqtgraph import TextItem, mkBrush, mkPen, ImageItem, PolyLineROI +from pyqtgraph.parametertree import ParameterTree, Parameter +from qtpy.QtCore import Qt, Signal +from qtpy.QtGui import QFont +from qtpy.QtWidgets import * +from sklearn.cluster import KMeans +from sklearn.neighbors import NearestNeighbors +from sklearn.decomposition import PCA +from sklearn.preprocessing import StandardScaler, Normalizer +from umap.umap_ import UMAP +from xicam.BSISB.widgets.mapviewwidget import MapViewWidget, toHtml +from xicam.BSISB.widgets.spectraplotwidget import SpectraPlotWidget +from xicam.BSISB.widgets.uiwidget import MsgBox +from xicam.core import msg + + +class ClusteringParameters(ParameterTree): + sigParamChanged = Signal(object) + + def __init__(self): + super(ClusteringParameters, self).__init__() + + self.parameter = Parameter(name='params', type='group', + children=[{'name': "Embedding", + 'values': ['PCA', 'UMAP'], + 'value': 'UMAP', + 'type': 'list'}, + {'name': "Components", + 'value': 3, + 'type': 'int'}, + {'name': "Neighbors", + 'value': 15, + 'type': 'int'}, + {'name': "Min Dist", + 'value': 0.1, + 'type': 'float'}, + {'name': "Metric", + 'values': ['euclidean', 'manhattan', 'correlation'], + 'value': 'euclidean', + 'type': 'list'}, + {'name': "Normalization", + 'values': ['L2', 'L1', 'None'], + 'value': 'L2', + 'type': 'list'}, + {'name': "Wavenumber Range", + 'value': '400, 4000', + 'type': 'str'}, + {'name': "Clusters", + 'value': 3, + 'type': 'int'}, + {'name': "X Component", + 'values': [1, 2, 3], + 'value': 1, + 'type': 'list'}, + {'name': "Y Component", + 'values': [1, 2, 3], + 'value': 2, + 'type': 'list'} + ]) + self.setParameters(self.parameter, showTop=False) + self.setIndentation(0) + self.parameter.child('Normalization').hide() + # change Fonts + self.fontSize = 12 + font = QFont("Helvetica [Cronyx]", self.fontSize) + boldFont = QFont("Helvetica [Cronyx]", self.fontSize, QFont.Bold) + self.header().setFont(font) + for item in self.listAllItems(): + if hasattr(item, 'widget'): + item.setFont(0, boldFont) + item.widget.setFont(font) + item.displayLabel.setFont(font) + item.widget.setMaximumHeight(40) + + # connect signals + self.parameter.child('Embedding').sigValueChanged.connect(self.updateMethod) + self.parameter.child('Components').sigValueChanged.connect(self.setComponents) + for entry in ['Components', 'Clusters', 'X Component', 'Y Component']: + self.parameter.child(entry).sigValueChanged.connect(partial(self.updateClusterParams, entry)) + + def updateMethod(self): + """ + Toggle parameter menu based on embedding method + :return: None + """ + if self.parameter["Embedding"] == 'UMAP': + self.parameter.child('Neighbors').show() + self.parameter.child('Min Dist').show() + self.parameter.child('Metric').show() + self.parameter.child('Normalization').hide() + else: + self.parameter.child('Neighbors').hide() + self.parameter.child('Min Dist').hide() + self.parameter.child('Metric').hide() + self.parameter.child('Normalization').show() + + def setComponents(self): + N = self.parameter['Components'] + for entry in ['X Component', 'Y Component']: + param = self.parameter.child(entry) + param.setLimits(list(range(1, N + 1))) + + def updateClusterParams(self, name): + self.sigParamChanged.emit(name) + + +class ScatterPlotWidget(SpectraPlotWidget): + sigScatterClicked = Signal(object) + sigScatterRawInd = Signal(object) + + def __init__(self): + super(ScatterPlotWidget, self).__init__(invertX=False, linePos=0) + # self.scene().sigMouseClicked.connect(self.setCrossPos) + self.scene().sigMouseMoved.connect(self.moveCrossPos) + self.line.hide() + self.scatterData = None + self.selectedPixels = None + self.selPx_rc2ind = None + self.selPx_ind2rc = None + self.rc2ind = None + self.ind2rc = None + self.nbr = None + + def setCrossPos(self, event): + pos = event.pos() + if (self.getViewBox().sceneBoundingRect().contains(pos)) and (self.scatterData is not None): + mousePoint = self.getViewBox().mapToView(pos) + x, y = mousePoint.x(), mousePoint.y() + _, ind = self.nbr.kneighbors(np.array([[x, y]])) + self.addItem(self.cross) + self.cross.setData(self.scatterData[ind[0], 0], self.scatterData[ind[0], 1]) + if self.selectedPixels is None: + self.sigScatterClicked.emit(ind[0, 0]) + self.sigScatterRawInd.emit(ind[0, 0]) + else: + row, col = self.selPx_ind2rc[ind[0, 0]] + raw_ind = self.rc2ind[(row, col)] + self.sigScatterClicked.emit(ind[0, 0]) + self.sigScatterRawInd.emit(raw_ind) + + def moveCrossPos(self, pos): + if (self.getViewBox().sceneBoundingRect().contains(pos)) and (self.scatterData is not None): + mousePoint = self.getViewBox().mapSceneToView(pos) + x, y = mousePoint.x(), mousePoint.y() + _, ind = self.nbr.kneighbors(np.array([[x, y]])) + self.addItem(self.cross) + self.cross.setData(self.scatterData[ind[0], 0], self.scatterData[ind[0], 1]) + if self.selectedPixels is None: + self.sigScatterClicked.emit(ind[0, 0]) + self.sigScatterRawInd.emit(ind[0, 0]) + else: + row, col = self.selPx_ind2rc[ind[0, 0]] + raw_ind = self.rc2ind[(row, col)] + self.sigScatterClicked.emit(ind[0, 0]) + self.sigScatterRawInd.emit(raw_ind) + + def clickFromImage(self, ind): + if self.selectedPixels is None: + self.cross.setData([self.scatterData[ind, 0]], [self.scatterData[ind, 1]]) + self.sigScatterClicked.emit(ind) + elif self.ind2rc[ind] in self.selPx_rc2ind: + row, col = self.ind2rc[ind] + ind = self.selPx_rc2ind[(row, col)] + self.cross.setData([self.scatterData[ind, 0]], [self.scatterData[ind, 1]]) + self.sigScatterClicked.emit(ind) + + def getNN(self): + msg.showMessage('Training NearestNeighbors model in scatter plot.') + self.nbr = NearestNeighbors(n_neighbors=1, algorithm='auto').fit(self.scatterData) + msg.showMessage('NearestNeighbors model training is finished.') + + +class ClusterSpectraWidget(SpectraPlotWidget): + def __init__(self): + super(ClusterSpectraWidget, self).__init__(txtPosRatio=0.35) + self._x = None + self.ymax, self.zmax = 0, 100 + self._plots = [] + + def getEnergy(self): + if self._y is not None: + x_val = self.line.value() + idx = val2ind(x_val, self._x) + x_val = self._x[idx] + y_val = self._y[idx] + txt_html = toHtml(f'X = {x_val: .2f}, Y = {y_val: .4f}') + self.txt.setHtml(txt_html) + self.cross.setData([x_val], [y_val]) + + def setColors(self, colorLUT): + self.colorLUT = colorLUT.copy() + # self.colorLUT[0, :] = np.ones(3) * 255 + + def plotClusterSpectra(self): + if self._data is not None: + if self._plots: + self.clearAll() + self._plots = [] + self.ymax = 0 + self.plotItem.addLegend(offset=(1, 1)) + self.nSpectra = len(self._data) + for i in range(self.nSpectra): + name = 'Cluster #' + str(i + 1) + tmp = self.plot(self.wavenumbers, self._data[i], pen=mkPen(self.colorLUT[i + 1], width=2), name=name) + tmp.curve.setClickable(True) + tmp.curve.sigClicked.connect(partial(self.curveHighLight, i)) + self._plots.append(tmp) + self.addItem(self.txt) + + def curveHighLight(self, k): + for i in range(self.nSpectra): + if i == k: + self._plots[i].setPen(mkPen(self.colorLUT[k + 1], width=6)) + self._plots[i].setZValue(50) + else: + self._plots[i].setPen(mkPen(self.colorLUT[i + 1], width=2)) + self._plots[i].setZValue(0) + self._x, self._y = self._plots[k].getData() + ymin, ymax = np.min(self._y), np.max(self._y) + self.getViewBox().setYRange(ymin, ymax, padding=0.1) + r = self.txtPosRatio + self.txt.setPos(r * self._x[-1] + (1 - r) * self._x[0], ymax) + self.getEnergy() + + def plot(self, x, y, *args, **kwargs): + # set up infinity line and get its position + plot_item = self.plotItem.plot(x, y, *args, **kwargs) + self.addItem(self.line) + self.addItem(self.cross) + x_val = self.line.value() + idx = val2ind(x_val, x) + x_val = x[idx] + y_val = y[idx] + txt_html = toHtml(f'X = {x_val: .2f}, Y = {y_val: .4f}') + self.txt = TextItem(html=txt_html, anchor=(0, 0)) + self.txt.setZValue(self.zmax - 1) + self.cross.setData([x_val], [y_val]) + self.cross.setZValue(self.zmax) + ymax = np.max(y) + if ymax > self.ymax: + self.ymax = ymax + self._x, self._y = x, y + r = self.txtPosRatio + self.txt.setPos(r * x[-1] + (1 - r) * x[0], self.ymax) + return plot_item + + +class ClusteringWidget(QSplitter): + def __init__(self, headermodel, selectionmodel): + super(ClusteringWidget, self).__init__() + self.headermodel = headermodel + self.selectionmodel = selectionmodel + # init some values + self.selectMapidx = 0 + self.embedding = None + self.labels = None + self.mean_spectra = None + + # split between cluster image and scatter plot + self.image_and_scatter = QSplitter() + # split between image&scatter and spec plot, vertical split + self.leftsplitter = QSplitter() + self.leftsplitter.setOrientation(Qt.Vertical) + # split between params, buttons and map list, vertical split + self.rightsplitter = QSplitter() + self.rightsplitter.setOrientation(Qt.Vertical) + + self.clusterImage = MapViewWidget() + self.clusterScatterPlot = ScatterPlotWidget() + self.rawSpecPlot = SpectraPlotWidget() + self.clusterMeanPlot = ClusterSpectraWidget() + + # ParameterTree + self.parametertree = ClusteringParameters() + self.parameter = self.parametertree.parameter + + # buttons layout + self.buttons = QWidget() + self.buttonlayout = QGridLayout() + self.buttons.setLayout(self.buttonlayout) + # set up buttons + self.fontSize = 12 + font = QFont("Helvetica [Cronyx]", self.fontSize) + self.computeBtn = QPushButton() + self.computeBtn.setText('Compute clusters') + self.computeBtn.setFont(font) + self.saveBtn = QPushButton() + self.saveBtn.setText('Save clusters') + self.saveBtn.setFont(font) + # add all buttons + self.buttonlayout.addWidget(self.computeBtn) + self.buttonlayout.addWidget(self.saveBtn) + + # Headers listview + self.headerlistview = QListView() + self.headerlistview.setModel(headermodel) + self.headerlistview.setSelectionModel(selectionmodel) # This might do weird things in the map view? + self.headerlistview.setSelectionMode(QListView.SingleSelection) + # add title to list view + self.mapListWidget = QWidget() + self.listLayout = QVBoxLayout() + self.mapListWidget.setLayout(self.listLayout) + mapListTitle = QLabel('Maps list') + mapListTitle.setFont(font) + self.listLayout.addWidget(mapListTitle) + self.listLayout.addWidget(self.headerlistview) + + # assemble widgets + self.image_and_scatter.addWidget(self.clusterImage) + self.image_and_scatter.addWidget(self.clusterScatterPlot) + self.leftsplitter.addWidget(self.image_and_scatter) + self.leftsplitter.addWidget(self.rawSpecPlot) + self.leftsplitter.addWidget(self.clusterMeanPlot) + self.leftsplitter.setSizes([200, 50, 50]) + self.rightsplitter.addWidget(self.parametertree) + self.rightsplitter.addWidget(self.buttons) + self.rightsplitter.addWidget(self.mapListWidget) + self.rightsplitter.setSizes([300, 50, 50]) + self.addWidget(self.leftsplitter) + self.addWidget(self.rightsplitter) + self.setSizes([500, 100]) + + # setup ROI item + sideLen = 10 + self.roi = PolyLineROI(positions=[[0, 0], [sideLen, 0], [sideLen, sideLen], [0, sideLen]], closed=True) + self.roi.hide() + self.roiInitState = self.roi.getState() + # set up mask item + self.maskItem = ImageItem(np.ones((1, 1)), axisOrder="row-major", autoLevels=True, opacity=0.3) + self.maskItem.hide() + # set up select mask item + self.selectMaskItem = ImageItem(np.ones((1, 1)), axisOrder="row-major", autoLevels=True, opacity=0.3, + lut=np.array([[0, 0, 0], [255, 0, 0]])) + self.selectMaskItem.hide() + self.clusterImage.view.addItem(self.roi) + self.clusterImage.view.addItem(self.maskItem) + self.clusterImage.view.addItem(self.selectMaskItem) + + # Connect signals + self.computeBtn.clicked.connect(self.computeEmbedding) + self.saveBtn.clicked.connect(self.saveCluster) + self.clusterImage.sigShowSpectra.connect(self.rawSpecPlot.showSpectra) + self.clusterImage.sigShowSpectra.connect(self.clusterScatterPlot.clickFromImage) + self.clusterScatterPlot.sigScatterRawInd.connect(self.rawSpecPlot.showSpectra) + self.clusterScatterPlot.sigScatterClicked.connect(self.showClusterMean) + self.clusterScatterPlot.sigScatterRawInd.connect(self.setImageCross) + self.parametertree.sigParamChanged.connect(self.updateClusterParams) + self.selectionmodel.selectionChanged.connect(self.updateMap) + self.selectionmodel.selectionChanged.connect(self.updateRoiMask) + + def computeEmbedding(self): + # get current map idx + if not self.isMapOpen(): + return + msg.showMessage('Compute embedding.') + # Select wavenumber region + wavROIList = [] + for entry in self.parameter['Wavenumber Range'].split(','): + try: + wavROIList.append(val2ind(int(entry), self.wavenumbers)) + except: + continue + if len(wavROIList) % 2 == 0: + wavROIList = sorted(wavROIList) + wavROIidx = [] + for i in range(len(wavROIList) // 2): + wavROIidx += list(range(wavROIList[2 * i], wavROIList[2 * i + 1] + 1)) + else: + msg.logMessage('"Wavenumber Range" values must be in pairs', msg.ERROR) + MsgBox('Clustering computation aborted.', 'error') + return + self.wavenumbers_select = self.wavenumbers[wavROIidx] + self.N_w = len(self.wavenumbers_select) + # get current dataset + if self.selectedPixels is None: + n_spectra = len(self.data) + self.dataset = np.zeros((n_spectra, self.N_w)) + for i in range(n_spectra): + self.dataset[i, :] = self.data[i][wavROIidx] + else: + n_spectra = len(self.selectedPixels) + self.dataset = np.zeros((n_spectra, self.N_w)) + for i in range(n_spectra): # i: ith selected pixel + row_col = tuple(self.selectedPixels[i]) + self.dataset[i, :] = self.data[self.rc2ind[row_col]][wavROIidx] + # get parameters and compute embedding + n_components = self.parameter['Components'] + if self.parameter['Embedding'] == 'UMAP': + n_neighbors = self.parameter['Neighbors'] + metric = self.parameter['Metric'] + min_dist = np.clip(self.parameter['Min Dist'], 0, 1) + self.umap = UMAP(n_neighbors=n_neighbors, + min_dist=min_dist, + n_components=n_components, + metric=metric, + random_state=0) + self.embedding = self.umap.fit_transform(self.dataset) + elif self.parameter['Embedding'] == 'PCA': + # normalize and mean center + if self.parameter['Normalization'] == 'L1': # normalize + data_norm = Normalizer(norm='l1').fit_transform(self.dataset) + elif self.parameter['Normalization'] == 'L2': + data_norm = Normalizer(norm='l2').fit_transform(self.dataset) + else: + data_norm = self.dataset + # subtract mean + data_centered = StandardScaler(with_std=False).fit_transform(data_norm) + # Do PCA + self.PCA = PCA(n_components=n_components) + self.PCA.fit(data_centered) + self.embedding = self.PCA.transform(data_centered) + # save embedding to standardModelItem + self.item.embedding = self.embedding + # update cluster map + self.computeCluster() + + def computeCluster(self): + # check if embeddings exist + if self.embedding is None: + return + msg.showMessage('Compute clusters.') + # get num of clusters + n_clusters = self.parameter['Clusters'] + # set colorLUT + self.colorLUT = cm.get_cmap('viridis', n_clusters + 1).colors[:, :3] * 255 + # compute cluster + cluster_object = KMeans(n_clusters=n_clusters, random_state=0).fit(self.embedding) + self.labels = cluster_object.labels_ + 1 + # update cluster image + if self.selectedPixels is None: # full map + self.cluster_map = self.labels.reshape(self.imgShape[0], self.imgShape[1]) + elif self.selectedPixels.size == 0: + self.cluster_map = np.zeros((self.imgShape[0], self.imgShape[1]), dtype=int) + else: + self.cluster_map = np.zeros((self.imgShape[0], self.imgShape[1]), dtype=int) + self.cluster_map[self.selectedPixels[:, 0], self.selectedPixels[:, 1]] = self.labels + self.cluster_map = np.flipud(self.cluster_map) + self.clusterImage.setImage(self.cluster_map, levels=[0, n_clusters]) + # self.clusterImage.setImage(self.cluster_map) + self.clusterImage._image = self.cluster_map + self.clusterImage.rc2ind = self.rc2ind + self.clusterImage.row, self.clusterImage.col = self.imgShape[0], self.imgShape[1] + self.clusterImage.txt.setPos(self.clusterImage.col, 0) + self.clusterImage.cross.show() + # update cluster mean + mean_spectra = [] + self.dfGroups = [] + if self.selectedPixels is None: + n_spectra = len(self.data) + self.dataList = np.zeros((n_spectra, len(self.wavenumbers))) + dataIdx = np.arange(n_spectra) + for i in range(n_spectra): + self.dataList[i] = self.data[i] + else: + n_spectra = len(self.selectedPixels) + self.dataList = np.zeros((n_spectra, len(self.wavenumbers))) + dataIdx = np.zeros(n_spectra, dtype=int) + for i in range(n_spectra): # i: ith selected pixel + row_col = tuple(self.selectedPixels[i]) + dataIdx[i] = self.rc2ind[row_col] + self.dataList[i] = self.data[dataIdx[i]] + + for ii in range(1, n_clusters + 1): + sel = (self.labels == ii) + # save each group spectra to a dataFrame + self.dfGroups.append(pd.DataFrame(self.dataList[sel], columns=self.wavenumbers.tolist(), + index=dataIdx[sel])) + this_mean = np.mean(self.dataset[sel, :], axis=0) + mean_spectra.append(this_mean) + self.mean_spectra = np.vstack(mean_spectra) + self.clusterMeanPlot.setColors(self.colorLUT) + self.clusterMeanPlot._data = self.mean_spectra + self.clusterMeanPlot.wavenumbers = self.wavenumbers_select + self.clusterMeanPlot.plotClusterSpectra() + # update scatter plot + self.updateScatterPlot() + + def saveCluster(self): + if hasattr(self, 'cluster_map') and hasattr(self, 'mean_spectra'): + filePath = self.pathList[self.selectMapidx] + # get dirname and old filename + dirName = os.path.dirname(filePath) + oldFileName = os.path.basename(filePath) + n_clusters = self.parameter['Clusters'] + for i in range(n_clusters): + # save dataFrames to csv file + csvName = oldFileName[:-3] + f'_cluster{i+1}.csv' + newFilePath = os.path.join(dirName, csvName) + self.dfGroups[i].to_csv(newFilePath) + MsgBox(f'Cluster spectra groups were successfully saved at: {newFilePath}!') + + def updateScatterPlot(self): + if (self.embedding is None) or (self.labels is None): + return + # get scatter x, y values + self.clusterScatterPlot.scatterData = self.embedding[:, + [self.parameter['X Component'] - 1, self.parameter['Y Component'] - 1]] + # get colormapings + brushes = [mkBrush(self.colorLUT[x, :]) for x in self.labels] + # make plots + if hasattr(self, 'scatterPlot'): + self.clusterScatterPlot.plotItem.clearPlots() + self.scatterPlot = self.clusterScatterPlot.plotItem.plot(self.clusterScatterPlot.scatterData, pen=None, + symbol='o', symbolBrush=brushes) + self.clusterScatterPlot.getViewBox().autoRange(padding=0.1) + self.clusterScatterPlot.getNN() + + def updateClusterParams(self, name): + if name == 'Components': + self.computeEmbedding() + elif name == 'Clusters': + self.computeCluster() + elif name in ['X Component', 'Y Component']: + self.updateScatterPlot() + + def updateMap(self): + # get current map idx + if not self.selectionmodel.selectedIndexes(): # no map is open + return + else: + self.selectMapidx = self.selectionmodel.selectedIndexes()[0].row() + # get current item + self.item = self.headermodel.item(self.selectMapidx) + if hasattr(self.item, 'embedding'): + # compute embedding + self.computeEmbedding() + else: + # reset custer image and plots + self.cleanUp() + + def showClusterMean(self, i): + if self.mean_spectra is None: + return + self.clusterMeanPlot.curveHighLight(self.labels[i] - 1) + + def setImageCross(self, ind): + row, col = self.ind2rc[ind] + # update cross + self.clusterImage.cross.setData([col + 0.5], [self.imgShape[0] - row - 0.5]) + # update text + self.clusterImage.txt.setHtml(toHtml(f'Point: #{ind}', size=8) + + toHtml(f'X: {col}', size=8) + + toHtml(f'Y: {row}', size=8) + + toHtml(f'Val: {self.clusterImage._image[self.imgShape[0] - row - 1, col] :d}', + size=8)) + + def cleanUp(self): + if self.selectionmodel.hasSelection(): + self.selectMapIdx = self.selectionmodel.selectedIndexes()[0].row() + elif self.headermodel.rowCount() > 0: + self.selectMapIdx = 0 + else: + return + + if hasattr(self, 'imgShapes') and (self.selectMapIdx < len(self.imgShapes)): + # self.clusterImage.clear() + img = np.zeros((self.imgShapes[self.selectMapIdx][0], self.imgShapes[self.selectMapIdx][1])) + self.clusterImage.setImage(img=img) + if hasattr(self, 'scatterPlot'): + self.clusterScatterPlot.plotItem.clearPlots() + self.clusterScatterPlot.scatterData = None + self.rawSpecPlot.clearAll() + self.rawSpecPlot._data = None + self.clusterMeanPlot.clearAll() + self.clusterMeanPlot._data = None + + def updateRoiMask(self): + if self.selectionmodel.hasSelection(): + self.selectMapIdx = self.selectionmodel.selectedIndexes()[0].row() + elif self.headermodel.rowCount() > 0: + self.selectMapIdx = 0 + else: + return + # update roi + try: + roiState = self.headermodel.item(self.selectMapIdx).roiState + if roiState[0]: # roi on + self.roi.show() + else: + self.roi.hide() + # update roi state + self.roi.blockSignals(True) + self.roi.setState(roiState[1]) + self.roi.blockSignals(False) + except Exception: + self.roi.hide() + # update automask + try: + maskState = self.headermodel.item(self.selectMapIdx).maskState + self.maskItem.setImage(maskState[1]) + if maskState[0]: # automask on + self.maskItem.show() + else: + self.maskItem.hide() + except Exception: + pass + # update selectMask + try: + selectMaskState = self.headermodel.item(self.selectMapIdx).selectState + self.selectMaskItem.setImage(selectMaskState[1]) + if selectMaskState[0]: # selectmask on + self.selectMaskItem.show() + else: + self.selectMaskItem.hide() + except Exception: + pass + + def setHeader(self, field: str): + self.headers = [self.headermodel.item(i).header for i in range(self.headermodel.rowCount())] + self.field = field + self.wavenumberList = [] + self.imgShapes = [] + self.rc2indList = [] + self.ind2rcList = [] + self.pathList = [] + self.dataSets = [] + + # get wavenumbers, imgShapes, rc2ind + for header in self.headers: + dataEvent = next(header.events(fields=[field])) + self.wavenumberList.append(dataEvent['wavenumbers']) + self.imgShapes.append(dataEvent['imgShape']) + self.rc2indList.append(dataEvent['rc_index']) + self.ind2rcList.append(dataEvent['index_rc']) + self.pathList.append(dataEvent['path']) + # get raw spectra + data = None + try: # spectra datasets + data = header.meta_array('spectra') + except IndexError: + msg.logMessage('Header object contained no frames with field ''{field}''.', msg.ERROR) + if data is not None: + self.dataSets.append(data) + self.cleanUp() + + def isMapOpen(self): + if not self.selectionmodel.selectedIndexes(): # no map is open + return False + else: + self.selectMapidx = self.selectionmodel.selectedIndexes()[0].row() + # get current data + self.item = self.headermodel.item(self.selectMapidx) + self.selectedPixels = self.item.selectedPixels + self.clusterScatterPlot.selectedPixels = self.selectedPixels + self.currentHeader = self.headers[self.selectMapidx] + self.wavenumbers = self.wavenumberList[self.selectMapidx] + self.rc2ind = self.rc2indList[self.selectMapidx] + self.ind2rc = self.ind2rcList[self.selectMapidx] + self.clusterScatterPlot.ind2rc = self.ind2rc + self.clusterScatterPlot.rc2ind = self.rc2ind + self.imgShape = self.imgShapes[self.selectMapidx] + self.data = self.dataSets[self.selectMapidx] + self.rawSpecPlot.setHeader(self.currentHeader, 'spectra') + if self.selectedPixels is not None: + self.clusterScatterPlot.selPx_rc2ind = {tuple(self.selectedPixels[i]): i for i in range(len(self.selectedPixels))} + self.clusterScatterPlot.selPx_ind2rc = {i: tuple(self.selectedPixels[i]) for i in range(len(self.selectedPixels))} + + return True diff --git a/xicam/BSISB/widgets/factorizationwidget.py b/xicam/BSISB/widgets/factorizationwidget.py index 84f45dd..e41e2ed 100644 --- a/xicam/BSISB/widgets/factorizationwidget.py +++ b/xicam/BSISB/widgets/factorizationwidget.py @@ -1,19 +1,20 @@ -from qtpy.QtWidgets import QSplitter, QGridLayout, QWidget, QListView +from qtpy.QtWidgets import * from xicam.core import msg -from xicam.gui.widgets.imageviewmixins import BetterButtons -from pyqtgraph import PlotWidget, mkPen +from pyqtgraph import PlotWidget, PlotDataItem, TextItem, mkPen, InfiniteLine, ImageItem, PolyLineROI from qtpy.QtCore import Qt, QItemSelectionModel -from qtpy.QtGui import QStandardItemModel +from qtpy.QtGui import QStandardItemModel, QFont from functools import partial from qtpy.QtCore import Signal -from sklearn.decomposition import PCA, NMF +from pymcr.mcr import McrAR +from sklearn.decomposition import PCA, NMF, FastICA from sklearn.preprocessing import StandardScaler, Normalizer import numpy as np import pandas as pd -import pyqtgraph as pg import matplotlib.pyplot as plt import seaborn as sns from xicam.BSISB.widgets.uiwidget import MsgBox +from xicam.BSISB.widgets.imshowwidget import SlimImageView +from xicam.BSISB.widgets.spectraplotwidget import SpectraPlotWidget from lbl_ir.data_objects.ir_map import val2ind from lbl_ir.tasks.preprocessing import data_prep from lbl_ir.tasks.NMF.multi_set_analyses import aggregate_data @@ -23,7 +24,6 @@ class FactorizationParameters(ParameterTree): - sigPCA = Signal(object) def __init__(self, headermodel: QStandardItemModel, selectionmodel: QItemSelectionModel): super(FactorizationParameters, self).__init__() @@ -31,11 +31,13 @@ def __init__(self, headermodel: QStandardItemModel, selectionmodel: QItemSelecti self.selectionmodel = selectionmodel self.parameter = Parameter(name='params', type='group', - children=[{'name': "# of Components", + children=[{'name': "Method", + 'values': ['PCA', 'NMF', 'MCR'], + 'value': 'PCA', + 'type': 'list'}, + {'name': "Components", 'value': 4, 'type': 'int'}, - {'name': "Calculate", - 'type': 'action'}, {'name': "Map 1 Component", 'values': [1, 2, 3, 4], 'value': 1, @@ -52,238 +54,105 @@ def __init__(self, headermodel: QStandardItemModel, selectionmodel: QItemSelecti 'values': [1, 2, 3, 4], 'value': 4, 'type': 'list'}, - {'name': "Wavenumber ROI", + {'name': "Wavenumber Range", 'value': '800,1800', 'type': 'str'}, {'name': "Normalization", 'values': ['L2', 'L1', 'None'], 'value': 'L2', 'type': 'list'}, - {'name': "Save results", - 'type': 'action'} + {'name': "C regressor", + 'values': ['OLS', 'NNLS'], + 'value': 'OLS', + 'type': 'list'} ]) self.setParameters(self.parameter, showTop=False) self.setIndentation(0) - - self.parameter.child('Calculate').sigActivated.connect(self.calculate) - self.parameter.child('Save results').sigActivated.connect(self.saveResults) - self.parameter.child('# of Components').sigValueChanged.connect(self.setNumComponents) - - def setHeader(self, wavenumbers, imgShapes, rc2indList, ind2rcList, field: str): - # get all headers selected - # headers = [self.headermodel.itemFromIndex(i).header for i in self.selectionmodel.selectedRows()] - self.headers = [self.headermodel.item(i).header for i in range(self.headermodel.rowCount())] - - self.field = field - self.wavenumbers = wavenumbers - self.N_w = len(self.wavenumbers) - self.imgShapes = imgShapes - self.rc2indList = rc2indList - self.ind2rcList = ind2rcList - self._dataSets = [] - - if field == 'spectra': # PCA workflow - for header in self.headers: - data = None - try: - data = header.meta_array(self.field) - except IndexError: - msg.logMessage('Header object contained no frames with field ''{field}''.', msg.ERROR) - - if data is not None: - self._dataSets.append(data) - elif field == 'volume': # NMF workflow - self.parameter.child('Normalization').setValue('None') - for header in self.headers: - volumeEvent = next(header.events(fields=['volume'])) - # readin filepath - path = volumeEvent['path'] - self._dataSets.append(path) + #constants + self.method = 'PCA' + self.field = 'spectra' + # change Fonts + self.fontSize = 12 + font = QFont("Helvetica [Cronyx]", self.fontSize) + boldFont = QFont("Helvetica [Cronyx]", self.fontSize, QFont.Bold) + self.header().setFont(font) + for item in self.listAllItems(): + if hasattr(item, 'widget'): + item.setFont(0, boldFont) + item.widget.setFont(font) + item.displayLabel.setFont(font) + item.widget.setMaximumHeight(40) + # connect signals + self.parameter.child('Components').sigValueChanged.connect(self.setNumComponents) + self.parameter.child('Method').sigValueChanged.connect(self.setMethod) + self.parameter.child('C regressor').hide() + + def setMethod(self): + if self.parameter['Method'] == 'PCA': + self.parameter.child('Normalization').setToDefault() + self.parameter.child('Normalization').show() + self.parameter.child('C regressor').hide() + elif self.parameter['Method'] == 'NMF': + self.parameter.child('Normalization').hide() + self.parameter.child('C regressor').hide() + elif self.parameter['Method'] == 'MCR': + self.parameter.child('Normalization').hide() + self.parameter.child('C regressor').show() def setNumComponents(self): - N = self.parameter['# of Components'] + N = self.parameter['Components'] for i in range(4): param = self.parameter.child(f'Map {i + 1} Component') param.setLimits(list(range(1, N + 1))) - def calculate(self): - - N = self.parameter['# of Components'] - - if hasattr(self, '_dataSets'): - wavROIList = [] - for entry in self.parameter['Wavenumber ROI'].split(','): - try: - wavROIList.append(val2ind(int(entry), self.wavenumbers)) - except: - continue - # Select wavenumber region - if len(wavROIList) % 2 == 0: - wavROIList = sorted(wavROIList) - wavROIidx = [] - for i in range(len(wavROIList) // 2): - wavROIidx += list(range(wavROIList[2 * i], wavROIList[2 * i + 1] + 1)) - else: - msg.logMessage('"Wavenumber ROI" values must be in pairs', msg.ERROR) - MsgBox('Factorization computation aborted.', 'error') - return - - self.wavenumbers_select = self.wavenumbers[wavROIidx] - # get map ROI selected region - self.selectedPixelsList = [self.headermodel.item(i).selectedPixels for i in - range(self.headermodel.rowCount())] - self.df_row_idx = [] # row index for dataframe data_fac - - print('Start computing factorization ...') - self.dataRowSplit = [0] # remember the starting/end row positions of each dataset - if self.field == 'spectra': # PCA workflow - self.N_w = len(self.wavenumbers_select) - self._allData = np.empty((0, self.N_w)) - print(self.imgShapes) - - for i, data in enumerate(self._dataSets): # i: map idx - if self.selectedPixelsList[i] is None: - n_spectra = len(data) - tmp = np.zeros((n_spectra, self.N_w)) - for j in range(n_spectra): - tmp[j, :] = data[j][wavROIidx] - self.df_row_idx.append((self.ind2rcList[i][j], j)) - else: - n_spectra = len(self.selectedPixelsList[i]) - tmp = np.zeros((n_spectra, self.N_w)) - for j in range(n_spectra): # j: jth selected pixel - row_col = tuple(self.selectedPixelsList[i][j]) - tmp[j, :] = data[self.rc2indList[i][row_col]][wavROIidx] - self.df_row_idx.append((row_col, self.rc2indList[i][row_col])) - - self.dataRowSplit.append(self.dataRowSplit[-1] + n_spectra) - self._allData = np.append(self._allData, tmp, axis=0) - - # define pop up plots labels - self.fac_method_name = 'PCA' - self.data_fac_name = 'data_PCA' - - if len(self._allData) > 0: - # normalize and mean center - if self.parameter['Normalization'] == 'L1':# normalize - data_norm = Normalizer(norm='l1').fit_transform(self._allData) - elif self.parameter['Normalization'] == 'L2': - data_norm = Normalizer(norm='l2').fit_transform(self._allData) - else: - data_norm = self._allData - #subtract mean - data_centered = StandardScaler(with_std=False).fit_transform(data_norm) - # Do PCA - self.PCA = PCA(n_components=N) - self.PCA.fit(data_centered) - self.data_PCA = self.PCA.transform(data_centered) - # pop up plots - self.popup_plots() - else: - msg.logMessage('The data matrix is empty. No PCA is performed.', msg.ERROR) - MsgBox('The data matrix is empty. No PCA is performed.', 'error') - self.PCA, self.data_PCA = None, None - # emit PCA and transformed data : data_PCA - self.sigPCA.emit((self.wavenumbers_select, self.PCA, self.data_PCA, self.dataRowSplit)) - - elif self.field == 'volume': # NMF workflow - data_files = [] - wav_masks = [] - row_idx = np.array([], dtype='int') - self.allDataRowSplit = [0] # row split for complete datasets - print(self.imgShapes) - - for i, file in enumerate(self._dataSets): - ir_data, fmt = read_map.read_all_formats(file) - n_spectra = ir_data.data.shape[0] - self.allDataRowSplit.append(self.allDataRowSplit[-1] + n_spectra) - data_files.append(ir_data) - ds = data_prep.data_prepper(ir_data) - wav_masks.append(ds.decent_bands) - # row selection - if self.selectedPixelsList[i] is None: - row_idx = np.append(row_idx, np.arange(self.allDataRowSplit[-2], self.allDataRowSplit[-1])) - for k, v in self.rc2indList[i].items(): - self.df_row_idx.append((k, v)) - else: - n_spectra = len(self.selectedPixelsList[i]) - for j in range(n_spectra): - row_col = tuple(self.selectedPixelsList[i][j]) - row_idx = np.append(row_idx, self.allDataRowSplit[-2] + - self.rc2indList[i][row_col]) - self.df_row_idx.append((row_col, self.rc2indList[i][row_col])) - - self.dataRowSplit.append(self.dataRowSplit[-1] + n_spectra) # row split for ROI selected rows - - # define pop up plots labels - self.fac_method_name = 'NMF' - self.data_fac_name = 'data_NMF' - - if len(self.df_row_idx) > 0: - # aggregate datasets - ir_data_agg = aggregate_data(self._dataSets, data_files, wav_masks) - col_idx = list(set(wavROIidx) & set(ir_data_agg.master_wmask)) - self.wavenumbers_select = self.wavenumbers[col_idx] - ir_data_agg.data = ir_data_agg.data[:, col_idx] - ir_data_agg.data = ir_data_agg.data[row_idx, :] - # perform NMF - NMF_obj = NMF(n_components=N) - self.data_NMF = NMF_obj.fit_transform(ir_data_agg.data) - self.NMF = NMF_obj - # pop up plots - self.popup_plots() - else: - msg.logMessage('The data matrix is empty. No NMF is performed.', msg.ERROR) - MsgBox('The data matrix is empty. No NMF is performed.', 'error') - self.NMF, self.data_NMF = None, None - # emit NMF and transformed data : data_NMF - self.sigPCA.emit((self.wavenumbers_select, self.NMF, self.data_NMF, self.dataRowSplit)) - - def popup_plots(self): - labels = [] - for i in range(getattr(self, self.fac_method_name).components_.shape[0]): - labels.append(self.fac_method_name + str(i + 1)) - plt.plot(self.wavenumbers_select, getattr(self, self.fac_method_name).components_[i, :], '-', - label=labels[i]) - loadings_legend = plt.legend(loc='best') - plt.setp(loadings_legend, draggable=True) - plt.xlim([max(self.wavenumbers_select), min(self.wavenumbers_select)]) - - groupLabel = np.zeros((self.dataRowSplit[-1], 1)) - for i in range(len(self.dataRowSplit) - 1): - groupLabel[self.dataRowSplit[i]:self.dataRowSplit[i + 1]] = int(i) - - df_scores = pd.DataFrame(np.append(getattr(self, self.data_fac_name), groupLabel, axis=1), - columns=labels + ['Group label']) - grid = sns.pairplot(df_scores, vars=labels, hue="Group label") - # change legend properties - legend_labels = [] - for i in range(self.headermodel.rowCount()): - if (self.selectedPixelsList[i] is None) or (self.selectedPixelsList[i].size > 0): - legend_labels.append(self.headermodel.item(i).data(0)) - for t, l in zip(grid._legend.texts, legend_labels): t.set_text(l) - plt.setp(grid._legend.get_texts(), fontsize=14) - plt.setp(grid._legend.get_title(), fontsize=14) - plt.setp(grid._legend, bbox_to_anchor=(0.2, 0.95), frame_on=True, draggable=True) - plt.setp(grid._legend.get_frame(), edgecolor='k', linewidth=1, alpha=1) - plt.show() - - def saveResults(self): - if (hasattr(self, 'PCA') and self.PCA is not None) or (hasattr(self, 'NMF') and self.NMF is not None): - name = self.fac_method_name - df_fac_components = pd.DataFrame(getattr(self, name).components_, columns=self.wavenumbers_select) - df_data_fac = pd.DataFrame(getattr(self, self.data_fac_name), index=self.df_row_idx) - df_fac_components.to_csv(name + '_components.csv') - df_data_fac.to_csv(name + '_data.csv') - np.savetxt(name + '_mapRowSplit.csv', np.array(self.dataRowSplit), fmt='%d', delimiter=',') - MsgBox(name + ' components successfully saved!') - else: - MsgBox('No factorization components available.') - +class ComponentPlotWidget(SpectraPlotWidget): + def __init__(self, *args, **kwargs): + super(ComponentPlotWidget, self).__init__(linePos=800, txtPosRatio=0.35, *args, **kwargs) + self.cross = PlotDataItem([800], [0], symbolBrush=(255, 255, 255), symbolPen=(255, 255, 255), + symbol='+', symbolSize=25) + self._x, self._y = None, None + self.ymax, self.zmax = 0, 100 + + def getEnergy(self): + if self._y is not None: + x_val = self.line.value() + idx = val2ind(x_val, self._x) + x_val = self._x[idx] + y_val = self._y[idx] + txt_html = f'
\ + X = {x_val: .2f}, Y = {y_val: .4f}
' + self.txt.setHtml(txt_html) + self.cross.setData([x_val], [y_val]) + + def plot(self, x, y, *args, **kwargs): + # set up infinity line and get its position + plot_item = self.plotItem.plot(x, y, *args, **kwargs) + self.addItem(self.line) + self.addItem(self.cross) + x_val = self.line.value() + idx = val2ind(x_val, x) + x_val = x[idx] + y_val = y[idx] + txt_html = f'
\ + X = {x_val: .2f}, Y = {y_val: .4f}
' + self.txt.setHtml(txt_html) + self.txt.setZValue(self.zmax - 1) + self.cross.setData([x_val], [y_val]) + self.cross.setZValue(self.zmax) + ymax = max(y) + if ymax > self.ymax: + self.ymax = ymax + self._x, self._y = x, y + r = self.txtPosRatio + self.txt.setPos(r * x[-1] + (1 - r) * x[0], self.ymax) + self.addItem(self.txt) + return plot_item class FactorizationWidget(QSplitter): + sigPCA = Signal(object) + def __init__(self, headermodel, selectionmodel): super(FactorizationWidget, self).__init__() self.headermodel = headermodel @@ -299,16 +168,16 @@ def __init__(self, headermodel, selectionmodel): self.gridwidget.setLayout(self.gridlayout) self.display = QSplitter() - self.componentSpectra = PlotWidget() + # self.componentSpectra = PlotWidget() + self.componentSpectra = ComponentPlotWidget() self._plotLegends = self.componentSpectra.addLegend() - self._colors = ['r', 'g', 'b', 'y', 'c', 'm', 'w'] # color for plots - self.componentSpectra.getViewBox().invertX(True) + self._colors = ['r', 'g', 'm', 'y', 'c', 'b', 'w'] # color for plots # self.spectraROI = PlotWidget() - self.NWimage = BetterButtons() - self.NEimage = BetterButtons() - self.SWimage = BetterButtons() - self.SEimage = BetterButtons() + self.NWimage = SlimImageView() + self.NEimage = SlimImageView() + self.SWimage = SlimImageView() + self.SEimage = SlimImageView() # setup ROI item sideLen = 10 self.roiList = [] @@ -321,32 +190,33 @@ def __init__(self, headermodel, selectionmodel): getattr(self, self._imageDict[i]).view.invertY(True) getattr(self, self._imageDict[i]).imageItem.setOpts(axisOrder="row-major") # set up roi item - roi = pg.PolyLineROI(positions=[[0, 0], [sideLen, 0], [sideLen, sideLen], [0, sideLen]], closed=True) + roi = PolyLineROI(positions=[[0, 0], [sideLen, 0], [sideLen, sideLen], [0, sideLen]], closed=True) roi.hide() self.roiInitState = roi.getState() self.roiList.append(roi) # set up mask item - maskItem = pg.ImageItem(np.ones((1,1)), axisOrder="row-major", autoLevels=True, opacity=0.3) + maskItem = ImageItem(np.ones((1,1)), axisOrder="row-major", autoLevels=True, opacity=0.3) maskItem.hide() self.maskList.append(maskItem) # set up select mask item - selectMaskItem = pg.ImageItem(np.ones((1, 1)), axisOrder="row-major", autoLevels=True, opacity=0.3, + selectMaskItem = ImageItem(np.ones((1, 1)), axisOrder="row-major", autoLevels=True, opacity=0.3, lut = np.array([[0, 0, 0], [255, 0, 0]])) selectMaskItem.hide() self.selectMaskList.append(selectMaskItem) + # set up image title + getattr(self, self._imageDict[i]).imageTitle = TextItem() getattr(self, self._imageDict[i]).view.addItem(roi) getattr(self, self._imageDict[i]).view.addItem(maskItem) getattr(self, self._imageDict[i]).view.addItem(selectMaskItem) + getattr(self, self._imageDict[i]).view.addItem(getattr(self, self._imageDict[i]).imageTitle) self.parametertree = FactorizationParameters(headermodel, selectionmodel) self.parameter = self.parametertree.parameter - self.parametertree.sigPCA.connect(self.showComponents) for i in range(4): self.parameter.child(f'Map {i + 1} Component').sigValueChanged.connect( partial(self.updateComponents, i)) self.addWidget(self.display) - self.rightsplitter.addWidget(self.parametertree) self.addWidget(self.rightsplitter) self.display.addWidget(self.gridwidget) self.display.addWidget(self.componentSpectra) @@ -359,12 +229,49 @@ def __init__(self, headermodel, selectionmodel): self.setOrientation(Qt.Horizontal) self.display.setOrientation(Qt.Vertical) + # buttons layout + self.buttons = QWidget() + self.buttonlayout = QGridLayout() + self.buttons.setLayout(self.buttonlayout) + # set up buttons + self.fontSize = 12 + font = QFont("Helvetica [Cronyx]", self.fontSize) + self.computeBtn = QPushButton() + self.computeBtn.setText('Decompose') + self.computeBtn.setFont(font) + self.saveBtn = QPushButton() + self.saveBtn.setText('Save Results') + self.saveBtn.setFont(font) + # add all buttons + self.buttonlayout.addWidget(self.computeBtn) + self.buttonlayout.addWidget(self.saveBtn) + # Headers listview self.headerlistview = QListView() self.headerlistview.setModel(headermodel) self.headerlistview.setSelectionModel(selectionmodel) # This might do weird things in the map view? - self.rightsplitter.addWidget(self.headerlistview) self.headerlistview.setSelectionMode(QListView.SingleSelection) + # add title to list view + self.fontSize = 12 + font = QFont("Helvetica [Cronyx]", self.fontSize) + self.mapListWidget = QWidget() + self.listLayout = QVBoxLayout() + self.mapListWidget.setLayout(self.listLayout) + mapListTitle = QLabel('Maps list') + mapListTitle.setFont(font) + self.listLayout.addWidget(mapListTitle) + self.listLayout.addWidget(self.headerlistview) + + # adjust right splitter + self.rightsplitter.addWidget(self.parametertree) + self.rightsplitter.addWidget(self.buttons) + self.rightsplitter.addWidget(self.mapListWidget) + self.rightsplitter.setSizes([300, 50, 50]) + + #connect signals + self.computeBtn.clicked.connect(self.calculate) + self.saveBtn.clicked.connect(self.saveResults) + self.sigPCA.connect(self.showComponents) def updateRoiMask(self): if self.selectionmodel.hasSelection(): @@ -385,17 +292,15 @@ def updateRoiMask(self): self.roiList[i].blockSignals(True) self.roiList[i].setState(roiState[1]) self.roiList[i].blockSignals(False) - except Exception: for i in range(4): self.roiList[i].hide() - # self.roiList[i].setState(self.roiInitState) - # update mask + # update automask try: maskState = self.headermodel.item(self.selectMapIdx).maskState for i in range(4): self.maskList[i].setImage(maskState[1]) - if maskState[0]: # roi on + if maskState[0]: # automask on self.maskList[i].show() else: self.maskList[i].hide() @@ -406,14 +311,13 @@ def updateRoiMask(self): selectMaskState = self.headermodel.item(self.selectMapIdx).selectState for i in range(4): self.selectMaskList[i].setImage(selectMaskState[1]) - if selectMaskState[0]: # roi on + if selectMaskState[0]: # selectmask on self.selectMaskList[i].show() else: self.selectMaskList[i].hide() except Exception: pass - def updateComponents(self, i): # i is imageview/window number # component_index is the PCA component index @@ -426,10 +330,7 @@ def updateComponents(self, i): # update PCA components if hasattr(self, '_plots'): # update plots - if self.field == 'spectra': - name = 'PCA' + str(component_index) - elif self.field == 'volume': - name = 'NMF' + str(component_index) + name = self.parameter['Method'] + str(component_index) self._plots[i].setData(self.wavenumbers, self._fac.components_[component_index - 1, :], name=name) # update legend label sample, label = self._plotLegends.items[i] @@ -461,7 +362,8 @@ def showComponents(self, fac_obj): # get map ROI selected region self.selectedPixelsList = [self.headermodel.item(i).selectedPixels for i in range(self.headermodel.rowCount())] # clear plots and legends - self.componentSpectra.clear() + self.componentSpectra.getViewBox().clear() + self.componentSpectra.ymax = 0 for sample, label in self._plotLegends.items[:]: self._plotLegends.removeItem(label.text) @@ -471,16 +373,18 @@ def showComponents(self, fac_obj): self._plots = [] for i in range(4): component_index = self.parameter[f'Map {i + 1} Component'] - if self.field == 'spectra': - name = 'PCA' + str(component_index) - elif self.field == 'volume': - name = 'NMF' + str(component_index) + name = self.parameter['Method'] + str(component_index) # show loading plots tmp = self.componentSpectra.plot(self.wavenumbers, self._fac.components_[component_index - 1, :], name=name, pen=mkPen(self._colors[i], width=2)) + tmp.curve.setClickable(True) + tmp.curve.sigClicked.connect(partial(self.curveHighLight, i)) self._plots.append(tmp) # show score plots self.drawMap(component_index, i) + # update the last image and loading plots as a recalculation complete signal + N = self.parameter['Components'] + self.parameter.child(f'Map 4 Component').setValue(N) # clear maps else: tab_idx = self.headermodel.rowCount() - 1 @@ -489,10 +393,6 @@ def showComponents(self, fac_obj): img = np.zeros((self.imgShapes[tab_idx][0], self.imgShapes[tab_idx][1])) getattr(self, self._imageDict[i]).setImage(img=img) - # update the last image and loading plots as a recalculation complete signal - N = self.parameter['# of Components'] - self.parameter.child(f'Map 4 Component').setValue(N) - def drawMap(self, component_index, i): # i is imageview/window number data_slice = self._data_fac[self._dataRowSplit[self.selectMapIdx]:self._dataRowSplit[self.selectMapIdx + 1], @@ -508,6 +408,22 @@ def drawMap(self, component_index, i): 1]] = data_slice img = np.flipud(img) getattr(self, self._imageDict[i]).setImage(img=img) + # set imageTitle + imageTitle = getattr(self, self._imageDict[i]).imageTitle + title = self.parameter['Method'] + str(component_index) + imageTitle.setHtml(f'
{title}
') + imageTitle.setPos(0, -5) + + def curveHighLight(self, k): + for i in range(4): + if i == k: + self._plots[i].setPen(mkPen(self._colors[k], width=6)) + self._plots[i].setZValue(50) + else: + self._plots[i].setPen(mkPen(self._colors[i], width=2)) + self._plots[i].setZValue(0) + self.componentSpectra._x, self.componentSpectra._y = self._plots[k].getData() + self.componentSpectra.getEnergy() def setHeader(self, field: str): @@ -517,16 +433,30 @@ def setHeader(self, field: str): self.imgShapes = [] self.rc2indList = [] self.ind2rcList = [] + self._dataSets = {'spectra': [], 'volume': []} # get wavenumbers, imgShapes for header in self.headers: dataEvent = next(header.events(fields=[field])) self.wavenumbers = dataEvent['wavenumbers'] + self.N_w = len(self.wavenumbers) wavenum_align.append( - (round(self.wavenumbers[0]), len(self.wavenumbers))) # append (first wavenum value, wavenum length) + (round(self.wavenumbers[0]), self.N_w)) # append (first wavenum value, wavenum length) self.imgShapes.append(dataEvent['imgShape']) self.rc2indList.append(dataEvent['rc_index']) self.ind2rcList.append(dataEvent['index_rc']) + # load data + data = None + try: # spectra datasets + data = header.meta_array('spectra') + except IndexError: + msg.logMessage('Header object contained no frames with field ''{field}''.', msg.ERROR) + if data is not None: + self._dataSets['spectra'].append(data) + # NMF path sets + volumeEvent = next(header.events(fields=['volume'])) + path = volumeEvent['path'] # readin filepath + self._dataSets['volume'].append(path) # init maps if len(self.imgShapes) > 0: @@ -536,7 +466,215 @@ def setHeader(self, field: str): MsgBox('Length of wavenumber arrays of displayed maps are not equal. \n' 'Perform PCA or NMF on these maps will lead to error.','warn') - self.parametertree.setHeader(self.wavenumbers, self.imgShapes, self.rc2indList, self.ind2rcList, field=field) + # self.parametertree.setHeader(self.wavenumbers, self.imgShapes, self.rc2indList, self.ind2rcList) + + def calculate(self): + + N = self.parameter['Components'] + #set decompose method + if self.parameter['Method'] == 'PCA': + self.method = 'PCA' + self.field = 'spectra' + elif self.parameter['Method'] == 'NMF': + self.method = 'NMF' + self.field = 'volume' + elif self.parameter['Method'] == 'MCR': + self.method = 'MCR' + self.field = 'spectra' + + if hasattr(self, '_dataSets'): + wavROIList = [] + for entry in self.parameter['Wavenumber Range'].split(','): + try: + wavROIList.append(val2ind(int(entry), self.wavenumbers)) + except: + continue + # Select wavenumber region + if len(wavROIList) % 2 == 0: + wavROIList = sorted(wavROIList) + wavROIidx = [] + for i in range(len(wavROIList) // 2): + wavROIidx += list(range(wavROIList[2 * i], wavROIList[2 * i + 1] + 1)) + else: + msg.logMessage('"Wavenumber Range" values must be in pairs', msg.ERROR) + MsgBox('Factorization computation aborted.', 'error') + return + + self.wavenumbers_select = self.wavenumbers[wavROIidx] + # get map ROI selected region + self.selectedPixelsList = [self.headermodel.item(i).selectedPixels for i in + range(self.headermodel.rowCount())] + self.df_row_idx = [] # row index for dataframe data_fac + + msg.showMessage('Start computing', self.method + '. Image shape:', str(self.imgShapes)) + self.dataRowSplit = [0] # remember the starting/end row positions of each dataset + if self.field == 'spectra': # PCA workflow + self.N_w = len(self.wavenumbers_select) + self._allData = np.empty((0, self.N_w)) + + for i, data in enumerate(self._dataSets['spectra']): # i: map idx + if self.selectedPixelsList[i] is None: + n_spectra = len(data) + tmp = np.zeros((n_spectra, self.N_w)) + for j in range(n_spectra): + tmp[j, :] = data[j][wavROIidx] + self.df_row_idx.append((self.ind2rcList[i][j], j)) + else: + n_spectra = len(self.selectedPixelsList[i]) + tmp = np.zeros((n_spectra, self.N_w)) + for j in range(n_spectra): # j: jth selected pixel + row_col = tuple(self.selectedPixelsList[i][j]) + tmp[j, :] = data[self.rc2indList[i][row_col]][wavROIidx] + self.df_row_idx.append((row_col, self.rc2indList[i][row_col])) + + self.dataRowSplit.append(self.dataRowSplit[-1] + n_spectra) + self._allData = np.append(self._allData, tmp, axis=0) + + if len(self._allData) > 0: + if self.method == 'PCA': + self.data_fac_name = 'data_PCA' # define pop up plots labels + # normalize and mean center + if self.parameter['Normalization'] == 'L1':# normalize + data_norm = Normalizer(norm='l1').fit_transform(self._allData) + elif self.parameter['Normalization'] == 'L2': + data_norm = Normalizer(norm='l2').fit_transform(self._allData) + else: + data_norm = self._allData + #subtract mean + data_centered = StandardScaler(with_std=False).fit_transform(data_norm) + # Do PCA + self.PCA = PCA(n_components=N) + self.PCA.fit(data_centered) + self.data_PCA = self.PCA.transform(data_centered) + # pop up plots + self.popup_plots() + elif self.method == 'MCR': + self.data_fac_name = 'data_MCR' # define pop up plots labels + # Do ICA to find initial estimate of ST matrix + self.ICA = FastICA(n_components=N) + self.ICA.fit(self._allData) + # Do MCR + self.MCR = McrAR(max_iter=100, c_regr=self.parameter['C regressor'], st_regr='NNLS', + tol_err_change=1e-6, tol_increase=0.5) + self.MCR.fit(self._allData, ST=self.ICA.components_) + self.MCR.components_ = self.MCR.ST_opt_ + self.data_MCR = self.MCR.C_opt_ + #test ICA + # self.MCR = self.ICA + # self.data_MCR = self.ICA.transform(self._allData) + # pop up plots + self.popup_plots() + else: + msg.logMessage('The data matrix is empty. No PCA is performed.', msg.ERROR) + MsgBox('The data matrix is empty. No PCA is performed.', 'error') + self.PCA, self.data_PCA = None, None + self.MCR, self.data_MCR = None, None + # emit PCA and transformed data + if self.method == 'PCA': + self.sigPCA.emit((self.wavenumbers_select, self.PCA, self.data_PCA, self.dataRowSplit)) + elif self.method == 'MCR': + self.sigPCA.emit((self.wavenumbers_select, self.MCR, self.data_MCR, self.dataRowSplit)) + + elif self.field == 'volume': # NMF workflow + data_files = [] + wav_masks = [] + row_idx = np.array([], dtype='int') + self.allDataRowSplit = [0] # row split for complete datasets + + for i, file in enumerate(self._dataSets['volume']): + ir_data, fmt = read_map.read_all_formats(file) + n_spectra = ir_data.data.shape[0] + self.allDataRowSplit.append(self.allDataRowSplit[-1] + n_spectra) + data_files.append(ir_data) + ds = data_prep.data_prepper(ir_data) + wav_masks.append(ds.decent_bands) + # row selection + if self.selectedPixelsList[i] is None: + row_idx = np.append(row_idx, np.arange(self.allDataRowSplit[-2], self.allDataRowSplit[-1])) + for k, v in self.rc2indList[i].items(): + self.df_row_idx.append((k, v)) + else: + n_spectra = len(self.selectedPixelsList[i]) + for j in range(n_spectra): + row_col = tuple(self.selectedPixelsList[i][j]) + row_idx = np.append(row_idx, self.allDataRowSplit[-2] + + self.rc2indList[i][row_col]) + self.df_row_idx.append((row_col, self.rc2indList[i][row_col])) + + self.dataRowSplit.append(self.dataRowSplit[-1] + n_spectra) # row split for ROI selected rows + + # define pop up plots labels + self.data_fac_name = 'data_NMF' + + if len(self.df_row_idx) > 0: + # aggregate datasets + ir_data_agg = aggregate_data(self._dataSets['volume'], data_files, wav_masks) + col_idx = list(set(wavROIidx) & set(ir_data_agg.master_wmask)) + self.wavenumbers_select = self.wavenumbers[col_idx] + ir_data_agg.data = ir_data_agg.data[:, col_idx] + ir_data_agg.data = ir_data_agg.data[row_idx, :] + # perform NMF + self.NMF = NMF(n_components=N) + self.data_NMF = self.NMF.fit_transform(ir_data_agg.data) + # pop up plots + self.popup_plots() + else: + msg.logMessage('The data matrix is empty. No NMF is performed.', msg.ERROR) + MsgBox('The data matrix is empty. No NMF is performed.', 'error') + self.NMF, self.data_NMF = None, None + # emit NMF and transformed data : data_NMF + self.sigPCA.emit((self.wavenumbers_select, self.NMF, self.data_NMF, self.dataRowSplit)) + + def popup_plots(self): + # component variance ratio plot + if self.method == 'PCA': + plt.plot(getattr(self, self.method).explained_variance_ratio_, 'o-b') + ax = plt.gca() + ax.set_ylabel('Explained variance ratio', fontsize=16) + ax.set_xlabel('Component number', fontsize=16) + ax.set_xticks(np.arange(self.parameter['Components'])) + # loadings plot + plt.figure() + labels = [] + for i in range(getattr(self, self.method).components_.shape[0]): + labels.append(self.method + str(i + 1)) + plt.plot(self.wavenumbers_select, getattr(self, self.method).components_[i, :], '-', + label=labels[i]) + loadings_legend = plt.legend(loc='best') + plt.setp(loadings_legend, draggable=True) + plt.xlim([max(self.wavenumbers_select), min(self.wavenumbers_select)]) + ax = plt.gca() + ax.set_xlabel('Wavenumber$(cm^{-1})$', fontsize=16) + # matrix plot + groupLabel = np.zeros((self.dataRowSplit[-1], 1)) + for i in range(len(self.dataRowSplit) - 1): + groupLabel[self.dataRowSplit[i]:self.dataRowSplit[i + 1]] = int(i) + df_scores = pd.DataFrame(np.append(getattr(self, self.data_fac_name), groupLabel, axis=1), + columns=labels + ['Group label']) + grid = sns.pairplot(df_scores, vars=labels, hue="Group label") + # change legend properties + legend_labels = [] + for i in range(self.headermodel.rowCount()): + if (self.selectedPixelsList[i] is None) or (self.selectedPixelsList[i].size > 0): + legend_labels.append(self.headermodel.item(i).data(0)) + for t, l in zip(grid._legend.texts, legend_labels): t.set_text(l) + plt.setp(grid._legend.get_texts(), fontsize=14) + plt.setp(grid._legend.get_title(), fontsize=14) + plt.setp(grid._legend, bbox_to_anchor=(0.2, 0.95), frame_on=True, draggable=True) + plt.setp(grid._legend.get_frame(), edgecolor='k', linewidth=1, alpha=1) + plt.show() + def saveResults(self): + if (hasattr(self, 'PCA') and self.PCA is not None) or (hasattr(self, 'NMF') and self.NMF is not None)\ + or (hasattr(self, 'MCR') and self.MCR is not None): + name = self.method + df_fac_components = pd.DataFrame(getattr(self, name).components_, columns=self.wavenumbers_select) + df_data_fac = pd.DataFrame(getattr(self, self.data_fac_name), index=self.df_row_idx) + df_fac_components.to_csv(name + '_components.csv') + df_data_fac.to_csv(name + '_data.csv') + np.savetxt(name + '_mapRowSplit.csv', np.array(self.dataRowSplit), fmt='%d', delimiter=',') + MsgBox(name + ' components successfully saved!') + else: + MsgBox('No factorization components available.') diff --git a/xicam/BSISB/widgets/imshowwidget.py b/xicam/BSISB/widgets/imshowwidget.py new file mode 100644 index 0000000..3f590fd --- /dev/null +++ b/xicam/BSISB/widgets/imshowwidget.py @@ -0,0 +1,33 @@ +from xicam.gui.widgets.imageviewmixins import BetterButtons +import numpy as np + + +class SlimImageView(BetterButtons): + def __init__(self, invertY=True): + super(SlimImageView, self).__init__() + # Shrink LUT + self.getHistogramWidget().setMinimumWidth(1) + # set up layout + # self.ui.gridLayout.addWidget(self.resetAxesBtn, 2, 2, 1, 1) + # self.ui.gridLayout.addWidget(self.resetLUTBtn, 3, 2, 1, 1) + # self.ui.gridLayout.addWidget(self.ui.graphicsView, 0, 0, 4, 1) + # set up colorbar + self.setPredefinedGradient("viridis") + self.view.invertY(invertY) + self.imageItem.setOpts(axisOrder="row-major") + # Setup late signal + self.sigTimeChangeFinished = self.timeLine.sigPositionChangeFinished + + def quickMinMax(self, data): + """ + Estimate the min/max values of *data* by subsampling. MODIFIED TO USE THE 99TH PERCENTILE instead of max. + """ + if data is None: + return 0, 0 + ax = np.argmax(data.shape) + sl = [slice(None)] * data.ndim + sl[ax] = slice(None, None, max(1, int(data.size // 1e4))) + data = data[sl] + return (np.nanmin(data), np.nanpercentile(np.where(data < np.nanmax(data), data, np.nanmin(data)), 99)) + + diff --git a/xicam/BSISB/widgets/mapconvertwidget.py b/xicam/BSISB/widgets/mapconvertwidget.py index 3753a54..5b10d24 100644 --- a/xicam/BSISB/widgets/mapconvertwidget.py +++ b/xicam/BSISB/widgets/mapconvertwidget.py @@ -7,9 +7,9 @@ from xicam.BSISB.widgets.uiwidget import MsgBox, YesNoDialog, uiGetFile, uiGetDir, uiSaveFile from xicam.BSISB.widgets.mapviewwidget import MapViewWidget from xicam.BSISB.widgets.spectraplotwidget import SpectraPlotWidget - from lbl_ir.data_objects import ir_map from lbl_ir.io_tools.read_omnic import read_and_convert +from lbl_ir.io_tools.read_numpy import read_npy class mapToH5(QSplitter): def __init__(self): @@ -34,8 +34,10 @@ def __init__(self): self.info.layout().addWidget(QLabel('Status Info:')) self.info.layout().addWidget(self.infoBox) # add tool bar buttons - self.openBtn = QToolButton() - self.openBtn.setText('Open Map') + self.openMapBtn = QToolButton() + self.openMapBtn.setText('Open Map') + self.openNpyBtn = QToolButton() + self.openNpyBtn.setText('Open Npy') self.saveBtn = QToolButton() self.saveBtn.setText('Save HDF5') self.batchBtn = QToolButton() @@ -47,7 +49,8 @@ def __init__(self): self.T2AConvert.setText('Auto T->A') self.T2AConvert.setChecked(True) # Assemble widgets - self.toollayout.addWidget(self.openBtn) + self.toollayout.addWidget(self.openMapBtn) + self.toollayout.addWidget(self.openNpyBtn) self.toollayout.addWidget(self.saveBtn) self.toollayout.addWidget(self.batchBtn) self.toollayout.addWidget(QLabel('Sample Name:')) @@ -67,19 +70,31 @@ def __init__(self): # Connect signals self.imageview.sigShowSpectra.connect(self.spectra.showSpectra) self.spectra.sigEnergyChanged.connect(self.imageview.setEnergy) - self.openBtn.clicked.connect(self.openBtnClicked) + self.openMapBtn.clicked.connect(self.openBtnClicked) + self.openNpyBtn.clicked.connect(self.openNpy) self.saveBtn.clicked.connect(self.saveBtnClicked) self.batchBtn.clicked.connect(self.batchBtnClicked) # Constants self.path = os.path.dirname(sys.path[1]) self.minYLimit = 5 self.epsilon = 1e-10 + self.fileFormat = 'map' + + def openNpy(self): + self.fileFormat = 'npy' + self.T2AConvert.setChecked(False) + self.openBtnClicked() + self.fileFormat = 'map' def openBtnClicked(self): # open omnic map file - self.filePath, self.fileName, canceled = uiGetFile('Open map file', self.path, "Omnic Map Files (*.map)") + if self.fileFormat == 'map': + self.filePath, self.fileName, canceled = uiGetFile('Open map file', self.path, "Omnic Map Files (*.map)") + elif self.fileFormat == 'npy': + self.filePath, self.fileName, canceled = uiGetFile('Open npy file', self.path, "Numpy array Files (*.npy)") + if canceled: - self.infoBox.setText('Open map canceled.') + self.infoBox.setText('Open file canceled.') return # set sample_id if self.sampleName.text() == 'None': @@ -88,31 +103,35 @@ def openBtnClicked(self): sample_info = ir_map.sample_info(sample_id=self.sampleName.text()) #try to open omnic map file try: - self.irMap = read_and_convert(self.filePath + self.fileName, sample_info=sample_info) + if self.fileFormat == 'map': + self.irMap = read_and_convert(self.filePath + self.fileName, sample_info=sample_info) + elif self.fileFormat == 'npy': + self.irMap = read_npy(self.filePath + self.fileName, sample_info=sample_info) except Exception as error: self.infoBox.setText(error.args[0] + f'\nFailed to open file: {self.fileName}.') else: - # check whether to perform T->A conversion - spec0 = self.irMap.imageCube[0, 0, :] - maxSpecY = np.max(spec0) - if (not self.T2AConvert.isChecked()) and (maxSpecY >= self.minYLimit): - userMsg = YesNoDialog(f'max(Y) of the first spectrum is greater than {self.minYLimit}, \ - while the "Auto T->A" box is not checked. \nPlease make sure data format is in absorbance.\ - \nDo you want to perform "Auto T->A" conversion?') - userChoice = userMsg.choice() - if userChoice == QMessageBox.Yes: - self.T2AConvert.setChecked(True) + if self.fileFormat == 'map': + # check whether to perform T->A conversion + spec0 = self.irMap.imageCube[0, 0, :] + maxSpecY = np.max(spec0) + if (not self.T2AConvert.isChecked()) and (maxSpecY >= self.minYLimit): + userMsg = YesNoDialog(f'max(Y) of the first spectrum is greater than {self.minYLimit}, \ + while the "Auto T->A" box is not checked. \nPlease make sure data format is in absorbance.\ + \nDo you want to perform "Auto T->A" conversion?') + userChoice = userMsg.choice() + if userChoice == QMessageBox.Yes: + self.T2AConvert.setChecked(True) + self.irMap.imageCube = -np.log10(self.irMap.imageCube / 100 + self.epsilon) + self.irMap.data = -np.log10(self.irMap.data / 100 + self.epsilon) + self.infoBox.setText(f'User chooses to perform T->A conversion in {self.fileName}.') + else: + self.infoBox.setText(f'User chooses not to perform T->A conversion in {self.fileName}.') + elif maxSpecY >= self.minYLimit: self.irMap.imageCube = -np.log10(self.irMap.imageCube / 100 + self.epsilon) self.irMap.data = -np.log10(self.irMap.data / 100 + self.epsilon) - self.infoBox.setText(f'User chooses to perform T->A conversion in {self.fileName}.') + self.infoBox.setText(f'T->A conversion is performed in {self.fileName}.') else: - self.infoBox.setText(f'User chooses not to perform T->A conversion in {self.fileName}.') - elif maxSpecY >= self.minYLimit: - self.irMap.imageCube = -np.log10(self.irMap.imageCube / 100 + self.epsilon) - self.irMap.data = -np.log10(self.irMap.data / 100 + self.epsilon) - self.infoBox.setText(f'T->A conversion is performed in {self.fileName}.') - else: - self.infoBox.setText(f"{self.fileName}'s datatype is absorbance. \nT->A conversion is not performed.") + self.infoBox.setText(f"{self.fileName}'s datatype is absorbance. \nT->A conversion is not performed.") self.dataCube = np.moveaxis(np.flipud(self.irMap.imageCube), -1, 0) # set up required data/properties in self.imageview @@ -129,6 +148,8 @@ def updateImage(self, row, col, wavenumbers, rc2ind, dataCube): self.imageview.row, self.imageview.col = row, col self.imageview.wavenumbers = wavenumbers self.imageview.rc2ind = rc2ind + self.imageview._data = dataCube + self.imageview._image = self.imageview._data[0] self.imageview.setImage(img=dataCube) def saveBtnClicked(self): diff --git a/xicam/BSISB/widgets/mapviewwidget.py b/xicam/BSISB/widgets/mapviewwidget.py index 1051bec..7ac93c7 100644 --- a/xicam/BSISB/widgets/mapviewwidget.py +++ b/xicam/BSISB/widgets/mapviewwidget.py @@ -1,12 +1,15 @@ import numpy as np -from xicam.gui.widgets.dynimageview import DynImageView +from xicam.BSISB.widgets.imshowwidget import SlimImageView from xicam.core import msg from xicam.core.data import NonDBHeader from pyqtgraph import ArrowItem, TextItem, PlotDataItem from qtpy.QtCore import Signal from lbl_ir.data_objects.ir_map import val2ind -class MapViewWidget(DynImageView): +def toHtml(txt, size=12): + return f'
{txt}
' + +class MapViewWidget(SlimImageView): sigShowSpectra = Signal(int) def __init__(self, *args, **kwargs): @@ -14,12 +17,8 @@ def __init__(self, *args, **kwargs): # self.scene.sigMouseMoved.connect(self.showSpectra) self.scene.sigMouseClicked.connect(self.showSpectra) self.view.invertY(True) - # add arrow - # self.arrow = ArrowItem(angle=60, headLen=15, tipAngle=45, baseAngle=30, brush = (200, 80, 20)) - # self.arrow.setPos(0, 0) self.cross = PlotDataItem([0], [0], symbolBrush=(200, 0, 0), symbolPen=(200, 0, 0), symbol='+', symbolSize=16) - self.view.addItem(self.cross) self.cross.hide() #add txt @@ -29,12 +28,11 @@ def __init__(self, *args, **kwargs): def setEnergy(self, lineobject): E = lineobject.value() # map E to index - i = val2ind(E, self.wavenumbers) - # print('E:', E, 'wav:', self.wavenumbers[i]) - self.setCurrentIndex(i) + idx = val2ind(E, self.wavenumbers) + self._image = self._data[idx] + self.setCurrentIndex(idx) def showSpectra(self, event): - pos = event.pos() if self.view.sceneBoundingRect().contains(pos): # Note, when axes are added, you must get the view with self.view.getViewBox() mousePoint = self.view.mapSceneToView(pos) @@ -44,14 +42,15 @@ def showSpectra(self, event): ind = self.rc2ind[(y,x)] self.sigShowSpectra.emit(ind) # print(x, y, ind, x + y * self.n_col) - - #update arrow + #update crosshair self.cross.setData([x + 0.5], [self.row - y - 0.5]) self.cross.show() # update text - self.txt.setHtml(f'
X: {x}
\ -
Y: {y}
\ -
Point: #{ind}
') + self.txt.setHtml(toHtml(f'Point: #{ind}', size=8) + + toHtml(f'X: {x}', size=8) + + toHtml(f'Y: {y}', size=8) + + toHtml(f'Val: {self._image[self.row - y -1, x]: .4f}', size=8) + ) except Exception: self.cross.hide() @@ -77,6 +76,7 @@ def setHeader(self, header: NonDBHeader, field: str, *args, **kwargs): # kwargs['transform'] = QTransform(1, 0, 0, -1, 0, data.shape[-2]) self.setImage(img=data, *args, **kwargs) self._data = data + self._image = self._data[0] def updateImage(self, autoHistogramRange=True): super(MapViewWidget, self).updateImage(autoHistogramRange) diff --git a/xicam/BSISB/widgets/preprocesswidget.py b/xicam/BSISB/widgets/preprocesswidget.py new file mode 100644 index 0000000..402f651 --- /dev/null +++ b/xicam/BSISB/widgets/preprocesswidget.py @@ -0,0 +1,679 @@ +import os +import numpy as np +import pandas as pd +from scipy.interpolate import interp1d +from functools import partial +from qtpy.QtWidgets import * +from qtpy.QtCore import Qt, QItemSelectionModel, Signal +from qtpy.QtGui import QStandardItemModel, QStandardItem, QFont +from pyqtgraph.parametertree import ParameterTree, Parameter +from xicam.core import msg +from lbl_ir.data_objects.ir_map import ir_map, val2ind +from lbl_ir.tasks.preprocessing.EMSC import Kohler_zero +from xicam.BSISB.widgets.spectraplotwidget import baselinePlotWidget +from xicam.BSISB.widgets.uiwidget import MsgBox, YesNoDialog + + +class Preprocessor: + def __init__(self, wavenumbers, spectrum): + self.wavenumbers = wavenumbers + self.spectrum = spectrum + self.preprocess_method = None + self.interp_method = None + self.wav_anchor = None + + def parse_anchors(self, anchors): + """ + parse anchor points str to real valued wavenumbers and confine energy range + :return: None + """ + anchor_idx = [] + for entry in anchors.split(','): + try: + anchor_idx.append(val2ind(int(entry.strip()), self.wavenumbers)) + except: + continue + anchor_idx = sorted(anchor_idx) + + self.energy = self.wavenumbers[anchor_idx[0]: anchor_idx[-1] + 1] + self.specTrim = self.spectrum[anchor_idx[0]: anchor_idx[-1] + 1] + + self.wav_anchor = self.wavenumbers[anchor_idx] + self.spec_anchor = self.spectrum[anchor_idx] + return None + + def isBaseFitOK(self, anchors, kind, w_regions): + """ + Check is there is enough anchor points to fit higher order baseline + :param anchors: anchor points + :param kind: fitting method + :return: decide if there is enough anchor points for baseline fit + """ + # parse anchor points + self.parse_anchors(anchors) + if self.preprocess_method == 'rubberband': + # decide if num of anchor points is enough for 'quadratic' or 'cubic' fit + if len(self.wav_anchor) < 2: + MsgBox('Baseline fitting needs at least 2 anchor points.\n' + + 'Please add more "anchor points" to correctly fit the baseline.', type='error') + return False + elif len(self.wav_anchor) < 3 and kind == 'quadratic': + MsgBox('Quadratic baseline needs more than 2 anchor points.\n' + + 'Please add more "anchor points" to correctly fit the baseline.', type='error') + return False + elif len(self.wav_anchor) < 4 and kind == 'cubic': + MsgBox('Cubic baseline needs more than 3 anchor points.\n' + + 'Please add more "anchor points" to correctly fit the baseline.', type='error') + return False + else: + return True + elif self.preprocess_method == 'kohler': + try: # read w_regions + self.w_regions = eval(w_regions) + except : + MsgBox('Fitting regions format is not correct.\n' + + 'Please consult default values.', type='error') + if self.w_regions is not None: + return True + else: + return False + + + def rubber_band(self, anchors, kind='linear', w_regions=None): + """ + Calculate rubberBaseline, debased spectrum and 2nd, 4th order derivative of the spectrum + :param anchors: rubberband anchor points + :param kind: spline fit curve order + :return: baseline fit success (bool) + """ + self.preprocess_method = 'rubberband' + self.interp_method = kind + # get rubberBaseline and debased spectrum + if not self.isBaseFitOK(anchors, kind, w_regions): + return False + + f = interp1d(self.wav_anchor, self.spec_anchor, kind=kind) + self.rubberBaseline = f(self.energy) + self.rubberDebased = self.specTrim - self.rubberBaseline + + # get 2nd order derivatives + self.get_derivative() + return True + + def kohler(self, anchors=None, kind=None, w_regions=None): + self.preprocess_method = 'kohler' + # get kohler baseline and debased spectrum + if not self.isBaseFitOK(anchors, kind, w_regions): + return False + + self.kohlerDebased, self.kohlerBaseline = Kohler_zero(self.energy, self.specTrim, self.w_regions) + # get 2nd order derivatives + self.get_derivative() + return True + + def get_derivative(self): + """ + Calculate 2nd and 4th order derivative of the spectrum + :param n: the derivative order + :return: None + """ + dx = self.energy[1] - self.energy[0] + if self.preprocess_method == 'rubberband': + self.deriv2_rubber = self.nthOrderGradient(dx, self.specTrim, n=2) + elif self.preprocess_method == 'kohler': + self.deriv2_kohler = self.nthOrderGradient(dx, self.kohlerDebased, n=2) + return None + + def nthOrderGradient(self, dx, y, n=1): + """ + Calculate nth order derivative of array y + :param dx: spacing of x + :param y: array + :param n: order of derivative + :return: nth order derivative + """ + for i in range(n): + z = np.gradient(y, dx) + y = z + y = np.where(np.abs(y) == np.inf, 0, y) # fix infinity values + return y + + +class PreprocessParameters(ParameterTree): + sigParamChanged = Signal(object) + + def __init__(self): + super(PreprocessParameters, self).__init__() + + self.parameter = Parameter(name='params', type='group', + children=[{'name': "Preprocess method", + 'values': ['Kohler_EMSC','Rubberband'], + 'value': 'Kohler_EMSC', + 'type': 'list'}, + {'name': "Anchor points", + 'value': '400, 4000', + 'type': 'str'}, + {'name': "Fitting regions", + 'value': '[(650, 750),(1780, 2680),(3680, 4000)]', + 'type': 'str'}, + {'name': "Interp method", + 'value': 'linear', + 'values': ['linear', 'quadratic', 'cubic'], + 'type': 'list'} + ]) + self.setParameters(self.parameter, showTop=False) + self.setIndentation(0) + self.parameter.child('Interp method').hide() + self.parameter.child('Anchor points').hide() + + # change Fonts + self.fontSize = 12 + font = QFont("Helvetica [Cronyx]", self.fontSize) + boldFont = QFont("Helvetica [Cronyx]", self.fontSize, QFont.Bold) + self.header().setFont(font) + for item in self.listAllItems(): + if hasattr(item, 'widget'): + item.setFont(0, boldFont) + item.widget.setFont(font) + item.displayLabel.setFont(font) + item.widget.setMaximumHeight(40) + # init params dict + self.argMap = {"Anchor points": 'anchors', + "Interp method": 'kind', + "Fitting regions": 'w_regions' + } + # set self.processArgs to default value + self.processArgs = {} + for child in self.parameter.childs: + if child.name() == "Anchor points": + self.processArgs['anchors'] = '400, 4000' + elif child.name() == "Interp method": + self.processArgs['kind'] = 'linear' + elif child.name() == "Fitting regions": + self.processArgs['w_regions'] = '[(650, 750), (1780, 2680), (3680, 4000)]' + + # connect signals + self.parameter.child('Preprocess method').sigValueChanged.connect(self.updateMethod) + for name in self.argMap.keys(): + self.parameter.child(name).sigValueChanged.connect(partial(self.updateParam, name)) + + def updateParam(self, name): + """ + get latest parameter values + :param name: parameter name + :return: None + """ + self.processArgs[self.argMap[name]] = self.parameter[name] + self.sigParamChanged.emit(self.processArgs) + + def updateMethod(self): + """ + Toggle parameter menu based on fit method + :return: + """ + if self.parameter["Preprocess method"] == 'Kohler_EMSC': + self.parameter.child('Fitting regions').show() + self.parameter.child('Interp method').hide() + self.parameter.child('Anchor points').hide() + else: + self.parameter.child('Fitting regions').hide() + self.parameter.child('Interp method').show() + self.parameter.child('Anchor points').show() + + +class PreprocessWidget(QSplitter): + def __init__(self, headermodel, selectionmodel): + super(PreprocessWidget, self).__init__() + self.headermodel = headermodel + self.mapselectmodel = selectionmodel + self.selectMapidx = 0 + self.resultDict = {} + self.isBatchProcessOn = False + self.out = None + self.dfDict = None + self.reportList = ['preprocess_method', 'wav_anchor', 'interp_method', 'w_regions'] + self.arrayList = ['kohlerDebased', 'kohlerBaseline', 'rubberDebased', 'deriv2_kohler', 'deriv2_rubber'] + self.mousePosList = [] + + # split between spectrum parameters and viewwindow, vertical split + self.params_and_specview = QSplitter() + self.params_and_specview.setOrientation(Qt.Vertical) + # split between buttons and parameters + self.buttons_and_params = QSplitter() + self.buttons_and_params.setOrientation(Qt.Horizontal) + # split between speclist and report + self.speclist_and_report = QSplitter() + self.speclist_and_report.setOrientation(Qt.Vertical) + + # buttons layout + self.buttons = QWidget() + self.buttonlayout = QGridLayout() + self.buttons.setLayout(self.buttonlayout) + # set up buttons + self.fontSize = 12 + font = QFont("Helvetica [Cronyx]", self.fontSize) + self.loadBtn = QPushButton() + self.loadBtn.setText('Load spectra') + self.loadBtn.setFont(font) + self.removeBtn = QPushButton() + self.removeBtn.setText('Remove spectrum') + self.removeBtn.setFont(font) + self.normBox = QComboBox() + self.normBox.addItems(['Raw spectrum', + 'Kohler EMSC baseline', + 'Rubberband baseline', + 'Kohler EMSC + 2nd derivative', + 'Rubberband + 2nd derivative', + ]) + self.normBox.setFont(font) + self.batchBtn = QPushButton() + self.batchBtn.setText('Batch process') + self.batchBtn.setFont(font) + self.saveResultBox = QComboBox() + self.saveResultBox.addItems(['Save kohler', + 'Save kohler baseline', + 'Save rubberband', + 'Save kohler 2nd derivative', + 'Save rubberband 2nd derivative', + 'Save all', + ]) + self.saveResultBox.setFont(font) + # add all buttons + self.buttonlayout.addWidget(self.loadBtn) + self.buttonlayout.addWidget(self.removeBtn) + self.buttonlayout.addWidget(self.normBox) + self.buttonlayout.addWidget(self.batchBtn) + self.buttonlayout.addWidget(self.saveResultBox) + # define report + self.reportWidget = QWidget() + self.reportWidget.setLayout(QVBoxLayout()) + self.infoBox = QTextEdit() + reportTitle = QLabel('Preprocess results') + reportTitle.setFont(font) + self.reportWidget.layout().addWidget(reportTitle) + self.reportWidget.layout().addWidget(self.infoBox) + # spectrum list view + self.specItemModel = QStandardItemModel() + self.specSelectModel = QItemSelectionModel(self.specItemModel) + self.speclistview = QListView() + self.speclistview.setModel(self.specItemModel) + self.speclistview.setSelectionModel(self.specSelectModel) + # add title to list view + self.specListWidget = QWidget() + self.listLayout = QVBoxLayout() + self.specListWidget.setLayout(self.listLayout) + specListTitle = QLabel('Spectrum List') + specListTitle.setFont(font) + self.listLayout.addWidget(specListTitle) + self.listLayout.addWidget(self.speclistview) + + # spectrum plot + self.rawSpectra = baselinePlotWidget() + self.resultSpectra = baselinePlotWidget() + # ParameterTree + self.parametertree = PreprocessParameters() + self.parameter = self.parametertree.parameter + self.processArgs = self.parametertree.processArgs + self.argMap = self.parametertree.argMap + + # assemble widgets + self.buttons_and_params.addWidget(self.parametertree) + self.buttons_and_params.addWidget(self.buttons) + self.buttons_and_params.setSizes([1000, 100]) + self.params_and_specview.addWidget(self.buttons_and_params) + self.params_and_specview.addWidget(self.rawSpectra) + self.params_and_specview.addWidget(self.resultSpectra) + self.params_and_specview.setSizes([150, 50, 50]) + self.speclist_and_report.addWidget(self.specListWidget) + self.speclist_and_report.addWidget(self.reportWidget) + self.speclist_and_report.setSizes([150, 100]) + self.addWidget(self.params_and_specview) + self.addWidget(self.speclist_and_report) + self.setSizes([1000, 200]) + + # Connect signals + self.loadBtn.clicked.connect(self.loadData) + self.removeBtn.clicked.connect(self.removeSpec) + self.batchBtn.clicked.connect(self.batchProcess) + self.saveResultBox.currentIndexChanged.connect(self.saveResults) + self.specSelectModel.selectionChanged.connect(self.updateSpecPlot) + self.normBox.currentIndexChanged.connect(self.updateSpecPlot) + self.parametertree.sigParamChanged.connect(self.updateSpecPlot) + self.rawSpectra.scene().sigMouseClicked.connect(self.setAnchors) + self.parameter.child('Preprocess method').sigValueChanged.connect(self.updateMethod) + + def setHeader(self, field: str): + self.headers = [self.headermodel.item(i).header for i in range(self.headermodel.rowCount())] + self.field = field + self.wavenumberList = [] + self.rc2indList = [] + self.ind2rcList = [] + self.pathList = [] + self.dataSets = [] + + # get wavenumbers, rc2ind + for header in self.headers: + dataEvent = next(header.events(fields=[field])) + self.wavenumberList.append(dataEvent['wavenumbers']) + self.rc2indList.append(dataEvent['rc_index']) + self.ind2rcList.append(dataEvent['index_rc']) + self.pathList.append(dataEvent['path']) + # get raw spectra + data = None + try: # spectra datasets + data = header.meta_array('spectra') + except IndexError: + msg.logMessage('Header object contained no frames with field ''{field}''.', msg.ERROR) + if data is not None: + self.dataSets.append(data) + + def isMapOpen(self): + if not self.mapselectmodel.selectedIndexes(): # no map is open + return False + else: + self.selectMapidx = self.mapselectmodel.selectedIndexes()[0].row() + return True + + def setAnchors(self, event): + # get current map idx and selected spectrum idx + specidx = self.getCurrentSpecid() + plotChoice = self.normBox.currentIndex() + if (not self.isMapOpen()) or (self.specItemModel.rowCount() == 0) or (specidx is None) or (plotChoice not in [2, 4]): + return + + pos = event.pos() + button = event.button() + parser = Preprocessor(self.wavenumberList[self.selectMapidx], self.dataSets[self.selectMapidx][specidx]) + parser.parse_anchors(self.parameter['Anchor points']) + anchor_low, anchor_high = parser.wav_anchor[0], parser.wav_anchor[-1] + if self.rawSpectra.getViewBox().sceneBoundingRect().contains(pos): + mousePoint = self.rawSpectra.getViewBox().mapToView(pos) + x = mousePoint.x() + if anchor_low < x < anchor_high: + if button == Qt.LeftButton:# left click, add point to mousePosList + self.mousePosList.append(x) + elif (button == Qt.MidButton) and self.mousePosList :# right click, remove last point from mousePosList + self.mousePosList.pop() + # set anchors list + anchors = [anchor_low] + sorted(self.mousePosList) + [anchor_high] + txt = ', '.join([str(int(round(x))) for x in anchors]) + self.parameter.child('Anchor points').setValue(txt) + + def getCurrentSpecid(self): + # get selected spectrum idx + specidx = None # default value + if self.specSelectModel.selectedIndexes(): + selectedSpecRow = self.specSelectModel.selectedIndexes()[0].row() + currentSpecItem = self.specItemModel.item(selectedSpecRow) + specidx = currentSpecItem.idx + return specidx + + def updateMethod(self): + if self.parameter["Preprocess method"] == 'Kohler_EMSC': + self.normBox.setCurrentIndex(1) + else: + self.normBox.setCurrentIndex(2) + + def updateSpecPlot(self): + # get current map idx and selected spectrum idx + specidx = self.getCurrentSpecid() + if not self.isMapOpen(): + return + elif self.specItemModel.rowCount() == 0: + MsgBox('No spectrum is loaded.\nPlease click "Load spectra" to import data.') + return + elif specidx is None: + return + + # get plotchoice + plotChoice = self.normBox.currentIndex() + + # create Preprocessor object + self.out = Preprocessor(self.wavenumberList[self.selectMapidx], self.dataSets[self.selectMapidx][specidx]) + baselineOK = self.out.rubber_band(**self.processArgs) and self.out.kohler(**self.processArgs) + + if not baselineOK: + return + + # make results report + if plotChoice != 0: + self.getReport(self.out, plotChoice) + + # if not batch processing, show plots + if not self.isBatchProcessOn: + # clean up plots + self.rawSpectra.clearAll() + self.resultSpectra.clearAll() + if plotChoice == 0: # plot raw spectrum + self.infoBox.setText('') # clear txt + self.rawSpectra.plotBase(self.out, plotType='raw') + elif plotChoice == 1: # plot raw, kohler + self.rawSpectra.plotBase(self.out, plotType='kohler_base') + self.resultSpectra.plotBase(self.out, plotType='kohler') + elif plotChoice == 2: # plot raw, rubberband + self.rawSpectra.plotBase(self.out, plotType='rubber_base') + self.resultSpectra.plotBase(self.out, plotType='rubberband') + elif plotChoice == 3: # plot raw, kohler 2nd derivative + self.rawSpectra.plotBase(self.out, plotType='kohler_base') + self.resultSpectra.plotBase(self.out, plotType='deriv2_kohler') + elif plotChoice == 4: # plot raw, rubberband 2nd derivative + self.rawSpectra.plotBase(self.out, plotType='rubber_base') + self.resultSpectra.plotBase(self.out, plotType='deriv2_rubberband') + + if plotChoice in [1, 3]: + self.parameter.child('Preprocess method').setValue('Kohler_EMSC', blockSignal=self.updateMethod) + elif plotChoice in [2, 4]: + self.parameter.child('Preprocess method').setValue('Rubberband', blockSignal=self.updateMethod) + + def getReport(self, output, plotChoice): + resultTxt = '' + # get baseline results + reportList = self.reportList.copy() + if plotChoice in [2, 4]: + reportList = self.reportList[:-1] + output.preprocess_method = 'rubberband' + elif plotChoice in [1, 3]: + reportList = [self.reportList[0], self.reportList[-1]] + output.preprocess_method = 'kohler' + + for item in dir(output): + if item in reportList: + if item == 'wav_anchor': + val = getattr(output, item) + printFormat = ('{:.2f}, ' * len(val))[:-1] + resultTxt += item + ': ' + printFormat.format(*val) + '\n' + else: + resultTxt += item + ': ' + str(getattr(output, item)) + '\n' + if (item in self.arrayList) or (item in self.reportList): + self.resultDict[item] = getattr(output, item) + + # send text to report info box + self.infoBox.setText(resultTxt) + + def loadData(self): + # get current map idx + if not self.isMapOpen(): + return + # pass the selected map data to plotwidget + self.rawSpectra.setHeader(self.headers[self.selectMapidx], 'spectra') + currentMapItem = self.headermodel.item(self.selectMapidx) + rc2ind = self.rc2indList[self.selectMapidx] + # get current map name + mapName = currentMapItem.data(0) + # get current selected pixels + pixelCoord = currentMapItem.selectedPixels + # get selected specIds + spectraIds = [] + if currentMapItem.selectedPixels is None: # select all + spectraIds = list(range(len(rc2ind))) + else: + for i in range(len(pixelCoord)): + row_col = tuple(pixelCoord[i]) + spectraIds.append(rc2ind[row_col]) + spectraIds = sorted(spectraIds) + # add specitem model + self.specItemModel.clear() + for idx in spectraIds: + item = QStandardItem(mapName + '# ' + str(idx)) + item.idx = idx + self.specItemModel.appendRow(item) + + def removeSpec(self): + # get current selectedSpecRow + if self.specSelectModel.selectedIndexes(): + selectedSpecRow = self.specSelectModel.selectedIndexes()[0].row() + self.specSelectModel.blockSignals(True) + self.specItemModel.removeRow(selectedSpecRow) + self.specSelectModel.blockSignals(False) + # clean up plots + self.rawSpectra.clearAll() + self.resultSpectra.clearAll() + self.infoBox.setText('') + + def cleanUp(self): + self.specItemModel.clear() + self.rawSpectra.clearAll() + self.resultSpectra.clearAll() + self.infoBox.setText('') + self.mousePosList = [] + self.normBox.setCurrentIndex(0) + + def batchProcess(self): + # get current map idx + if not self.isMapOpen(): + return + elif self.specItemModel.rowCount() == 0: + MsgBox('No spectrum is loaded.\nPlease click "Load spectra" to import data.') + return + # check if baseline fit OK + if self.out is None: + self.out = Preprocessor(self.wavenumberList[self.selectMapidx], self.dataSets[self.selectMapidx][0]) + + # get plotchoice + plotChoice = self.normBox.currentIndex() + if plotChoice != 0: + # calculate rubberband and kohler baseline + baselineOK = self.out.rubber_band(**self.processArgs) and self.out.kohler(**self.processArgs) + else: + MsgBox('Plot type is "Raw spectrum".\nPlease change plot type to "Kohler" or "Rubberband".') + return + if not baselineOK: + return + + # notice to user + userMsg = YesNoDialog(f'Ready to batch process selected spectra.\nDo you want to continue?') + userChoice = userMsg.choice() + if userChoice == QMessageBox.No: # user choose to stop + return + + self.isBatchProcessOn = True + + # init resultSetsDict, paramsDict + self.resultSetsDict = {} + self.paramsDict = {} + self.paramsDict['specID'] = [] + self.paramsDict['row_column'] = [] + ind2rc = self.ind2rcList[self.selectMapidx] + energy = self.out.energy + n_energy = len(energy) + for item in self.arrayList: + self.resultSetsDict[item] = np.empty((0, n_energy)) + for item in self.reportList: + self.paramsDict[item] = [] + # batch process begins + n_spectra = self.specItemModel.rowCount() + for i in range(n_spectra): + msg.showMessage(f'Processing {i + 1}/{n_spectra} spectra') + # select each spec and collect results + self.specSelectModel.select(self.specItemModel.index(i, 0), QItemSelectionModel.ClearAndSelect) + # get spec idx + currentSpecItem = self.specItemModel.item(i) + self.paramsDict['specID'].append(currentSpecItem.idx) + self.paramsDict['row_column'].append(ind2rc[currentSpecItem.idx]) + # append all results into a single array/list + for item in self.arrayList: + self.resultSetsDict[item] = np.append(self.resultSetsDict[item], self.resultDict[item].reshape(1, -1), + axis=0) + for item in self.reportList: + self.paramsDict[item].append(self.resultDict[item]) + + # result collection completed. convert paramsDict to df + self.dfDict = {} + self.dfDict['param'] = pd.DataFrame(self.paramsDict).set_index('specID') + for item in self.arrayList: + # convert resultSetsDict to df + self.dfDict[item] = pd.DataFrame(self.resultSetsDict[item], columns=energy.tolist(), + index=self.paramsDict['specID']) + + # batch process completed + self.isBatchProcessOn = False + msg.showMessage(f'Batch processing is completed! Saving results to csv files.') + # save df to files + self.saveResults() + + def saveResults(self): + if self.dfDict is None: + return + filePath = self.pathList[self.selectMapidx] + energy = self.out.energy + saveDataChoice = self.saveResultBox.currentIndex() + if saveDataChoice != 5: # save a single result + saveDataType = self.arrayList[saveDataChoice] + dirName, csvName, h5Name = self.saveToFiles(energy, self.dfDict, filePath, saveDataType) + if h5Name is None: + MsgBox(f'Processed data was saved as csv file at: \n{dirName + csvName}') + else: + MsgBox( + f'Processed data was saved as: \n\ncsv file at: {dirName + csvName} and \n\nHDF5 file at: {dirName + h5Name}') + else: # save all results + csvList = [] + h5List = [] + for saveDataType in self.arrayList: + dirName, csvName, h5Name = self.saveToFiles(energy, self.dfDict, filePath, saveDataType) + csvList.append(csvName) + h5List.append(h5Name) + + allcsvName = (', ').join(csvList) + if h5Name is None: + MsgBox(f'Processed data was saved as csv files at: \n{dirName + allcsvName}') + else: + allh5Name = (', ').join(h5List) + MsgBox( + f'Processed data was saved as: \n\ncsv files at: {dirName + allcsvName} and \n\nHDF5 files at: {dirName + allh5Name}') + + # save parameter + xlsName = csvName[:-4] + '_param.xlsx' + self.dfDict['param'].to_excel(dirName + xlsName) + + def saveToFiles(self, energy, dfDict, filePath, saveDataType): + + ind2rc = self.ind2rcList[self.selectMapidx] + n_spectra = self.specItemModel.rowCount() + + # get dirname and old filename + dirName = os.path.dirname(filePath) + dirName += '/' + oldFileName = os.path.basename(filePath) + + # save dataFrames to csv file + csvName = oldFileName[:-3] + '_' + saveDataType + '.csv' + dfDict[saveDataType].to_csv(dirName + csvName) + + # if a full map is processed, also save results to a h5 file + h5Name = None + if n_spectra == len(ind2rc): + fullMap = ir_map(filename=filePath) + fullMap.add_image_cube() + fullMap.wavenumbers = energy + fullMap.N_w = len(energy) + fullMap.data = np.zeros((fullMap.data.shape[0], fullMap.N_w)) + fullMap.imageCube = np.zeros((fullMap.imageCube.shape[0], fullMap.imageCube.shape[1], fullMap.N_w)) + for i in self.paramsDict['specID']: + fullMap.data[i, :] = self.resultSetsDict[saveDataType][i, :] + row, col = ind2rc[i] + fullMap.imageCube[row, col, :] = fullMap.data[i, :] = self.resultSetsDict[saveDataType][i, :] + # save data as hdf5 + h5Name = oldFileName[:-3] + '_' + saveDataType + '.h5' + fullMap.write_as_hdf5(dirName + h5Name) + + return dirName, csvName, h5Name diff --git a/xicam/BSISB/widgets/spectramaproiwidget.py b/xicam/BSISB/widgets/spectramaproiwidget.py new file mode 100644 index 0000000..2f427e6 --- /dev/null +++ b/xicam/BSISB/widgets/spectramaproiwidget.py @@ -0,0 +1,298 @@ +import os +from qtpy.QtCore import * +from qtpy.QtGui import * +from qtpy.QtWidgets import * +import pickle +import pyqtgraph as pg +from pyqtgraph.parametertree import ParameterTree, Parameter +import numpy as np +from xicam.core.data import NonDBHeader +from xicam.gui.widgets.imageviewmixins import BetterButtons +from xicam.BSISB.widgets.uiwidget import MsgBox, uiSaveFile, uiGetFile +from xicam.BSISB.widgets.mapviewwidget import MapViewWidget +from xicam.BSISB.widgets.spectraplotwidget import SpectraPlotWidget + +class MapView(QSplitter): + sigRoiPixels = Signal(object) + sigRoiState = Signal(object) + sigAutoMaskState = Signal(object) + sigSelectMaskState = Signal(object) + + def __init__(self, header: NonDBHeader = None, stream: str = 'primary', field: str = 'primary' ): + """ + A widget to display imageCube like dataset with ROI buttons + :param header: Xi-cam datahandler header + :param field: header's field param + """ + super(MapView, self).__init__() + # layout set up + self.setOrientation(Qt.Vertical) + self.imageview = MapViewWidget() + self.spectra = SpectraPlotWidget() + self.spectraSplitter = QSplitter() + self.spectraSplitter.addWidget(self.spectra) + # self.spectraSplitter.insertWidget(1, BetterButtons()) # add a 2D spectrum window + # self.spectra.getViewBox().setXRange(0, 4000) # set xrange + + self.imageview_and_toolbar = QSplitter() + self.imageview_and_toolbar.setOrientation(Qt.Horizontal) + self.toolbar_and_param = QSplitter() + self.toolbar_and_param.setOrientation(Qt.Vertical) + #define tool bar + self.toolBar = QWidget() + self.gridlayout = QGridLayout() + self.toolBar.setLayout(self.gridlayout) + #add tool bar buttons + self.roiBtn = QToolButton() + self.roiBtn.setText('Manual ROI') + self.roiBtn.setCheckable(True) + self.roiMeanBtn = QToolButton() + self.roiMeanBtn.setText('ROI Mean') + self.autoMaskBtn = QToolButton() + self.autoMaskBtn.setText('Auto ROI') + self.autoMaskBtn.setCheckable(True) + self.selectMaskBtn = QToolButton() + self.selectMaskBtn.setText('Mark Select') + self.selectMaskBtn.setCheckable(True) + self.saveRoiBtn = QToolButton() + self.saveRoiBtn.setText('Save ROI') + self.saveRoiBtn.setCheckable(False) + self.loadRoiBtn = QToolButton() + self.loadRoiBtn.setText('Load ROI') + self.loadRoiBtn.setCheckable(False) + self.gridlayout.addWidget(self.roiBtn, 0, 0, 1, 1) + self.gridlayout.addWidget(self.autoMaskBtn, 0, 1, 1, 1) + self.gridlayout.addWidget(self.selectMaskBtn, 1, 0, 1, 1) + self.gridlayout.addWidget(self.roiMeanBtn, 1, 1, 1, 1) + self.gridlayout.addWidget(self.saveRoiBtn, 2, 0, 1, 1) + self.gridlayout.addWidget(self.loadRoiBtn, 2, 1, 1, 1) + + self.parameterTree = ParameterTree() + self.parameter = Parameter(name='Threshhold', type='group', + children=[{'name': 'Amide II', + 'value': 0, + 'type': 'float'}, + {'name': "ROI type", + 'values': ['+', '-'], + 'value': '+', + 'type': 'list'}, + ]) + self.parameter.child('Amide II').setOpts(step=0.1) + self.parameterTree.setParameters(self.parameter, showTop=False) + self.parameterTree.setHeaderLabels(['Params','Value']) + self.parameterTree.setIndentation(0) + + # Assemble widgets + self.toolbar_and_param.addWidget(self.toolBar) + self.toolbar_and_param.addWidget(self.parameterTree) + self.toolbar_and_param.setSizes([1000, 1]) #adjust initial splitter size + self.imageview_and_toolbar.addWidget(self.toolbar_and_param) + self.imageview_and_toolbar.addWidget(self.imageview) + self.imageview_and_toolbar.setSizes([1, 1000])#adjust initial splitter size + self.addWidget(self.imageview_and_toolbar) + self.addWidget(self.spectraSplitter) + self.setSizes([1000, 1000]) # adjust initial splitter size + + # readin header + self.imageview.setHeader(header, field='image') + self.spectra.setHeader(header, field='spectra') + self.header = header + + #setup ROI item + sideLen = 10 + self.roi = pg.PolyLineROI(positions=[[0, 0], [sideLen, 0], [sideLen, sideLen], [0, sideLen]], closed=True) + self.imageview.view.addItem(self.roi) + self.roiInitState = self.roi.getState() + self.roi.hide() + + #constants + self.path = os.path.expanduser('~/') + self.pixSelection = {'ROI': None, 'Mask': None} # init pixel selection dict + + # Connect signals + self.imageview.sigShowSpectra.connect(self.spectra.showSpectra) + self.spectra.sigEnergyChanged.connect(self.imageview.setEnergy) + self.roiBtn.clicked.connect(self.roiBtnClicked) + self.roi.sigRegionChangeFinished.connect(self.roiSelectPixel) + self.roi.sigRegionChangeFinished.connect(self.showSelectMask) + self.sigRoiPixels.connect(self.spectra.getSelectedPixels) + self.roiMeanBtn.clicked.connect(self.spectra.showMeanSpectra) + self.autoMaskBtn.clicked.connect(self.showAutoMask) + self.selectMaskBtn.clicked.connect(self.showSelectMask) + self.saveRoiBtn.clicked.connect(self.saveRoi) + self.loadRoiBtn.clicked.connect(self.loadRoi) + self.parameter.child('Amide II').sigValueChanged.connect(self.showAutoMask) + self.parameter.child('Amide II').sigValueChanged.connect(self.intersectSelection) + self.parameter.child('ROI type').sigValueChanged.connect(self.intersectSelection) + + def roiBtnClicked(self): + self.roiSelectPixel() + if self.roiBtn.isChecked(): + self.imageview.cross.hide() + self.roi.show() + self.sigRoiState.emit((True, self.roi.getState())) + else: + self.roi.hide() + self.roi.setState(self.roiInitState) + self.sigRoiState.emit((False, self.roi.getState())) + + def saveRoi(self): + parameterDict = {name: self.parameter[name] for name in self.parameter.names.keys()} + roiStates = {'roiBtn': self.roiBtn.isChecked(), 'maskBtn': self.autoMaskBtn.isChecked(), + 'roiState': self.roi.getState(), 'parameter': parameterDict} + filePath, fileName, canceled = uiSaveFile('Save ROI state', self.path, "Pickle Files (*.pkl)") + if not canceled: + with open(filePath + fileName, 'wb') as f: + pickle.dump(roiStates, f) + MsgBox(f'ROI state file was saved! \nFile Location: {filePath + fileName}') + + def loadRoi(self): + filePath, fileName, canceled = uiGetFile('Open ROI state file', self.path, "Pickle Files (*.pkl)") + if not canceled: + with open(filePath + fileName, 'rb') as f: + roiStates = pickle.load(f) + self.roiBtn.setChecked(roiStates['roiBtn']) + self.roi.setState(roiStates['roiState']) + if roiStates['roiBtn']: + self.roi.show() + self.autoMaskBtn.setChecked(roiStates['maskBtn']) + self.selectMaskBtn.setChecked(True) + self.showSelectMask(True) + for k, v in roiStates['parameter'].items(): + self.parameter[k] = v + MsgBox(f'ROI states were loaded from: \n{filePath + fileName}') + else: + return + + def roiMove(self, roi): + roiState = roi.getState() + self.roi.setState(roiState) + + def getImgShape(self, imgShape, rc2ind): + self.row, self.col = imgShape[0], imgShape[1] + self.rc2ind = rc2ind + # determine whether spectra data is sparse + if len(rc2ind) == self.row * self.col: + self.isDenseImage = True + else: + self.isDenseImage = False + #set up X,Y grid + x = np.linspace(0, self.col - 1, self.col) + y = np.linspace(self.row - 1, 0, self.row) + self.X, self.Y = np.meshgrid(x, y) + if self.isDenseImage: + self.fullMap = list(zip(self.Y.ravel(), self.X.ravel())) + else: + self.fullMap = list(rc2ind.keys()) + # setup automask item + self.autoMask = np.ones((self.row, self.col)) + self.autoMaskItem = pg.ImageItem(self.autoMask, axisOrder="row-major", autoLevels=True, opacity=0.3) + self.imageview.view.addItem(self.autoMaskItem) + self.autoMaskItem.hide() + # setup selctmask item to mark selected pixels + self.selectMask = np.ones((self.row, self.col)) + self.selectMaskItem = pg.ImageItem(self.selectMask, axisOrder="row-major", autoLevels=True, opacity=0.3, + lut = np.array([[0, 0, 0], [255, 0, 0]])) + self.imageview.view.addItem(self.selectMaskItem) + self.selectMaskItem.hide() + + def roiSelectPixel(self): + if self.roiBtn.isChecked(): + #get x,y positions list + xPos = self.roi.getArrayRegion(self.X, self.imageview.imageItem) + xPos = np.round(xPos[xPos > 0]) + yPos = self.roi.getArrayRegion(self.Y, self.imageview.imageItem) + yPos = np.round(yPos[yPos > 0]) + + # extract x,y coordinate from selected region + selectedPixels = list(zip(yPos, xPos)) + self.intersectSelection('ROI', selectedPixels) + self.sigRoiState.emit((True, self.roi.getState())) + else: + self.intersectSelection('ROI', None) # no ROI, select all pixels + self.sigRoiState.emit((False, self.roi.getState())) + + def showSelectMask(self, signalReceived): + if self.selectMaskBtn.isChecked(): + # show roi and autoMask + if self.roiBtn.isChecked(): + self.roi.show() + self.sigRoiState.emit((True, self.roi.getState())) + if self.autoMaskBtn.isChecked(): + self.autoMaskItem.show() + self.sigAutoMaskState.emit((True, self.autoMask)) + # update and show mask + self.selectMaskItem.setImage(self.selectMask) + self.selectMaskItem.show() + self.sigSelectMaskState.emit((True, self.selectMask)) + else: + self.selectMaskItem.hide() + self.sigSelectMaskState.emit((False, self.selectMask)) + if signalReceived == False: + self.roi.hide() + self.autoMaskItem.hide() + self.sigRoiState.emit((False, self.roi.getState())) + self.sigAutoMaskState.emit((False, self.autoMask)) + + + def showAutoMask(self): + if self.autoMaskBtn.isChecked(): + # update and show mask + self.autoMask = self.imageview.makeMask([self.parameter['Amide II']]) + self.autoMaskItem.setImage(self.autoMask) + self.autoMaskItem.show() + # select pixels + mask = self.autoMask.astype(np.bool) + selectedPixels = list(zip(self.Y[mask], self.X[mask])) + self.intersectSelection('Mask', selectedPixels) + self.sigAutoMaskState.emit((True, self.autoMask)) + else: + self.autoMaskItem.hide() + self.autoMask[:, :] = 1 + self.intersectSelection('Mask', None) # no mask, select all pixels + self.sigAutoMaskState.emit((False, self.autoMask)) + + def intersectSelection(self, selector, selectedPixels): + # update pixel selection dict + if (selector == 'ROI') or (selector == 'Mask'): + self.pixSelection[selector] = selectedPixels + # reverse ROI selection + if (self.parameter['ROI type'] == '-') and (self.pixSelection['ROI'] is not None): + roi_copy = self.pixSelection['ROI'] + reverseROI = set(self.fullMap) - set(self.pixSelection['ROI']) + self.pixSelection['ROI'] = list(reverseROI) + + if (self.pixSelection['ROI'] is None) and (self.pixSelection['Mask'] is None): + if self.isDenseImage: + self.sigRoiPixels.emit(None) # no ROI, select all pixels + self.selectMask = np.ones((self.row, self.col)) + else: + allSelected = np.array(list(self.rc2ind.keys()), dtype='int') + self.sigRoiPixels.emit(allSelected) # no ROI, select all sparse rc2ind + self.selectMask = np.zeros((self.row, self.col)) + self.selectMask[allSelected[:, 0], allSelected[:, 1]] = 1 + self.selectMask = np.flipud(self.selectMask) + return + elif self.pixSelection['ROI'] is None: + allSelected = set(self.pixSelection['Mask']) #de-duplication of pixels + elif self.pixSelection['Mask'] is None: + allSelected = set(self.pixSelection['ROI']) #de-duplication of pixels + else: + allSelected = set(self.pixSelection['ROI']) & set(self.pixSelection['Mask']) + + if self.isDenseImage: + allSelected = np.array(list(allSelected), dtype='int') # convert to array + else: + allSelected &= set(self.rc2ind.keys()) + allSelected = np.array(list(allSelected), dtype='int') + + self.selectMask = np.zeros((self.row, self.col)) + if len(allSelected) > 0: + self.selectMask[allSelected[:, 0], allSelected[:, 1]] = 1 + self.selectMask = np.flipud(self.selectMask) + self.sigRoiPixels.emit(allSelected) + # show SelectMask + self.showSelectMask(selector) + #recover ROI selection + if (self.parameter['ROI type'] == '-') and (self.pixSelection['ROI'] is not None): + self.pixSelection['ROI'] = roi_copy \ No newline at end of file diff --git a/xicam/BSISB/widgets/spectraplotwidget.py b/xicam/BSISB/widgets/spectraplotwidget.py index 9a27aa1..152aee4 100644 --- a/xicam/BSISB/widgets/spectraplotwidget.py +++ b/xicam/BSISB/widgets/spectraplotwidget.py @@ -1,28 +1,44 @@ -from pyqtgraph import PlotWidget, TextItem, PlotDataItem +from qtpy.QtCore import Qt +from pyqtgraph import PlotWidget, TextItem, PlotDataItem, mkPen from xicam.core import msg from xicam.core.data import NonDBHeader import numpy as np from pyqtgraph import InfiniteLine from qtpy.QtCore import Signal from lbl_ir.data_objects.ir_map import val2ind +from xicam.BSISB.widgets.mapviewwidget import toHtml class SpectraPlotWidget(PlotWidget): sigEnergyChanged = Signal(object) - def __init__(self, *args, **kwargs): + def __init__(self, linePos=650, txtPosRatio=0.35, invertX=True, *args, **kwargs): + """ + A widget to display a 1D spectrum + :param linePos: the initial position of the InfiniteLine + :param txtPosRatio: a coefficient that determines the relative position of the textItem + :param invertX: whether to invert X-axis + """ super(SpectraPlotWidget, self).__init__(*args, **kwargs) self._data = None + assert (txtPosRatio >= 0) and (txtPosRatio <= 1), 'Please set txtPosRatio value between 0 and 1.' + self.txtPosRatio = txtPosRatio self.positionmap = dict() self.wavenumbers = None - self._meanSpec = True # whether current spectrum is a mean spectrum + self._meanSpec = True # whether current spectrum is a mean spectrum self.line = InfiniteLine(movable=True) self.line.setPen((255, 255, 0, 200)) - self.line.setZValue(100) + self.line.setValue(linePos) self.line.sigPositionChanged.connect(self.sigEnergyChanged) self.line.sigPositionChanged.connect(self.getEnergy) self.addItem(self.line) - self.getViewBox().invertX(True) + self.cross = PlotDataItem([linePos], [0], symbolBrush=(255, 0, 0), symbolPen=(255, 0, 0), symbol='+', + symbolSize=20) + self.cross.setZValue(100) + self.addItem(self.cross) + self.txt = TextItem() + self.getViewBox().invertX(invertX) + self.spectrumInd = 0 self.selectedPixels = None self._y = None @@ -33,15 +49,13 @@ def getEnergy(self, lineobject): x_val = self.wavenumbers[idx] y_val = self._y[idx] if not self._meanSpec: - txt_html = f'
\ - Spectrum #{self.spectrumInd}
' + txt_html = toHtml(f'Spectrum #{self.spectrumInd}') else: - txt_html = f'
\ - {self._mean_title}
' + txt_html = toHtml(f'{self._mean_title}') - txt_html += f'
\ - X = {x_val: .2f}, Y = {y_val: .4f}
' + txt_html += toHtml(f'X = {x_val: .2f}, Y = {y_val: .4f}') self.txt.setHtml(txt_html) + self.cross.setData([x_val], [y_val]) def setHeader(self, header: NonDBHeader, field: str, *args, **kwargs): self.header = header @@ -56,15 +70,15 @@ def setHeader(self, header: NonDBHeader, field: str, *args, **kwargs): try: data = header.meta_array(field) except IndexError: - msg.logMessage('Header object contained no frames with field ''{field}''.', msg.ERROR) + msg.logMessage(f'Header object contained no frames with field {field}.', msg.ERROR) if data is not None: # kwargs['transform'] = QTransform(1, 0, 0, -1, 0, data.shape[-2]) self._data = data def showSpectra(self, i=0): - if self._data is not None: - self.clear() + if (self._data is not None) and (i < len(self._data)): + self.getViewBox().clear() self._meanSpec = False self.spectrumInd = i self.plot(self.wavenumbers, self._data[i]) @@ -73,9 +87,16 @@ def getSelectedPixels(self, selectedPixels): self.selectedPixels = selectedPixels # print(selectedPixels) + def clearAll(self): + # remove legend + _legend = self.plotItem.legend + if (_legend is not None) and (_legend.scene() is not None): + _legend.scene().removeItem(_legend) + self.getViewBox().clear() + def showMeanSpectra(self): self._meanSpec = True - self.clear() + self.getViewBox().clear() if self.selectedPixels is not None: n_spectra = len(self.selectedPixels) tmp = np.zeros((n_spectra, self.N_w)) @@ -99,6 +120,44 @@ def plot(self, x, y, *args, **kwargs): # set up infinity line and get its position self.plotItem.plot(x, y, *args, **kwargs) self.addItem(self.line) + self.addItem(self.cross) + x_val = self.line.value() + if x_val == 0: + y_val = 0 + else: + idx = val2ind(x_val, self.wavenumbers) + x_val = self.wavenumbers[idx] + y_val = y[idx] + + if not self._meanSpec: + txt_html = toHtml(f'Spectrum #{self.spectrumInd}') + else: + txt_html = toHtml(f'{self._mean_title}') + + txt_html += toHtml(f'X = {x_val: .2f}, Y = {y_val: .4f}') + self.txt.setHtml(txt_html) + ymax = np.max(y) + self._y = y + r = self.txtPosRatio + self.txt.setPos(r * x[-1] + (1 - r) * x[0], ymax) + self.cross.setData([x_val], [y_val]) + self.addItem(self.txt) + +class baselinePlotWidget(SpectraPlotWidget): + def __init__(self): + super(baselinePlotWidget, self).__init__() + self.line.setValue(800) + self.txt = TextItem('', anchor=(0, 0)) + self.cross = PlotDataItem([800], [0], symbolBrush=(255, 255, 0), symbolPen=(255, 255, 0), + symbol='+',symbolSize=20) + self.line.sigPositionChanged.connect(self.getMu) + self._mu = None + + def plot(self, x, y, *args, **kwargs): + # set up infinity line and get its position + self.plotItem.plot(x, y, *args, **kwargs) + self.addItem(self.line) + self.addItem(self.cross) x_val = self.line.value() if x_val == 0: y_val = 0 @@ -117,7 +176,94 @@ def plot(self, x, y, *args, **kwargs): txt_html += f'
\ X = {x_val: .2f}, Y = {y_val: .4f}
' self.txt = TextItem(html=txt_html, anchor=(0, 0)) - ymax = max(y) + ymax = np.max(y) self._y = y - self.txt.setPos(1500, 0.95 * ymax) + r = self.txtPosRatio + self.txt.setPos(r * x[-1] + (1 - r) * x[0], 0.95 * ymax) + self.cross.setData([x_val], [y_val]) + self.addItem(self.txt) + + def getMu(self): + if self._mu is not None: + x_val = self.line.value() + if x_val == 0: + y_val = 0 + else: + idx = val2ind(x_val, self._x) + x_val = self._x[idx] + y_val = self._mu[idx] + txt_html = f'
\ + X = {x_val: .2f}, Y = {y_val: .4f}
' + self.txt.setHtml(txt_html) + self.cross.setData([x_val], [y_val]) + + def addDataCursor(self, x, y): + self.addItem(self.line) + self.addItem(self.cross) + ymax = np.max(y) + self.txt.setText('') + r = self.txtPosRatio + self.txt.setPos(r * x[-1] + (1 - r) * x[0], 0.95 * ymax) self.addItem(self.txt) + self.getMu() + + def plotBase(self, dataGroup, plotType='raw'): + """ + make plots for Larch Group object + :param dataGroup: Larch Group object + :return: + """ + # add legend + self.plotItem.addLegend(offset=(-1, -1)) + x = self._x = dataGroup.energy # self._x, self._mu for getEnergy + y = self._mu = dataGroup.specTrim + n = len(x) # array length + self._y = None # disable getEnergy func + if plotType == 'raw': + self.plotItem.plot(x, y, name='Raw', pen=mkPen('w', width=2)) + elif plotType == 'rubber_base': + self.plotItem.plot(x, y, name='Raw', pen=mkPen('w', width=2)) + self.plotItem.plot(x, dataGroup.rubberBaseline, name='Rubberband baseline', pen=mkPen('g', style=Qt.DotLine, width=2)) + self.plotItem.plot(dataGroup.wav_anchor, dataGroup.spec_anchor, symbol='o', symbolPen='r', symbolBrush=0.5) + elif plotType == 'kohler_base': + self.plotItem.plot(x, y, name='Raw', pen=mkPen('w', width=2)) + self.plotItem.plot(x, dataGroup.kohlerBaseline, name='Kohler EMSC baseline', pen=mkPen('g', style=Qt.DotLine, width=2)) + elif plotType == 'rubberband': + y = self._mu = dataGroup.rubberDebased + self.plotItem.plot(x, y, name='Rubberband debased', pen=mkPen('r', width=2)) + elif plotType == 'kohler': + y = self._mu = dataGroup.kohlerDebased + self.plotItem.plot(x, y, name='Kohler EMSC debased', pen=mkPen('r', width=2)) + elif plotType == 'deriv2_rubberband': + y = self._mu = dataGroup.rubberDebased # for data cursor + scale, offset = self.alignTwoCurve(dataGroup.rubberDebased[n//4:n*3//4], dataGroup.deriv2_rubber[n//4:n*3//4]) + deriv2Scaled = dataGroup.deriv2_rubber * scale + offset + ymin, ymax = np.min(y), np.max(y) + self.getViewBox().setYRange(ymin, ymax, padding=0.1) + self.plotItem.plot(x, y, name='Rubberband debased', pen=mkPen('r', width=2)) + self.plotItem.plot(x, deriv2Scaled, name='2nd derivative (scaled, Rubberband)', pen=mkPen('g', width=2)) + elif plotType == 'deriv2_kohler': + y = self._mu = dataGroup.kohlerDebased + scale, offset = self.alignTwoCurve(dataGroup.kohlerDebased[n//4:n*3//4], dataGroup.deriv2_kohler[n//4:n*3//4]) + deriv2Scaled = dataGroup.deriv2_kohler * scale + offset + ymin, ymax = np.min(y), np.max(y) + self.getViewBox().setYRange(ymin, ymax, padding=0.1) + self.plotItem.plot(x, y, name='Rubberband debased', pen=mkPen('r', width=2)) + self.plotItem.plot(x, deriv2Scaled, name='2nd derivative (scaled, Kohler)', pen=mkPen('g', width=2)) + # add infinityline, cross + self.addDataCursor(x, y) + + def alignTwoCurve(self, y1, y2): + """ + Align the scale of y2 to that of y1 + :param y1: the main curve + :param y2: the curve to be aligned + :return: + scale: scale factor + offset: y offset + """ + y1Range, y2Range = np.max(y1) - np.min(y1), np.max(y2) - np.min(y2) + scale = y1Range / y2Range + y = y2 * scale + offset = np.max(y1) - np.max(y) + return scale, offset \ No newline at end of file