EmbVision(CQ版) チュートリアル（７）

回帰分析と曲線検出 - MATLAB編 -

新潟大学 村松　正吾，高橋　勇希

Copyright (c), All rights reserved, 2014-2025, Shogo MURAMATSU and Yuki TAKAHASHI

Contents

Part6 | | メニュー | | Part8

概要

本演習では、System object の応用例として、 多項式回帰を利用した画像からの曲線検出を実装する。

準備として、開いている全ての Figure を close 関数で 閉じておく。

close all

回帰分析

回帰分析では、ある変数 $y$ とある変数 $x$ の間を

$$ y = f(x) $$

のように関係づける未知の関数 $f(\cdot)$ を仮定する。 回帰分析は、このような未知の関数 $f(\cdot)$ を観測可能なデータから推論する 統計的手法である。

変数 $x$ を説明変数とよび、変数 $y$ を目的変数とよぶ。 また、この未知関数の推論をフィッティングとよぶ。

以下では、回帰分析の一例としてデータの多項式近似を試みる。 まず、多項式回帰を試すためのデータを用意しよう。

load census

loadコマンドで読み込まれた変数 cdate, popは いずれも21次元ベクトルであり、その内容は以下のとおり。

plot関数で、cdateをX軸、popをY軸にしてX-Yプロットを 表示しよう。

plot(cdate,pop,'o')
xlabel('year')
ylabel('population')

プロットされたグラフを見ると二次関数のような増え方をしている。

MATLABは多項式近似を行う polyfit関数を提供している。 $x$ を説明変数の観測データ， $y$ を目的変数の観測データとして同じ次元の ベクトルで与えられたとすれば、コマンド

>> p = polyfit(x,y,n)

は、未知の関数 $f(x)$ を $n$ 次多項式

$$ p(x) = p_1x^n+p_2x^{n-1}+\cdots+p_nx+p_{n+1} $$

で近似する。ただし、 $p_1,p_2,\cdots,p_{n+1}$ は未知パラメータで、 データとの整合性が最小自乗誤差の意味で最適になるように選択される。

polyfit　関数をデータに適用して近似多項式を求めてみよう。

n = 2;
p = polyfit(cdate,pop,n);

づづけて、得られた近似多項式を polyval 関数により 評価しよう。 linspace 関数で、年度 cdate の最小値と 最大値の間を100点等間隔にサンプルする。

x = linspace(min(cdate),max(cdate),100);
y = polyval(p, x);
hold on
plot(x,y)

[ トップ ]

曲線検出

System object の応用例として、曲線を含む二値画像から多項式近似により 曲線を近似するモジュールを作成する。

以下のテストクラスを定義する。

classdef CurveDetectionSystemTestCase < matlab.unittest.TestCase
    methods(Test)
        % Test methods
        function testDefaultDegree(testCase)
            % 期待値　
            degExpctd = 3;
            % ターゲットクラスのインスタンス化
            obj = CurveDetectionSystem();
            % プロパティ Degree の取得
            degActual = obj.Degree;
            % プロパティ Degree の検証
            testCase.verifyEqual(degActual,degExpctd)
        end
        function testCoefs(testCase)
           % 準備
           xmax = 4;
           x = 1:xmax; % 説明変数
           y = x.^2;   % 目的変数
           BW = zeros(xmax^2,xmax); % 座標を二値画像に変換
           for idx = 1:length(x)
               BW(y(idx),x(idx)) = 1;
           end
           deg = 3;  % 次数
           % 期待値の設定
           coefsExpctd = polyfit(x,y,deg); % 多項式係数の期待値
           y = round(polyval(coefsExpctd,x));   % 目的変数
           lineExpctd = zeros(size(BW),'like',BW); % 抽出される曲線の期待値
           for idx = 1:length(x)
               lineExpctd(y(idx),x(idx)) = 1;
           end
           % ターゲットのインスタンス化
           obj = CurveDetectionSystem();
           % 実行結果（実現値の取得）
           [lineActual,coefsActual] = obj.step(BW);
           % 配列の値の検証
           testCase.verifyEqual(coefsActual,coefsExpctd,'AbsTol',1e-9);
           testCase.verifyEqual(lineActual,lineExpctd,'AbsTol',1e-9);
      end
        function testSetDegree(testCase)
            % 期待値
            degExpctd = 2;
            % ターゲットクラスのインスタンス化
            obj = PolyfitSystem('Degree',degExpctd);
            % プロパティー Degree の取得
            degActual = obj.Degree;
            % プロパティー Degree の検証
            testCase.verifyEqual(degActual,degExpctd)
        end
    end
end

上記のテストクラスを満足するターゲットクラスとして以下のSystem objectを定義しよう。

classdef CurveDetectionSystem < matlab.System
% 二値画像から近似曲線画像に変換
    properties (Nontunable)
        Degree = 3;
    end
    methods
        %コンストラクタ
        function obj = CurveDetectionSystem(varargin)
            setProperties(obj,nargin,varargin{:})
        end
    end
    methods(Access = protected)
        function [Line,p] = stepImpl(obj,BW)
            [y,x] = find(BW); % 検出点を座標に変換
            p = polyfit(x,y,obj.Degree); % 多項式近似

             line = zeros(size(bw),'like',bw);
             xline = 1:size(bw,2); % 説明変数を画像の端から端とする
             yline = uint8(polyval(p,xline));    % 多項式から目的変数を求める画素に収まるように整数型にしておく
             yline = min(max(yline,1),size(BW,1));   % 目的変数が画像内に収まるように変更する。
             ind = sub2ind(size(BW),yline,xline);  % 座標から画像位置に変換
             line(ind) = 1;    % 線を描画
         end
     end
 end
% sz = [288 352];
% BW = zeros(sz);
% ndata = 21;
% x = randi(sz(2),ndata,1);
% deg = 3;
% c = 1e-6*randn(deg,1);
% y = round(x.^(deg:-1:1)*c+sz(1)/2);
% for idx = 1:ndata
%     BW(y(idx),x(idx)) = 1;
% end
% p = polyfit(x,y,deg);
% obj = PolyfitSystem('Degree',deg);
% obj.step(BW)

[ トップ ]

曲線描画

[ トップ ]

演習課題

課題7-1. XXX

...

課題7-2. XXX

...

(処理例）

Part6 | | メニュー | | トップ | | Part8


Published with MATLAB® R2024b