algorithm - MATLABの画像処理アルゴリズム

Question

このホワイトペーパーで説明されているアルゴリズムを実装しようとしています。

一時的なスペクトルバンドでのバイオスペックル画像の分解

アルゴリズムの説明は次のとおりです。

Nサンプリング周波数で連続したスペックル画像のシーケンスを記録しましたfs。このようにして、N画像を通してピクセルがどのように進化するかを観察することができました。その進化は時系列として扱うことができ、次の方法で処理できます。すべてのピクセルの進化に対応する各信号がフィルターバンクへの入力として使用されました。強度値は、オブジェクトの反射率または照明の局所的な違いを最小限に抑えるために、以前は時間平均値で除算されていました。適切に分析できる最大周波数は、サンプリング定理とサンプリング周波数の半分によって決定されますfs。後者は、CCDカメラ、画像のサイズ、およびフレームグラバーによって設定されます。フィルタのバンクの概要を図1に示します。

私たちの場合、10個の5°次のバターワースフィルターが使用されましたが、この数は必要な識別に応じて変えることができます。バンクは、MATLABソフトウェアを使用してコンピューターに実装されました。バターワースフィルターを選択したのは、その単純さに加えて、最大限にフラットであるためです。他のフィルター、無限インパルス応答、または有限インパルス応答を使用できます。

このフィルターバンクを使用して、各一時的なピクセル展開の各フィルターの10個の対応する信号が出力として取得されました。次に、各信号の平均エネルギーEbを計算しました。

ここpb(n)で、はフィルターのn番目の画像のフィルター処理されたピクセルの強度をbその平均値で割ったもので、Nは画像の総数です。このようにEnして、各ピクセルのエネルギー値が取得されました。各裾は、図1の周波数帯域の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の異なるバンドの結果画像です：

Answer 1

私が見つけることができるいくつかの問題：

• 信号をその平均で割るときは、それがゼロではなかったことを確認する必要があります。それ以外の場合、結果はになりますNaN。

• 著者（私はこの記事をフォローしています）は、ナイキスト周波数までの全範囲をカバーする周波数帯域を持つフィルターのバンクを使用しました。あなたはその半分をやっています。あなたが渡す正規化された周波数は、 （に対応する）butterまでずっと行く必要があります1fs/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

いくつかの変更を加え、上記のコードを少し簡略化しました。これにより、メモリフットプリントが削減されます。たとえば、結果を格納するために、単一の多次元行列の代わりにセル配列を使用しました。そうすれば、連続したメモリの1つの大きなブロックを割り当てません。また、各中間ステップで新しい変数を導入する代わりに、同じ変数を再利用しました...

Answer 2

この論文では、時系列の平均を差し引くことについては触れられていませんが、それが必要であると確信していますか？また、最後の画像からnew_maxとnew_minを1回だけ計算します。