MATLAB中的图像处理算法

我正试图实现本文中描述的算法：

临时光谱带中生物斑纹图像的分解

这是对算法的解释：

我们以采样频率fs记录了一系列连续的N个斑点图像。 用这种方法可以观察像素如何通过N图像进行演变。 这种演变可以看作是一个时间序列，可以按照以下方式进行处理：每个像素的演变对应的每个信号都被用作一组滤波器的输入。 之前将强度值除以它们的时间平均值以最小化对象的反射率或照明的局部差异。 可以充分分析的最大频率取决于采样定理和采样频率fs的一半。 后者由CCD相机，图像大小和图像采集卡设置。 过滤器组在图1中概述。

在我们的例子中，使用了10个5阶巴特沃斯滤波器，但是这个数量可以根据所需的辨别度而变化。 该银行使用MATLAB软件在计算机中实施。 我们选择黄油价值过滤器，因为除了简单之外，它最大限度地平坦。 可以使用其他滤波器，无限脉冲响应或有限脉冲响应。

通过这组滤波器，每个临时像素进化的每个滤波器的十个对应信号被获得作为输出。 然后计算每个信号的平均能量Eb：

其中pb(n)是滤波器b的第n个图像中滤波像素的强度除以其平均值， N是图像的总数。 以这种方式，获得每个像素的能量的En值，每个边缘属于图1中的一个频带。

通过这些值，可以构建10个活动对象的图像，每个图像显示在特定频带内存在多少时变斑点能量。 对结果中的灰色级别进行虚假颜色分配将有助于区分。

这里是我的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<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

但我没有得到与纸张相同的结果。 有谁知道为什么？ 或我出错的地方？

通过从@Amro获得帮助编辑并使用他的代码我最终以下图像：这是我的原始图像从72小时发芽扁豆（400图像，每秒5帧）：

这里是10个不同频段的结果图像：

我可以发现几个问题：

img = ( img - min(img(:)) ) ./ ( max(img(:)) - min(img(:)) );

编辑：

考虑到以上几点，我尝试以矢量化的方式重写代码。 由于您没有发布样本输入图像，我无法测试结果是否如预期的那样......另外，我不知道如何解释最终图像:)

%# 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)

（我使用这里的图像来获得上述结果）

编辑＃2

做了一些修改并简化了上面的代码。 这将减少内存占用。 例如，我使用单元阵列而不是单个多维矩阵来存储结果。 这样我们就不会分配一大块连续的内存。 我也重复使用了相同的变量，而不是在每个中间步骤引入新变量。

该论文没有提到减去时间序列的平均值，你确定这是必要的吗？ 此外，您只能从最后一张图像计算一次new_max和new_min。

上一篇: image processing algorithm in MATLAB

下一篇: Implementing the intelligent recursive algorithm in matlab