一尘不染

MATLAB中的图像处理算法

algorithm

我正在尝试实现本文描述的算法:

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

这是算法的说明:

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

过滤器库

在我们的案例中,使用了十个5°阶Butterworth滤波器,但是该数量可以根据所需的区分度而变化。该银行是使用MATLAB软件在计算机上实现的。我们选择Butter-
worth滤波器是因为它不仅简单,而且最大程度地平坦。可以使用其他滤波器,无限冲激响应或有限冲激响应。

通过该组滤波器,获得了每个临时像素演变的每个滤波器的十个对应信号作为输出。然后计算每个信号中的平均能量Eb:

能量方程

其中,pb(n)是第n个要过滤的图像中经过滤波的像素的强度b除以其平均值,并且N是图像总数。这样,En获得了每个像素的能量值,每个下摆都属于图1中的一个频带。

利用这些值,可以构建活动对象的十个图像,每个图像都显示在某个频带中存在多少随时间变化的散斑能量。错误地将颜色分配给结果中的灰度级将有助于区分。

这是我基于此的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获得帮助并使用他的代码进行 *编辑
,我得到了以下图像:这是我从72hrs发芽的扁豆中提取的原始图像(400张图像,每秒5帧):
在此处输入图片说明

这是10个不同波段的结果图像:

乐队1
乐队2
乐队3
乐队4
乐队5
乐队6
乐队7
BAnd8
乐队9
乐队10


阅读 323

收藏
2020-07-28

共1个答案

一尘不染

我可以发现几个问题:

  • 当您将信号除以均值时,需要检查它是否不为零。否则结果将是NaN

  • 作者(我正在关注本文)使用了一组滤波器,这些滤波器的频带覆盖了直到奈奎斯特频率的整个范围。您正在做一半。您传递给的归一化频率butter应该一直上升到1(对应于fs/2

  • 在计算每个滤波信号的能量时,我认为您不应该除以其均值(您之前已经考虑过了)。取而代之的是:E = sum(sig.^2);对每个滤波信号

  • 在最后的后处理步骤中,应将其归一化为[0,1]范围,然后应用中值滤波算法medfilt2。计算看起来不正确,应该是这样的:

    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

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

2020-07-28