This repository has been archived by the owner on Dec 16, 2019. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 10
/
im_analyze5.m
232 lines (192 loc) · 6.98 KB
/
im_analyze5.m
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
function im_analyze5()
% im_analyze4.m
% clear
% -------------- Program Variables --------------------------------%
im_dir = 'Trial 32_8\'; % Folder where images or video are stored
% im_dir = 'Nov 6\';
% im_list = 'All'; % Images or frames to read. 1st image sets reference plane
im_list = [10:30];
plots_on = 2; % Generate diagnostic plots
% ----------- Experimental Parameters -----------------------------%
L_0 = (28+4)*25.4; % Distance from object to camera & projector: mm
D = 4.75*25.4; % Distance from camera to projector: mm
% f_0 = 13.5/25.4; % Fundamental frequency of observed grating, gratings/mm
% pixel_pitch = 190/25.4; % pixels / mm
pixel_pitch = 6.59; % Trial 30 -- pixels / mm
% win_length = 40; % = 1/standard deviation: pixels^-1
% -------- Data analysis variables -----------------------------%
ROI = [200, 90, 1050, 990]; % X Start (1st column), Y start (1st row),
% % X End (last column), Y end (last row)
% ROI = 'off';
rotate_on = 0; % 0 if gratings are horizontal (default), 1 rotates images
high_pass_wid = 80; % Width of high pass gaussian filter
gauss_filter_ratio = 1/3; % Ratio of standard deviation to wavenumber
% ------------ Read images ------------------------------%
files = dir(im_dir);
if strcmpi(files(3).name(end-2:end),'avi')
if ischar(im_list)
im_list = 1:length(files); % Create a list of all frames in video to import
end
num_images = length(im_list);
readerobj = mmreader(strcat(im_dir, files(3).name));
im.name = readerobj.Name;
for i = 1:num_images
temp = read(readerobj,im_list(i)); % Read frame
im(i).image = temp(:,:,1); % for B&W images only use R of RGB
end
else
if ischar(im_list)
im_list = 1:(length(files)-2); % Create a list of all images in directory
end
num_images = length(im_list);
for i = 1:num_images
im(i).image = imread(strcat(im_dir, files(im_list(i)+2).name));
im(i).name = files(im_list(i)+2).name;
end
end
% ------------- Region of Interest ------------------------%
if plots_on >= 1
subplot(1,2,1)
imshow(im(1).image, [0 mean(max(im(1).image))])
title('Reference image')
end
if ~ischar(ROI(1))
if plots_on
hold on
plot([ROI(1),ROI(3),ROI(3),ROI(1),ROI(1)],[ROI(2),ROI(2),...
ROI(4),ROI(4),ROI(2)], 'Color','red', 'LineWidth',1) % Rectangle over ROI
hold off
end
% Discard subset of image not in ROI
for i = 1:num_images
im(i).image = im(i).image(ROI(2):ROI(4),ROI(1):ROI(3));
end
end
if rotate_on
for i = 1:num_images
im(i).image = im(i).image';
end
end
% ----------- Analysis & processing ------------------------%
[N c] = size(im(1).image); % N is the length of the column vectors = fft length
% NFFT = 2^nextpow2(N); % Next power of 2
mid_x = floor(c/2);
mid_y = floor(N/2);
x_p = 0:c-1; % X Pixels
y_p = 0:N-1; % Y Pixels
y_range = floor(.45*N:.55*N);
mich_contrast = double(max(im(1).image(y_range,mid_x)) - min(im(1).image(y_range,mid_x)))...
/ double(max(im(1).image(y_range,mid_x)) + min(im(1).image(y_range,mid_x)));
if plots_on >= 1
subplot(1,2,1)
hold on
plot([ROI(1)+mid_x,ROI(1)+mid_x],[ROI(2)+y_range(1),ROI(2)+y_range(end)],...
'Color','blue', 'LineWidth',1);
hold off
subplot(1,2,2);
plot(y_range,im(1).image(y_range,mid_x))
title(sprintf('Michelson Contrast = %1.2f', mich_contrast))
xlabel('Y (pixels)')
ylabel('Intensity')
% xlim([1 2])
end
% % Remove average
% im_1 = im_1 - mean(mean(im_1));
% im_2 = im_2 - mean(mean(im_2));
% ---------- Take Fourier Transform -----------
filt = fspecial('disk',3); % Initialize averaging filter
for i = 1:num_images
% im(i).image = imfilter(im(i).image, filt); % Filter image first
im(i).fft = fft(imfilter(double(im(i).image),filt));
% im(i).image = fft(im(i).image); % Fourier transform columns
end
% --------- Plot power spectrum ---------------
f0 = (-N/2:N/2-1);
% % Plot single-sided amplitude spectrum.
% f = 1/2*linspace(0,1,NFFT/2+1);
% plot(f,2*abs(im(1).fft(1:NFFT/2+1,mid_x)))
% title('Single-Sided Amplitude Spectrum of y(t)')
% xlabel('Frequency (Hz)')
% ylabel('|Y(f)|')
if plots_on >= 3
% 2D Plot
mesh(x_p,f0,fftshift(abs(im(1).fft)))
end
% ---------- Window (filter) data ----------
% High-pass filter (remove lowest frequencies)
win = 1-gausswin(N,high_pass_wid);
if plots_on >= 2
plot(fftshift(abs(im(1).fft(:,mid_x))))
hold on
plot(win*max(abs(im(1).fft(:,mid_x))))
hold off
end
for i = 1:num_images
im(i).fft = ifftshift(fftshift(im(i).fft).*repmat(win,1,c));
end
% im_1 = ifftshift(fftshift(im_1).*repmat(win,1,c));
% im_2 = ifftshift(fftshift(im_2).*repmat(win,1,c));
% Identify where to filter using reference image (this happens once)
[C I] = max(abs(im(1).fft(:,mid_x)));
% shift_length = I;
shift_length = I-round(N/2);
win_length = (I-1)*gauss_filter_ratio; % Standard deviation of window
win = circshift(gausswin(N,(N/2)/(win_length)),shift_length); % Generate and shift gaussian window
win = win/max(win);
if plots_on >= 1
plot(abs(im(1).fft(:,mid_x))/N)
hold on
plot(win*max(abs(im(1).fft(:,mid_x)))/N)
title(sprintf('Gaussian Filter with Sigma = %1.2f', win_length))
xlim([0 3*I])
hold off
end
% Apply window to fft data
for i = 1:num_images
im(i).fft = im(i).fft.*repmat(win,1,c);
end
% ----- Transform back to real space -------
filt = fspecial('disk',5); % Initialize averaging filter
for i = 1:num_images
im(i).phase = ifft(im(i).fft);
end
f_0 = (I/N)*pixel_pitch; % Gratings / pixel * pixel/ mm
% ---- Calculate delta_phi (change of phase with respect to reference)---
for i = 2:num_images
im(i).phase = angle(im(i).phase.*conj(im(1).phase));
im(i).phase = imfilter(im(i).phase, filt); % Image filter
im(i).phase = unwrap(im(i).phase); % Remove jumps along columns
im(i).phase = unwrap(im(i).phase')'; % Check for jumps along rows too
im(i).h = L_0 * im(i).phase./(im(i).phase-2*pi*f_0*D); % Compute height from phase
if rotate_on==1
im(i).h = im(i).h'; % flip back to original orientation
end
figure(1)
mesh(x_p/pixel_pitch, y_p/pixel_pitch, im(i).h)
title(sprintf('Image %g of %g', i-1, num_images-1))
xlabel('X (mm)')
ylabel('Y (mm)')
zlabel('Z (mm)')
zlim([-4 4])
figure(2)
subplot(2,2,1)
plot(f0,fftshift(abs(im(i).fft(:,mid_x))))
title('Filtered Frequency')
subplot(2,2,2)
plot(y_p,im(i).phase(:,mid_x))
title('Unwrapped Phase')
ylabel('Radians')
subplot(2,2,3)
plot(y_p/pixel_pitch,im(i).h(:,mid_x))
title('Middle X Slice')
xlabel('Y (mm)')
subplot(2,2,4)
plot(x_p/pixel_pitch,im(i).h(mid_y,:))
title('Middle Y')
ylabel('Height (mm)')
xlabel('X (mm)')
end
% zlim([0 .3])
% figure
% plot(y_p/pixel_pitch,im.h(mid_x,:))
end