colfiltを使った演算について

colfiltを使って、任意のブロックサイズの分散を1画素ごとずらしながら求めようとしています。

しかし、エラー：reshape要素数を変更してはなりません。と出てしまいます。

プログラムは以下の様に書きました。

img=imread('1111.bmp');
img2=imread('2222.bmp');

fun=@(block_struct) std2(block_struct)^2;

varimg=colfilt(img,[5 5],'sliding',fun);%imgの分散
varimg2=colfilt(img2,[5 5],'sliding',fun);%img2の分散

この時何が問題で、プログラムにエラーが出ているのでしょうか？

ご教授お願い致します。

Answer 1

[Shunichi Kusano]

colfiltのドキュメントには下記のように書いています。

• まず、colfilt は A の m 行 n 列の各近傍を、関数 im2col を使用して一時行列の列に再配置します。
• 次に、colfilt は関数 fun をこの一時行列に適用します。fun は、一時行列の各列の単一値を含む行ベクトルを返さなければなりません (sum などの列圧縮関数は適切なタイプの出力を返します)。
• 最後に、colfilt は fun で返されるベクトルを、関数 reshape を使用して A と同じサイズの行列に形状を変更します。

どういうことかというと例えば入力する画像が100x100の大きさだったとして、colfilt(A,[5 5],'sliding',fun)を実行すると、

1. まずim2colにより、大きさ25x10000の配列ができます。これは、行方向に各ブロックの値が1列に並び替えられてできたものです。slidingを指定しているので、となり合う列同士で値の重複は結構あります。
2. そのあとこの配列に例えばsumみたいな列圧縮関数（何も引数を指定しないと列ごと行方向に和をとって一つの値が計算される）を処理すると、結果1x10000の配列ができます。
3. これをreshapeで100x100に形を直すことで、スライディングブロックごとに和を取った結果の配列ができます。

ということで、エラーが出る理由はstd2関数を使っているからです。std2を使うと2次元配列の全ての値を使って標準偏差が計算されるので、上記で言うと25x10000の配列が1x1の大きさに圧縮されてしまいます。そのあと100x100にreshapeできなくなってしまうのでエラーが出ます。std2ではなくstdを使えば大丈夫のはずです。

[朋貴 熊田]

ご回答ありがとうございます。早速、熟考します!!

[朋貴 熊田]

加えて、質問なのですが2枚の画像を使って共分散を求めたいとも思っています。

以前、blockproc関数を使って分割画像で求めた時は、2枚の画像を3次元方向に重ね合わせて、

stacks=cut(3,img,img2);
fun3=@(block_struct) corr2(block_struct.data(:,:,1),block_struct.data(:,:,2));
corrimg=blockproc(stacks,[bl,br],fun3);

上記のように関数を渡し行っていました。

今回のcolfilt関数も同じ様に出来ますでしょうか？

[Shunichi Kusano]

corr2は相互相関係数を計算するもので、共分散とは異なります。欲しいのはどちらでしょうか？

いずれにしてもcolfiltは3次元配列を入れられないようです。また、似たようなことができる関数にnlfilterというものもありますが、こちらも3次元配列はアウトでした。

なのでここは原点に返ってim2colを使った方法が考えられます。2枚の画像をそれぞれim2colで各ブロックの値を列方向に並べ替えた配列を作ります。それぞれの列に対して共分散を計算して、reshapeで戻せば完成、という流れです。

下記は、共分散を計算するコートです。もし相互相関係数が欲しい場合は適当に改変してください。

img = imread('peppers.png');
img1 = img(:,:,1);
img2 = img(:,:,2);
% im2colでcolfilt内部で一時作成される行列を計算
colimg1 = im2col(im2double(img1),[3 3],'sliding');
colimg2 = im2col(im2double(img2),[3 3],'sliding');
% 各列でペアを作り共分散を計算する
colimg1_centered = colimg1 - mean(colimg1);
colimg2_centered = colimg2 - mean(colimg2);
c = sum(colimg1_centered .* colimg2_centered) / (size(colimg1,1)-1);
c = reshape(c,size(img,1:2)-2);
% 可視化
figure;imagesc(c);

im2col使うと画像の端のパディングが行われないので、できた画像は元の画像より少し小さくなってしまいます。これを防ぎたければ事前にpadarray関数で端っこのデータを作っておいてからim2colに入れれば大丈夫です。

・・・かなり手間ですね。nlfilterが3次元配列に対応するように開発に要望を出しておきますね。

[朋貴 熊田]

ご回答ありがとうございます。

示したコードの方は私が間違えていいたので、共分散のcov関数で大丈夫です。ありがとうございます。

colfiltは3次元を入れられないのですね...

来年のMatlabで対応されることを待ち望んでおります。

ご教授ありがとうございました。