diff --git a/check.py b/check.py index f9bc884..fff81c5 100644 --- a/check.py +++ b/check.py @@ -7,6 +7,7 @@ import matplotlib.pyplot as plt from optifik.analysis import * +from optifik.fft import Prominence_from_fft from optifik.auto import auto from optifik import io @@ -35,20 +36,31 @@ def play(): FILE_NAME = '000004310.xy' #TEST# spectrum_file = os.path.join(DATA_FOLDER, FILE_NAME) + + + spectrum_file = 'tests/spectraVictor1/T5469.xy' #TEST# + lambdas, intensities = io.load_spectrum(spectrum_file) plot_spectrum(lambdas, intensities, title='Raw') - + lambdas, intensities = io.load_spectrum(spectrum_file, lambda_min=450) plot_spectrum(lambdas, intensities, title='Raw, cropped') - + smoothed_intensities = smooth_intensities(intensities) plot_spectrum(lambdas, smoothed_intensities, title='Smoothed') - + prominence = 0.02 peaks_min, peaks_max = finds_peak(lambdas, smoothed_intensities, min_peak_prominence=prominence, plot=True) + indice = 1.324188 + 3102.060378 / (lambdas**2) + p, s, w = Prominence_from_fft(lambdas, smoothed_intensities, indice, plot=True) + + print(p) + peaks_min, peaks_max = finds_peak(lambdas, smoothed_intensities, + min_peak_prominence=p, + plot=True) def check_basic(): @@ -71,28 +83,28 @@ def check_basic(): def check_SV1(): - - + + DATA_FOLDER = os.path.join('tests', 'spectraVictor1') - + import yaml yaml_file = os.path.join(DATA_FOLDER, 'known_value.yaml') with open(yaml_file, "r") as yaml_file: thickness_dict = yaml.safe_load(yaml_file) - - + + for fn, val in thickness_dict.items(): #auto(DATA_FOLDER, fn) - + spectre_file = os.path.join(DATA_FOLDER, fn) lambdas, raw_intensities = load_spectrum(spectre_file, lambda_min=450) - + ##### Affichage du spectre lissé ##### - + #smoothed_intensities, intensities, lambdas = Data_Smoothed(spectre_file) - + smoothed_intensities = smooth_intensities(raw_intensities) @@ -109,22 +121,22 @@ def check_SV1(): peaks_min, peaks_max = finds_peak(lambdas, smoothed_intensities, min_peak_prominence=prominence, plot=False) - + result = thickness_from_minmax(lambdas, smoothed_intensities, refractive_index=indice, min_peak_prominence=prominence) print(f'thickness: {result.thickness:.2f} nm') - - + + print(f'expected: {val}') print('#-' * 10) - + if __name__ == '__main__': - check_basic() - check_SV1() - play() \ No newline at end of file + #check_basic() + #check_SV1() + play() diff --git a/optifik/analysis.py b/optifik/analysis.py index 3956f7c..1512662 100644 --- a/optifik/analysis.py +++ b/optifik/analysis.py @@ -1,6 +1,3 @@ - -import pandas as pd - from scipy.signal import savgol_filter from scipy.signal import find_peaks @@ -13,54 +10,56 @@ plt.rcParams.update({ 'legend.fontsize': 23, }) -from .io import load_spectrum -from .fft import * -from .minmax import * +def plot_spectrum(wavelengths, intensities, title=''): - - -def plot_spectrum(lambdas, intensities, title=''): - - plt.figure(figsize=(10, 6),dpi =600) - plt.plot(lambdas, intensities, 'o-', markersize=2) + plt.figure(figsize=(10, 6), dpi=300) + plt.plot(wavelengths, intensities, 'o-', markersize=2) plt.xlabel(r'$\lambda$ (nm)') plt.ylabel(r'$I^*$') plt.title(title) - plt.tight_layout() + plt.tight_layout() plt.show() - -def finds_peak(lambdas, intensities, min_peak_prominence, min_peak_distance=10, plot=None): +def finds_peak(wavelengths, intensities, min_peak_prominence, min_peak_distance=10, plot=None): """ - Charge un fichier .xy et affiche les données avec les extrema détectés (minima et maxima). + Detect minima and maxima Parameters ---------- - filename : str - Chemin vers le fichier .xy (2 colonnes : lambda et intensité). + wavelengths : array + Wavelength values in nm. + intensities : array + Intensity values. min_peak_prominence : float - Importance minimale des pics. - min_peak_distance : float - Distance minimale entre les pics. + min prominence for scipy find_peak. + min_peak_distance : int, optional + min peak distance for scipy find_peak. The default is 10. + plot : bool, optional + Display a curve, useful for checking or debuging. The default is None. + + Returns + ------- + (peaks_min, peaks_max) + """ - + peaks_max, _ = find_peaks(intensities, prominence=min_peak_prominence, distance=min_peak_distance) peaks_min, _ = find_peaks(-intensities, prominence=min_peak_prominence, distance=min_peak_distance) - + if plot: - plt.figure(figsize=(10, 6),dpi =600) - plt.plot(lambdas, intensities, 'o-', markersize=2, label="Smoothed data") - plt.plot(lambdas[peaks_max], intensities[peaks_max], 'ro') - plt.plot(lambdas[peaks_min], intensities[peaks_min], 'ro') + plt.figure(figsize=(10, 6), dpi=300) + plt.plot(wavelengths, intensities, 'o-', markersize=2, label="Smoothed data") + plt.plot(wavelengths[peaks_max], intensities[peaks_max], 'ro') + plt.plot(wavelengths[peaks_min], intensities[peaks_min], 'ro') plt.xlabel(r'$\lambda$ (nm)') plt.ylabel(r'$I^*$') plt.legend() - plt.tight_layout() + plt.tight_layout() plt.show() @@ -68,8 +67,23 @@ def finds_peak(lambdas, intensities, min_peak_prominence, min_peak_distance=10, -def smooth_intensities(intensities): - WIN_SIZE = 11 - smoothed_intensities = savgol_filter(intensities, WIN_SIZE, 3) +def smooth_intensities(intensities, window_size=11): + """ + Return a smoothed intensities array with a Savitzky-Golay filter. + + Parameters + ---------- + intensities : ndarray + Intensity values + window_size : int, optional + The length of the filter window. The default is 11. + + Returns + ------- + smoothed_intensities + + """ + polynom_order = 3 + smoothed_intensities = savgol_filter(intensities, window_size, 3) return smoothed_intensities diff --git a/optifik/auto.py b/optifik/auto.py index 436e052..c56249c 100644 --- a/optifik/auto.py +++ b/optifik/auto.py @@ -24,10 +24,10 @@ def auto(spectrum_file, plot=None): ##### Determination de la prominence associé ##### - prominence = Prominence_from_fft(lambdas=lambdas, - intensities=smoothed_intensities, - refractive_index=indice, - plot=plot) + prominence, signal, wavelength = Prominence_from_fft(lambdas, + smoothed_intensities, + refractive_index=indice, + plot=plot) prominence = 0.03 ##### Find Peak ##### @@ -48,7 +48,7 @@ def auto(spectrum_file, plot=None): result = thickness_from_fft(lambdas, smoothed_intensities, refractive_index=indice, plot=plot) - + print(f'thickness: {result.thickness:.2f} nm') @@ -58,7 +58,7 @@ def auto(spectrum_file, plot=None): refractive_index=indice, min_peak_prominence=prominence, plot=plot) - + print(f'thickness: {result.thickness:.2f} nm') if total_extrema <= 4 and total_extrema >= 2: #& 2peak minimum: @@ -67,7 +67,7 @@ def auto(spectrum_file, plot=None): refractive_index=indice, min_peak_prominence=prominence, plot=plot) - + print(f'thickness: {result.thickness:.2f} nm') if total_extrema <= 4 and len(peaks_max) == 1 and len(peaks_min) == 0 : #dans l'ordre zéro ! diff --git a/optifik/fft.py b/optifik/fft.py index d0ccf9f..45a2229 100644 --- a/optifik/fft.py +++ b/optifik/fft.py @@ -17,7 +17,7 @@ from .utils import OptimizeResult -def thickness_from_fft(lambdas, intensities, +def thickness_from_fft(wavelengths, intensities, refractive_index, num_half_space=None, plot=None): @@ -26,7 +26,7 @@ def thickness_from_fft(lambdas, intensities, Parameters ---------- - lambdas : array + wavelengths : array Wavelength values in nm. intensities : array Intensity values. @@ -34,7 +34,7 @@ def thickness_from_fft(lambdas, intensities, Value of the refractive index of the medium. num_half_space : scalar, optional Number of points to compute FFT's half space. - If `None`, default corresponds to `10*len(lambdas)`. + If `None`, default corresponds to `10*len(wavelengths)`. debug : boolean, optional Show plot of the transformed signal and the peak detection. @@ -44,32 +44,32 @@ def thickness_from_fft(lambdas, intensities, The attribute `thickness` gives the thickness value in nm. """ if num_half_space is None: - num_half_space = 10 * len(lambdas) + num_half_space = 10 * len(wavelengths) # FFT requires evenly spaced data. # So, we interpolate the signal. - # Interpolate to get a linear increase of 1 / lambdas. - inverse_lambdas_times_n = refractive_index / lambdas - f = interp1d(inverse_lambdas_times_n, intensities) + # Interpolate to get a linear increase of 1 / wavelengths. + inverse_wavelengths_times_n = refractive_index / wavelengths + f = interp1d(inverse_wavelengths_times_n, intensities) - inverse_lambdas_linspace = np.linspace(inverse_lambdas_times_n[0], - inverse_lambdas_times_n[-1], + inverse_wavelengths_linspace = np.linspace(inverse_wavelengths_times_n[0], + inverse_wavelengths_times_n[-1], 2*num_half_space) - intensities_linspace = f(inverse_lambdas_linspace) + intensities_linspace = f(inverse_wavelengths_linspace) # Perform FFT - density = (inverse_lambdas_times_n[-1]-inverse_lambdas_times_n[0]) / (2*num_half_space) - inverse_lambdas_fft = fftfreq(2*num_half_space, density) + density = (inverse_wavelengths_times_n[-1]-inverse_wavelengths_times_n[0]) / (2*num_half_space) + inverse_wavelengths_fft = fftfreq(2*num_half_space, density) intensities_fft = fft(intensities_linspace) # The FFT is symetrical over [0:N] and [N:2N]. # Keep over [N:2N], ie for positive freq. intensities_fft = intensities_fft[num_half_space:2*num_half_space] - inverse_lambdas_fft = inverse_lambdas_fft[num_half_space:2*num_half_space] + inverse_wavelengths_fft = inverse_wavelengths_fft[num_half_space:2*num_half_space] idx_max_fft = np.argmax(abs(intensities_fft)) - freq_max = inverse_lambdas_fft[idx_max_fft] + freq_max = inverse_wavelengths_fft[idx_max_fft] thickness_fft = freq_max / 2. @@ -77,7 +77,7 @@ def thickness_from_fft(lambdas, intensities, plt.figure(figsize=(10, 6),dpi =600) if plot: - plt.loglog(inverse_lambdas_fft, np.abs(intensities_fft)) + plt.loglog(inverse_wavelengths_fft, np.abs(intensities_fft)) plt.loglog(freq_max, np.abs(intensities_fft[idx_max_fft]), 'o') plt.xlabel('Frequency') plt.ylabel(r'FFT($I^*$)') @@ -86,67 +86,90 @@ def thickness_from_fft(lambdas, intensities, return OptimizeResult(thickness=thickness_fft,) -def Prominence_from_fft(lambdas, intensities, refractive_index, num_half_space=None, plot=True): +def Prominence_from_fft(wavelengths, intensities, refractive_index, + num_half_space=None, plot=None): if num_half_space is None: - num_half_space = 10 * len(lambdas) + num_half_space = len(wavelengths) - # Interpolation pour que les données soient uniformément espacées - inverse_lambdas_times_n = refractive_index / lambdas - f = interp1d(inverse_lambdas_times_n, intensities) + # # # 1. Spectre original + # if plot: + # plt.figure(figsize=(10, 6), dpi=150) + # plt.plot(1/wavelengths, intensities, label='Spectre original') + # plt.xlabel('1/Longueur d\'onde (nm-1)') + # plt.ylabel('Intensité') + # plt.legend() + # plt.show() - inverse_lambdas_linspace = np.linspace(inverse_lambdas_times_n[0], - inverse_lambdas_times_n[-1], - 2*num_half_space) - intensities_linspace = f(inverse_lambdas_linspace) - # FFT - density = (inverse_lambdas_times_n[-1] - inverse_lambdas_times_n[0]) / (2*num_half_space) - freqs = fftfreq(2*num_half_space, density) - fft_vals = fft(intensities_linspace) + # 2. Conversion lambda → k = n(λ) / λ + k_vals = refractive_index / wavelengths + f_interp = interp1d(k_vals, intensities, kind='linear', fill_value="extrapolate") - # On conserve uniquement les fréquences positives - freqs = freqs[num_half_space:] - fft_vals = fft_vals[num_half_space:] + # 3. Axe k uniforme + interpolation + k_linspace = np.linspace(k_vals[0], k_vals[-1], 2 * num_half_space) + intensities_k = f_interp(k_linspace) - # Trouver le pic principal - abs_fft = np.abs(fft_vals) - idx_max = np.argmax(abs_fft) - F_max = freqs[idx_max] + # 4. FFT + delta_k = (k_vals[-1] - k_vals[0]) / (2 * num_half_space) + fft_vals = fft(intensities_k) + freqs = fftfreq(2 * num_half_space, delta_k) + + # 5. Pic FFT + freqs_pos = freqs[freqs > 0] + abs_fft_pos = np.abs(fft_vals[freqs > 0]) + idx_max = np.argmax(abs_fft_pos) + F_max = freqs_pos[idx_max] if plot: - print(f"F_max detected at: {F_max:.4f}") - plt.figure(figsize=(10, 6),dpi = 600) - plt.plot(freqs, abs_fft, label='|FFT|') - plt.axvline(F_max, color='r', linestyle='--', label='F_max') - plt.xlabel('Fréquence') - plt.ylabel('Amplitude FFT') - plt.yscale('log') + plt.figure(figsize=(10, 6), dpi=150) + plt.plot(freqs_pos, abs_fft_pos, label='|FFT|') + plt.axvline(F_max, color='r', linestyle='--', label='Pic principal') + plt.xlabel('Distance optique [nm]') + plt.ylabel(r'FFT($I^*$)') plt.xscale('log') + plt.yscale('log') plt.legend() plt.show() - # Filtrage : on garde les composantes au-dessus de 10 * F_max - cutoff = 10 * F_max - mask = freqs >= cutoff - fft_filtered = np.zeros_like(fft_vals) - fft_filtered[mask] = fft_vals[mask] + # 6. Filtrage (garde hautes fréquences) + cutoff_HF = 2 * F_max - fft_full = np.zeros(2 * num_half_space, dtype=complex) - fft_full[num_half_space:] = fft_filtered # fréquences positives - fft_full[:num_half_space] = np.conj(fft_filtered[::-1]) - # IFFT - signal_filtered = np.real(ifft(fft_full)) + mask_HF = np.abs(freqs) >= cutoff_HF + fft_filtered_HF = np.zeros_like(fft_vals, dtype=complex) + fft_filtered_HF[mask_HF] = fft_vals[mask_HF] - # Max amplitude après filtrage - max_amplitude = np.max(np.abs(signal_filtered)) + # 7. Filtrage (garde basses fréquences) + cutoff_BF = 10 * F_max + mask_BF = np.abs(freqs) <= cutoff_BF + fft_filtered_BF = np.zeros_like(fft_vals, dtype=complex) + fft_filtered_BF[mask_BF] = fft_vals[mask_BF] + + # 8. Reconstruction + signal_filtered_HF = np.real(ifft(fft_filtered_HF)) + signal_filtered_BF = np.real(ifft(fft_filtered_BF)) + lambda_reconstructed = np.interp(k_linspace, k_vals[::-1], wavelengths[::-1]) + + # Masque dans la plage [550, 700] nm + mask_Cam_Sensitivity = (lambda_reconstructed >= 550) & (lambda_reconstructed <= 700) + + # 9. Affichage reconstruction if plot: - plt.figure(figsize=(10, 6),dpi = 600) - plt.plot(signal_filtered, label='Signal filtered') - plt.xlabel('Échantillons') - plt.ylabel('Amplitude') + plt.figure(figsize=(10, 6), dpi=150) + plt.plot(lambda_reconstructed, intensities_k, '--', label='Original interpolé') + plt.plot(lambda_reconstructed, signal_filtered_HF,'--', color='gray') + + plt.plot(lambda_reconstructed[mask_Cam_Sensitivity], signal_filtered_HF[mask_Cam_Sensitivity], + color='orange', label='Spectre filtré HF') + plt.plot(lambda_reconstructed, signal_filtered_BF, + color='red', label='Spectre filtré BF') + + plt.xlabel('Wavelength (nm)') + plt.ylabel('Intensity') plt.legend() plt.show() - print(f"Amplitude Mal filtered : {max_amplitude:.4f}") - return max_amplitude + max_amplitude = np.max(np.abs(signal_filtered_HF[mask_Cam_Sensitivity])) + + return max_amplitude, signal_filtered_BF, lambda_reconstructed + diff --git a/optifik/io.py b/optifik/io.py index 8a9bc12..bcbb35c 100644 --- a/optifik/io.py +++ b/optifik/io.py @@ -1,7 +1,9 @@ import numpy as np -def load_spectrum(spectrum_path, lambda_min=0, lambda_max=np.inf, +def load_spectrum(spectrum_path, + lambda_min=0, + lambda_max=np.inf, delimiter=','): """ Load a spectrum file. diff --git a/optifik/minmax.py b/optifik/minmax.py index 41b7b58..4085332 100644 --- a/optifik/minmax.py +++ b/optifik/minmax.py @@ -15,21 +15,21 @@ plt.rcParams.update({ from .utils import OptimizeResult - -def thickness_from_minmax(lambdas, - intensities, + +def thickness_from_minmax(wavelengths, + intensities, refractive_index, - min_peak_prominence, + min_peak_prominence, min_peak_distance=10, - method='linreg', + method='linreg', plot=None): - - """ + + """ Return the thickness from a min-max detection. - + Parameters ---------- - lambdas : array + wavelengths : array Wavelength values in nm. intensities : array Intensity values. @@ -44,12 +44,12 @@ def thickness_from_minmax(lambdas, for Randon Sampling Consensus. debug : boolean, optional Show plots of peak detection and lin regression. - + Returns ------- results : Instance of `OptimizeResult` class. The attribute `thickness` gives the thickness value in nm. - + Notes ----- For more details about `min_peak_prominence` and `min_peak_distance`, @@ -63,11 +63,12 @@ def thickness_from_minmax(lambdas, peaks.sort() k_values = np.arange(len(peaks)) - n_over_lambda = refractive_index[peaks][::-1] / lambdas[peaks][::-1] + n_over_lambda = refractive_index[peaks][::-1] / wavelengths[peaks][::-1] if k_values.size < 2: # Can't fit if less than two points. - return np.nan + return OptimizeResult(thickness=np.nan) + if isinstance(refractive_index, np.ndarray): refractive_index = refractive_index[peaks][::-1] @@ -77,7 +78,7 @@ def thickness_from_minmax(lambdas, min_samples = 2 # Scikit-image data = np.column_stack([k_values, n_over_lambda]) - model_robust, inliers = ransac(data, LineModelND, + model_robust, inliers = ransac(data, LineModelND, min_samples=min_samples, residual_threshold=residual_threshold, max_trials=100) @@ -85,7 +86,7 @@ def thickness_from_minmax(lambdas, thickness_minmax = 1 / slope / 4 # Scikit-learn - #X, y = k_values.reshape(-1, 1), 1/lambdas[peaks][::-1] + #X, y = k_values.reshape(-1, 1), 1/wavelengths[peaks][::-1] ## Fit line using all data #lr = linear_model.LinearRegression() diff --git a/optifik/scheludko.py b/optifik/scheludko.py index 3b49b96..5c02e08 100644 --- a/optifik/scheludko.py +++ b/optifik/scheludko.py @@ -1,6 +1,4 @@ import numpy as np - - from scipy.optimize import curve_fit import matplotlib.pyplot as plt @@ -16,20 +14,21 @@ from .io import load_spectrum from .utils import OptimizeResult from .analysis import finds_peak -def thickness_scheludko_at_order(wavelengths, - intensity, - interference_order, - refractive_index, + +def thickness_scheludko_at_order(wavelengths, + intensities, + interference_order, + refractive_index, Imin=None): """ Compute the film thickness for a given interference order. Parameters ---------- - wavelengths : TYPE - DESCRIPTION. - intensity : TYPE - DESCRIPTION. + wavelengths : array + Wavelength values in nm. + intensities : array + Intensity values. interference_order: TYPE DESCRIPTION. refractive_index : TYPE @@ -44,27 +43,27 @@ def thickness_scheludko_at_order(wavelengths, """ if Imin is None: - Imin = np.min(intensity) + Imin = np.min(intensities) n = refractive_index m = interference_order - I = (np.asarray(intensity) - Imin) / (np.max(intensity) - Imin) - + I = (np.asarray(intensities) - Imin) / (np.max(intensities) - Imin) + prefactor = wavelengths / (2 * np.pi * n) argument = np.sqrt(I / (1 + (1 - I) * (n**2 - 1)**2 / (4 * n**2))) - + if m % 2 == 0: - term1 = (m / 2) * np.pi + term1 = (m / 2) * np.pi else: term1 = ((m+1) / 2) * np.pi - + term2 = (-1)**m * np.arcsin(argument) - + return prefactor * (term1 + term2) -def Delta(wavelengths, thickness, interference_order, refractive_index): + """ Calculates the Delta value for arrays of wavelengths, thicknesses h and r_indexs n. @@ -78,6 +77,29 @@ def Delta(wavelengths, thickness, interference_order, refractive_index): - delta: ndarray of corresponding Δ values """ + +def Delta(wavelengths, thickness, interference_order, refractive_index): + """ + + + Parameters + ---------- + wavelengths : array + Wavelength values in nm. + thickness : array_like (or float) + Film thickness. + interference_order : int + Interference order. + refractive_index : array_like (or float) + Refractive index. + + Returns + ------- + TYPE + DESCRIPTION. + + """ + # ensure that the entries are numpy arrays wavelengths = np.asarray(wavelengths) h = np.asarray(thickness) @@ -103,36 +125,50 @@ def Delta(wavelengths, thickness, interference_order, refractive_index): return (A * (1 + alpha)) / (1 + A * alpha) + def Delta_fit(xdata, thickness, interference_order): """ - Wrapper pour curve_fit : on fixe m, et lambda & n seront passés comme "xdata" - """ - lambd, n = xdata - return Delta(lambd, thickness, interference_order, n) - - - - -def thickness_from_scheludko(lambdas, - smoothed_intensities, - refractive_index, - min_peak_prominence, - plot=None): - """ - + Wrapper on Delta() for curve_fit. Parameters ---------- - lambdas : TYPE - DESCRIPTION. - raw_intensities : TYPE - DESCRIPTION. - smoothed_intensities : TYPE - DESCRIPTION. - refractive_index : TYPE - DESCRIPTION. - plot : TYPE, optional - DESCRIPTION. The default is None. + xdata : tuple + (wavelengths, n) + thickness : array_like (or float) + Film thickness. + interference_order : int + Interference order. + + Returns + ------- + ndarray + Delta values. + + """ + lambdas, n = xdata + return Delta(lambdas, thickness, interference_order, n) + + + + +def thickness_from_scheludko(wavelengths, + intensities, + refractive_index, + min_peak_prominence, + plot=None): + """ + + + Parameters + ---------- + wavelengths : array + Wavelength values in nm. + intensities : array + Intensity values. + refractive_index : scalar, optional + Value of the refractive index of the medium. + plot : bool, optional + Display a curve, useful for checking or debuging. The default is None. Returns ------- @@ -142,149 +178,153 @@ def thickness_from_scheludko(lambdas, """ max_tested_order = 12 r_index = refractive_index - - peaks_min, peaks_max = finds_peak(lambdas, smoothed_intensities, - min_peak_prominence=min_peak_prominence, - plot=False) - - - lambda_min = lambdas[peaks_min[-1]] - lambda_max = lambdas[peaks_max[-1]] - - # On s'assure que lambda_min < lambda_max + + peaks_min, peaks_max = finds_peak(wavelengths, intensities, + min_peak_prominence=min_peak_prominence, + plot=False) + + # Get the last oscillation peaks + lambda_min = wavelengths[peaks_min[-1]] + lambda_max = wavelengths[peaks_max[-1]] + + # Order them lambda_start = min(lambda_min, lambda_max) - lambda_end = max(lambda_min, lambda_max) - - # On crée le masque pour ne garder que les lambdas entre les deux extrema - mask = (lambdas >= lambda_start) & (lambdas <= lambda_end) - lambdas_masked = lambdas[mask] + lambda_stop = max(lambda_min, lambda_max) + + # mask input data + mask = (wavelengths >= lambda_start) & (wavelengths <= lambda_stop) + wavelengths_masked = wavelengths[mask] r_index_masked = r_index[mask] - intensities_masked = smoothed_intensities[mask] + intensities_masked = intensities[mask] + + min_ecart = np.inf best_m = None meilleure_h = None - + if plot: plt.figure(figsize=(10, 6),dpi =600) plt.ylabel(r'$h$ ($\mathrm{{nm}}$)') plt.xlabel(r'$\lambda$ ($ \mathrm{nm} $)') - - + + for m in range(0, max_tested_order+1): - h_values = thickness_scheludko_at_order(lambdas_masked, intensities_masked, m, r_index_masked) - + h_values = thickness_scheludko_at_order(wavelengths_masked, intensities_masked, m, r_index_masked) + if plot: - plt.plot(lambdas_masked, h_values,'.', markersize=3, label=f"Épaisseur du film (Scheludko, m={m})") + plt.plot(wavelengths_masked, h_values,'.', markersize=3, label=f"Épaisseur du film (Scheludko, m={m})") ecart = np.max(h_values)-np.min(h_values) - + print(f"Écart pour m={m} : {ecart:.3f} nm") - + if ecart < min_ecart: min_ecart = ecart best_m = m meilleure_h = h_values - - if plot: - plt.figure(figsize=(10, 6), dpi=600) - + + DeltaVrai = (intensities_masked -np.min(intensities_masked))/(np.max(intensities_masked) -np.min(intensities_masked)) #DeltaVrai = (intensities_raw_masked -np.min(intensities_raw_masked))/(np.max(intensities_raw_masked) -np.min(intensities_raw_masked)) - - DeltaScheludko = Delta(lambdas_masked, np.mean(meilleure_h), best_m, r_index_masked) + + DeltaScheludko = Delta(wavelengths_masked, np.mean(meilleure_h), best_m, r_index_masked) #print(np.mean(meilleure_h),np.std(meilleure_h)) + if plot: - plt.plot(lambdas_masked,DeltaVrai,'bo-', markersize=2,label=r'$\mathrm{{Smoothed}}\ \mathrm{{Data}}$') - plt.plot(lambdas_masked,DeltaScheludko,'go-', markersize=2,label = rf'$\mathrm{{Scheludko}}\ (h = {np.mean(meilleure_h):.1f} \pm {np.std(meilleure_h):.1f}\ \mathrm{{nm}})$') - - - xdata = (lambdas_masked, r_index_masked) + plt.figure(figsize=(10, 6), dpi=600) + plt.plot(wavelengths_masked, DeltaVrai, + 'bo-', markersize=2, label=r'$\mathrm{{Smoothed}}\ \mathrm{{Data}}$') + plt.plot(wavelengths_masked, DeltaScheludko, + 'go-', markersize=2, label = rf'$\mathrm{{Scheludko}}\ (h = {np.mean(meilleure_h):.1f} \pm {np.std(meilleure_h):.1f}\ \mathrm{{nm}})$') + + + xdata = (wavelengths_masked, r_index_masked) popt, pcov = curve_fit(lambda x, h: Delta_fit(x, h, m), xdata, DeltaVrai, p0=[np.mean(meilleure_h)]) fitted_h = popt[0] - - + + if plot: - plt.plot(lambdas_masked, Delta(lambdas_masked, fitted_h, best_m, r_index_masked ), 'ro-',markersize=2, label=rf'$\mathrm{{Fit}}\ (h = {fitted_h:.1f}\pm {np.sqrt(pcov[0][0]):.1f} \ \mathrm{{nm}})$') + plt.plot(wavelengths_masked, Delta(wavelengths_masked, fitted_h, best_m, r_index_masked ), 'ro-',markersize=2, label=rf'$\mathrm{{Fit}}\ (h = {fitted_h:.1f}\pm {np.sqrt(pcov[0][0]):.1f} \ \mathrm{{nm}})$') plt.legend() plt.ylabel(r'$\Delta$') plt.xlabel(r'$\lambda$ ($ \mathrm{nm} $)') - - + + return OptimizeResult(thickness=fitted_h ,) -def thickness_for_order0(lambdas, - smoothed_intensities, - refractive_index, +def thickness_for_order0(wavelengths, + intensities, + refractive_index, min_peak_prominence, plot=None): - - + + File_I_min = 'tests/spectre_trou/000043641.xy' r_index = refractive_index - - peaks_min, peaks_max = finds_peak(lambdas, smoothed_intensities, + + peaks_min, peaks_max = finds_peak(wavelengths, intensities, min_peak_prominence=min_peak_prominence, plot=False) - - - - - lambdas_I_min, intensities_I_min = load_spectrum(File_I_min, lambda_min=450) - - lambda_unique = lambdas[peaks_max[0]] - # On crée le masque pour ne garder que les lambdas superieures a lambdas unique - mask = lambdas >= lambda_unique - lambdas_masked = lambdas[mask] + + + wavelengths_I_min, intensities_I_min = load_spectrum(File_I_min, lambda_min=450) + + lambda_unique = wavelengths[peaks_max[0]] + + + # On crée le masque pour ne garder que les wavelengths superieures a wavelengths unique + mask = wavelengths >= lambda_unique + wavelengths_masked = wavelengths[mask] r_index_masked = r_index[mask] - intensities_masked = smoothed_intensities[mask] + intensities_masked = intensities[mask] intensities_I_min_masked =intensities_I_min[mask] min_ecart = np.inf best_m = None meilleure_h = None - + m = 0 - h_values = thickness_scheludko_at_order(lambdas_masked, - intensities_masked, - 0, - r_index_masked, + h_values = thickness_scheludko_at_order(wavelengths_masked, + intensities_masked, + 0, + r_index_masked, Imin=intensities_I_min_masked) - + if plot: plt.figure(figsize=(10, 6), dpi=600) - plt.plot(lambdas_masked, h_values, label=r"Épaisseur du film (Scheludko, m=0)") - + plt.plot(wavelengths_masked, h_values, label=r"Épaisseur du film (Scheludko, m=0)") + ecart = np.max(h_values) - np.min(h_values) best_m = m meilleure_h = h_values - - - + + + DeltaVrai = (intensities_masked -np.min(intensities_I_min_masked))/(np.max(intensities_masked) -np.min(intensities_I_min_masked)) #DeltaVrai = (intensities_masked -np.min(intensities_masked))/(np.max(intensities_masked) -np.min(intensities_masked)) - DeltaScheludko = Delta(lambdas_masked, np.mean(meilleure_h), best_m, r_index_masked) + DeltaScheludko = Delta(wavelengths_masked, np.mean(meilleure_h), best_m, r_index_masked) #print(np.mean(meilleure_h),np.std(meilleure_h)) - - + + if plot: plt.figure(figsize=(10, 6), dpi=600) - plt.plot(lambdas_masked,DeltaVrai,'bo-', markersize=2,label=r'$\mathrm{{Raw}}\ \mathrm{{Data}}$') - plt.plot(lambdas_masked,DeltaScheludko,'ro-', markersize=2,label = rf'$\mathrm{{Scheludko}}\ (h = {np.mean(meilleure_h):.1f} \pm {np.std(meilleure_h):.1f}\ \mathrm{{nm}})$') + plt.plot(wavelengths_masked,DeltaVrai,'bo-', markersize=2,label=r'$\mathrm{{Raw}}\ \mathrm{{Data}}$') + plt.plot(wavelengths_masked,DeltaScheludko,'ro-', markersize=2,label = rf'$\mathrm{{Scheludko}}\ (h = {np.mean(meilleure_h):.1f} \pm {np.std(meilleure_h):.1f}\ \mathrm{{nm}})$') - xdata = (lambdas_masked, r_index_masked) + xdata = (wavelengths_masked, r_index_masked) popt, pcov = curve_fit(lambda x, h: Delta_fit(x, h, m), xdata, DeltaVrai, p0=[np.mean(meilleure_h)]) fitted_h = popt[0] - + if plot: - plt.plot(lambdas_masked, Delta(lambdas_masked, fitted_h, best_m, r_index_masked ), 'go-',markersize=2, label=rf'$\mathrm{{Fit}}\ (h = {fitted_h:.1f}\pm {np.sqrt(pcov[0][0]):.1f} \ \mathrm{{nm}})$') + plt.plot(wavelengths_masked, Delta(wavelengths_masked, fitted_h, best_m, r_index_masked ), 'go-',markersize=2, label=rf'$\mathrm{{Fit}}\ (h = {fitted_h:.1f}\pm {np.sqrt(pcov[0][0]):.1f} \ \mathrm{{nm}})$') plt.legend() plt.ylabel(r'$\Delta$') plt.xlabel(r'$\lambda$ (nm)') - return OptimizeResult(thickness=fitted_h ,) \ No newline at end of file + return OptimizeResult(thickness=fitted_h ,) diff --git a/tests/test_SV1.py b/tests/test_SV1.py index 241d2fb..88be384 100644 --- a/tests/test_SV1.py +++ b/tests/test_SV1.py @@ -1,12 +1,13 @@ import os.path import numpy as np -from numpy.testing import assert_allclose +from numpy.testing import assert_allclose, assert_equal import pytest from optifik.minmax import thickness_from_minmax from optifik.io import load_spectrum from optifik.analysis import smooth_intensities from optifik.analysis import finds_peak +from optifik.fft import Prominence_from_fft import yaml @@ -26,13 +27,27 @@ def test_minmax(spectrum_path, expected): lambdas, raw_intensities = load_spectrum(spectrum_path, lambda_min=450) smoothed_intensities = smooth_intensities(raw_intensities) + assert_equal(len(lambdas), len(smoothed_intensities)) + indice = 1.324188 + 3102.060378 / (lambdas**2) prominence = 0.02 + #prominence, s, w = Prominence_from_fft(lambdas, smoothed_intensities, indice) + #prominence *= 10 + thickness_minmax = thickness_from_minmax(lambdas, smoothed_intensities, refractive_index=indice, min_peak_prominence=prominence) + + #prominence, s, w = Prominence_from_fft(lambdas, smoothed_intensities, indice) + #print(f'Prom: {prominence}') + #indice = 1.324188 + 3102.060378 / (w**2) + #thickness_minmax = thickness_from_minmax(w, + # s, + # refractive_index=indice, + # min_peak_prominence=prominence) + result = thickness_minmax.thickness assert_allclose(result, expected, rtol=1e-1) diff --git a/tests/test_basic.py b/tests/test_basic.py index 10dba36..3038c3a 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -8,7 +8,7 @@ from optifik.minmax import thickness_from_minmax from optifik.scheludko import thickness_from_scheludko from optifik.scheludko import thickness_for_order0 from optifik.analysis import smooth_intensities -from optifik.analysis import Prominence_from_fft +from optifik.fft import Prominence_from_fft from optifik.io import load_spectrum @@ -39,10 +39,10 @@ def test_minmax_ransac(): smoothed_intensities = smooth_intensities(raw_intensities) indice = 1.324188 + 3102.060378 / (lambdas**2) - prominence = Prominence_from_fft(lambdas=lambdas, - intensities=smoothed_intensities, - refractive_index=indice, - plot=False) + prominence, signal, wavelength = Prominence_from_fft(lambdas, + smoothed_intensities, + refractive_index=indice, + plot=False) thickness_minmax = thickness_from_minmax(lambdas, @@ -66,7 +66,10 @@ def test_scheludko_4peaks(): smoothed_intensities = smooth_intensities(raw_intensities) indice = 1.324188 + 3102.060378 / (lambdas**2) - prominence = Prominence_from_fft(lambdas=lambdas, intensities=smoothed_intensities, refractive_index=indice, plot=False) + prominence, signal, wavelength = Prominence_from_fft(lambdas, + smoothed_intensities, + refractive_index=indice, + plot=False) result = thickness_from_scheludko(lambdas, smoothed_intensities,