画像処理ブログ

画像処理エンジニアのおっさんが興味あることや調べたことを徒然のままに垂れ流すブログ

2017/04/08 関東CV勉強会参加報告 6巻第2章 幾何学的推定のための最適化手法:最小化を超えて

画像処理の勉強のため、関東CV勉強会に参加してきました。へたれなので聴講だけですが。
コンピュータビジョン勉強会

この勉強会は、コンピュータビジョン最先端ガイドシリーズを教科書に、画像処理のトピックについて発表形式で勉強を進める会です。
教科書の内容を扱う他、CVPR等学会の読み会、トレンドとなっている内容を扱っていらっしゃいます。
connpassで参加登録ができます。

コンピュータビジョン勉強会@関東 - connpass


今回は幾何学的推定のための最適化手法というトピックでした。
書かれたのは岡山大学の金谷教授です。

今回の内容については、近い(というかほぼ同じ)内容がWebで公開されています。
http://www.iim.cs.tut.ac.jp/~kanatani/papers/hypertutor.pdf

幾何学的問題の最適化

画像処理の3次元復元や自己位置同定、AR等にも関係の深いトピックになります。
例えば画像中の楕円を検出する際には、楕円のエッジ上にある特徴点から、楕円の式を最適化します。
ステレオカメラの両方のカメラに写った物体からそれぞれ特徴点を抽出し、ステレオカメラのカメラ間の関係を表す基礎行列を求める際にも最適化が使われます。

勉強会の発表資料はこちらで公開されています。
www.slideshare.net

最小化による方法

勉強会の発表資料はこちらで公開されています。
CVIM最先端ガイド6 幾何学的推定のための最適化手法 3.5 - 3.8:embed:cite]

最小化によらない方法

勉強会の発表資料はこちらで公開されています。
speakerdeck.com
大変温かみのある資料でございました。

超精度くりこみ法
勉強会の発表資料はこちらで公開されています。

www.slideshare.net

KLT法解説 まとめ

KLT(Kanade-Lucas-Tomasi Feature Tracker)法はCarnegie Mellon UniversityのKanade Takeo先生らが中心となって開発されたImage Registration, Feature Trackingの手法です。 手法を勉強して、大変感銘を受けたので、ブログ記事にまとめました。 こちらの記事はその索引になります。

まずは最初に発表された論文の式を追いました。

Bruce D. Lucas and Takeo Kanade. An Iterative Image Registration Technique with an Application to Stereo Vision. International Joint Conference on Artificial Intelligence, pages 674-679, 1981. https://cecas.clemson.edu/~stb/klt/lucas_bruce_d_1981_1.pdf

m-yoshi-1700.hatenablog.com

m-yoshi-1700.hatenablog.com

m-yoshi-1700.hatenablog.com

m-yoshi-1700.hatenablog.com

式の導出が終わったので、実際にC++で動かしてみました。

m-yoshi-1700.hatenablog.com

KLT: Kanade-Lucas-Tomasi Feature Trackerをリアルタイムで動かして特徴点追跡をするまで

こちらの記事では、KLT法(KLT: Kanade-Lucas-Tomasi Feature Tracker)をmac、あるいはlinux上で読み込み、リアルタイムで特徴点抽出、追跡をするまでを説明します。
筆者のPC環境は、一昔前のMac book air(Core 2 Duo)です。言語はC++のみを使います。gccのバージョンは4.2.1です。画像の読み込み部分にOpenCVを使っています。


この技術についての投稿は、一つの索引用の記事にまとめられています。
全体を俯瞰する場合はそちらを御覧ください。


m-yoshi-1700.hatenablog.com

ライブラリをダウンロード

以下のサイトからライブラリをダウンロードして展開します。
https://cecas.clemson.edu/~stb/klt/installation.html

ライブラリのビルド

ライブラリをmakeします

cd klt
make

make後、エラーがなく、libklt.aという静的リンクライブラリが生成されていれば完了です。

読み込み側のプログラム

以下の様にプログラムを作成します。
非常に汚いプログラムで申し訳ありません...

sample.cpp:

#include <iostream>
#include <vector>
#include <opencv2/opencv.hpp>

#include "klt.h"

using namespace std;
using namespace cv;


Rect box;
bool dragging = false;
bool reset = false;

bool track = false;

typedef struct {
  int xdiff;
  int ydiff;
} TrackResult;

void my_mouse_callback(int event, int x, int y, int flags, void* param){
  switch (event){
  case cv::EVENT_MOUSEMOVE:
    if (dragging){
      box.width = x - box.x;
      box.height = y - box.y;
    }
    break;

  case cv::EVENT_LBUTTONDOWN:
    track = false;
    reset = false;
    dragging = true;
    box = cv::Rect(x, y, 0, 0);
    break;

  case cv::EVENT_LBUTTONUP:
    dragging = false;
    reset = true;
    // track = true;
    if (box.width < 0){
      box.x += box.width;
      box.width *= -1;
    }
    if (box.height < 0){
      box.y += box.height;
      box.height *= -1;
    }
    break;
  }
}

int main(int argc, char** argv)
{
  cv::VideoCapture cap(0);
  if(!cap.isOpened())
    return -1;

  int nFeatures = 20;

  KLT_TrackingContext tc;// =
  KLT_FeatureList fl;// = KLTCreateFeatureList(nFeatures);

  Mat previousFrame;

  cv::namedWindow("window");
  cv::setMouseCallback("window", my_mouse_callback);

  // box = Rect(0, 0, 640, 480);
  box = Rect(0, 0, 320, 240);
  reset = true;
  
  while(1){
    cv::Mat frame;
    cap >> frame;

    cv::resize(frame, frame, cv::Size(), 0.5, 0.5);
    
    Mat grayImage;
    cv::cvtColor(frame, grayImage, CV_BGR2GRAY);

    vector<Point> points;

    vector<Point> beforeTrack;
    vector<TrackResult> afterTrack;

    if(reset){
      tc = KLTCreateTrackingContext();

      tc->sequentialMode = TRUE;
      tc->writeInternalImages = FALSE;
      // tc->affineConsistencyCheck = -1;

      tc->mindist = 20;
      tc->window_width  = 3;
      tc->window_height = 3;
      KLTChangeTCPyramid(tc, 10);
      KLTUpdateTCBorder(tc);
  
      KLTPrintTrackingContext(tc);
      
      fl = KLTCreateFeatureList(nFeatures);
      KLTSelectGoodFeatures(tc, grayImage(box).data, box.width, box.height, fl);

      reset = false;
      track = true;

      previousFrame = grayImage.clone();
      continue;
    }
    if (track){
      beforeTrack.clear();
      for(int i = 0; i < nFeatures; ++i){
	KLT_Feature feature = fl->feature[i];
	Point point(feature->x, feature->y);

	beforeTrack.push_back(point);
      }
      
      KLTTrackFeatures(tc,
		       (previousFrame(box)).data,
		       (grayImage(box)).data,
		       box.width,
		       box.height,
		       fl);

      afterTrack.clear();
      for(int i = 0; i < nFeatures; ++i){
	KLT_Feature feature = fl->feature[i];
	if (feature->val < 0)
	  continue;

	TrackResult result;
	result.xdiff = feature->x - beforeTrack[i].x;
	result.ydiff = feature->y - beforeTrack[i].y;
	if (result.xdiff == 0 && result.ydiff == 0)
	  continue;
	afterTrack.push_back(result);
      }

      KLTReplaceLostFeatures(tc, (grayImage(box)).data, box.width, box.height, fl);
      points.reserve(fl->nFeatures);
    }

    if (track){
      for(int i = 0; i < fl->nFeatures; ++i){
	KLT_Feature feature = fl->feature[i];
	Point center(box.x + feature->x, box.y + feature->y);

	if (center.x < 0 || center.y < 0|| center.x > grayImage.cols || center.y > grayImage.rows)
	  continue;

	points.push_back(center);
      
	circle(frame, center, 3, Scalar(0, 255, 0), -1);
      }

      double xsum = 0;
      double ysum = 0;
      for(int i = 0; i < afterTrack.size(); ++i){
	xsum += afterTrack[i].xdiff;
	ysum += afterTrack[i].ydiff;
      }

      xsum /= afterTrack.size();
      ysum /= afterTrack.size();

      TrackResult diff;

      diff.xdiff = (int)xsum;
      diff.ydiff = (int)ysum;

	diff.xdiff = 0;
	diff.ydiff = 0;

      box.x += diff.xdiff;
      box.y += diff.ydiff;

      for(int i = 0; i < nFeatures; ++i){
	fl->feature[i]->x += diff.xdiff;
	fl->feature[i]->y += diff.ydiff;
      }

    }

    previousFrame = grayImage.clone();

    cv::imshow("window", frame);

    int key = cv::waitKey(1);
    if(key == 113)
      break;
    else if(key == 115)
      cv::imwrite("img.png", frame);
  }
  cv::destroyAllWindows();
  return 0;
}

cmake用のCMakeLists.txtも用意します。
CMakeLists.txt

cmake_minimum_required(VERSION 2.8)

set(cmake_build_type debug)

project(cv_shokyu)

find_package(opencv)
find_package(boost)

include_directories("${PROJECT_SOURCE_DIR}")
include_directories("/usr/local/include")
include_directories(${OpenCV_INCLUDE_DIRS})
include_directories(${Boost_INCLUDE_DIRS})
include_directories("${PROJECT_SOURCE_DIR}/klt")

link_directories("/usr/local/lib")

add_library(klt
    STATIC
    IMPORTED
)

set_target_properties(klt
	PROPERTIES
	IMPORTED_LOCATION ${CMAKE_SOURCE_DIR}/klt/libklt.a)

add_executable(sample sample.cpp)

target_link_libraries(sample klt ${OpenCV_LIBRARIES} ${Boost_LIBRARIES})

以下のコマンドでビルドします。

cmake .
make

これで準備は完了です。
実行ファイルsampleを実行すると、PCに接続されたカメラを使って、リアルタイムに特徴点の検出、追跡をします。

https://youtu.be/ZLj8A0yOpqo

KLT: Kanade-Lucas-Tomasi Feature Trackerをリアルタイムで動かして特徴点追跡をするまで

LK法の解説 2次元への拡張

前回までに1次元の場合のLK法の式の導出、その一般化と実装をしてパフォーマンスの検証をしました。


この技術についての投稿は、一つの索引用の記事にまとめられています。
全体を俯瞰する場合はそちらを御覧ください。


m-yoshi-1700.hatenablog.com


こちらの記事では、画像などの2次元の情報についてLK法で解を求める方法を解説します。

2次元への拡張

まず、これまでに求めたLK法の式は以下になります。

{ \displaystyle
\begin{eqnarray}
h_{0} &=& 0\\\
h\_{k + 1} &=& h\_{k} +
\left[ \frac {\sum\_{x} w(x) F'(x + h\_{k})[G(x) - F(x + h\_{k})]} {\sum\_{x} w(x) F'(x + h\_{k})^{2}} \right]
\end{eqnarray}
\tag{10}
}

この式を2次元に拡張する上で問題になるのは、Fの導関数を求めている部分です。

2次元の信号について、E2ノルムは以下の様に定義されます。

{ \displaystyle
E = \sum_{x \in R} \left[ F(\textbf{x} + \textbf{h}) - G(\textbf{x}) \right]^2
}

ここで、xhはn次元のベクトルです。

前回の式(8)の様に線形近似をしていきます。

{ \displaystyle
F(\textbf{x} + \textbf{h})
\approx F(\textbf{x}) + \textbf{h}
\frac {\delta}{\delta \textbf{x}}
F(\textbf{x})
}

偏微分をして誤差を最小化します。

{ \displaystyle
\begin{eqnarray}
0 &=& \frac {\delta E} {\delta h}\\\

&\approx& \frac {\delta} {\delta \textbf{h}}
\sum_{\textbf{x}}
\left[ F(\textbf{x}) + \textbf{h}
\frac {\delta}{\delta \textbf{x}}
F(\textbf{x}) - G(\textbf{x})
\right]^{2}\\\

&=&
\sum_{\textbf{x}}
2 \frac {\delta F}{\delta \textbf{x}}
\left[ F(\textbf{x}) + \textbf{h}
\frac {\delta}{\delta \textbf{x}}
F(\textbf{x}) - G(\textbf{x})
\right]
\end{eqnarray}
}

式を変形して以下のhを得ます。

{ \displaystyle
\textbf{h} \approx
\left[ \sum_{\textbf{x}} ( \frac {\delta F} {\delta \textbf{x}} )^{\mathrm{T}}
\left[ G(\textbf{x}) - F(\textbf{x}) \right]
\right]
\left[ \sum_{\textbf{x}} ( \frac {\delta F} {\delta \textbf{x}} )^{\mathrm{T}} ( \frac {\delta F} {\delta \textbf{x}} ) \right]^{-1}
}

これでLK法の式の導出は完了です。

次回は特徴点の抽出についてまとめます。

LK法のパフォーマンス

こちらの記事では、1次元の場合についてLK法を実装し、パフォーマンスについて解説します。

この技術についての投稿は、一つの索引用の記事にまとめられています。 全体を俯瞰する場合はそちらを御覧ください。

m-yoshi-1700.hatenablog.com

前々回からLK法の式の導出をしてきました。

ここまでで求めた式は以下になります。

{ \displaystyle
\label{eq10}
h_{0} = 0\\
h_{k + 1} = h_{k} +
\left[ \frac {\sum_{x} w(x) F'(x + h_{k})[G(x) - F(x + h_{k})]} {\sum_{x} w(x) F'(x + h_{k})^{2}} \right]
\tag{10}
}

こちらの式で実際のhに収束する条件と、収束するための速度について考えます。

LK法の実装

ここで、F(x)とG(x)を以下の様に置きます。

{ \displaystyle
F(x) = \sin x\\
G(x) = F(x + h) = sin(x + h)
}

こちらをpythonで実装します。 まず、F(x)とG(x)をグラフに描画してみます。

#!/usr/bin/env python
#-*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

# hの真値
h_g = 0.3


def f(x):
    return np.sin(x)


def g(x):

    return np.sin(x + h_g)


def main():
    x = np.arange(0, 2 * np.pi, 0.01)

    f_val = f(x)

    g_val = g(x)

    plt.plot(x, f_val, label="f(x)")
    plt.plot(x, g_val, label="g(x)")
    plt.legend()
    plt.show()

if __name__ == "__main__":
    main()

こちらを実行すると、次のようなグラフが表示されます。

f:id:m-yoshi-1700:20170319122604p:plain

Fが正弦波で、Gがそれをhだけずらしたものになっています。

ここから、式\eqref{eq10}を用いてhを求めます。

def weight(x):
    diff = np.absolute(g_1st_deriv(x) - f_1st_deriv(x))
    if diff == 0:
        return 0
    return (1 / diff)

ERROR = 10e-20

def estimate_h(f, g, x):
    itr_count = 0
    h_estimated = [0]
    while True:
        cum = 0
        sigma_count = 0
        for sub_x in x:
            cum += ((g(sub_x) - f(sub_x + h_estimated[-1])) *
                    weight(sub_x) * f_1st_deriv(sub_x + h_estimated[-1]))
            sigma_count += weight(sub_x) * f_1st_deriv(sub_x + h_estimated[-1]) * f_1st_deriv(sub_x + h_estimated[-1])
        h_estimated.append(h_estimated[-1] + cum / sigma_count)

        if np.absolute(g(sub_x) - f(sub_x + h_estimated[-1])) < ERROR:
            break

        itr_count += 1
        if 100 == itr_count:
            break

こちらの関数を使うと、3回のイテレーションで目的のhを求められました。

LK法で解が求まる条件

LKで解が求められる限界について考えます。

上で例に挙げたFの場合、

{ \displaystyle
h \lt | \pi\ |
}

となる範囲でのみ解が正しく求められます。

試しに、先のプログラムでhを3.2に設定すると、値が発散して解が求められないことが確認できると思います。

LK法で解を求める場合、LPFを用いて低周波にしてからかけることで、解が収束する可能性が高くなります。

そこで、低周波の信号で解を求めた後に、その解を初期値として高周波の信号で解を求めるという戦略が考えられます。

LK法のパフォーマンス

LK法では、重みw(x)を設定することで、より正確に近似ができるということでした。

より正確な近似ができるということは、より少ない回数で解に近づけるということになります。こちらを実際にプログラムを動かして確かめます。

まず、重み関数w(x)を以下の様に変更します。

def weight(x):
    return 1

解が求めやすいと違いがわからないので、hを3、ERRORを10e-50に変更します。

すると、w(x)を1とした時のイテレーションは7回、w(x)を式\eqref{eq10}とした時のイテレーションは6回となり、たしかに計算回数が減っていることがわかります。

今回の記事は以上になります。次回は2次元へ拡張をする予定です。

LK法の式の導出 一般化

前回の記事では、LK法の概要と、1次元の信号についてのLK法の式の導出を解説しました。

m-yoshi-1700.hatenablog.com

この記事では、一般化したLK法の式を導出します。

この技術についての投稿は、一つの索引用の記事にまとめられています。 全体を俯瞰する場合はそちらを御覧ください。

m-yoshi-1700.hatenablog.com

最終的な目的は、次のLK法の式を導出することです。

LK法の定義

画像を与える関数をF(x)とする、その関数をhだけ動かしたものをG(x)とする。 hの初期値をh0と置いた時、hは以下の式で求められる。

{ \displaystyle
\def\vector#1{\mbox{\boldmath $#1$}}
\textbf{h} \approx
\left[ \sum_{\textbf{x}} ( \frac {\delta F} {\delta \textbf{x}} )^{\mathrm{T}}
\left[ G(\textbf{x}) - F(\textbf{x}) \right]
\right]
\left[ \sum_{\textbf{x}} ( \frac {\delta F} {\delta \textbf{x}} )^{\mathrm{T}} ( \frac {\delta F} {\delta \textbf{x}} ) \right] ^{-1}
}

一般化

前回の記事で導出した式は、1次元の信号について考えたもので、一般化されているわけではありません。

一般化が必要な理由は以下になります。

  • 2次元の導関数は1次元とは定義が違う。
  • F'(x)=0となる点でhが定義されていない。(F'(x)はhの近似式の分母になっています。)

まず、式(1)を次のように変形します。

{ \displaystyle
F'(\textbf{x}) \approx
\frac {F(x + h) - F(x)} {h}
= \frac {G(x) - F(x)} {h}
\tag{1}
}

{ \displaystyle
F(x + h) \approx F(x) + hF'(x)
\tag{8}
}

※ 式番号は論文と対応させています。

ここで、L2 normを最小化することを考えます。

{ \displaystyle
E = \sum_{x} \left[ F(x + h) - G(x) \right]^{2}
}

誤差を最小化するhを以下の様に置きます。

{ \displaystyle
0 = \frac {\delta E} {\delta h}\\
   \approx \frac {\delta} {\delta h} \sum_{x} \left[ F(x) + hF'(x) - G(x) \right]^{2}\\
   = \sum_{x} 2F'(x) \left[ F(x) + hF'(x) - G(x) \right]
}

この式を変形し、

{ \displaystyle
h \approx \frac {\sum_{x} F'(x) \left[ G(x) - F(x) \right]} {\sum_{x} F'(x)^{2}}
\tag{9}
}

この式は基本的には前の記事の式(6)と同じ意味で、w(x) = F'(x)2と置いたものになります。

また、この式は、前の式の式(3)と比較して、F'(x) = 0となる箇所があった場合にも対応しています。

式(3)の場合は、どこか1箇所でもF'(x)=0となる箇所があると成り立たなくなるのに対し、この式は、全てF'(x) = 0となるxの範囲以外は成立します。

式(9)を用いて、前の式(6)を書き直すと次のようになります。

{ \displaystyle
h_{0} = 0\\
h_{k + 1} = h_{k} +
\left[ \frac {\sum_{x} w(x) F'(x + h_{k})[G(x) - F(x + h_{k})]} {\sum_{x} w(x) F'(x + h_{k})^{2}} \right]
\tag{10}
}

こちらで導出は以上になります。

次回は1次元の信号についてのLK法で、パフォーマンスについて解説します。

LK法の式の導出

勉強会で学んだことを書きます。

この技術についての投稿は、一つの索引用の記事にまとめられています。 全体を俯瞰する場合はそちらを御覧ください。

m-yoshi-1700.hatenablog.com

LK法

画像の位置合わせ(Image Registration)の手法 CMUのBruce D. Lucas and Takeo Kanadeが開発

応用例

特徴

  • テンプレートマッチに比べて高速
  • Newton-Raphson法を用いた球解

Webサイト

KLT: An Implementation of the Kanade-Lucas-Tomasi Feature Tracker

論文

Bruce D. Lucas and Takeo Kanade. An Iterative Image Registration Technique with an Application to Stereo Vision. International Joint Conference on Artificial Intelligence, pages 674–679, 1981.

https://cecas.clemson.edu/~stb/klt/lucas_bruce_d_1981_1.pdf

内容の解説

画像を与える関数をF(x)とする、その関数をhだけ動かしたものをG(x)とする。 hの初期値をh0と置いた時、hは以下の式で求められる。

{ \displaystyle
\def\vector#1{\mbox{\boldmath $#1$}}
\textbf{h} \approx
\left[ \sum_{\textbf{x}} ( \frac {\delta F} {\delta \textbf{x}} )^{\mathrm{T}}
\left[ G(\textbf{x}) - F(\textbf{x}) \right]
\right]
\left[ \sum_{\textbf{x}} ( \frac {\delta F} {\delta \textbf{x}} )^{\mathrm{T}} ( \frac {\delta F} {\delta \textbf{x}} ) \right] ^{-1}
}

この式がLK法の基本になります。 この式を以下で導出していきます。

Image Registration

まずは問題の定義から。 画像をF(x)で与えられる関数とします。 ※ xは画像上の座標、F(x)は輝度値を返します。 対象画像をF(x)、テンプレートをG(x)とした時に、 F(x + h)とG(x)の違いを最小化するようなhを求めたいというのがImage Registrationの問題になります。

f:id:m-yoshi-1700:20170305165203p:plain

F(x + h)とG(x)の違いは以下の様に定義します。

  • L1ノルム

{ \displaystyle
L_{1} norm = \sum_{\textbf{x} \in R} | F(\textbf{x} + \textbf{h}) - G(\textbf{x}) |\
}

  • L2ノルム

{ \displaystyle
L_{2} norm = \left( \sum_{\textbf{x} \in R} \left[ F(\textbf{x} + \textbf{h}) - G(\textbf{x}) \right]^{2} \right) ^{\frac {1} {2}}
}

  • Negative of normalized correlation

{ \displaystyle
negative of normalized correlation = \frac
{- \sum_{\textbf{x} \in R} F(\textbf{x} + \textbf{h}) G(\textbf{x}) }
{
\left( \sum_{\textbf{x} \in R} \left[ F(\textbf{x} + \textbf{h})\right]^{2} \right) ^{\frac {1} {2}}
\left( \sum_{\textbf{x} \in R} \left[ G(\textbf{x})\right]^{2} \right) ^{\frac {1} {2}}
}
}

1次元の場合を考える

まずは画像(2次元)でなく、1次元の関数で考えてみます。 対象の関数をF(x)、そこからhだけ平行移動した関数をG(x)と置きます。

{ \displaystyle
G(x) = F(x + h)
}

ここで、hが十分に小さいと仮定して、F(x)の一次の導関数を考える。

{ \displaystyle
F'(\textbf{x}) \approx
\frac {F(x + h) - F(x)} {h}
= \frac {G(x) - F(x)} {h}
\tag{1}
}

この式を変形する。

{ \displaystyle
h \approx
\frac {G(x) - F(x)} {F'(\textbf{x})}
\tag{2}
}

hの式でのhの近似式はxに依存する、そこで一定の範囲について積分をする。

{ \displaystyle
h \approx
\left[ \sum_{x} \frac {G(x) - F(x)} {F'(\textbf{x})} \right]
\left[ \sum_{x} 1 \right] ^{-1}
\tag{3}
}

さて、ここでどのようなxがhを求めるのに有効かを考える。 xに対してF(x)の変化が線形であることが望ましいので、F'‘(x)の値が小さいほどhを求めるのに適したxと言える。 そこで、F’‘(x)の値を重みとして先ほどの近似の式に用いる。

{ \displaystyle
F^{\prime\prime}(x) \approx
\frac {G'(x) - F'(x)} {h}
\tag{4}
}

ここから、重みをF'‘(x)の逆数と置くと、

{ \displaystyle
w(x) =
\frac {1} { | G'(x) - F'(x) | }
\tag{5}
\label{5}
}

この重みを用いて、hの式を書き直します。

{ \displaystyle
h \approx
\left[ \sum_{x} \frac {w(x) [G(x) - F(x)]} {F'(\textbf{x})} \right]
\left[ \sum_{x} w(x) \right] ^{-1}
\tag{6}
\label{6}
}

ここで、w(x)は式\eqref{5}で定義されたものです。

いやー、はてなブログ便利ですね、式の参照もリンクにしてくれる。

式\eqref{6}を繰り返すことで真値に近づきます。

{ \displaystyle
h_{0} = 0\\
h_{k + 1} = h_{k} +
\left[ \sum_{x} \frac {w(x) [G(x) - F(x)]} {F'(\textbf{x})} \right]
\left[ \sum_{x} w(x) \right] ^{-1}
\tag{7}
}

これで1次元の場合のLK方の式の導出は完了です。 本日はここまで、一般化と2次元(画像)への拡張は次の記事で。