在MATLAB图像处理算法 [英] image processing algorithm in MATLAB
问题描述
我想实现本文中所描述的算法：
I am trying to implement an algorithm described in this paper:
下面是该算法的说明：
我们记录了 N A序列连续斑点图像采样
频率 FS 。以这种方式，有可能观察如何像素
演变通过 N 图像。进化可以被视为一个时间
系列，并且可以通过以下方式进行处理：每个信号
对应于每个像素的演变被用作输入到一个
滤波器的银行。强度值是previously除以他们的
时间的平均值，以减少反射率的地方差异或
照明的对象。这可能是最大频率
充分分析由采样定理和同父异母的决定
采样频率 FS 。后者是由CCD照相机设置，则
图像的尺寸，和所述图像采集。过滤器的银行
在图中概述。 1。
We recorded a sequence of
Nsuccessive speckle images with a sampling frequencyfs. In this way it was possible to observe how a pixel evolves through theNimages. 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 frequencyfs. 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.
在我们的例子中，十5°阶Butterworth滤波器 根据所需被使用，但这个数目可以变化 歧视。该银行是在利用MATLAB计算机实现 软件。我们选择了黄油价值滤波器，因为，除了其 简单，它是最平坦的。其他滤波器，无限脉冲 反应，或有限脉冲响应可以使用。
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.
通过这种手段 滤波器，每个滤波器的每个滤波器的10对应的信号的银行， 临时像素进化得到作为输出。平均能量为EB 每个信号中，然后计算：
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:
其中， PB（N）是滤波像素的第n个图像的亮度
对于过滤器 B 以其平均值除以和 N 是总数
图像。以这种方式，恩获得的能量为每个像素值，
各下摆属于图1中的频带中的一个。 1。
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 code的基础上：
and here is my MATLAB code base on that :
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<new_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
但我没有得到同样的结果如纸。有任何人有任何想法，为什么？或者我已经错了？
but I am not getting the same result as paper. has anybody has any idea why? or where I have gone wrong?
修改
从@Amro获得帮助，并使用他的code我endup用下面的图片：
这是我从72小时原始图像发芽扁豆（400幅图像，以每秒5帧）：
EDIT by getting help from @Amro and using his code I endup with the following images: here is my Original Image from 72hrs germinated Lentil (400 images, with 5 frame per second):
这里是10个不同频段的结果图片：
here is the results images for 10 different band :
推荐答案
几个问题我可以当场：
-
在它的平均分的信号，你需要检查，这是不是零。否则，结果将是
NaN的。
作者（我下面这篇文章）使用银行过滤器具有覆盖整个范围可达奈奎斯特频率的频段。您正在做的一半。标准化后的频率传递给黄油应该去一路攀升至 1 （对应于 FS / 2 ）
the authors (I am following this article) used a bank of filters with frequency bands covering the entire range up to the Nyquist frequency. You are doing half of that. The normalized frequencies you pass to butter should go all the way up to 1 (corresponds to fs/2)
在计算每个滤波信号的能量，我想你不应该由它的平均分（你之前已经占到了这一点）。而不是简单地做： E =总和（SIG ^ 2）; 每个滤波后的信号的
When computing the energy of each filtered signal, I think you should not divide by its mean (you have already accounted for that before). Instead simply do: E = sum(sig.^2); for each of the filtered signals
在最后的后处理步骤，你应该正常化的范围[0,1]，然后应用中值滤波算法 medfilt2 。该计算不看的权利，应该是这样的：
In the last post-processing step, you should normalize to the range [0,1], and then apply the median filtering algorithm medfilt2. The computation doesn't look right, it should be something like:
img = ( img - min(img(:)) ) ./ ( max(img(:)) - min(img(:)) );
通过记住以上几点，我试图重写code以量化的方式。既然你没有张贴样品输入图像，我无法测试，如果结果如预期......加上我不知道如何跨preT最终图像呢：）
With the above points in mind, I tried to rewrite the code in a vectorized way. Since you didn't post sample input images, I can't test if the result is as expected... Plus I am not sure how to interpret the final images anyway :)
%# read biospeckle images
fnames = dir( fullfile('folder','myimages*.jpg') );
fnames = {fnames.name};
N = numel(fnames); %# number of images
Fs = 1; %# sampling frequency in Hz
sz = [209 278]; %# image sizes
T = zeros([sz N],'uint8'); %# store all images
for i=1:N
T(:,:,i) = imread( fullfile('folder',fnames{i}) );
end
%# timeseries corresponding to every pixel
T = reshape(T, [prod(sz) N])'; %# columns are the signals
T = double(T); %# work with double class
%# normalize signals before filtering (avoid division by zero)
mn = mean(T,1);
T = bsxfun(@rdivide, T, mn+(mn==0)); %# divide by temporal mean
%# bank of filters
numBanks = 10;
order = 5; % butterworth filter order
fCutoff = linspace(0, Fs/2, numBanks+1)'; % lower/upper cutoff freqs
W = [fCutoff(1:end-1) fCutoff(2:end)] ./ (Fs/2); % normalized frequency bands
W(1,1) = W(1,1) + 1e-5; % adjust first freq
W(end,end) = W(end,end) - 1e-5; % adjust last freq
%# filter signals using the bank of filters
Tf = cell(numBanks,1); %# filtered signals using each filter
for i=1:numBanks
[b,a] = butter(order, W(i,:)); %# bandpass filter
Tf{i} = filter(b,a,T); %# apply filter to all signals
end
clear T %# cleanup unnecessary stuff
%# compute average energy in each signal across frequency bands
Tf = cellfun(@(x)sum(x.^2,1), Tf, 'Uniform',false);
%# normalize each to [0,1], and build corresponding images
Tf = cellfun(@(x)reshape((x-min(x))./range(x),sz), Tf, 'Uniform',false);
%# show images
for i=1:numBanks
subplot(4,3,i), imshow(Tf{i})
title( sprintf('%g - %g Hz',W(i,:).*Fs/2) )
end
colormap(gray)
（我使用图像这里对于上述结果）的
(I used the image from here for the above result)
做了一些修改，简化了上述codeA位。这将减少内存占用。比如我用，而不是一个单一的多维矩阵单元阵列来存储结果。这样，我们不分配的连续内存的一个大块。我还重用，而不是引入新的在每一个中间步骤相同的变量...
Made some changes and simplified the above code a bit. This shall reduce memory footprint. For example I used cell array instead of a single multidimensional matrix to store the result. That way we don't allocate one big block of contiguous memory. I also reused same variables instead of introducing new ones at each intermediate step...
这篇关于在MATLAB图像处理算法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
