feat: Add UTSR module for 3D app
chore: Update UTSR module with necessary functions for simulation, data processing, and visualization.
This commit is contained in:
58
3D_app/UTSR.py
Normal file
58
3D_app/UTSR.py
Normal file
@ -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)
|
||||
122
3D_app/utsr/DatNDTRead.py
Normal file
122
3D_app/utsr/DatNDTRead.py
Normal file
@ -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
|
||||
171
3D_app/utsr/GenerateModelFilter.py
Normal file
171
3D_app/utsr/GenerateModelFilter.py
Normal file
@ -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)
|
||||
77
3D_app/utsr/ModS2.py
Normal file
77
3D_app/utsr/ModS2.py
Normal file
@ -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)
|
||||
89
3D_app/utsr/ModS2T.py
Normal file
89
3D_app/utsr/ModS2T.py
Normal file
@ -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)
|
||||
137
3D_app/utsr/PlotResults.py
Normal file
137
3D_app/utsr/PlotResults.py
Normal file
@ -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()
|
||||
239
3D_app/utsr/SimulNDT.py
Normal file
239
3D_app/utsr/SimulNDT.py
Normal file
@ -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
|
||||
|
||||
260
3D_app/utsr/UTSRModS2.py
Normal file
260
3D_app/utsr/UTSRModS2.py
Normal file
@ -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
|
||||
}
|
||||
Reference in New Issue
Block a user