Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

動いている振子の長さを検出

この例では、振子をセグメント化し、多数のイメージ フレームから振子の中心を計算し、中心座標を円の方程式に当てはめることによって、動いている振子の長さを計算する方法を説明します。

手順 1: イメージの読み込み

MAT ファイルから振子のデータを読み込みます。このファイルには、2 つの変数が含まれています。変数 frames には、ビデオが 4 次元数値配列として格納されます。変数 rect には、振子の範囲を制限する関心領域が格納されます。

load pendulum;

implay を使用してビデオをプレビューします。各フレームの上半分で振子が揺れているのがわかります。

implay(frames);

手順 2: 振子が揺れている領域の選択

セグメント化の効率を上げるために、振子が揺れている領域のみが含まれる新しい一連のフレームを作成します。

最初に、プリロードされた変数 rect としてトリミング領域を指定して、imcrop を 1 つのフレーム上で実行します。次に、トリミングしたすべてのフレームを格納するための配列を、最初にトリミングしたフレームのサイズに基づいて初期化します。最後に、すべてのフレームをループ処理し、同じトリミング領域を使用して各フレームをトリミングし、初期化した配列に結果を保存します。

first_frame = frames(:,:,:,1);
first_region = imcrop(first_frame,rect);
nFrames = size(frames,4);
frame_regions = repmat(uint8(0), [size(first_region) nFrames]);
for count = 1:nFrames
  frame_regions(:,:,:,count) = imcrop(frames(:,:,:,count),rect);
end
imshow(frame_regions(:,:,:,1))

手順 3: 各フレームでの振子のセグメント化

振子は背景よりも暗くなっていることに注目してください。フレームをグレースケールに変換し、imbinarize を使用してしきい値処理し、関数 imopen および imclearborder を使用して背景構造を除去することで、各フレームで振子をセグメント化できます。

バイナリ配列を初期化して、セグメント化された振子のフレームを格納します。

seg_pend = false([size(first_region,1) size(first_region,2) nFrames]);
centroids = zeros(nFrames,2);
se_disk = strel("disk",3);

すべてのフレームをループ処理し、セグメント化と後処理を実行します。1 つのモンタージュに、元のフレームとセグメント化の結果を表示します。

for count = 1:nFrames
    fr = frame_regions(:,:,:,count);

    gfr = im2gray(fr);
    gfr = imcomplement(gfr);

    bw = imbinarize(gfr,0.7);  % Threshold is determined experimentally
    bw = imopen(bw,se_disk);
    bw = imclearborder(bw);
    seg_pend(:,:,count) = bw;

    montage({fr,labeloverlay(gfr,bw)});
    pause(0.2)

end

手順 4: 各フレームでセグメント化された振子の中心の検出

フレームによって振子の形状が異なっていることがわかります。必要なのは振子の中心のみなので、これは重大な問題ではありません。振子の中心を使用して、振子の長さを検出します。

regionprops を使用して、振子の中心を計算します。

pend_centers = zeros(nFrames,2);
for count = 1:nFrames
    property = regionprops(seg_pend(:,:,count),"Centroid");
    pend_centers(count,:) = property.Centroid;
end

plot を使用して、振子の中心を表示します。

x = pend_centers(:,1);
y = pend_centers(:,2);
figure
plot(x,y,"m.")
axis ij
axis equal
hold on;
xlabel("x");
ylabel("y");
title("Pendulum Centers");

手順 5: 振子の中心を介した円の近似による半径の計算

以下の円の基本方程式を書き換えます。

(x-xc)^2 + (y-yc)^2 = radius^2

ここで、(xc,yc) は中心であり、パラメーター abc で表すと、以下のようになります。

x^2 + y^2 + a*x + b*y + c = 0

ここで、a = -2*xcb = -2*ycc = xc^2 + yc^2 - radius^2 です。

最小二乗法を使用して、パラメーター ab、および c について解くことができます。上記の方程式を以下のように書き換えます。

a*x + b*y + c = -(x^2 + y^2)

以下のように書き換えることもできます。

[x y 1] * [a;b;c] = -x^2 - y^2.

バックスラッシュ (\) 演算子を使用してこの方程式を解きます。

円の半径は振子の長さです (ピクセル単位)。

abc = [x y ones(length(x),1)] \ -(x.^2 + y.^2);
a = abc(1);
b = abc(2);
c = abc(3);
xc = -a/2;
yc = -b/2;
circle_radius = sqrt((xc^2 + yc^2) - c);
pendulum_length = round(circle_radius)
pendulum_length =

   253

振子の中心のプロットの上に、円および円の中心を重ね合わせます。

circle_theta = pi/3:0.01:pi*2/3;
x_fit = circle_radius*cos(circle_theta)+xc;
y_fit = circle_radius*sin(circle_theta)+yc;

plot(x_fit,y_fit,"b-");
plot(xc,yc,"bx","LineWidth",2);
plot([xc x(1)],[yc y(1)],"b-");
text(xc-110,yc+100,sprintf("Pendulum length = %d pixels",pendulum_length));

参考

| | | | |