MATLABの質問

質問

a = imread('AM100-1.png');
b = imread('AM100-2.png');
a = im2uint16(a);
b = im2uint16(b);
c=b-a;
c2=a-b;
d=c/10;
d2=c2/10;
p1=a+d;
P1=p1-d2;
p2=P1+d;
P2=p2-d2;
p3=P2+d;
P3=p3-d2;
p4=P3+d;
P4=p4-d2;
p5=P4+d;
P5=p5-d2;

PP1 = label2rgb(P1, @jet);
PP2 = label2rgb(P2, @jet);
PP3 = label2rgb(P3, @jet);
PP4 = label2rgb(P4, @jet);
PP5 = label2rgb(P5, @jet);

[nRows nCols] = size(P1);
for iCol = 1:nCols
for iRow = 1:nRows
if P1(iRow,iCol)<17000; PP1(iRow,iCol,:)=0;
if P2(iRow,iCol)<17000; PP2(iRow,iCol,:)=0;
if P3(iRow,iCol)<17000; PP3(iRow,iCol,:)=0;
if P4(iRow,iCol)<17000; PP4(iRow,iCol,:)=0;
if P5(iRow,iCol)<17000; PP5(iRow,iCol,:)=0;

end
end
end
end
end
end
end
A = cat(4,PP1,PP2,PP3,PP4,PP5);

初心者ながら上のようなプログラムを作ったんですが、うまくまとめられないでしょうか？

回答No.6

sgwjn · ベストアンサー

#2です。

> P1が17000以上であっても、同座標のP2～P5が17000未満ならデータを0にしたい

やっぱりそうですよね。
質問にあるコードにインデントを付けると、

for iCol = 1:nCols
　　for iRow = 1:nRows
　　　　if P1(iRow,iCol)<17000;
　　　　　　PP1(iRow,iCol,:)=0;

　　　　　　if P2(iRow,iCol)<17000;
　　　　　　　　PP2(iRow,iCol,:)=0;

　　　　　　　　if P3(iRow,iCol)<17000;
　　　　　　　　　　PP3(iRow,iCol,:)=0;

　　　　　　　　　　if P4(iRow,iCol)<17000;
　　　　　　　　　　　　PP4(iRow,iCol,:)=0;

　　　　　　　　　　　　if P5(iRow,iCol)<17000;
　　　　　　　　　　　　　　PP5(iRow,iCol,:)=0;
　　　　　　　　　　　　end
　　　　　　　　　　end
　　　　　　　　end
　　　　　　end
　　　　end
　　end
end

となり、P5はP1～P4が成立しない限り更新できなくなっていますから。
では、それも踏まえて３例ほど提示します。

a = imread('AM100-1.png');
b = imread('AM100-2.png');
a = im2uint16(a);
b = im2uint16(b);

以上は共通部分とし、

■ 行数重視
N = 5
d = (b - a) / 10;
d2 = (a - b) / 10;
A = zeros(size(a, 1), size(a, 2), 3, N);
for i = 1:N
　　P = (a + i * d) - i * d2;
　　A(;, ;, ;, i) = (17000 <= repmat(P, [1, 1, 3])) .* label2rgb(P, @jet);
end

■ 速度重視
d = (b - a) / 10;
d2 = (a - b) / 10;
P1 = (a + d) - d2;
P2 = (P1 + d) - d2;
P3 = (P2 + d) - d2;
P4 = (P3 + d) - d2;
P5 = (P4 + d) - d2;

PP1 = (17000 <= repmat(P1, [1, 1, 3])) .* label2rgb(P1, @jet);
PP2 = (17000 <= repmat(P2, [1, 1, 3])) .* label2rgb(P2, @jet);
PP3 = (17000 <= repmat(P3, [1, 1, 3])) .* label2rgb(P3, @jet);
PP4 = (17000 <= repmat(P4, [1, 1, 3])) .* label2rgb(P4, @jet);
PP5 = (17000 <= repmat(P5, [1, 1, 3])) .* label2rgb(P5, @jet);

A = cat(4, PP1, PP2, PP3, PP4, PP5);

■ 可読性重視
d = (b - a) / 10;
d2 = (a - b) / 10;
P1 = (a + d) - d2;
P2 = (P1 + d) - d2;
P3 = (P2 + d) - d2;
P4 = (P3 + d) - d2;
P5 = (P4 + d) - d2;

PP1 = label2rgb(P1, @jet);
PP2 = label2rgb(P2, @jet);
PP3 = label2rgb(P3, @jet);
PP4 = label2rgb(P4, @jet);
PP5 = label2rgb(P5, @jet);

for i = 1:size(P1, 1)
　　PP1(i, find(P1(i, :) < 17000), :) = 0;
　　PP2(i, find(P2(i, :) < 17000), :) = 0;
　　PP3(i, find(P3(i, :) < 17000), :) = 0;
　　PP4(i, find(P4(i, :) < 17000), :) = 0;
　　PP5(i, find(P5(i, :) < 17000), :) = 0;
end

A = cat(4, PP1, PP2, PP3, PP4, PP5);

という感じでしょうか。
行数重視は、一見しただけではもはや処理内容が分からなくなっていますので、これはやり過ぎだと思います。
もしコードの最適化を行うにしても、速度重視あたりまでが妥当ではないでしょうか？

実行環境がないため、これらのコードは机上チェックしかしていません。動作するとは思いますが、もし問題があればおっしゃってください。

お礼

ありがとうございます。大変助かりました。
また何かあったらよろしくお願いします

回答No.5

sgwjn

#2です。

例として提示させて頂いたコードの意図が上手く伝わっていないようですので、少し補足させていただきます。

ご質問の意図は、単純にコードの行数が減ればそれでＯＫということでしょうか？
そうであれば、#1の方が提示されている反復処理を多用する方法が良いでしょう。P1～P5は構造も処理内容も似通っているのでかなり簡略化できると思います。

ただし、その場合かなり処理速度が低下するはずです。（もちろん画像のサイズにもよりますが）

MATLABは、反復処理や反復処理内の条件分岐が存在すると極端に処理速度が低下します。コードの可読性だけでなく、現実の使用を考えた場合、反復処理は極力使用せず行列演算を用いるべきです。

そういった意図もあり、#2は反復処理を減らすように変更を加えています。提示させて頂いた例は、完全に反復処理をなくすことも可能ですが、可読性とのバランスを考えて一つ残しています。

あと、一つ疑問なのですが、P1が17000以上であれば、同座標のP2～P5が17000未満でもデータを0に更新しないのは意図した動作でしょうか？

お礼

単純にコード行数を減らしたかったのですが、処理速度の低下も困るので２がいいですね

P1が17000以上であっても、同座標のP2～P5が17000未満ならデータを0にしたいのですが、なっていませんでしたか・・・

回答No.4

irija_bari

ANo. 1 です。
＞・一応動作しております
ということなので、まずは単純にコードをまとめたいだけということなのですね。

imread で読み込んだ画像 a は、a[row, column, color] となるので、
(ただし、グレイスケールのとき color=1, カラーのとき color = 1(Red)～3(Blue)）
複数の画像を 1 つの行列に収めるために P[row, column, color, index] という行列を
想定します。

P(index) = ((indexMax-index)/indexMax)*a + (index/indexMax)*b
(ただし、0≦index≦indexMax)

のように index を変化させて、それを P に代入していきます。
あとは、P の各画素が閾値に達しているかをチェックしていくことになります。

またまた疑問が生じたので挙げておきます。
（私は画像処理の経験が無いので的外れの質問である可能性があります。）
・元画像の色の諧調が 2^16 であるなら、a = im2uint16(a); を実行することに意味はあるのでしょうか？
・PP1 = label2rgb(P1, @jet); の意図が不明です。Color Map を変更したいのなら、colormap(jet); ではないでしょうか？
・for iCol = 1:nCols 以下の文ですが、ご質問の文章で提示されたコードだと画像 P1 の画素が
　閾値(17000)に達していないときだけ画像 P2 の画素の閾値をチェックすることになりますがそれで良いのでしょうか？

※読みやすくするために、インデント（タブ）を全角スペースに置き換えています。
※カラー画像の場合は "==== RGB の場合 ====" の部分のコメントを外して、"==== グレイスケールの場合 ===="
　の部分をコメント アウトしてください。
==== test02.m ====
clear;
% 変化の段階数
x = 5;

% 画像の読み込み
a = imread('AM100-1.png');
b = imread('AM100-2.png');
a = im2uint16(a);
b = im2uint16(b);

% 画像 a から画像 b への変化
% P = (1-x)*a + x*b （ただし、0≦x≦1）という意味の式
for i = 1:x+1
　　P(:,:,:,i) = a(:,:,:) * ((x-i+1)/x) + b(:,:,:) * ((i-1)/x);
end

% Color Map を変化させたいということなのか？（意図不明）
for i = 1:x+1
　　PP(:,:,:,i) = label2rgb(P(:,:,:,i), @jet);
end

% 画像の画素数の取得
[nRows nCols tmp1 tmp2] = size(P);

% 各画素が閾値に達しない場合 0 にする
for iCol = 1:nCols
　　for iRow = 1:nRows
　　　　for i = 1:x+1
　　　　　　% ==== RGB の場合 ====
　　　　　　% RGB を取得
　　　　　　%r = PP(iRow, iCol, 1, i);
　　　　　　%g = PP(iRow, iCol, 2, i);
　　　　　　%b = PP(iRow, iCol, 3, i);
　　　　　　%
　　　　　　% 輝度を計算
　　　　　　%if (r * 0.29891) + (g * 0.58661) + (b * 0.11448) < 17000
　　　　　　%　　PP(iRow, iCol, 1, i) = 0;
　　　　　　%　　PP(iRow, iCol, 2, i) = 0;
　　　　　　%　　PP(iRow, iCol, 3, i) = 0;
　　　　　　%endif
　　　　　　% ここまで ==== RGB の場合 ====

　　　　　　% ==== グレイスケールの場合 ====
　　　　　　grayScale = PP(iRow, iCol, 1, i);

　　　　　　% 輝度を計算
　　　　　　if grayScale < 17000
　　　　　　　　PP(iRow, iCol, 1, i) = 0;
　　　　　　endif
　　　　　　% ここまで ==== グレイスケールの場合 ====
　　　　end
　　end
end

% 画像の連結
A(:,:,:) = PP(:,:,:,1);
for i = 2:x+1
　　A = cat(2, A, PP(:,:,:,i));
end
==== end of test02.m ====

お礼

・すいません２５６でした
・colormap(jet);だとカラー画像として保存してくれないようなのですが、表示のみ？
・すいません、よく分らないのですが。それでいいとおもいます

回答No.3

irija_bari

ANo. 1 です。
まだ途中ですが、とりあえずそれっぽいソースを書いてみました。
（手元に Matlab が無いため、Matlab のクローンソフトである Octave を使用しています。）

以下のソースコードを「test01.m」などのファイル名で保存し、
matlab で test01 と入力することで実行します。

いくつか疑問があるので挙げます。
・画像の色の諧調は 2^16 (= 65536) なのでしょうか？
・画像はグレイスケールなのでしょうか、それともカラーなのでしょうか？
・もう一度お尋ねしますが、提示されたプログラムは動作しているのでしょうか？
　（エラーで止まる or 動作はするが思ったとおりの結果ではない or 思ったとおりに動作する）

※読みやすくするために、インデント（タブ）を全角スペースに置き換えています。
==== test01.m ====
clear;
% 変化数
x = 5;

% 画像の読み込み
a = imread('AM100-1.png');
b = imread('AM100-2.png');

% 画像 a から画像 b への変化
for i = 1:x+1
　　P(:,:,:,i) = a(:,:,:) * ((x-i+1)/x) + b(:,:,:) * ((i-1)/x);
end

[nRows nCols tmp1 tmp2] = size(P);
for iCol = 1:nCols
　　for iRow = 1:nRows
　　　　for i = 1:x+1
　　　　　　% RGB の場合
　　　　　　% RGB を取得
　　　　　　%r = P(iRow, iCol, 1, i);
　　　　　　%g = P(iRow, iCol, 2, i);
　　　　　　%b = P(iRow, iCol, 3, i);
　　　　　　%
　　　　　　% 輝度を計算
　　　　　　%if (r * 0.29891) + (g * 0.58661) + (b * 0.11448) < 17000
　　　　　　%　　P(iRow, iCol, 1, i) = 0;
　　　　　　%　　P(iRow, iCol, 2, i) = 0;
　　　　　　%　　P(iRow, iCol, 3, i) = 0;
　　　　　　%endif

　　　　　　% グレイスケールの場合
　　　　　　grayScale = P(iRow, iCol, 1, i);

　　　　　　% 輝度を計算
　　　　　　if grayScale < 17000
　　　　　　　　P(iRow, iCol, 1, i) = 0;
　　　　　　endif
　　　　end
　　end
end
==== end of test01.m ====

お礼

・画像の色の諧調は65536です
・画像はグレイスケールです（できればそこからカラーにしたいのですが・・）
・一応動作しております

回答No.2

sgwjn

単純に整理するだけなら、こんなところでしょうか。

a = imread('AM100-1.png');
b = imread('AM100-2.png');
a = im2uint16(a);
b = im2uint16(b);

d = (b - a) / 10;
d2 = (a - b) / 10;
P1 = (a + b) - d2;
P2 = (P1 + d) - d2;
P3 = (P2 + d) - d2;
P4 = (P3 + d) - d2;
P5 = (P4 + d) - d2;

PP1 = label2rgb(P1, @jet);
PP2 = label2rgb(P2, @jet);
PP3 = label2rgb(P3, @jet);
PP4 = label2rgb(P4, @jet);
PP5 = label2rgb(P5, @jet);

for i = 1:size(P1, 1)
　　PP1(i, find(P1(i, :) < 17000), :) = 0;
　　PP2(i, find((P1(i, :) < 17000) & (P2(i, :) < 17000)), :) = 0;
　　PP3(i, find((P1(i, :) < 17000) & (P2(i, :) < 17000) & (P3(i, :) < 17000)), :) = 0;
　　PP4(i, find((P1(i, :) < 17000) & (P2(i, :) < 17000) & (P3(i, :) < 17000) & (P4(i, :) < 17000)), :) = 0;
　　PP5(i, find((P1(i, :) < 17000) & (P2(i, :) < 17000) & (P3(i, :) < 17000) & (P4(i, :) < 17000) & (P5(i, :) < 17000)), :) = 0;
end

A = cat(4, PP1, PP2, PP3, PP4, PP5);

お礼

ありがとうございます参考になりました

回答No.1

irija_bari

回答される方がいらっしゃらないようなので。。。

・何を意図したプログラムなのか明記しましょう。
・動作はするのでしょうか？デバグなのか？それともコードを簡潔にしたいのか？

お礼

2つの熱画像A,B（ここでは強度イメージ）の差分を取りその差分をXで割りX回たすことで徐々にAの画像をBにちかづけて、
更にcolormapをjetにして強度イメージの輝度が17000以下の所は輝度０にしてそれの変化を動画のようにするというもです