1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
|
# -*- coding: utf-8 -*-
"""
Created on Sat Sep 15 12:21:45 2018
@author: gcg
1d grid method for determining pixel displacements.
Version 0
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage.interpolation import map_coordinates
from scipy import signal
from scipy.interpolate import interp1d
class grid_method_1d():
def __init__(self, X0, Y0, ROI, carrier_frequency=None):
self.X0 = X0
self.Y0 = Y0
#print "ROI = ", ROI
self.ROI = ROI # ROI is only used for carrier frequency identification and unwrapping
self.lambda0 = self.freq_from_fft()
self.w0 = 2 * np.pi / self.lambda0
self.treatBorder = True
print("measured sample wavelength is: ", self.lambda0)
self.phi0, self.border = self.compute_phases_LSA(Y0)
def run_analysis(self, Y1):
self.Y1 = Y1
self.phi1, self.border = self.compute_phases_LSA(Y1)
#self.unwrap_by_similarity()
self.unwrap_by_constrained_displacement()
# only take care of the phases for now
# the idea is to perform temporal unwrapping later on
#self.phi0[:self.ROI[0]] = self.phi0[self.ROI[0]]
#self.phi0[self.ROI[1]:] = self.phi0[self.ROI[1]]
#self.phi1[:self.ROI[0]] = self.phi1[self.ROI[0]]
#self.phi1[self.ROI[1]:] = self.phi1[self.ROI[1]]
self.u = self.compute_displacements()
self.strain = self.compute_strain()
self.Y1_pulled_back = self.shift_to_reference_frame(self.Y1, self.u)
def unwrap_by_similarity(self):
# compute native displacement directly from phases
u = -( self.phi1 - self.phi0 ) * self.lambda0/ (2*np.pi)
# using the native displacement, shift Y1 onto Y0 and compute autocorrelation
Y1_pulled_back = self.shift_to_reference_frame(self.Y1, u)
score0 = np.dot(Y1_pulled_back[self.ROI[0]:self.ROI[1]], self.Y0[self.ROI[0]:self.ROI[1]])
# add 2 Pi and repeat
u = -( self.phi1 + 2*np.pi - self.phi0 ) * self.lambda0/ (2*np.pi)
Y1_pulled_back = self.shift_to_reference_frame(self.Y1, u)
score_plus = np.dot(Y1_pulled_back[self.ROI[0]:self.ROI[1]], self.Y0[self.ROI[0]:self.ROI[1]])
# subtract 2 Pi and repeat
u = -( self.phi1 - 2*np.pi - self.phi0 ) * self.lambda0/ (2*np.pi)
Y1_pulled_back = self.shift_to_reference_frame(self.Y1, u)
score_minus = np.dot(Y1_pulled_back[self.ROI[0]:self.ROI[1]], self.Y0[self.ROI[0]:self.ROI[1]])
scores = [score0, score_plus, score_minus]
shifts = [0, 2*np.pi, -2*np.pi]
idx = np.argmax(scores)
shift = shifts[idx]
print( "scores", scores - np.min(scores))
print( "adding %f to phase 1" % (shift))
self.phi1 += shift
def unwrap_by_constrained_displacement(self):
""" assume that the displacement must be small compared to wavelength """
self.phi0[:self.ROI[0]] = self.phi0[self.ROI[0]]
self.phi0[self.ROI[1]:] = self.phi0[self.ROI[1]]
self.phi1[:self.ROI[0]] = self.phi1[self.ROI[0]]
self.phi1[self.ROI[1]:] = self.phi1[self.ROI[1]]
u = -( self.phi1 - self.phi0 ) * self.lambda0/ (2*np.pi)
u_avg = np.mean(u[self.ROI[0]:self.ROI[1]])
if u_avg > 0.5 * self.lambda0:
#u -= this.lambda0
self.phi1 += 2.0 * np.pi
elif u_avg < -0.5 * self.lambda0:
#u += this.lambda0
self.phi1 -= 2.0 * np.pi
#average_displacement = np.mean(u[p["ROI start"]:p["ROI stop"]])
# print "****************************************"
# print "... detected displacement exceeding carrier wavelength=%f" % (self.lambda0)
# print "u extreme = %f" % (u_extreme)
# print "proposed shift: %f" % (shift * self.lambda0/ (2*np.pi))
# print "****************************************"
#self.phi1 += shift
def compute_strain(self):
#N = 5
#window = np.ones(N) / (1.0 * N)
window = signal.gaussian(np.ceil(8*self.lambda0), std=self.lambda0) # Gaussian window
norm = np.sqrt(2*np.pi*self.lambda0**2)
window /= norm
convolution = signal.fftconvolve(self.u, window, mode="same")
strain = np.gradient(convolution)
return strain
def shift_to_reference_frame(self, sig, shift):
""" shift the input signal to back to the new frame of reference
which is given by original coordinates + shift """
# make sure the new frame of reference does not extend beyond
# the boundaries of the original coordinates becuase otherwise the
# interpolation will fail.
X1 = self.X0 + shift
X1 = np.where(X1 > self.X0[-1], self.X0[-1], X1)
X1 = np.where(X1 < self.X0[0], self.X0[0], X1)
f = interp1d(self.X0, sig)
Y1 = f(X1)
#Y1 = map_coordinates(sig, [self.X0+shift],order=3) # shift u into the initial coordinate frame
return Y1
def compute_displacements(self):
u = -( self.phi1 - self.phi0 ) * self.lambda0/ (2*np.pi)
MaxIter = 100
for n in np.arange(MaxIter):
interpPHI1 = self.shift_to_reference_frame(self.phi1,u)
newu = -( interpPHI1 - self.phi0 ) * self.lambda0/ (2*np.pi)
diff = np.abs(u - newu).sum()
#print("iteration %d, residual is %g" % (n, diff))
u = newu
if diff < 1.0e-3: break
return u
def compute_phases_LSA(self, sig):
""" Local Spectrum Analysis (LSA) of the image phases by convolving the
product of input signal and carrier frequency with a filtering window of compact support
(hence, local analysis). This is the implementaion of Eq. 4.30 in:
M. Grédiac, F. Sur, and B. Blaysat. The grid method for in-plane
displacement and strain measurement: a review and analysis. Strain, 2016.
"""
window = signal.gaussian(np.ceil(8*self.lambda0), std=self.lambda0) # Gaussian window
#window = signal.triang(np.ceil(4*self.lambda0))
kernel = sig*np.exp(-1j * self.w0 * self.X0)
convolution = signal.fftconvolve(kernel, window, mode="same")
phases = np.angle(convolution)
# take care of regions where the autocorrelation is not well defined
# because window and kernel do not fully overlap
if self.treatBorder:
border = len(window)
phaseLeft = phases[border]
phaseRight = phases[-border]
phases[0:border] = phaseLeft
phases[-border:] = phaseRight
self.valid = np.ones(len(self.X0))
self.valid[0:border] = 0
self.valid[-border:] = 0
else:
self.valid = np.ones(len(self.X0))
border = 0
return np.unwrap(phases), border
#return phases, border
def freq_from_fft(self):
"""
Estimate frequency from peak of FFT
"""
s = signal.detrend(self.Y0[self.ROI[0]:self.ROI[1]])
# used Blackman window as the signal is non-periodic and do FFT
black= np.blackman(len(s))
Y = np.fft.fft(s*black)
N = len(Y)//2 + 1
X = np.linspace(0, 0.5, N, endpoint=True)
self.periods = 1.0/(X[1:])
self.power = np.abs(Y[1:N])
maxIdx = np.argmax(self.power)
maxPeriod = self.periods[maxIdx]
return maxPeriod
if __name__ == "__main__":
# generate reference 1d signal
carrier_wavelength = 5
w0 = 2.0 * np.pi / (carrier_wavelength)
print( "carrier_frequency is", w0)
X0 = np.arange(1024)
Y0 = 10*np.cos(X0 * w0)
#Y0 = signal.square(X0 * w0)
kernel = signal.gaussian(2*carrier_wavelength, std=4)
blurred = signal.fftconvolve(Y0, kernel, mode='same')
#Y0 = blurred
# gaussian displacement field
center = 512
width = 32
deformation = 0.095 * np.exp(-(X0-center)**2/(2.0*(width**2))) # check if this integrates to unity
start = 256
length = 512
deformation = np.zeros(len(X0))
for i in range(len(X0)):
deformation[i] = 0.0
if start < i < start + length:
deformation[i] = (i - start) * 0.002
# new approach for displacement field: scale X0
from scipy.integrate import cumtrapz
u_prescribed = cumtrapz(deformation, initial = 0.0)
u_max = np.max(u_prescribed)
print("maximum prescribed displacement is: ", u_max, np.argmax(u_prescribed))
X1 = X0 + u_prescribed
f = interp1d(X1, Y0)
Y1 = f(X0)
# run grid method
this = grid_method_1d(X0, Y0, Y1, carrier_wavelength)
# plot results
plt.subplot(3, 1, 1)
plt.title('1D grid method test')
plt.ylabel('image amplitude')
plt.plot(this.X0, this.Y0, label="Y0")
plt.plot(this.X0, this.Y1, label="Y1")
plt.legend()
plt.subplot(3, 1, 2)
plt.ylabel('displacements')
plt.plot(u_prescribed, label="prescribed displacement")
plt.plot(this.X0, this.u, label="measured displacement")
plt.legend()
plt.subplot(3, 1, 3)
plt.ylabel('strains')
strain_prescribed = np.gradient(u_prescribed, X0)
max_strain = np.max(np.max(np.abs(strain_prescribed)))
strain_error = (this.strain - (strain_prescribed)) / max_strain
plt.plot(strain_prescribed, label = "prescribed strain")
plt.plot(this.X0, this.strain, "r--", label = "measured strain")
plt.ylim(-max_strain,max_strain)
plt.plot(X0, strain_error, label="relative strain err")
plt.legend()
plt.show()
|