我正试图实现本文中描述的算法:
We recorded a sequence of N successive speckle images with a sampling
frequency fs. In this way it was possible to observe how a pixel
evolves through the N images. That evolution can be treated as a time
series and can be processed in the following way: Each signal
corresponding to the evolution of every pixel was used as input to a
bank of filters. The intensity values were previously divided by their
temporal mean value to minimize local differences in reflectivity or
illumination of the object. The maximum frequency that can be
adequately analyzed is determined by the sampling theorem and s half
of sampling frequency fs. The latter is set by the CCD camera, the
size of the image, and the frame grabber. The bank of filters is
outlined in Fig. 1.
In our case, ten 5° order Butterworth filters
were used, but this number can be varied according to the required
discrimination. The bank was implemented in a computer using MATLAB
software. We chose the Butter-worth filter because, in addition to its
simplicity, it is maximally flat. Other filters, an infinite impulse
response, or a finite impulse response could be used.
By means of this
bank of filters, ten corresponding signals of each filter of each
temporary pixel evolution were obtained as output. Average energy Eb
in each signal was then calculated:
where pb(n) is the intensity of the filtered pixel in the nth image
for filter b divided by its mean value and N is the total number of
images. In this way, En values of energy for each pixel were obtained,
each of hem belonging to one of the frequency bands in Fig. 1.
With these values it is possible to build ten images of the active object,
each one of which shows how much energy of time-varying speckle there
is in a certain frequency band. False color assignment to the gray
levels in the results would help in discrimination.
这里是我的MATLAB代码库:
for i=1:520
for j=1:368
ts = [];
for k=1:600
ts = [ts D{k}(i,j)]; %%% kth image pixel i,j --- ts is time series
end
ts = double(ts);
temp = mean(ts);
if (temp==0)
for l=1:10
filtImag1{l}(i,j)=0;
end
continue;
end
ts = ts-temp;
ts = ts/temp;
N = 5; % filter order
W = [0.0 0.10;0.10 0.20;0.20 0.30;0.30 0.40;0.40 0.50;0.50 0.60 ;0.60 0.70;0.70 0.80 ;0.80 0.90;0.90 1.0];
[B,A]=butter(N,0.10,'low');
ts_f(1,:) = filter(B,A,ts);
N1 = 5;
for ind = 2:9
Wn = W(ind,:);
[B,A] = butter(N1,Wn);
ts_f(ind,:) = filter(B,A,ts);
end
[B,A]=butter(N,0.90,'high');
ts_f(10,:) = filter(B,A,ts);
for ind=1:10
%Following Paper Suggestion
filtImag1{ind}(i,j) =sum(ts_f(ind,:).^2);
end
end
end
for i=1:10
figure,imshow(filtImag1{i});
colorbar
end
pre_max = max(filtImag1{1}(:));
for i=1:10
new_max = max(filtImag1{i}(:));
if (pre_max
pre_max=max(filtImag1{i}(:));
end
end
new_max = pre_max;
pre_min = min(filtImag1{1}(:));
for i=1:10
new_min = min(filtImag1{i}(:));
if (pre_min>new_min)
pre_min = min(filtImag1{i}(:));
end
end
new_min = pre_min;
%normalize
for i=1:10
temp_imag = filtImag1{i}(:,:);
x=isnan(temp_imag);
temp_imag(x)=0;
t_max = max(max(temp_imag));
t_min = min(min(temp_imag));
temp_imag = (double(temp_imag-t_min)).*((double(new_max)-double(new_min))/double(t_max-t_min))+(double(new_min));
%median filter
%temp_imag = medfilt2(temp_imag);
imag_test2{i}(:,:) = temp_imag;
end
for i=1:10
figure,imshow(imag_test2{i});
colorbar
end
for i=1:10
A=imag_test2{i}(:,:);
B=A/max(max(A));
B=histeq(A);
figure,imshow(B);
colorbar
imag_test2{i}(:,:)=B;
end但我没有得到与纸张相同的结果。有谁知道为什么?或我出错的地方?
编辑
通过从@Amro获得帮助并使用他的代码,我最终获得了以下图像:
这里是我的原始图像从72小时发芽扁豆(400图像,每秒5帧):
这里是10个不同频段的结果图像: