%%%% ロボットビジョン(5T) sample06_moment.m  2015/5/28
clear;
Iin = imread('image03.jpg');
I = double(Iin);

% 画像の2値化-------------------------------------------------
threshold = 180;  % 2値化のしきい値

[n,m] = size(I);
Ibin = zeros(n,m); % nxmのゼロ行列（2値化画像の格納先）
for j = 1 : n
    for i = 1 : m
        if I(j,i) <= threshold;  % しきい値処理
            Ibin(j,i) = 1;
        end;
    end;
end;

% 0次モーメントおよび1次モーメント------------------------------
M00 = 0;
M10 = 0;
M01 = 0;
for j = 1 : n
    for i = 1 : m
        % 0次モーメント
        M00 = M00 + Ibin(j,i);
        % 1次モーメント
        M10 = M10 + i * Ibin(j,i);
        M01 = M01 + j * Ibin(j,i);
    end;
end;

% 0次モーメントを用いて，面積占有率を計算
M00ratio = 100 * M00/(n*m);

% 1次モーメントと0次モーメントを用いて，重心を計算
CGx = M10/M00;
CGy = M01/M00;

% 重心位置の目安を表示(重心が交点になるような濃度50%の十字線)
IbinCG = Ibin;
IbinCG(round(CGy),:) = 0.5; % 注）round()により小数点以下は切り捨てている
IbinCG(:,round(CGx)) = 0.5;


figure(1);
subplot(1,3,1); %画像を1x3並べて表示
imagesc(I); colormap(gray);

subplot(1,3,2);
imagesc(Ibin); colormap(gray);
title(['Area(pixels)=',num2str(M00), ', Area(%)=',num2str(M00ratio)]); %面積と占有率を表示

subplot(1,3,3);
imagesc(IbinCG); colormap(gray);
title(['CGx=',num2str(CGx), ', CGy=',num2str(CGy)]); %重心の座標を表示

%%%% EOF