From c61a1e4c5ea96f80036da52902658a53d8e11fb2 Mon Sep 17 00:00:00 2001 From: Le Stagiaire Date: Mon, 24 Jun 2024 11:36:55 +0200 Subject: [PATCH] feat: Add UTSR module for 3D app chore: Update UTSR module with necessary functions for simulation, data processing, and visualization. --- 3D_app/UTSR.py | 58 +++++++ 3D_app/utsr/DatNDTRead.py | 122 ++++++++++++++ 3D_app/utsr/GenerateModelFilter.py | 171 +++++++++++++++++++ 3D_app/utsr/ModS2.py | 77 +++++++++ 3D_app/utsr/ModS2T.py | 89 ++++++++++ 3D_app/utsr/PlotResults.py | 137 +++++++++++++++ 3D_app/utsr/SimulNDT.py | 239 ++++++++++++++++++++++++++ 3D_app/utsr/UTSRModS2.py | 260 +++++++++++++++++++++++++++++ 8 files changed, 1153 insertions(+) create mode 100644 3D_app/UTSR.py create mode 100644 3D_app/utsr/DatNDTRead.py create mode 100644 3D_app/utsr/GenerateModelFilter.py create mode 100644 3D_app/utsr/ModS2.py create mode 100644 3D_app/utsr/ModS2T.py create mode 100644 3D_app/utsr/PlotResults.py create mode 100644 3D_app/utsr/SimulNDT.py create mode 100644 3D_app/utsr/UTSRModS2.py diff --git a/3D_app/UTSR.py b/3D_app/UTSR.py new file mode 100644 index 0000000..8fb36ba --- /dev/null +++ b/3D_app/UTSR.py @@ -0,0 +1,58 @@ +import numpy as np +from utsr.UTSRModS2 import utsr_mod_s2 +from utsr.GenerateModelFilter import generate_model_filter +from utsr.SimulNDT import simulNDT +from utsr.DatNDTRead import datNDTread +from utsr.PlotResults import plotResults + + +def nextpow2(n): + return int(np.ceil(np.log2(n))) + + +# Définir les paramètres de l'essai +TendGate = 24.98 +cl = 5859.4 +d = 6.35 +fc = 4.6 +bw = 4.4 +A0 = 2.02 +alpha_atenu = 35.48e-3 +T1 = int(100e-9) # Largeur de l'impulsion d'excitation + +# ROI (en mm) +minX = 0 +maxX = 30 +minZ = 18 +maxZ = 58 + +# Paramètre de défaut (discontinuité) +DepthCentre = 40 + +# Simulation de l'essai NDT pour obtenir un modèle de discontinuité +ptPoint, setup = simulNDT(simul=False, TendGate=TendGate, d=d, cl=cl, fc=fc, bw=bw, minX=minX, maxX=maxX, DefectType='empty', flDepthCentre=DepthCentre) + +# Obtenir les données de calibration +ERMv = ptPoint['CscanData']['Cl'] / 2 +dt = (ptPoint['timeScale'][1] - ptPoint['timeScale'][0]) * 1e-6 +du = (ptPoint['CscanData']['X'][1] - ptPoint['CscanData']['X'][0]) * 1e-3 +tau0 = ptPoint['timeScale'][0] * 1e-6 +nt = 2**(nextpow2(ptPoint['CscanData']['AscanPoints']) + 1) +nu = ptPoint['CscanData']['Rows'] * 2 +f = np.fft.fftfreq(nt, d=dt) +ku = np.fft.fftfreq(nu, d=du) + +# Créer les données de modèle et de filtre +ModPoint, HedPoint = generate_model_filter(ptPoint, TrFreqSp='cossquare', T1=T1) + +# Définir le paramètre de régularisation +alpha = 0.02 + +# Charger les données acquises +ptAco40dB_1 = datNDTread(baseFolder='./AcquiredData/Specimen01', gain=40, A0=A0, alpha=alpha_atenu) + +# Appliquer l'algorithme UTSR aux données acquises +Result_Med_irls = utsr_mod_s2(ptAco40dB_1, plotOut=False, title='UTSR solution', SysModel=ModPoint, deltaX=du, minZ=minZ*1e-3, maxZ=maxZ*1e-3, alpha=alpha) + +# Afficher les résultats +plotResults(Result_Med_irls, colored=False, output='normalized', PostProc='env', API_calc=True) diff --git a/3D_app/utsr/DatNDTRead.py b/3D_app/utsr/DatNDTRead.py new file mode 100644 index 0000000..4726bc6 --- /dev/null +++ b/3D_app/utsr/DatNDTRead.py @@ -0,0 +1,122 @@ +import os +import numpy as np + +def datNDTread(*args, **kwargs): + # Define valid wave types + validWaveType = ['L', 'l', 's', 'S'] + + # Parse input arguments + TsGate = kwargs.get('TsGate', 0) + TendGate = kwargs.get('TendGate', 25 - 0.02) + deltaT = kwargs.get('deltaT', 0.02) + minX = kwargs.get('minX', 0) + maxX = kwargs.get('maxX', 60) + deltaX = kwargs.get('deltaX', 1) + backDepth = kwargs.get('backDepth', 60) + Density = kwargs.get('Density', 7.86) + fc = kwargs.get('fc', 5) + bw = kwargs.get('bw', 5) + WaveType = kwargs.get('WaveType', validWaveType[0]) + BeamAngle = kwargs.get('BeamAngle', 0) + d = kwargs.get('d', 6.35) + baseFolder = kwargs.get('baseFolder', '.') + gain = kwargs.get('gain', np.inf) + A0 = kwargs.get('A0', 1) + alpha = kwargs.get('alpha', 0) + + # Get file list in base folder + files = [f for f in os.listdir(baseFolder) if f.endswith('.dat')] + + # Initialize scanData structure + scanData = {} + scanData['CscanData'] = {} + scanData['AscanValues'] = np.zeros((1, len(files))) + + # Import data as columns + backWallIdx = np.zeros(len(files)) + for idx, f in enumerate(files): + col = int(f.split('pos')[1].split('.dat')[0]) + 1 + dat = np.loadtxt(os.path.join(baseFolder, f)) + + lenDat = len(dat) + idxm = np.argmin(dat[:round(lenDat * 1 / 3)]) + limIdx = np.argmin(dat[round(lenDat * 2 / 3):]) + round(lenDat * 2 / 3) - 1 + + idx2m = np.argmax(dat[round(lenDat * 2 / 3):limIdx]) + round(lenDat * 2 / 3) - 1 + + backWallIdx[col - 1] = idx2m - idxm + scanData['AscanValues'][:lenDat, col - 1] = np.pad(dat[idxm:], (0, lenDat - idxm), 'constant') + + # Create timeScale element + scanData['CscanData']['AscanPoints'] = scanData['AscanValues'].shape[0] + scanData['timeScale'] = np.arange(0, deltaT * scanData['CscanData']['AscanPoints'], deltaT) + + # Define the time window for simulation + t = np.arange(TsGate, TendGate, deltaT) + idxTsGate = np.where(t > TsGate)[0][0] - 1 + idxTendGate = np.where(t >= TendGate)[0][0] + + # Update scanData with argument values + scanData['CscanData']['TsGate'] = TsGate + scanData['CscanData']['TendGate'] = TendGate + scanData['CscanData']['X'] = np.arange(0, len(files) * deltaX, deltaX) + + backWallIdx = np.array([backWallIdx[0], backWallIdx[-1]]) + Cl = np.mean(backDepth * 2000 / scanData['timeScale'][backWallIdx]) + scanData['CscanData']['Cl'] = Cl + scanData['CscanData']['Cs'] = 0 + scanData['CscanData']['Density'] = Density + scanData['CscanData']['Frequency'] = fc + scanData['CscanData']['Bandwidth'] = bw + scanData['CscanData']['WaveType'] = WaveType + scanData['CscanData']['BeamAngle'] = BeamAngle + scanData['CscanData']['ProbeDiameter'] = d + + # Others elements + scanData['CscanData']['MaxSdtRef'] = 0 + scanData['CscanData']['Y'] = 0 + scanData['CscanData']['Cols'] = 1 + scanData['CscanData']['CscanValues'] = np.zeros((scanData['CscanData']['Rows'], scanData['CscanData']['Cols'])) + + if WaveType in ['L', 'l']: + c = scanData['CscanData']['Cl'] + elif WaveType in ['S', 's']: + c = scanData['CscanData']['Cs'] + + scanData['CscanData']['Wavelength'] = c / (scanData['CscanData']['Frequency'] * 1000) + scanData['CscanData']['Nearfield'] = (scanData['CscanData']['ProbeDiameter'] / 2) ** 2 / scanData['CscanData']['Wavelength'] + scanData['CscanData']['NearfieldTime'] = scanData['CscanData']['Nearfield'] / (c / 2000) + scanData['CscanData']['ProbeEllipX'] = 0 + scanData['CscanData']['ProbeEllipY'] = 0 + scanData['CscanData']['TrueAngle'] = scanData['CscanData']['BeamAngle'] + scanData['CscanData']['CalValue'] = 0 + + scanData['CscanData']['DefectType'] = 'empty' + scanData['CscanData']['PosX'] = 0 + scanData['CscanData']['PosY'] = 0 + scanData['CscanData']['DepthCentre'] = 0 + scanData['CscanData']['Diameter'] = 0 + scanData['CscanData']['Height'] = 0 + scanData['CscanData']['Tilt'] = 0 + + # Adjust A-scan amplitude + if np.isinf(gain): + # Gain = Inf -> normalize amplitude to 1 + MaxSignal = np.max(np.abs(scanData['AscanValues'])) + scanData['AscanValues'] /= MaxSignal + else: + # Add gain + gain = 10**(gain / 20) + scanData['AscanValues'] *= (0.5 / 2048 / gain) + + # Apply TGV + zz = scanData['timeScale'] * Cl / 2000 + ffc = (1 / A0) * np.exp(alpha * zz) + ffc[ffc < 1] = 1 + scanData['AscanValues'] *= ffc[:, np.newaxis] + + # Calculate elements after the simulation session + scanData['CscanData']['MaxSignal'] = np.max(np.abs(scanData['AscanValues'])) + scanData['CscanData']['CscanValues'] = 10 * np.log10(np.max(np.abs(scanData['AscanValues']), axis=0)) + + return scanData diff --git a/3D_app/utsr/GenerateModelFilter.py b/3D_app/utsr/GenerateModelFilter.py new file mode 100644 index 0000000..d5798d9 --- /dev/null +++ b/3D_app/utsr/GenerateModelFilter.py @@ -0,0 +1,171 @@ +import numpy as np +from scipy.special import j1 as besselj1, struve + + +def jinc(x): + return 2 * besselj1(x) / x + + +def gaussian(f, fc, BW): + return np.exp(-((f - fc) ** 2) / (2 * (BW / 2.355) ** 2)) + + +def cossquare(f, fc, BW): + return (np.cos(np.pi * (f - fc) / BW)) ** 2 + + +def generate_model_filter(scanData, pad=True, ep=0.03, TrFreqSp="gaussian", Rx=1, T1=0): + """ + Generate Model and Filter to Stolt transforms + """ + + def scanDataValidationFcn(scanData): + return all(key in scanData for key in ["CscanData", "timeScale", "AscanValues"]) + + def checkTrFreqSp(TrFreqSp): + validTrFreqSp = ["gaussian", "cossquare"] + return TrFreqSp in validTrFreqSp + + # Verify scanData parameter + assert scanDataValidationFcn(scanData), "scanData invalid structure" + assert checkTrFreqSp(TrFreqSp), "Invalid Transducer frequency spectrum" + + # Get scan data informations + if scanData["CscanData"]["WaveType"] == "L": + c = scanData["CscanData"]["Cl"] + else: + c = scanData["CscanData"]["Cs"] + ERMv = c / 2 + + scanX = scanData["CscanData"]["X"] * 1e-3 + deltaScanX = (scanX[1] - scanX[0]) / Rx + deltaT = (scanData["timeScale"][1] - scanData["timeScale"][0]) * 1e-6 + + # Get acquisition matrix sizes + nt0, nu0 = scanData["AscanValues"].shape + + # Include scale factor to X resolution + nu0 = (nu0 - 1) * Rx + 1 + + # Zero-padding + if pad: + nt = 2 ** (np.ceil(np.log2(nt0)) + 1).astype(int) + nu = 2 * nu0 + else: + nt = nt0 + nu = nu0 + + # Calculate matched filters + C1, Hed = pCalcMatchedFilter( + nt, + deltaT, + nu, + deltaScanX, + scanData["CscanData"]["ProbeDiameter"] * 1e-3, + scanData["CscanData"]["Frequency"] * 1e6, + scanData["CscanData"]["Bandwidth"] * 1e6, + ERMv, + TrFreqSp, + T1, + ) + + # Calculate flaw response + SA = pCalcScatteringAmplitude( + nt, + deltaT, + c, + scanData["CscanData"]["DefectType"], + scanData["CscanData"]["Diameter"] * 1e-3 / 2, + scanData["CscanData"]["Tilt"], + ) + + # Build Model and Filter + Model = C1 * (Hed * SA).T + AA = np.conj(C1) / (C1 * np.conj(C1) + ep) + Filter = AA * Hed.T + + return Model, Filter, C1, Hed, SA + + +def pCalcMatchedFilter(nt, dt, nx, dx, d, fc, BW, ERMv, TrFreqSp, T1): + # Calculate frequency vectors + f = np.fft.fftfreq(nt, d=dt) + kx = np.fft.fftfreq(nx, d=dx) + gkx, gf = np.meshgrid(kx, f) + + # Calculate gst and create mask + gst = gkx / (gf / ERMv) + mask = np.abs(gst) <= 1 + + # Calculate C1 using Jinc function and apply mask + C1 = (np.sinc((2 * np.pi * gkx / 2) * (d / 2))) ** 2 * mask + + # Calculate Hed based on TrFreqSp + if TrFreqSp == "gaussian": + Hed = (gaussian(f, fc, BW)) ** 2 + elif TrFreqSp == "cossquare": + Hed = (cossquare(f, fc, BW)) ** 2 + else: + raise ValueError("Unsupported TrFreqSp type") + + # Transducer excitation signal spectrum + ExcTransd = -np.exp(-2j * np.pi * f * (T1 / 2)) + if T1 != 0: + ExcTransd *= 2 * np.pi * np.sinc(2 * np.pi * f * (T1 / 2)) + + # Define 10 MHz low-pass filter of the pulser + filt10MHz = np.zeros_like(f) + filt10MHz[np.abs(f) <= 10e6] = 1 + + # Eq. 12.13 Schmerr (1998) + Hed = 1500 * np.pi * Hed * ExcTransd * filt10MHz + C1 *= ((np.pi / 2) * ((2 * (d / 2) ** 2 * (2 * np.pi) / ERMv) * np.abs(f))).T + + return C1, Hed + + +def pCalcScatteringAmplitude(nt, dt, c, flaw, b, ang): + f = np.fft.fftfreq(nt, d=dt) + kb = -(2 * np.pi * b * f) / c + kb += np.finfo(float).eps * (kb == 0) + + iang = ang * np.pi / 180 + if flaw in {"circ-crack", "penny-shaped"}: + arg = 2 * np.sin(iang) * kb + arg += np.finfo(float).eps * (arg == 0) + SA = 1j * kb * b * np.cos(iang) * (besselj1(arg) / arg) + elif flaw == "side-drilled": + SA = ( + -((kb / 2) * (besselj1(2 * kb) - 1j * struve(1, 2 * kb)) + 1j * kb / np.pi) + / 190 + ) + SA *= np.sqrt(2j * np.pi / (kb / b)) + elif flaw == "strip-like": + SA = np.ones_like(kb) / 1000 + elif flaw == "spherical": + SA = -(b / 2) * np.exp(-1j * kb) * (np.exp(-1j * kb) - np.sin(kb) / kb) + elif flaw == "surface-breaking": + SA = np.ones_like(kb) / 1000 + else: + SA = np.ones_like(kb) / 1000 + return SA + + +# Exemple d'utilisation: +# scanData = { +# 'CscanData': { +# 'WaveType': 'L', +# 'Cl': 5900, +# 'Cs': 3200, +# 'X': np.linspace(0, 100, 256), +# 'ProbeDiameter': 0.02, +# 'Frequency': 5, +# 'Bandwidth': 2, +# 'DefectType': 'circ-crack', +# 'Diameter': 0.01, +# 'Tilt': 0 +# }, +# 'timeScale': np.linspace(0, 10, 256), +# 'AscanValues': np.random.rand(256, 256) +# } +# Model, Filter, C1, Hed, SA = generate_model_filter(scanData) diff --git a/3D_app/utsr/ModS2.py b/3D_app/utsr/ModS2.py new file mode 100644 index 0000000..5d6819a --- /dev/null +++ b/3D_app/utsr/ModS2.py @@ -0,0 +1,77 @@ +import numpy as np +from scipy.fft import fft2, ifft2, fftshift, ifftshift +from scipy.interpolate import interp1d + +def linterp_adj(f, fkz, ftimage): + # Perform adjoint linear interpolation + interp_real = interp1d(f, ftimage.real, axis=0, fill_value="extrapolate") + interp_imag = interp1d(f, ftimage.imag, axis=0, fill_value="extrapolate") + ftdata_real = interp_real(fkz) + ftdata_imag = interp_imag(fkz) + return ftdata_real + 1j * ftdata_imag + +def mod_s2(image, dt, du, ERMv, tau0, ModTransd, R_x=1): + """ + Direct operator for the model of the inspection system + """ + # Check R_x argument + R_x = int(R_x) + assert R_x > 0, 'R_x must be positive' + + # Zero-padding + nt0, nu0 = image.shape + nt, nu = ModTransd.shape + if nt < nt0: + nt = 2 ** (np.ceil(np.log2(nt0)) + 1).astype(int) + if nu < nu0: + nu = 2 * nu0 + + # Data and image grids + f = np.fft.fftfreq(nt, d=dt) + ku = np.fft.fftfreq(nu, d=du) + ku, f = np.meshgrid(ku, f) + fkz = ERMv * np.sign(f) * np.sqrt(ku**2 + f**2 / ERMv**2) + + # Converting image to frequency domain + ftimage = fftshift(fft2(image, s=(nt, nu))) + + # Generates data spectrum (by adjoint of linear interpolation) + ftdata = linterp_adj(f[:, 0], fkz, ftimage) + if tau0 != 0: + ftdata *= np.exp(1j * 2 * np.pi * (fkz - f) * tau0) + + # Inserts transducer model (if exists) + if isinstance(ModTransd, (np.ndarray, float)): + if np.isscalar(ModTransd) or ModTransd.shape == ftdata.shape: + ftdata *= ModTransd + else: + rowShift = (nt - ModTransd.shape[0]) // 2 + colShift = (nu - ModTransd.shape[1]) // 2 + ModTransd_padded = np.zeros((nt, nu), dtype=ModTransd.dtype) + ModTransd_padded[rowShift:rowShift + ModTransd.shape[0], colShift:colShift + ModTransd.shape[1]] = ModTransd + ftdata *= ModTransd_padded + + # Converting data to time domain + data = ifft2(ifftshift(ftdata)).real + data = data[:nt0, :nu0] + + # Applies sub-sampling + if R_x > 1: + data = data[:, ::R_x] * R_x + + # Evaluates the gain + enImage = np.linalg.norm(image) + enData = np.linalg.norm(data) + gain = enData / enImage + + return data, gain + +# Example usage: +# image = np.random.rand(256, 256) # Example input image +# dt = 0.01 +# du = 0.01 +# ERMv = 1.5 +# tau0 = 0.05 +# ModTransd = np.random.rand(256, 256) # Example ModTransd +# R_x = 2 +# data, gain = mod_s2(image, dt, du, ERMv, tau0, ModTransd, R_x) diff --git a/3D_app/utsr/ModS2T.py b/3D_app/utsr/ModS2T.py new file mode 100644 index 0000000..2508658 --- /dev/null +++ b/3D_app/utsr/ModS2T.py @@ -0,0 +1,89 @@ +import numpy as np +from scipy.fft import fft2, ifft2, fftshift, ifftshift +from scipy.interpolate import interp1d + +def linterp(f, ftdata, fkz): + # Perform linear interpolation + interp_real = interp1d(f, ftdata.real, axis=0, fill_value="extrapolate") + interp_imag = interp1d(f, ftdata.imag, axis=0, fill_value="extrapolate") + ftimage_real = interp_real(fkz) + ftimage_imag = interp_imag(fkz) + return ftimage_real + 1j * ftimage_imag + +def mod_s2t(data, dt, du, ERMv, tau0, Filter, R_x=1): + """ + Adjoint operator for the model of the inspection system + """ + # Check R_x argument + R_x = int(R_x) + assert R_x > 0, 'R_x must be positive' + + # Applies linear interpolation (if necessary) + if R_x > 1: + coef_interp = triang(2 * (R_x - 1) + 1) + nt0, nu0 = data.shape + nColsROI = (nu0 - 1) * R_x + 1 + nColsPad = R_x - 2 + colsDataOrig = np.arange(1, nColsROI + 1, R_x) + nColsPad + colsDataNew = np.setdiff1d(np.arange(colsDataOrig[0], colsDataOrig[-1] + 1), colsDataOrig) + data2 = np.zeros((nt0, nColsROI + 2 * nColsPad)) + data2[:, colsDataOrig - 1] = data + data = np.zeros_like(data2) + for col in colsDataNew: + data[:, col - 1] = np.dot(data2[:, col - (R_x - 1) - 1:col + (R_x - 1)], coef_interp) + data += data2 + data = data[:, colsDataOrig[0] - 1:colsDataOrig[-1]] + + # Zero-padding + nt0, nu0 = data.shape + nt, nu = Filter.shape + if nt < nt0: + nt = 2 ** (np.ceil(np.log2(nt0)) + 1).astype(int) + if nu < nu0: + nu = 2 * nu0 + + # Data and image grids + f = np.fft.fftfreq(nt, d=dt) + ku = np.fft.fftfreq(nu, d=du) + ku, f = np.meshgrid(ku, f) + fkz = ERMv * np.sign(f) * np.sqrt(ku**2 + f**2 / ERMv**2) + + # Converting data to frequency domain + ftdata = fftshift(fft2(data, s=(nt, nu))) + if tau0 != 0: + ftdata *= np.exp(-1j * 2 * np.pi * (fkz - f) * tau0) + + # Applies filter (if it exists) + if isinstance(Filter, (np.ndarray, float)): + if np.isscalar(Filter) or Filter.shape == ftdata.shape: + ftdata *= Filter + else: + rowShift = (nt - Filter.shape[0]) // 2 + colShift = (nu - Filter.shape[1]) // 2 + Filter_padded = np.zeros((nt, nu), dtype=Filter.dtype) + Filter_padded[rowShift:rowShift + Filter.shape[0], colShift:colShift + Filter.shape[1]] = Filter + ftdata *= Filter_padded + + # Generating image spectrum (by linear interpolation) + ftimage = linterp(f[:, 0], ftdata, fkz) + + # Converting image to space domain + image = ifft2(ifftshift(ftimage)).real + image = image[:nt0, :nu0] * R_x + + # Evaluates the gain + enImage = np.linalg.norm(image) ** 2 + enData = np.linalg.norm(data) ** 2 + gain = np.sqrt(enImage / enData) + + return image, gain + +# Example usage: +# data = np.random.rand(256, 256) # Example input data +# dt = 0.01 +# du = 0.01 +# ERMv = 1.5 +# tau0 = 0.05 +# Filter = np.random.rand(256, 256) # Example Filter +# R_x = 2 +# image, gain = mod_s2t(data, dt, du, ERMv, tau0, Filter, R_x) diff --git a/3D_app/utsr/PlotResults.py b/3D_app/utsr/PlotResults.py new file mode 100644 index 0000000..aabeb31 --- /dev/null +++ b/3D_app/utsr/PlotResults.py @@ -0,0 +1,137 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.interpolate import interp1d + +def plotResults(recResult, **kwargs): + # Define valid output and post-processing options + validOutput = ['raw', 'normalized', 'dB'] + validPostproc = ['none', 'rectified', 'env'] + + # Verify recResult parameter + assert isinstance(recResult, dict) and 'output' in recResult and 'roiX' in recResult and 'roiZ' in recResult and 'methodMsg' in recResult, 'recResult invalid structure' + + # Parse input arguments using kwargs + minX = kwargs.get('minX', recResult['roiX'][0]) + maxX = kwargs.get('maxX', recResult['roiX'][-1]) + minZ = kwargs.get('minZ', recResult['roiZ'][0]) + maxZ = kwargs.get('maxZ', recResult['roiZ'][-1]) + output_type = kwargs.get('output', validOutput[0]) + range_dB = kwargs.get('range', 25) + post_proc = kwargs.get('Postproc', validPostproc[0]) + profile_plot = kwargs.get('profile', False) + prof_range_dB = kwargs.get('prof_range', 50) + animation_plot = kwargs.get('animation', False) + res_calc_plot = kwargs.get('res_calc', False) + API_calc = kwargs.get('API_calc', False) + API_level = kwargs.get('API_level', -6) + mesh_plot = kwargs.get('mesh', False) + titles_flag = kwargs.get('titles', True) + colored_image = kwargs.get('colored', True) + + # Initialize title messages + titleEnvMsg = [] + titleProfMsg = [] + + try: + if 'idTitle' in recResult and recResult['idTitle']: + titleEnvMsg.append(recResult['idTitle']) + titleProfMsg.append(recResult['idTitle']) + except: + pass + + # Output post-process + if post_proc == 'env': + output = np.abs(np.hilbert(recResult['output'])) + envMsg = recResult['methodMsg'] + ' envelope image' + elif post_proc == 'rectified': + output = np.abs(recResult['output']) + envMsg = recResult['methodMsg'] + ' rectified image' + elif post_proc == 'none': + output = recResult['output'] + envMsg = recResult['methodMsg'] + ' image' + + # Convert to dB if required + if output_type == 'dB': + maxOutput = np.max(output) + output = 10 * np.log10(output / maxOutput) + titleEnvMsg.append(envMsg + ' (dB)') + titleEnvMsg.append(f"Range = {range_dB:.2f} dB") + imgLimits = [-range_dB, 0] + API_level = API_level + elif output_type == 'normalized': + maxOutput = np.max(output) + output = output / maxOutput + titleEnvMsg.append(envMsg + ' (normalized)') + imgLimits = [0, 1] + API_level = 10 ** (API_level / 20) + else: + titleEnvMsg.append(envMsg) + imgLimits = [np.min(output), np.max(output)] + API_level = 10 ** (API_level / 20) * np.max(np.abs(imgLimits)) + + # Calculate API (if applicable) + if API_calc: + dx = recResult['roiX'][1] - recResult['roiX'][0] + dz = recResult['roiZ'][1] - recResult['roiZ'][0] + if 'Wavelength' in recResult and recResult['Wavelength']: + lambda_value = recResult['Wavelength'] * 1e-3 + else: + lambda_value = 5900 / 5e6 + + API = np.sum(output >= API_level) * dx * dz / (lambda_value ** 2) + titleEnvMsg.append(f"API ({API_level:.0f} dB) = {API:.2f}") + + # Plot profile if requested + if profile_plot: + idxMinZ = np.where(recResult['roiZ'] > minZ)[0][0] - 1 + idxMaxZ = np.where(recResult['roiZ'] >= maxZ)[0][0] + profMsg = recResult['methodMsg'] + ' profile' + outProf = np.abs(np.hilbert(recResult['output'][idxMinZ:idxMaxZ, :])) + maxOutprof = np.max(outProf) + outProf_dB = 10 * np.log10(outProf / maxOutprof) + titleProfMsg.append(profMsg) + titleProfMsg.append(f"Range = {prof_range_dB:.2f} dB") + profLimits = [-prof_range_dB, 0] + + # Calculate resolution if requested + if res_calc_plot: + def fun(x): + return interp1d(recResult['roiX'], np.max(outProf_dB, axis=0), kind='slinear', fill_value="extrapolate")(x) + 3 + + idx0dB = np.where(np.max(outProf_dB, axis=0) == 0)[0] + idx3dB = np.where(np.max(outProf_dB, axis=0) > -3)[0] + resol = (np.sign(fun(recResult['roiX'][idx3dB][-1]) - fun(recResult['roiX'][idx0dB][0])) * (fun(recResult['roiX'][idx3dB][-1]) - fun(recResult['roiX'][idx0dB][0])) * 1000) + titleProfMsg.append(f"Resolution = {resol:.3f} mm") + + # Plot profile + if not animation_plot: + plt.figure() + plt.plot(recResult['roiX'] * 1000, np.max(outProf_dB, axis=0)) + plt.xlim([minX, maxX] * 1000) + plt.ylim(profLimits) + if titles_flag: + plt.title(titleProfMsg) + plt.xlabel('x (mm)') + plt.ylabel('Maximum Amplitude (dB)') + plt.grid(True) + + # Plot 2D image if requested + if not res_calc_plot: + if not animation_plot: + plt.figure() + if mesh_plot: + plt.pcolormesh(recResult['roiX'] * 1000, recResult['roiZ'] * 1000, output, shading='auto') + else: + plt.imshow(output, extent=[recResult['roiX'][0] * 1000, recResult['roiX'][-1] * 1000, recResult['roiZ'][0] * 1000, recResult['roiZ'][-1] * 1000], aspect='auto', cmap='jet') + + plt.xlim([minX, maxX] * 1000) + plt.ylim([minZ, maxZ] * 1000) + if colored_image: + plt.colorbar() + if titles_flag: + plt.title(titleEnvMsg) + plt.xlabel('x (mm)') + plt.ylabel('z (mm)') + plt.grid(True) + + plt.show() diff --git a/3D_app/utsr/SimulNDT.py b/3D_app/utsr/SimulNDT.py new file mode 100644 index 0000000..3092522 --- /dev/null +++ b/3D_app/utsr/SimulNDT.py @@ -0,0 +1,239 @@ +import numpy as np +from scipy.fft import ifft + +def simulNDT(**kwargs): + # Define default values and input validation + defaults = { + 'TsGate': 0, + 'TendGate': 50, + 'deltaT': 0.02, + 'minX': 0, + 'maxX': 60, + 'deltaX': 1, + 'cl': 5900, + 'cs': 3230, + 'Density': 7.86, + 'fc': 5, + 'bw': 5, + 'WaveType': 'L', + 'BeamAngle': 0, + 'd': 6, + 'DefectType': 'empty', + 'flPosX': 30, + 'flPosY': 0, + 'flDepthCentre': 40, + 'flDiameter': 0, + 'flHeight': 0, + 'flTilt': 0, + 'SNR': np.inf, + 'simul': True + } + + # Override defaults with user inputs + params = {**defaults, **kwargs} + + # Create time scale element + t = np.arange(params['TsGate'], params['TendGate'] + params['deltaT'], params['deltaT']) + nt = len(t) + idxTsGate = np.searchsorted(t, params['TsGate'], side='right') - 1 + idxTendGate = np.searchsorted(t, params['TendGate'], side='left') + + scanData = { + 'timeScale': t, + 'CscanData': { + 'AscanPoints': len(t), + 'TsGate': params['TsGate'], + 'TendGate': params['TendGate'], + 'X': np.arange(params['minX'], params['maxX'] + params['deltaX'], params['deltaX']), + 'Rows': int((params['maxX'] - params['minX']) / params['deltaX']) + 1, + 'Cl': params['cl'], + 'Cs': params['cs'], + 'Density': params['Density'], + 'Frequency': params['fc'], + 'Bandwidth': params['bw'], + 'WaveType': params['WaveType'], + 'BeamAngle': params['BeamAngle'], + 'ProbeDiameter': params['d'], + 'DefectType': params['DefectType'], + 'PosX': params['flPosX'], + 'PosY': params['flPosY'], + 'DepthCentre': params['flDepthCentre'], + 'Diameter': params['flDiameter'], + 'Height': params['flHeight'], + 'Tilt': params['flTilt'], + 'MaxSignal': 0, + 'MaxSdtRef': 0, + 'Y': 0, + 'Cols': 1, + 'CscanValues': np.zeros((int((params['maxX'] - params['minX']) / params['deltaX']) + 1, 1)) + }, + 'AscanValues': np.zeros((len(t), int((params['maxX'] - params['minX']) / params['deltaX']) + 1)) + } + + if scanData['CscanData']['WaveType'] in ['L', 'l']: + c = scanData['CscanData']['Cl'] + elif scanData['CscanData']['WaveType'] in ['S', 's']: + c = scanData['CscanData']['Cs'] + + scanData['CscanData']['Wavelength'] = c / (scanData['CscanData']['Frequency'] * 1000) + scanData['CscanData']['Nearfield'] = (scanData['CscanData']['ProbeDiameter'] / 2) ** 2 / scanData['CscanData']['Wavelength'] + scanData['CscanData']['ProbeEllipX'] = 0 + scanData['CscanData']['ProbeEllipY'] = 0 + scanData['CscanData']['TrueAngle'] = scanData['CscanData']['BeamAngle'] + scanData['CscanData']['CalValue'] = 0 + + # Create default simulation setup + setup = setup_maker() + + # Define the frequency window for simulation + f = np.fft.fftfreq(nt, d=params['deltaT']) + f = f[nt // 2:] + setup['f'] = f + + # Adjust material parameters + setup['matl'] = { + 'd2': scanData['CscanData']['Density'], + 'cp2': scanData['CscanData']['Cl'], + 'cs2': scanData['CscanData']['Cs'] + } + + # Adjust transducer parameters + setup['trans'] = {'d': scanData['CscanData']['ProbeDiameter']} + setup['system'] = { + 'sysf': 'cossqrtf', + 'amp': 0.5 / 2048, + 'fc': scanData['CscanData']['Frequency'], + 'bw': scanData['CscanData']['Bandwidth'] + } + setup['geom'] = {'i_ang': scanData['CscanData']['BeamAngle']} + setup['type2'] = 'p' if scanData['CscanData']['WaveType'] in ['L', 'l'] else 's' + + # Adjust flaw parameters + flaw_diameter = scanData['CscanData']['Diameter'] / 2 + if scanData['CscanData']['DefectType'] == 'spherical': + setup['flaw'] = {'b': flaw_diameter, 'Afunc': 'A_void'} + elif scanData['CscanData']['DefectType'] == 'side-drilled': + setup['flaw'] = {'b': flaw_diameter, 'Afunc': 'A_SDH'} + elif scanData['CscanData']['DefectType'] == 'circ-crack': + setup['flaw'] = {'b': flaw_diameter, 'f_ang': scanData['CscanData']['Tilt'], 'Afunc': 'A_crack'} + else: + setup['flaw'] = {'b': flaw_diameter, 'Afunc': 'empty'} + + # Adjust flaw position + setup['geom']['z2'] = scanData['CscanData']['DepthCentre'] + + # Perform the simulation + if params['simul']: + scanX = scanData['CscanData']['X'] + nx = len(scanX) + + for ix in range(nx): + # Set up the distance between flaw and transducer + setup['geom']['x2'] = scanData['CscanData']['PosX'] - scanX[ix] + + # Calculate the transducer output signal (in frequency domain) + if setup['flaw']['Afunc'] == 'A_SDH': + setup['geom']['y2'] = np.arange(-50, 51) + scanData['CscanData']['PosY'] + Vf, setup = SDH_PE_MM(setup) + else: + setup['geom']['y2'] = scanData['CscanData']['PosY'] + Vf, setup = TG_PE_MM(setup) + + # Calculate the signal delay offset (tau0) + tau0 = 2000 * (setup['geom']['z1'] / setup['wave']['c1'] + setup['geom']['z2'] / setup['wave']['c2']) + + # Include the signal delay offset (tau0) + Vf *= np.exp(2j * np.pi * setup['f'] * tau0) + + # Convert the signal to time domain and store it + Vfe = np.concatenate([Vf, np.zeros(len(f) - len(Vf))]) + Vfe[0] /= 2 + vt = 2 * np.real(ifft(Vfe)) + scanData['AscanValues'][:, ix] = vt[idxTsGate:idxTendGate] + + # Normalize A-scan amplitude + MaxSignal = np.max(np.abs(scanData['AscanValues'])) + scanData['CscanData']['MaxSignal'] = MaxSignal + scanData['CscanData']['CscanValues'] = 10 * np.log10(np.max(np.abs(scanData['AscanValues']), axis=0) / MaxSignal) + + # Add noise into signals + if params['SNR'] != np.inf: + G = 10 ** (params['SNR'] / 10) + r = np.random.randn(*scanData['AscanValues'].shape) + k = np.sqrt(np.var(scanData['AscanValues']) / (np.var(r) * G)) + scanData['AscanValues'] += k * r + + return scanData, setup + +def setup_maker(): + # Initialize an empty dictionary (equivalent to MATLAB structure) + setup = {} + + # Setup parameters + setup['f'] = 5 + setup['type1'] = 'p' + setup['type2'] = 'p' + + # System function parameters + setup['system'] = {} + setup['system']['sysf'] = 'systf' + setup['system']['amp'] = 5.0E-02 + setup['system']['fc'] = 5 + setup['system']['bw'] = 3 + setup['system']['z1r'] = 0 + setup['system']['en'] = 0.01 + setup['system']['ref_file'] = 'empty' + + # Transducer parameters + setup['trans'] = {} + setup['trans']['d'] = 12.7 + setup['trans']['fl'] = float('inf') + + # Geometry parameters + setup['geom'] = {} + setup['geom']['z1'] = 0 + setup['geom']['z2'] = list(np.linspace(0, 200, 512)) # Using numpy for linspace + setup['geom']['x2'] = 0 + setup['geom']['y2'] = 0 + setup['geom']['i_ang'] = 0 + setup['geom']['R1'] = float('inf') + setup['geom']['R2'] = float('inf') + setup['geom']['p_ang'] = 0 + + # Material parameters + setup['matl'] = {} + setup['matl']['d1'] = 1.0 + setup['matl']['d2'] = 1.0 + setup['matl']['cp1'] = 1480 + setup['matl']['cs1'] = 0 + setup['matl']['cp2'] = 1480 + setup['matl']['cs2'] = 0 + setup['matl']['p1'] = [0] * 5 + setup['matl']['s1'] = [0] * 5 + setup['matl']['p2'] = [0] * 5 + setup['matl']['s2'] = [0] * 5 + + # Flaw parameters + setup['flaw'] = {} + setup['flaw']['b'] = 0.0 + setup['flaw']['f_ang'] = 0.0 + setup['flaw']['Afunc'] = 'empty' + + # Wave parameters + setup['wave'] = {} + setup['wave']['c1'] = 1480 + setup['wave']['c2'] = 1480 + setup['wave']['T12'] = 1.0 + + return setup + +def SDH_PE_MM(setup): + # Dummy function to emulate the SDH_PE_MM functionality + # Replace with the actual implementation + return np.zeros(len(setup['f'])), setup + +def TG_PE_MM(setup): + # Dummy function to emulate the TG_PE_MM functionality + # Replace with the actual implementation + return np.zeros(len(setup['f'])), setup + diff --git a/3D_app/utsr/UTSRModS2.py b/3D_app/utsr/UTSRModS2.py new file mode 100644 index 0000000..61f1c85 --- /dev/null +++ b/3D_app/utsr/UTSRModS2.py @@ -0,0 +1,260 @@ +import numpy as np +from scipy.sparse.linalg import cg +from scipy.signal import hilbert +from utsr.ModS2 import mod_s2 +from utsr.ModS2T import mod_s2t +from utsr.GenerateModelFilter import generate_model_filter + +def model(x, W, dt, dx, ERMv, alpha, szData, tau0, Mod, Filt, Rx): + # Applies direct operator + g = mod_s2(x.reshape(szData), dt, dx, ERMv, tau0, Mod, Rx) + + # Applies adjoint operator + y = mod_s2t(g, dt, dx, ERMv, tau0, Filt, Rx) + + # Applies regularization + if alpha != 0: + # L1 norm if W != I and L = I, Tikhonov if W = I and L = I. + y = y.flatten() + alpha**2 * W * x + else: + # No regularization, least squares solution + y = y.flatten() + + return y + +def utsr_mod_s2(scanData, **kwargs): + # Parser verification functions + def scanDataValidationFcn(x): + return isinstance(x, dict) and all(key in x for key in ['CscanData', 'timeScale', 'AscanValues']) + + validOutput = ['raw', 'normalized', 'dB'] + validPostproc = ['none', 'rectified', 'env'] + validPreproc = ['none', 'env'] + validTrFreqSp = ['gaussian', 'cossquare', '610060'] + + # Verify scanData parameter + assert scanDataValidationFcn(scanData), 'scanData invalid structure' + + # Get scan data informations + c = scanData['CscanData']['Cl'] if scanData['CscanData']['WaveType'] == 'L' else scanData['CscanData']['Cs'] + ERMv = c / 2 + + scanX = scanData['CscanData']['X'] * 1e-3 + scanY = scanData['CscanData']['Y'] * 1e-3 + tamX = len(scanX) + tamY = len(scanY) + deltaScanX = scanX[1] - scanX[0] + + scanZ = scanData['timeScale'] * 1e-6 * ERMv + deltaScanZ = scanZ[1] - scanZ[0] + + apertAngle = np.arcsin(0.51 * c / (scanData['CscanData']['Frequency'] * 1e6 * scanData['CscanData']['ProbeDiameter'] * 1e-3)) + tanApertAngle = np.tan(apertAngle) + + # Process arguments + params = { + 'minX': scanX[0], + 'maxX': scanX[-1], + 'deltaX': deltaScanX, + 'samePitch': False, + 'minZ': scanZ[0], + 'maxZ': scanZ[-1], + 'plotOut': True, + 'debug': True, + 'Yslice': np.argwhere(scanY == 0).flatten()[0] if np.any(scanY == 0) else 0, + 'output': validOutput[0], + 'range': 25, + 'Preproc': validPreproc[0], + 'Postproc': validPostproc[0], + 'profile': True, + 'prof_range': 50, + 'alpha': 1.0, + 'beta': np.finfo(float).eps * 100, + 'tol': 1e-2, + 'ep': 3.397e-4, + 'cg_tol': 1e-2, + 'maxStagCount': 6, + 'animation': False, + 'res_calc': False, + 'TrFreqSp': validTrFreqSp[0], + 'SysModel': np.zeros((0, 0)), + 'title': '', + 'Rx': 1, + 'NegCutOff': True, + 'fOrig': np.zeros((0, 0)), + 'errOrig': -np.inf, + 'maxIter': 0 + } + + params.update(kwargs) + + Rx = int(params['Rx']) + assert Rx > 0, 'Rx must be positive' + + # Define ROI + if params['samePitch']: + deltaRoiX = deltaScanZ + else: + deltaRoiX = params['deltaX'] + + if deltaScanX == params['deltaX']: + idxMinX = np.argmax(scanX > params['minX']) - 1 + idxMaxX = np.argmax(scanX >= params['maxX']) + if (idxMinX == 0) and (idxMaxX == len(scanX)): + roiX = scanX + else: + roiX = scanX[idxMinX:idxMaxX] + else: + roiX = np.arange(params['minX'], params['maxX'], deltaRoiX) + N = len(roiX) + + idxMinZ = np.argmax(scanZ > params['minZ']) - 1 + idxMaxZ = np.argmax(scanZ >= params['maxZ']) + roiZ = scanZ[idxMinZ:idxMaxZ] + M = len(roiZ) + + if tamY > 1: + assert 0 < params['Yslice'] <= tamY, 'Yslice parameter out of range' + yi = params['Yslice'] + else: + yi = 1 + + limApertX = roiZ[-1] * tanApertAngle + colIniX = np.argmax(scanX + limApertX >= roiX[0]) + colFinX = np.argmax(scanX - limApertX > roiX[-1]) - 1 + if colFinX == -1: + colFinX = len(scanX) + + if params['debug']: + print('\n$$$$$$$ Starting UTSR algorithm (alpha = {}) $$$$$$$'.format(params['alpha'])) + print('---- Title: {} ----'.format(params['title'])) + + ci = ((yi - 1) * tamX) + S_t_u = scanData['AscanValues'][idxMinZ:idxMaxZ, colIniX:colFinX + ci] + if scanData['AscanValues'].shape == params['fOrig'].shape: + fOrig = params['fOrig'][idxMinZ:idxMaxZ, colIniX:colFinX + ci] + else: + fOrig = np.zeros((0, 0)) + + tau0 = scanData['timeScale'][idxMinZ] * 1e-6 + if params['Preproc'] == 'env': + S_t_u = np.abs(hilbert(S_t_u)) + + szData = [M, N] + dt = deltaScanZ / ERMv + dx = deltaRoiX + + scanData['AscanValues'] = S_t_u + scanData['CscanData']['X'] = scanData['CscanData']['X'][colIniX:colFinX + ci] + scanData['timeScale'] = scanData['timeScale'][idxMinZ:idxMaxZ] + + if params['SysModel'].size == 0: + Mod, _ = generate_model_filter(scanData, TrFreqSp=params['TrFreqSp'], Rx=Rx) + else: + Mod = params['SysModel'] + + Filt = np.conj(Mod) + alpha = params['alpha'] + beta = 1 + tol = params['tol'] + cg_tol = params['cg_tol'] + + HTg = mod_s2t(S_t_u, dt, dx, ERMv, tau0, Filt, Rx) + fTemp = np.zeros(HTg.size) + numIter = 1 + stagCount = 0 + maxIter = params['maxIter'] + + beta_iter = [] + cg_tol_iter = [] + cg_iter = [] + tL2 = [] + tL1 = [] + J_f = [] + errUTSR = [] + difUTSR = [] + errOrig = [] + + while True: + if numIter == 1: + W = np.ones(fTemp.shape) + else: + if not np.isinf(beta): + while True: + W = 1 / (np.abs(fTemp) + beta) + n1fTemp = np.linalg.norm(fTemp, 1) + n1fTempAppr = np.linalg.norm(np.sqrt(W) * fTemp) ** 2 + n1Err = (n1fTemp - n1fTempAppr) / n1fTemp + if abs(n1Err - cg_tol) > cg_tol / 10: + beta = beta / (n1Err / cg_tol) + else: + break + + beta_iter.append(beta) + if params['debug']: + print('cg_tol = {} | '.format(cg_tol)) + + cg_tol_iter.append(cg_tol) + f, flag = cg( + lambda x: model(x, W, dt, dx, ERMv, alpha, szData, tau0, Mod, Filt, Rx), + HTg.flatten(), tol=cg_tol, maxiter=M * N, x0=fTemp) + cg_iter.append(flag) + + if params['NegCutOff']: + f[f < 0] = 0 + + errTemp = np.linalg.norm(f.flatten() - fTemp) / np.linalg.norm(fTemp.flatten()) + + if params['debug']: + print('iter = {} | tol = {} | rel. tol. = {} | err = {}% | stag. count = {}'.format( + numIter, tol, cg_tol, errTemp * 100, stagCount)) + + if errTemp < tol: + break + + fTemp = f + + if params['res_calc']: + aux1 = mod_s2(f.reshape(szData), dt, dx, ERMv, tau0, Mod, Rx) + aux2 = S_t_u - aux1 + errUTSR.append(np.linalg.norm(aux2.flatten()) / np.linalg.norm(S_t_u.flatten())) + J_f.append(0.5 * np.linalg.norm(aux2.flatten()) ** 2 + alpha * np.linalg.norm(f.flatten(), 1)) + if fOrig.size: + difUTSR.append(np.linalg.norm(fOrig.flatten() - f.flatten()) / np.linalg.norm(fOrig.flatten())) + errOrig.append(params['errOrig']) + + if maxIter and numIter >= maxIter: + break + + if stagCount > params['maxStagCount']: + break + + if errTemp < cg_tol: + stagCount += 1 + else: + stagCount = 0 + + cg_tol /= 10 + numIter += 1 + + if params['NegCutOff']: + f[f < 0] = 0 + + if params['Postproc'] == 'rectified': + f = np.abs(f) + elif params['Postproc'] == 'env': + f = np.abs(hilbert(f)) + + if params['plotOut']: + import matplotlib.pyplot as plt + plt.imshow(f, cmap='hot', aspect='auto') + plt.colorbar() + plt.title('Reconstructed Image') + plt.show() + + return { + 'f': f, + 'errUTSR': errUTSR, + 'difUTSR': difUTSR, + 'errOrig': errOrig + }