Visualise 2D maps from the Horiba Evolution confocal Raman microscope/spectrometer

Horiba_Labspec_Raman_map_analysis.py

#!/usr/bin/python3  
#-*- coding: utf-8 -*-
# Visualise 2D maps from the Horiba Evolution confocal Raman microscope/spectrometer (needs converting .L6M to .TXT)

## Import common moduli
import sys, os, time, collections
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
#from scipy.constants import c, hbar, pi

## Use LaTeX (optional)
# matplotlib.rc('text', usetex=True)
# matplotlib.rc('font', size=8)
# matplotlib.rc('text.latex', preamble = r'\usepackage{amsmath}, \usepackage{yfonts}, \usepackage{txfonts}, \usepackage{lmodern},')


## Load data
wl = np.genfromtxt(sys.argv[1], unpack=True, max_rows=1)
raw = np.genfromtxt(sys.argv[1], unpack=True, skip_header=1)
x, y, rawI = np.unique(raw[0]), np.unique(raw[1]), raw[2:]
print(x, y, rawI.shape)



## Raman spectra preprocessing

def retouch_outliers(y, sigma_criterion=2, peaks=True, notches=False, iterations=3):
    kernel = [.5,0,.5] # averaging of neighbors #kernel = np.exp(-np.linspace(-2,2,5)**2) ## Gaussian
    for n in range(iterations):
        conv = np.convolve(y, kernel, mode='same')
        norm = np.convolve(np.ones_like(y), kernel, mode='same')
        smooth = conv/norm    # find the average value of neighbors
        rms_noise = np.average((y[1:]-y[:-1])**2)**.5   # estimate what the average noise is (rms derivative)
        if peaks and notches:
            outlier_mask = (np.abs(y-smooth) > rms_noise*sigma_criterion)    # find all points with difference from average less than 3sigma
        elif peaks:
            outlier_mask = ((y-smooth) > rms_noise*sigma_criterion)    # find all points with difference from average less than 3sigma
        elif notches:
            outlier_mask = ((y-smooth) < rms_noise*sigma_criterion)    # find all points with difference from average less than 3sigma
        #y[outlier_mask] = np.roll(y,1)[outlier_mask] # smooth[outlier_mask+1]
        y[outlier_mask] = smooth[outlier_mask]
    #print(np.sum(outlier_mask), end=' ',flush=True)
    #print(rms_noise, end=', ',flush=True)
    return y

rawIr = np.apply_along_axis(retouch_outliers, 0, rawI.copy())

def smooth(y, width=10):
    kernel = 2**(-np.linspace(-2, 2, width)**2) # truncated Gaussian
    conv = np.convolve(y, kernel, mode='same')
    norm = np.convolve(np.ones_like(y), kernel, mode='same')
    return conv/norm

rawIr = np.apply_along_axis(smooth, 0, rawIr)

def rm_bg(y, iter=50, coef=0.75, blurpx=250):
    """ subtracts smooth slowly varying background, keeping peaks and similar features,
    (almost) never resulting in negative values """
    def edge_safe_convolve(arr,ker): 
        return np.convolve(np.pad(arr,len(ker),mode='edge'),ker,mode='same')[len(ker):-len(ker)]

    y0 = y[:]
    ker = 2**-np.linspace(-2,2,blurpx)**2; 
    for i in range(iter):
      y = edge_safe_convolve(y,ker/np.sum(ker))
      y = y - ( np.abs(y-y0) + y - y0)*coef
    return y0-y

rawI = np.apply_along_axis(rm_bg, 0, rawI)
rawIr = np.apply_along_axis(rm_bg, 0, rawIr)


# == Plotting with object model==
fig, (ax, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(20, 10))
I = rawI.reshape([len(wl), len(x), len(y)])
Ir = rawIr.reshape([len(wl), len(x), len(y)])

for yy in range(I.shape[2]):
    for xx in range(I.shape[1]):
        #ax.plot(wl,I[:,xx,yy], alpha=.2, c='r'); 
        #ax.plot(wl,I[:,xx,yy], alpha=1, c='r' ,lw=.1, marker='.', markersize=1); 
        ax.plot(wl,Ir[:,xx,yy], alpha=1, c='k',lw=.1) #, marker='.', markersize=1); 
ax.set_xlabel(u"Raman shift (cm⁻¹)"); 
ax.set_ylabel(u"Intensity (a.u.)"); 
ax.grid()
ax.legend(prop={'size':10}, loc='upper right')

ax2.imshow(I.sum(axis=0))
ax2.set_xlabel(u"Position x (μm)"); 
ax2.set_ylabel(u"Position y (μm)"); 

## ==== Outputting ====
## Simple axes
#ax.set_ylim((-0.1,1.1)); ax.set_yscale('linear')
#ax.set_xlim((-0.1,1.1)); ax.set_xscale('linear')

## Show or save
fig.savefig("output.png", bbox_inches='tight')

#np.savetxt(sys.argv[1]+"_corrected.dat", np.array([]).T)

fig.tight_layout()
fig.canvas.mpl_connect('close_event', quit);  
fig.show()