JP2007128237A - 最適平滑化スプラインによる極値検出方法及び極値検出プログラム - Google Patents

Abstract

【課題】離散データ点の局部のノイズの影響を過大に受けることなく極値を精度良く安定的に検出する。
【解決手段】データが入力される入力手段と、プログラムを格納する記憶手段と、プログラムを実行しデータを処理する演算手段と、演算手段の処理結果を出力する出力手段とを備える処理装置で実行される極値検出方法であって、入力手段から入力された離散データ群又は関数から、演算手段が最適平滑化スプラインを生成するステップと、演算手段が、最適平滑化スプライン又はその導関数あるいは偏導関数の各節点領域について各々極値を求めることにより、最適平滑化スプライン又はその導関数あるいは偏導関数の極値のすべてを求め、極値を出力手段へ出力するステップとを有する。
【選択図】図２

Description

本発明は、離散データ点や関数に関する極値検出方法及び極値検出プログラムに関し、特に、離散データ点や関数に関する最適平滑化スプラインを設計することにより、すべての極値を検出する極値検出方法及び極値検出プログラム、およびその極値検出方法及び極値検出プログラムを利用したデジタル画像のエッジ検出方法及びエッジ検出プログラムに関する。

離散データ点や関数に関する極値検出の問題は、工学を始めとする様々な分野で現れ、いわゆる「数理計画法」あるいは「最適化法」という分野の中で重要な問題として知られている。離散データ点や関数に関する極値検出の方法は、山登り法（山下り法）に代表されるような数値的な極値探索法が一般的であり、従って局所的な方法である。

しかしながら、このような局所的な方法による極値探索は、離散データ点の局部のノイズの影響を受けやすく、極値を精度良く安定的に検出することが困難である。したがって、離散データ点に関して大域的に評価することが可能な極値検出方法が求められている。

離散データ点に関して大域的に評価することが可能な極値検出方法に関しては、例えば、特許文献１に記載の方法がある。この方法によれば、１組のデータ点について大域的にサンプリングを行い、この１組のデータ点に基づいて単一のポリラインを生成し、このポリラインの極値を求める方法が開示されている。
特開２００５−２２２５５４号公報

特許文献１に記載の方法においては、曲線を作成する際に使用する「特定の点」や「分割点」の選択に、曲率（２次導関数）が用いられている。特に「特定の点」は，曲率が極値（極小値，極大値）をとる点あるいは０になる点として選ばれている。一方、「分割点」は曲率の積分値に基づいて選択されている。

しかしながら、曲率の極値や，曲率が０になる点の計算法は何ら具体的に示されていない。曲率の値を幾つかの点で計算し，得られた値から極大や極小，さらに０になる点を選んでいるのではないかと推測されるが、定かではない。

また、特許文献１に記載の方法は、曲線設計に関する高い技術を有していないユーザを対象として、ユーザの希望する曲線を描き、それを修正することを支援するものに過ぎず、離散データ点に関して大域的に評価することが可能ではない。

本発明は、このような極値に関する従来の問題も解決するもので、離散データ点や関数に関する最適平滑化スプラインを設計することにより、最適平滑化スプラインに対する極値を精度良く安定的に検出する方法を提供しようとするものである。また、すべての極値が求まることから、最大値や最小値も得られる方法を提供しようとするものである。

本発明は、データが入力される入力手段と、プログラムを格納する記憶手段と、プログラムを実行しデータを処理する演算手段と、演算手段の処理結果を出力する出力手段とを備える処理装置で実行される極値検出方法であって、入力手段から入力された平面における離散データ群又は関数から、演算手段が最適平滑化スプライン曲線を生成するステップと、演算手段が、前記最適平滑化スプライン曲線又はその導関数の各節点区間について各々極値を求めることにより、最適平滑化スプライン曲線又はその導関数の極値のすべてを求め、極値を出力手段へ出力するステップとを有することを特徴とする最適平滑化スプライン曲線による極値検出方法である。

また、本発明は、データが入力される入力手段と、プログラムを格納する記憶手段と、プログラムを実行しデータを処理する演算手段と、演算手段の処理結果を出力する出力手段とを備える処理装置で実行される極値検出方法であって、入力手段から入力された空間における離散データ群又は関数から演算手段が最適平滑化スプライン曲面を生成するステップと、演算手段が最適平滑化スプライン曲面又はその偏導関数の各節点領域について各々極値を求めることにより、最適平滑化スプライン曲面又はその偏導関数の極値のすべてを求め、極値を前記出力手段へ出力するステップとを有することを特徴とする最適平滑化スプライン曲面による極値検出方法である。

また、本発明は、ピクセルデータが入力される入力手段と、プログラムを格納する記憶手段と、プログラムを実行しピクセルデータを処理する演算手段と、演算手段の処理結果を出力する出力手段とを備える画像処理装置で実行されるデジタル画像のエッジ検出方法であって、入力手段から入力されたピクセルの座標及び輝度値からなるピクセルデータから、演算手段が最適平滑化スプライン曲面を生成するステップと、演算手段が最適平滑化スプライン曲面を第１の座標軸上の複数の座標で固定した複数の曲線からなる第１の曲線群を生成し、該第１の曲線群を構成する複数の曲線の１次導関数の各節点区間について各々極値を求めることにより、第１の曲線群の１次導関数の極値のすべてを求めるステップと、演算手段が最適平滑化スプライン曲面を第２の座標軸上の複数の座標で固定した複数の曲線からなる第２の曲線群を生成し、第２の曲線群を構成する複数の曲線の１次導関数の各節点区間について各々極値を求めることにより、第２の曲線群の１次導関数の極値のすべてを求めるステップと、演算手段が第１の曲線群の１次導関数の極値又は第２の曲線群の１次導関数の極値の絶対値が所定のしきい値より大きい箇所を、デジタル画像のエッジ部分と判定し、判定結果を出力手段へ出力するステップとを有することを特徴とするデジタル画像のエッジ検出方法である。

また、本発明は、データが入力される入力手段と、プログラムを格納する記憶手段と、プログラムを実行しデータを処理する演算手段と、演算手段の処理結果を出力する出力手段とを備える処理装置に、入力手段から入力された平面における離散データ群又は関数から、演算手段が最適平滑化スプライン曲線を生成するステップと、演算手段が最適平滑化スプライン曲線又はその導関数の各節点区間について各々極値を求めることにより、最適平滑化スプライン曲線又はその導関数の極値のすべてを求め、極値を出力手段に出力するステップとを実行させるための最適平滑化スプライン曲線による極値検出プログラムである。

また、本発明は、データが入力される入力手段と、プログラムを格納する記憶手段と、プログラムを実行しデータを処理する演算手段と、演算手段の処理結果を出力する出力手段とを備える処理装置に、入力手段から入力された空間における離散データ群又は関数から、演算手段が最適平滑化スプライン曲面を生成するステップと、演算手段が最適平滑化スプライン曲面又はその偏導関数の各節点領域について各々極値を求めることにより、最適平滑化スプライン曲面又はその偏導関数の極値のすべてを求め、極値を出力手段に出力するステップとを実行させるための最適平滑化スプライン曲面による極値検出プログラムである。

また、本発明は、ピクセルデータが入力される入力手段と、プログラムを格納する記憶手段と、プログラムを実行しピクセルデータを処理する演算手段と、演算手段の処理結果を出力する出力手段とを備える画像処理装置に、入力手段から入力されたピクセルの座標及び輝度値からなるピクセルデータから演算手段が最適平滑化スプライン曲面を生成するステップと、演算手段が最適平滑化スプライン曲面を第１の座標軸上の複数の座標で固定した複数の曲線からなる第１の曲線群を生成し、第１の曲線群を構成する複数の曲線の１次導関数の各節点区間について各々極値を求めることにより、第１の曲線群の１次導関数の極値のすべてを求めるステップと、演算手段が最適平滑化スプライン曲面を第２の座標軸上の複数の座標で固定した複数の曲線からなる第２の曲線群を生成し、第２の曲線群を構成する複数の曲線の１次導関数の各節点区間について各々極値を求めることにより、第２の曲線群の１次導関数の極値のすべてを求めるステップと、演算手段が第１の曲線群の１次導関数の極値又は第２の曲線群の１次導関数の極値の絶対値が所定のしきい値より大きい箇所を、デジタル画像のエッジ部分と判定し、判定結果を出力手段へ出力するステップとを実行させるためのデジタル画像のエッジ検出プログラムである。

本発明の最適平滑化スプラインによる極値検出方法によれば、設計されたスプラインを対象とするため、各節点ごとに定義された多項式の極値（極小値、極大値）を求めれば良い。これによりすべての極値が求まり、最小値や最大値を求めることが可能となる。

（第１実施形態）
以下、本発明の第１実施形態である最適平滑化スプライン曲線による極値検出方法及び極値検出プログラムについて、図を参照して詳細に説明をする。

よく知られるように、スプライン関数は様々な工学分野で用いられており、これまで幅広い研究がなされてきている。特に、最適制御理論に基づいた「動的スプライン」の理論が提案されており、そのようなアプローチは、Ｂスプライン関数の解析や、最適平滑化曲線の設計に利用されている。さらには、最適内挿および平滑化スプラインの枠組みが日本の書道で見られるような文字や文字列の設計および再設計を行うために利用されている。

本実施形態では、まず、与えられた離散データの集合に対して最適平滑化スプライン曲線の設計を行う。そのような曲線の設計には、正規化された一様なＢスプライン関数を基底関数として用いる。まず、最適平滑化曲線の設計について、簡潔な計算手順が得られる基本的な結果を示す。次に、最適曲線のすべての極値、すなわち極小値、極大値を計算するためのアルゴリズムを示す。その際、スプラインが区分的多項式であり、連続かつ微分可能であるということを利用するものである。スプライン曲線は、区分的多項式であり制御点によって表すことができるので、その極値を検出し計算するアルゴリズムは簡潔かつ数値計算に対しても利用しやすいものとなっている。

B_k(・)を、整数値の節点をもつk次の正規化された一様なＢスプライン関数とし、一定間隔でシフトし重みづけしたＢスプライン関数の和としてスプライン曲線x(t)を生成する。すなわち、ｋ次スプライン曲線x(t)は、

のように表現される。

ここで、mは、整数であって、予め定めうる設計パラメータである。τ_i∈Rは制御点と呼ばれる重み係数である。α(>0)は、予め定めうる設計パラメータであり、等間隔に配置された節点t_iの間隔（節点区間）を、t_i+1−t_i=1/αに従って調整するための定数である。

数式１によってスプライン曲線x(t)が構成される様子を図１に示す。各つり鐘状の波形は適宜シフトされたｋ次ＢスプラインB_k(・)に対応する制御点τ_i（「×」印で示す）によって重み付けされた関数を表し，スプライン曲線x(t)はこれらをすべて加え合わせて構成されている。制御点τ_i∈Rを適切に選ぶことによって，任意のｋ次スプライン曲線x(t)を区間[t₀,t_m]で設計することができる。

ｋ次Ｂ−スプライン関数B_k(t)は、次式

で定義されている。B_k(t)は整数値の節点をもち，j=0，…,kに対する各節点区間[j,j+1]で基底関数N_j，k(t−j)によって指定され，[0,k+1]以外の区間では恒等的に0である。

ここで基底関数N_j，k(t)は、

によって，初期関数N_0，0(t)=1から順次生成することができる。

実際上、k=3、すなわち３次のスプライン関数を選定した場合、基底関数N_j，3(t)は次式

のように、４つ（すなわちk+1個）の関数を含むことになる。

このことは、各tにおけるｋ次スプライン関数の値は、ｔが含まれる区間におけるk+1個の基底関数（従って対応するk+1個の制御点）の影響を受けるだけで、その範囲外の基底関数（従って対応する制御点）の影響を受けないことを意味している。

（最適化過程）
次に、離散データＤの集合が

によって与えられているとし、また、τ∈R^M(M=m+k)を以下のように定義される制御点ベクトルとする。

そのとき、最適平滑化スプライン曲線の設計の基本問題は、評価関数

を最小にするスプライン曲線x(t)もしくは等価的に制御点ベクトルτ∈R^Mを見つける問題となる。ただし、λ>0、w_i∈(0，1]∀iである。

この問題に対する最適解τは以下の線形方程式の解として得られる。

ここで、Q∈R^M×M は

によって定義されるグラミアンで、与えられる離散データＤとは無関係に予め計算しておくことができる。

またW∈R^N×N、d∈R^NおよびB∈R^M×Nは

である。

ただし数式９，１２におけるb(t)∈R^Mは、

なるベクトルである。

数式８は可解であり、従って解τは常に存在する。もちろんλQ+BWB^T>0のときに限り唯一となる。この解τを数式１に用いることによって、区間[t₀，t_m]におけるｋ次の最適平滑化スプラインx(t)が設計される。

このように、制御点ベクトルτに関する最適化解を求め、数式１を用いることによって、ｋ次最適平滑化スプライン曲線x(t)が生成される。

（スプライン曲線の極値）
次に、最適平滑化スプライン曲線x(t)の極値を検出する方法を説明する。具体的には、x⁽¹⁾(t)=0とx⁽²⁾(t)>0もしくはx⁽²⁾(t)<0を満たしている点t∈(t₀，t_m)を見つければ良い。

区間[t₀，t_m]で設計されたk次最適平滑化スプラインx(t)は、各区間[t_j，t_j+1)(k=0,1，2，・・・，m−1)で各々k次多項式である。従って、区間毎に極値を検出することによって、x(t)の最大値、最小値を含むすべての極値が求まる。

x(t)は区間[t_j，t_j+1)において、

となり、ｋ＋１個の制御点τ_j-k、τ_j-k+1、・・・、τ_jによって決まるｋ次多項式である。

ｔからδへの変数変換を次のように行う。

このとき変数ｔの区間[t_j，t_j+1)は変数δに対する単位区間[0，1)に正規化され、関数x(t)はδの関数

に変換される。x(δ)はδに関するｋ次多項式であり、これを

とする。

なお、以下に使用される数式に関して、δを変数とする関数xには、すべて‘ハット（＾）’がつくものであるが、明細書本文中では省略するものとする。

明らかにx(t)＝x(δ)、またx⁽¹⁾(t)=αx ⁽¹⁾(δ)、x⁽²⁾(t)=α²x⁽²⁾(δ)かつα＞０より、x(t)の区間[t_j，t_j+1)における極値はx(δ)の単位区間[0，1)における極値として得ることができる。x(δ)がδ＝δ_ｉ（０≦δ_ｉ＜１）で極値x(δ_ｉ)を持てばx(t)はｔ＝（１／α）δ_ｉ＋ｔ_ｊで極値x(t)＝x(δ_ｉ)を持つ。

x(δ)の単位区間[0，1)における極値の求め方は、後述する図３のステップ１０３〜ステップ１０５の通りである。その際、１次導関数x^(１)(δ)の根の計算は、既存の多項式の根の計算プログラムを用いることができる。

（３次スプライン）
実際に良く用いられている３次スプライン（すなわちｋ＝３）の場合には、以下の仮定１を置くことによって、上のステップはより詳細なアルゴリズムとして書くことができる。これはx^(１)(δ)が２次多項式となり、その根が陽的に表現できるからである。

ここでは、実用的な目的上、次の仮定を導入する。

（仮定１）各節点区間[t_j，t_j+1)のx(t)の極値の個数は、多くても1つである。ただし、j=0, 1，2，・・・，m−1

x(t)が、最適平滑化スプライン曲線であることに注意すると、節点区間[t_j，t_j+1)が小さければ仮定１は妥当なものである。

３次スプライン（すなわちｋ＝３）の場合には、多項式x(δ)は４つの重み係数τ_j-3，τ_j-2，τ_j-1，τ_jだけに依存し、δの変域[0，1)に対して次式のように表現することができる。

すなわち、係数p_j、q_j、r_j、s_jは、上記最適化過程ですでに得られている重み係数（ベクトル）τのうち、τ_j-3、τ_j-2、τ_j-1、τ_jの４つから唯一に決定される。

前述のように、x(δ)がδ=δ_e∈[0，1)で極値をもつことがわかれば、そのときx(t)は、t_e=t_j+δ_e(1/α)∈[t_j，t_j+1)で極値をもち、その値は、x(t_e)= x(δ_e)により得られることとなる。

δ_eを得るためには、数式１８のx(δ)の1次導関数、すなわち
が重要な役割をはたす。p_j≠0である場合には、x⁽¹⁾(δ)の2つの根は、

のように表すことができる。そのとき、x⁽¹⁾(0)=r_j/2およびx⁽¹⁾(1)=r_j+1/2となるので、極値が存在するかどうかはr_jとr_j+1の符号により調べることができる。具体的には、仮定１のもとでの幾何学的な観察により、以下の3つの場合が得られる。

(P1) r_j・r_j+1<0であれば、極値は[0，1)に存在する。

(P2) r_j・r_j+1=0であれば、極値はいくつかの付加的な条件のもとで[0，1)に存在する。

(P3) r_j・r_j+1>0であれば、極値は[0，1)には存在しない。

(P1)と(P2)の場合について、以下の結果が得られる。

（命題１）r_j・r_j+1<0と仮定する。そのとき、r_j>0かつr_j+1<0（r_j<0かつr_j+1>0）ならば、関数x(δ)はδ=δ_e∈(0，1)で極大(極小)をもつ。ただし、

（命題２）r_j・r_j+1=0と仮定するとき、以下の２つの場合が生じる。

(i)r_j=0の場合: q_j<0 (q_j>0)であれば、そのときx(δ)はδ=0で極大(極小)をもつ。q_j=0である場合には、極値は存在しない。

(ii)r_j≠0かつr_j+1=0の場合: r_j>0かつq_j+1>0(r_j<0かつq_j+1<0)であれば、そのときx(δ)はδ=r_j/p_j∈(0，1)で極大(極小)をもつ。また、r_j・q_j+1≦0であれば極値は存在しない。

したがって、３次最適平滑化スプライン関数x(t)の極値を検出し計算するアルゴリズムは、以下のようにまとめることができる。

重み係数τ_-3，τ_-2，τ_-1，・・・，τ_m-1が与えられ、r_jとr_j+1を、j=0, 1，2，・・・，m-1に対して計算し、そのときr_j・r_j+1<0およびr_j・r_j+1=0の場合には、それぞれ命題１および２を使う。δ=δ_eで極値x(δ)が見つかれば、t_e=t_j+(1/α)δ_eで極値x(t_e)を計算する。

図２は、本実施形態のｋ次最適平滑化スプライン曲線による極値検出方法及び極値検出プログラムに用いられる処理装置を示すブロック図である。本実施形態の最適平滑化スプラインによる極値検出方法に用いられる処理装置１は、ＣＰＵ５（演算手段）とメモリ装置４（ＲＡＭやＲＯＭ等：記憶手段）を有し、このメモリ装置４に格納されたプログラムを実行し、図３、４のフローチャートに記載した動作を実行することにより、図２のブロック図に示すように、ＣＰＵ５と、メモリ装置４とマウス、キーボード、スキャナ、データロガー等の入力装置６（入力手段）とバス３と、ディスプレイ装置７（出力手段）と、プリンタ装置８（出力手段）と、外部ネットワーク１０と接続するネットワーク装置９（入力手段・出力手段）を備えた処理装置として機能する。

図３は、本実施形態の一般のｋ次最適平滑化スプラインによる極値検出方法の処理の流れを示すフローチャートである。

ステップ１０１（図ではステップをＳと略す）は、離散データＤについて重み係数τ_iを求め入力する工程である。すなわち、ＣＰＵ５は、入力した離散データＤ（u_i，d_i），i=1，・・・，Nについて、最適な重み係数τ_i，i=−k, −k+1，・・・，m−1を決定し、最適平滑化スプライン曲線x(t)を設計する。

ステップ１０２では、ＣＰＵ５は、節点区間[t_j，t_j+1）での最適平滑化スプライン曲線x(t)を正規化して得られる、単位区間[0，1 ）でのk次多項式x(δ)を生成する。すなわち、x(δ)の各係数a₀，・・・，a_kの計算を行う。

x(δ)は、数式１６で定義されている。ここで、基底関数（多項式）N_0，k(δ)、N_1，k(δ)、・・・、N_h，k(δ)は、数式３の逐次式から計算機を用いて予め求めておくことができ、いずれもｋ次多項式である。いま、

とし、求めるx(δ)の多項式表現を数式１７とすると、その係数a₀，・・・，a_kは数式１７、数式２２によって次式から求められる。

ステップ１０３では、ＣＰＵ５は、ステップ１０２で求めた多項式x(δ)の１次導関数x⁽¹⁾(δ)の根を求める。

ステップ１０４では、ＣＰＵ５は、多項式x(δ)が単位区間[0，1)で実根を有するか否か判定する。実根を有する場合にはステップ１０５へ進み、実根を有しない場合にはステップ１０６へ進む。

ステップ１０５では、ＣＰＵ５は、x(δ)の２次導関数x ⁽²⁾(δ)へ、ステップ１０３で求めた多項式x(δ)の実根δ_eを代入し、x⁽²⁾(δ_e)＞0の場合には、t_e=t_j+(1/α)δ_eで極小値x(δ_e)をもつものとし、x⁽²⁾(δ_e)＜0の場合には、t_e=t_j+(1/α)δ_eで極大値x(δ_e)をもつものとする。これにより、最適平滑化スプライン関数x(t)の所定の節点区間[t_j，t_j+1)における、極大値または極小値が求まることとなる。

ステップ１０６では、jが、m−1より小さい場合にはステップ１０２へ戻り、jが、m−1の場合にはステップ１０７へ進み、区間[t₀，t_m]のすべての節点区間[t_j，t_j+1)についての極値検出手続きが終了するまで、ステップ１０２からステップ１０５を繰り返す。

ステップ１０７では、ＣＰＵ５は、最適平滑化スプライン関数x(t)の全区間[t₀，t_m]に存在するすべての極値から、すべての極値の大小を判別することにより、最適平滑化スプライン関数x(t)の最大値及び最小値を求める。そして、最適平滑化スプライン関数x(t)の極大値、極小値、最大値及び最小値を、ディスプレイ装置７やプリンタ装置８に出力する。

図４は、本実施形態の最適平滑化スプラインによる極値検出方法の処理の流れにおいて、３次スプライン（ｋ＝３）を用いた場合を示すフローチャートである。なお、図３のフローチャートとはステップ１０２〜１０４のみ相違するため、対応する部分を、ステップ１０２ａ及び１０３ａとして記載している。

ステップ１０２ａは、図３のフローチャートのステップ１０２に対応する部分である。ステップ１０２ａでは、ＣＰＵ５は、３次多項式であるx(δ)について、各係数p_j、q_j、r_j、s_jを算出する。上述したように、各係数p_j、q_j、r_j、s_jは、τ_j-3、τ_j-2、τ_j-1、τ_jの４つから唯一に決定される。

ステップ１０３ａは、図３のフローチャートのステップ１０３及び１０４に対応する部分である。ステップ１０３ａでは、ＣＰＵ５は、ステップ１０２ａで求めた多項式x(δ)の１次導関数x⁽¹⁾(δ)の根を求める。q_j ²−p_jr_j≧0の場合には図３のステップ１０５ヘ進み、q_j ²−p_jr_j＜0の場合には図３のステップ１０６へ進む。

なお、本実施形態においては、離散データＤが与えられた場合について極値をもとめたが、関数f(t)，t₀≦t≦t_mが入力データとして与えられた場合についても、同様な方法で極値を求めることが可能である。この場合にも、まずf(t)を近似する最適スプラインを以下のように設計する。

関数f(t)が与えられたとき、評価関数は、数式７の代わりに

を用いる。この最適解τは線形方程式

の解として得られる。ここで、Ｑは、数式９で定義される行列、Ｑ_０は

で定義される行列で、いずれもf(t)とは無関係に予め計算しておくことができる。右辺は、f(t)が与えられると適当な数値積分法を用いて計算できる。

一旦、最適なベクトルτが決まると、極値の検出等、後の扱いは離散データの場合と同様である。ただし、既存の最適化アルゴリズムとの併用によって以下のような有効性がある。一般に関数f(t)が与えられたとき、ある区間（一般性を失わずに[t₀，t_m]と考えてよい）での最大値や最小値を求める方法は繰り返しによる数値的探索法によるのが通常である。

その際、繰り返し式の初期値の取り方が問題となり、それによって多峰性関数の場合は到達する極値が異なる。最大値や最小値が得られる保証は全くない。本方法によれば、f(t)を近似する最適スプラインが得られ、その最大値や最小値が得られる。さらに対応するｔの値を初期値として数値的探索法を適用すればf(t)の最大値や最小値をより正確に求めることができる。

（数値実験例）
図５は、本発明の第１実施形態の最適平滑化スプライン曲線とその極値を示す図である。図５は、最適平滑化スプライン曲線x(t)（図中の実線で示された関数）から極値を抽出する性能を示す図である。図５は、k=3，α=1， t₀=1およびt_m=m=50の場合での結果を示している。ただし、三角形と逆三角形はそれぞれ検出された極大値および極小値を示す。ここで、データd_i（図中の＊印）は、関数f(t)= e^a(t−1)cos(b(t−1))+1（図中の点線で示された関数）をサンプリングすることにより生成されている。ただし、a=1/(m−1) log1/4， b=2π×3/(m−1)である。 データの個数はN=30とし、u_iは区間[t₀，t_m]=(1，50)においてランダムに配置されたデータ点であり、またデータに加えるノイズの大きさはσ=0.02とする。さらに、平滑化パラメータλを推定するためにいわゆるcross-validation手法を採用し、その最適値はλ^★=0.1585のように得られた。

図からわかるように、データノイズがかなり大きいにもかかわらず、元の関数であるf(t)のすべての極値がかなり精度よく検出されていることがわかる。特に極値を与えるｔの値はほぼ正確に得られていることがわかる。

なお、上記実施形態では、最適平滑化スプライン関数x(t)の極値を求めたが、同様な方法で、x(t)の第ｉ導関数x⁽ⁱ⁾(t)の極値を求めることが可能である。

区間[t₁，t_m]で設計された最適スプラインx(t)の第ｉ導関数x⁽ⁱ⁾(t)(i=1，2，・・・，k)も各節点区間[t_j，t_j+1) (k=0, 1，・・・，m−1)毎に考えれば良い。すなわち、節点区間[t_j，t_j+1 )においては数式１４と数式１５の変数変換によって、x⁽ⁱ⁾(t)=αⁱx⁽ⁱ⁾(δ)の関係が成立する。またα＞０よりx⁽ⁱ⁾(t)の節点区間[t_j，t_j+1 ) における極値は、x ⁽ⁱ⁾(δ)の単位区間[0，1 )における極値となる。この第ｉ導関数x ⁽ⁱ⁾(δ)は、数式１７と数式２３の結果を用いて、

によって計算できる。

x⁽ⁱ⁾(δ)（ただしi≦k−2）の極値の検出は、図３のステップ１０２における多項式x(δ)をx⁽ⁱ⁾(δ)に置き換え、これに伴いx⁽¹⁾(δ)， x⁽²⁾(δ)をそれぞれx⁽ⁱ⁺ⁱ⁾(δ)、x⁽ⁱ⁺²⁾(δ) に置き換えて、ステップ１０３以降に同様な方法を適用すれば良い。ここで、i≦k−2に限定されるのは、k−(i+2)次の多項式であるi+2次導関数が必要であるからで、k−(i+2)≧0が必要なためである。

以上説明したように、本実施形態の最適平滑化スプラインによる極値検出方法によれば、設計されたスプラインを対象とするため、各節点区間ごとに定義された多項式の極値（極小値、極大値）を求めれば良い。これによりすべての極値が求まり、最小値や最大値を求めることが可能となる。

また、基底関数としてＢスプライン関数を用いているため、設計されたスプラインの表現が簡潔であり、その導関数も容易に求めることができる。

そして、極値計算の対象は、離散データでも関数でも良い。いずれの場合にも最適平滑化スプラインが設計でき、その際、節点の間隔やＢスプラインの次数を調整することによって十分な精度で近似でき、ノイズの影響を受けにくい数値的に非常に安定した方法といえる。もちろん、この極値検出法は最適平滑化スプライン曲線に限らず、数式１で表現されるスプライン曲線であれば適用できるものである。

特に、３次Ｂスプライン関数を用いた場合には、各区間毎に３次多項式となり、極値は解析的に得ることができる。

（第２実施形態）
以下、本発明の第２実施形態である最適平滑化スプライン曲面による極値検出方法及び極値検出プログラムについて、図を参照して詳細に説明をする。

（最適平滑化スプライン曲面）
本発明の第１実施形態では、平面における離散データについて最適平滑化スプライン曲線を生成し極値の検出を行ったが、本実施形態は、空間における離散データを近似する最適平滑化スプライン曲面を生成し極値の検出を行うものである。また、本実施形態は、本発明の第１実施形態の最適平滑化スプライン曲線の極値検出方法を最適平滑化スプライン曲面に素直に拡張した方法を提供するものである。

図６は、本実施形態の最適平滑化スプライン曲面を概略的に示す図である。図６Ａは、空間における離散データから求められた最適な制御点（重み係数）の分布を示す図であり、図６Ｂは、空間における離散データに基づいて生成された最適平滑化スプライン曲面x(s，t)を示す図である。

（曲面の生成）
本実施形態の最適平滑化スプライン曲面x(s，t)は、ｋ次Bスプライン関数B_k(・)を用いて、数式１の自然な拡張である以下の数式２８で表すことが可能である。

ここで、α、β（＞0）は定数、m_１、m_２（＞2）は整数、s_i、t_jは以下のように等間隔に配置された節点を表す。

このとき、制御点τ_ｉ，ｊを適切に選ぶことによって、任意のｋ次スプライン曲面を（ｓ，ｔ）の領域

上で設計することができる。この制御点をまとめて行列τ∈Ｒ^M1×M2（ただしM₁=m₁+k、M₂=m₂+k）で表すことにする。

いま、空間内に離散データが

のように与えられたとする。このとき最適平滑化スプライン曲面を、下記の評価関数J(τ)が最小となるように設計する。

この問題に対する解、すなわち最適な制御点は次の線形方程式の解として得られる。

以下で数式３４において現れる各行列、ベクトルを示す。まず、Qは次式で与えられる M₁M₂×M₁M₂の行列である。

ここで、Q_l ^(ij)∈Ｒ^M1×M2は

によって定義される。ただし、ベクトルb₁(t)∈Ｒ^M1、b₂(t)∈Ｒ^M2 は、Bスプライン関数から構成されている。

数式３６に示す各行列、従ってQは与えられるデータとは独立にBスプラインから予め計算しておくことができる。

行列W∈Ｒ^N×Nは、数式３３に示す評価関数における重みから構成される対角行列である

行列B∈Ｒ^M1M2×N、 およびベクトルd∈Ｒ^Nは与えられるデータによって以下のように定められる。

以上によって、数式３４に現れる各行列、ベクトルが定まる。この方程式を解くことによって最適なτ、従って最適な制御点行列τが求められ、最適平滑化曲面x(s，t)が数式２８によって構成できる。

（曲面の極値）
次に、上記の方法により求めた最適平滑化曲面x(s，t)の極値を求める方法について説明する。

（対象とする曲面）
数式２８のｋ次スプライン曲面の極値を求める。ここでパラメータｋ、α、β、ｍ_１、ｍ_２は設定済みで，制御点τ_i，jは前述の方法等によってすでに与えられているとする。なお、(s，t) の定義域はS=[s₀，s_m1]×[t₀，t_m2]⊂R²である。

（極値）
曲面x(s，t)の極値は

をみたす点(s，t)∈Sとして得られる。

曲面x(s，t)のヘッセ行列が正定値となる場合が極小値、曲面x(s，t)のヘッセ行列が負定値となる場合が極大値を与える。x(s，t)の領域Sでの極値を求める問題は、節点間の部分領域S_κ，μ=[s_κ，s_κ+1)×[t_κ，t_κ+1 )における極値を求める問題を全領域κ=0，1，・・・，m₁−1；μ=0，1，・・・，m₂−1に渡って繰り返せば良い。

（部分領域の正規化）
領域S_κ，μに注目し、u=α(s−s_κ)，v=β(t−t_μ)なる変数変換を行うと、u，vの変域はいずれも[0，1)と正規化される。すなわち、(s，t)の領域S_κ，μは、κ、μにかかわらず、(u，v)の単位領域ε=[0，1)×[0，1)になる。このとき、S_κ，μにおける関数x(s，t)は、

と表すことができる。なお、以下に使用される数式に関して、ｕ、ｖを変数とする関数xには、すべて‘ハット（＾）’がつくものであるが、明細書本文中では省略するものとする。

さらにx(s，t)の偏導関数とx(u，v)の偏導関数の間の関係は

となる。α、β＞０より、x(s，t)の領域S_κ，μ上での極値はx(u，v)の領域ε上での極値と一致する。

（部分領域での極値）
数式４０に示す、曲面x(s，t)の勾配ベクトルが０、かつ、曲面x(s，t)のヘッセ行列が正定値となるか、又は、曲面x(s，t)のヘッセ行列が負定値となる点(u^＊，v^＊)∈εを検出すればよい。このときx(s，t)は、点(s^＊，t^＊)=((1/α) u^＊+s_κ，(1/β) v^＊+t_μ)で極値x(s^＊，t^＊)=x(u^＊，v^＊)をもつ。

（x(u，v)の極値の求め方）
（関数x(u，v)の表現）
x(u，v)はu，vについてそれぞれｋ次の２変数多項式であり、以下のように簡潔に表すことができる。

(k+1)×(k+1)行列Aは

で定義される。また、S，T_κ，μは次のようにして得られる(k+1)×(k+1)行列である。

行列Sは、その各行がBスプラインの基底関数N_i，k(t)，i=0，1，・・・，kの各係数から構成される行列である。すなわち、

としたとき、N_k(t)=Sh_k+1(t)となる行列Sであり、例えばk=3の場合、数式４よりSはつぎのようになる。

行列T_κ，μは、曲面x(u，v)の定義式に現れる(k+1)²個の制御点から次のように構成される。

すなわち、曲面x(s，t)を構成する全体の制御点行列τ（数式３１）のうちの、要素T_κ，μから左上の(k+1)×(k+1)のブロックとして得られるτの部分行列である。

（勾配ベクトルとヘッセ行列）
x(u，v)の勾配ベクトルおよびヘッセ行列は以下のように表現される。

ただしC_kは次のように定義される(k+1)×k行列である。

（極値の候補点）
曲面x(u，v)の勾配ベクトルが０となる点(u，v)∈εを求める。曲面x(u，v)の勾配ベクトルの各要素は、次のように与えられる。

ただし、

は次式で定義されている。

x_u(u，v)=0、x_v(u，v)=0が同時に成り立つためには、これらをuに関する多項式とみなしたときの終結式が０となる必要がある。すなわち、

ただしR(v)は以下の(2k−1)×(2k−1)行列である。

｜R(v)｜はvに関する多項式であり、この多項式の根を求め、0≦v＜1をみたす根をv=v^＊とする。なお、このような根がない場合にはx(u，v)は、単位領域εに極値をもたない。このv^＊に対応するuは、uに関する連立方程式x_u(u，v^＊)=0、x_v(u，v^＊)=0の解として求められる。具体的には、q_i（バー）=q_i(v^＊)、r_i（バー）=r_i(v^＊)として得られる２つの多項式

の共通因子f₀(u)を求める。その際、例えばユークリッド互助法を適用することができる。なお、f₁(u)、f₂(u)が共通因子をもつことは保証されており、もちろん、f₀(u)はuに関する次数１以上の多項式である。f₀(u)の根uが0≦u＜1をみたすとき、u=u^＊とおく。みたさないときはx(u，v)は領域εに極値をもたない。

（極値の判定）
得られた点(u^＊，v^＊)∈εが関数x(u，v)の極値の候補となる点であり、この点に対するヘッセ行列が正定（極小）あるいは負定（極大）の極値条件をみたすかどうか判定すればよい。なお、ヘッセ行列の計算の際、次のように定義される多項式

を導入し、さらにq(v)、r(v)を用いることによって得られる次式を用いるのがよい。

図７は、本実施形態の最適平滑化スプラインによる極値検出方法及び極値検出プログラムに用いられる処理装置を示すブロック図である。本実施形態の極値検出処理装置１１は、ＣＰＵ１５（演算手段）とメモリ装置１４（ＲＡＭやＲＯＭ等：記憶手段）を有し、このメモリ装置１４に格納されたプログラムを実行し、図８〜１１のフローチャートに記載した動作を実行することにより、図７のブロック図に示すように、ＣＰＵ１５（演算手段）と、メモリ装置１４（記憶手段）と、マウス、キーボード、スキャナ、データロガー等の入力装置１６（入力手段）と、バス１３と、ディスプレイ装置１７（出力手段）と、プリンタ装置１８（出力手段）と、外部ネットワーク２０と接続するネットワーク装置１９（入力手段・出力手段）を備えた処理装置として機能する。

図８は、本実施形態の最適平滑化スプラインによる極値検出方法及び極値検出プログラムの処理の概要を示すフローチャートである。図９は、図８のフローチャートのステップ２０３の処理の詳細を示すフローチャートであり、図１０は、図８のフローチャートのステップ２０６の処理の詳細を示すフローチャートであり、図１１は、図８のフローチャートのステップ２１１の処理の詳細を示すフローチャートである。

ステップ２０１では、極値検出処理装置１１に曲面の最適設計に用いられたパラメータκ、α、β、ｍ_１、ｍ_２および得られた最適制御点τ_i，j、i=−k，−k+1，・・・，m₁−1、j=−k，−k+1，・・・，m₂−1を入力する。

ステップ２０２では、ＣＰＵ１５は、(k+1)×(k+1)行列S、(k+1)×k行列C_kおよびk×(k−1)行列C_k-1を準備する。

以下、すべての(κ，μ)、κ=1，2，・・・，m₁−1、μ=1，2，・・・，m₂−1、の組み合わせに対して、以下の手順に従って領域S_κ，μでのx(s，t)の極値の検出、計算を行う。

ＣＰＵ１５は、(k+1)×(k+1)行列A（ステップ２０３ａ）、多項式q_i(v)，i=1，2，・・・，k、 多項式r_i(v)，i=1，2，・・・，k+1（ステップ２０３ｂ）、および終結式｜R(v)｜を求め（ステップ２０３ｃ）、vの多項式である終結式の根を求める（ステップ２０３）。この終結式｜R(v)｜の計算には、数式処理を援用することができる。0≦v＜1をみたす根vが存在しない場合には、次の(κ，μ)の組に移る（ステップ２０４）。

ＣＰＵ１５は、得られたv(0≦v＜1)をv^＊とおき（ステップ２０５）、さらに、q_i（バー）=q_i(v^＊)，i=1，2，・・・，k、およびr_i（バー）=r_i(v^＊)，i=1，2，・・・，k+1を計算する。それぞれq_i（バー）、r_i（バー）を係数とするuの多項式f₁(u)、f₂(u)を定め（ステップ２０６ａ）、例えば、ユークリッド互助法を適用するなどして、f₁(u)とf₂(u)の共通因子f₀(u)を求める（ステップ２０６ｂ）。f₀(u)の根u（ステップ２０６ｃ）が、0≦u＜1をみたすとき（ステップ２０７）、u=u^＊とおく（ステップ２０８）。そのような根uが存在しない場合には、次の(κ，μ)の組に移る。

ＣＰＵ１５は、点(u^＊，v^＊)におけるヘッセ行列を計算し（ステップ２０９）、点(u^＊，v^＊)におけるヘッセ行列が正定となる場合には極小、点(u^＊，v^＊)におけるヘッセ行列が負定となる場合には極大と判断し（ステップ２１０）、対応する極値x(u^＊，v^＊)を計算する（ステップ２１１ａ）。これら両ケース以外の場合は、次の(κ，μ)の組に移る。

ＣＰＵ１５は、x(s，t)の極値をとる点(s^＊，t^＊)を、(s^＊，t^＊)=((1/α)u^＊+s_κ，(1/β)v^＊+t_μ)で計算し（ステップ２１１ｂ）、その値x(s^＊，t^＊)=x(u^＊，v^＊)とともに出力する（ステップ２１１ｃ）。

（２変数関数に対する最適曲面の設計）
なお、本実施形態においては、離散データＤが与えられた場合について極値をもとめたが、関数f(s，t)が入力データとして与えられた場合についても、同様な方法で極値を求めることが可能である。関数f(s，t)が与えられたときのk次最適スプライン曲面x(s，t)の設計方法は以下の通りである。評価関数は、数式３３の代わりに次式を用いる。

この評価関数を最小にするτ(=vecτ)は、数式３４の代わりに線形方程式

を解けばよい。この右辺のベクトルは、数値積分を用いて求めることができる。なお、行列Q，Q₁ ⁽⁰⁰⁾，Q₂ ⁽⁰⁰⁾は数式３５、ベクトルb₁(s)、b₂(t)は数式３７で定義されている。

（数値実験例）
極値検出の数値例を示す。極値検出の対象となる最適平滑化スプライン曲面x(s，t)は、次式をサンプルして得られたデータに対して設計した。

ここでのパラメータの値を次の表に示す。

曲面の設計の際に用いた各パラメータは、k=3、α=β=1、 s₀=t₀=0、m₁=m₂=10 である。従って、３次スプライン曲面を領域S=[0，10]×[0，10]で設計した。設計に用いたデータは、上記d(s，t)をs，t方向とも0.5間隔の格子点でサンプリングして生成した。設計された曲面x(s，t)に上述の極値検出法を適用して得られた結果を図１２に示す。図では、検出された極大値（△印）と極小値（▽印）を、曲面x(s，t)の３次元表示に重ねて示した。また、st平面上には曲面の等高線表示と極値点とを示した。いずれの図からもすべての極値が検出されていることが確認できる。

以上説明したように、本実施形態の最適平滑化スプライン曲面による極値検出方法及び極値検出プログラムによれば、与えられたk次スプライン曲面x(s，t)のすべての極値（極小値、極大値）を求めることが可能となる。従って、最小値や最大値も求まったことになる。曲面x(s，t)の設計は、与えられた離散空間データD（数式３２）に基づく既述の方法でも良いし、あるいは与えられた関数f(s，t)を近似する曲面の設計によっても良い。後者の場合、設計されたスプライン曲面x(s，t)の最大値や最小値を求め、これらを与える点を初期値として、関数f(s，t)にさらに山登り法等の数値的探索法を適用することによって、f(s，t)の最大値や最小値がより正確、かつ系統だった方法で求められることになる。

なお、上記実施形態では、最適平滑化スプライン関数x(s，t)の極値を求めたが、同様な方法で、x(s，t) の任意の偏導関数の極値を求めることが可能である。

（偏導関数の極値）
いま記号の簡単化のためにx(s，t)とx(u，v)の偏導関数を以下のようにおく。

また、問題が意味をもつためにi，j≦k−2とする。このとき、数式４２の関係は、x_i，j(s，t)=αⁱβ^jx_i，j(u，v)となり、偏導関数x_i，j(s，t)の極値は、x_i，j(u，v)の極値から求めることができる。すなわち、x_i，j(u，v)が点(u^＊，v^＊)∈εで極値をとれば、x_i，j(s，t)は点(s^＊，t^＊)=((1/α)u^＊+s_κ，(1/β)v^＊+t_μ)で、極値x_i，j(s，t)=αⁱβ^jx_i，j(u^＊，v^＊)をもつことが分かる。

（関数x_i，j(u，v)の表現）
x_i，j(u，v)は，数式４３から

と表される。ここで、

である。このx_i，j(u，v)は、数式４３のx(u，v)と同形式であることから、その極値もx(u，v)の極値の検出、計算の場合と同様な手順に従って求めることができる。もちろん、偏導関数x_i，j(s，t)のすべての極値（極小値、極大値）を求めることができ、従って、その最小値や最大値も求まったことになる。

（第３実施形態）
以下、本発明の第３実施形態である最適平滑化スプラインによる極値検出方法及び極値検出プログラムについて、図を参照して詳細に説明をする。

本発明の第２実施形態では、空間における離散データを近似する最適平滑化スプライン曲面を生成し極値の検出を行ったが、本実施形態では、離散データとしてデジタル画像の各ピクセルの輝度値データについて極値の検出を行う。さらに、ここでは、最適平滑化スプライン曲面x(s，t)のs方向およびt方向の1次偏導関数の極値を検出することにより、デジタル画像のエッジ検出を行う。

（曲面の生成）
デジタル画像のエッジ検出に上記の極値検出法を適用した。対象としたのは２５６×２５６（pixel）のモノクロ画像であり、各ピクセルの輝度値をデータとして、まず、最適平滑化スプライン曲面x(s，t)を設計し、画像のエッジは、輝度値の変化率が急激に変化する部分（つまり、最適平滑化スプライン曲面x(s，t)のs方向およびt方向の1次偏導関数の極値が存在する部分）に存在するとの考えに基づき検出する。

本実施形態では、最適平滑化スプライン曲面x(s，t)をライン毎に（たとえばsをs=pに固定し）処理している。その場合、エッジは変数tに関するスプライン曲線x(p，t)の１次導関数の極値として求めればよい。

図６は、本実施形態の最適平滑化スプライン曲面を概略的に示す図である。図６Ａは、各ピクセルの位置（座標）及びその輝度の分布から求められた最適な制御点（重み係数）の分布を示す図であり、図６Ｂは、各ピクセルの空間における離散データに基づいて生成された最適平滑化スプライン曲面x(s，t)及びライン毎に生成された関数x(s，q)、x(p，t)を示す図である。

図１３は、本実施形態の最適平滑化スプラインによる極値検出方法及び極値検出プログラムを用いたデジタル画像のエッジ検出処理装置を示すブロック図である。

本実施形態の画像のエッジ検出処理装置２１は、ＣＰＵ２５（演算手段）とメモリ装置２４（ＲＡＭやＲＯＭ等：記憶手段）を有し、このメモリ装置２４に格納されたプログラムを実行し、図１４〜１６のフローチャートに記載した動作を実行することにより、図１３のブロック図に示すように、ＣＰＵ２５（演算手段）と、メモリ装置２４（記憶手段）と、マウス、キーボード、スキャナ、データロガー等の入力装置２６（入力手段）と、バス２３と、ディスプレイ装置２７（出力手段）と、プリンタ装置２８（出力手段）と、外部ネットワーク３０と接続するネットワーク装置２９（入力手段・出力手段）を備えた処理装置として機能する。

図１４は、本実施形態の最適平滑化スプラインによるデジタル画像のエッジ検出方法の処理の概要を示すフローチャートである。

ステップ３０１では、ピクセルデータ（i，j，f_i，j）、i= 1，2,・・・，I、j= 1，2,・・・，J、を、画像のエッジ検出処理装置２１に入力する工程である。ここで、f_i，jは、ピクセルの座標（i，j）におけるピクセルの輝度値を示す。ピクセルデータは、例えば、不図示の光学的読取手段にてデータとして取得されうるものである。また、ピクセルデータは、ネットワーク装置２９を介して不図示の外部装置から取得してもよい。そして、しきい値ρをメモリ装置２４に格納する。しきい値ρは、極値がエッジ部分かどうか判別するためのパラメータである。x⁽¹⁾(s，t)の極値（極大、極小）の絶対値がある程度以上の大きさの箇所をエッジとみなす。その程度は、画像の状態やエッジ検出条件により異なるため、エッジ検出においてしきい値ρを設定する。しきい値ρよりも極値の絶対値が大きければ、その箇所をエッジ部分としてみなし、それ以外の場合は、極値であってもエッジ部分とはみなさない。

ステップ３０２では、ＣＰＵ２５は、上記入力ピクセルデータ（i，j，f_i，j）に基づいて、最適平滑化スプライン曲面を生成する。具体的には、ＣＰＵ２５は、最適な重み係数τ_ij（i=−k, −k+1，・・・，I−1、j=−k, −k+1，・・・，J−1）を求め、数式２８に示す最適平滑化スプライン曲面を生成する。最適平滑化スプライン曲面を生成する方法については、本発明の第２実施形態で説明した通りである。その際のパラメータはα＝β＝１，s₀=t₀=１、m₁=I−１，m₂=Ｊ−１とする。

最適平滑化スプライン曲面を生成したのち、その曲面をライン毎に（たとえばsをs=qに固定し）処理する。図６のｓ軸方向に対する処理をするために、ステップ３０３へ進み、図６のｔ軸方向の処理をするためにステップ３０４へ進む。

ステップ３０３では、ＣＰＵ２５は、ステップ３０２で生成された最適平滑化スプライン曲面x(s，t)について、ｔ座標をq(q=1，2,・・・，J)で固定した複数の曲線x(s，q)（図６Ｂ参照）を生成し、曲線x(s，q)の１次導関数x⁽¹⁾(s，q)について極値を検出する（ステップ３０３ａ）。そして、ＣＰＵ２５は、この極値に基づいて、各々エッジの判定の処理を行う（ステップ３０３ｂ）。

ステップ３０４では、ＣＰＵ２５は、ステップ３０２で生成された最適平滑化スプライン曲面x(s，t)について、s座標をp(p=1，2,・・・，I)で固定した複数の曲線x(p，t)（図６Ｂ参照）を生成し、曲線x(p，t)の１次導関数x⁽¹⁾(p，t)について極値を検出する（ステップ３０４ａ）。そして、ＣＰＵ２５は、この極値に基づいて、各々エッジの判定の処理を行う（ステップ３０４ｂ）。

ステップ３０５では、ＣＰＵ２５は、ステップ３０３とステップ３０４で得られた判定により検出されたエッジ箇所を総和としてまとめる。すなわちｓ方向、ｔ方向のうち少なくとも一方でエッジと判定された箇所をエッジと判定する。

ステップ３０６では、ＣＰＵ２５は、ステップ３０５の結果に基づいて、入力画像に対するエッジ検出結果を画像に変換し、ディスプレイ装置２７やプリンタ装置２８に出力する。

図１５は、ステップ３０３における処理の詳細を示すフローチャートである。ステップ４０１〜４０３が、図１４のステップ３０３ａに相当し、ステップ４０４〜４０７が、図１４のステップ３０３ｂに相当する。

ステップ４０１は、ＣＰＵ２５に、図１４のステップ３０１で求めたピクセルデータ（i，j，f_i，j）について重み係数τ_i,j（i=−k, −k+1，・・・，I−1、j= −k, −k+1，・・・，J−1）を入力する工程である。

なお、以下に使用される数式に関して、δを変数とする関数xには、すべて‘ハット（＾）’がつくものであるが、明細書本文中では省略するものとする。

ステップ４０２では、ＣＰＵ２５は、ステップ３０２で生成された最適平滑化スプライン曲面x(s，t)を、ｔ座標をq(q= 1，2,・・・，J)で固定して生成された複数の曲線x(s，q)について、節点区間[s_j，s_j+1)において正規化を行う。そして、単位区間[0，1)における正規化された曲線x(δ_s，q)を生成し、１次導関数x⁽¹⁾(δ_s，q)を求め、多項式であるx⁽¹⁾(δ_s，q)について、各係数a₀，・・・，a_k-1の計算を行う。

ステップ４０３では、ＣＰＵ２５は、多項式x⁽¹⁾(δ_s，q)の極値点δ_s1、δ_s2、・・・、δ_sm (m≦k−1)の検出を行う。多項式x⁽¹⁾(δ_s，q)の極値点δ_s1、δ_s2、・・・、δ_smは、x⁽²⁾(δ_s，q)の根のうち、単位区間[0，1)にある実根を検出して求めることができる。

ステップ４０４、４０５では、ＣＰＵ２５は、すべての各極値点δ_s1、δ_s2、・・・、δ_smについて、x⁽¹⁾(δ_sl，q)・x⁽³⁾(δ_sl，q)＜0、かつ、｜x⁽¹⁾(δ_sl，q)｜＞ρ（0≦l≦m）か否かを判断する。x⁽¹⁾(δ_sl，q)・x⁽³⁾(δ_sl，q)＜0、かつ、｜x⁽¹⁾(δ_sl，q)｜＞ρ（0≦l≦m）の場合には、ステップ４０６へ進む。そうでない場合には、ステップ４０４に戻る。

ステップ４０６では、ＣＰＵ２５は、x⁽¹⁾(δ_sl，q)・x⁽³⁾(δ_sl，q)＜0、かつ、｜x⁽¹⁾(δ_sl，q)｜＞ρ（0≦l≦m）の場合には、点（p+(1/α)δ_sl ，q）においてエッジがあるものと判断する。これは、x⁽¹⁾(t,q)の極値の箇所が必ずしもエッジの位置に対応しないことがあり、x⁽¹⁾(δ_sl，q)・x⁽³⁾(δ_sl，q)＜0を満たす箇所が、実際のエッジに対応するからである。

ステップ４０７では、ＣＰＵ２５は、ｐ＜Ｉ−１の場合には、ステップ４０２〜ステップ４０６の各工程を繰り返し、ステップ４０８では、ｑ＜Ｊの場合には、ステップ４０２〜ステップ４０７の各工程を繰り返し、ステップ３０２で生成された最適平滑化スプライン曲面x(s，t)のｔ座標をq(q=1，2,・・・，J)で固定して生成された複数の曲線x(s，q)のすべてについてエッジを判断する。

図１６は、ステップ３０４における処理の詳細を示すフローチャートである。

ステップ５０１は、ＣＰＵ２５に、図１４のステップ３０１で求めたピクセルデータ（i，j，f_i，j）について重み係数τ_i,j（i=−k, −k+1，・・・，I−1、j=−k, −k+1，・・・，J−1）を入力する工程である。

ステップ５０２では、ＣＰＵ２５は、ステップ３０２で生成された最適平滑化スプライン曲面x(s，t)を、s座標をp(p=1，2,・・・，I)で固定して生成された複数の曲線x(p，t)について、区間[t_j，t_j+1)において正規化を行う。そして、単位区間[0，1)における正規化された曲線x(p，δ_t)を生成し、１次導関数x⁽¹⁾(p，δ_t)を求め、多項式であるx⁽¹⁾(p，δ_t)について、各係数b₀，・・・，b_k-1の計算を行う。

ステップ５０３では、ＣＰＵ２５は、多項式x⁽¹⁾(p，δ_t)の極値点δ_t1、δ_t2、・・・、δ_tn (n≦k−1)の検出を行う。多項式x⁽¹⁾(p，δ_t)の極値点δ_t1、δ_t2、・・・、δ_tnは、x⁽²⁾(p， δ_t)の根のうち、単位区間[0，1)にある実根を検出して求めることができる。

ステップ５０４、５０５では、ＣＰＵ２５は、すべての各極値点δ_t1、δ_t2、・・・、δ_tnについて、x⁽¹⁾(p，δ_tl)・x⁽³⁾(p，δ_tl)＜0、かつ、｜x⁽¹⁾(p，δ_tl)｜＞ρ（0≦l≦n）か否かを判断する。x⁽¹⁾(p，δ_tl)・x⁽³⁾(p，δ_tl)＜0、かつ、｜x⁽¹⁾(p，δ_tl)｜＞ρ（0≦l≦n）の場合には、ステップ５０６へ進む。そうでない場合には、ステップ５０４に戻る。

ステップ５０６では、ＣＰＵ２５は、x⁽¹⁾(p，δ_tl)・x⁽³⁾(p，δ_tl)＜0、かつ、｜x⁽¹⁾(p，δ_tl)｜＞ρ（0≦l≦n）の場合には、点（p，q+(1/α)δ_tl）においてエッジがあるものと判断する。

ステップ５０７では、ｑ＜Ｊ−１の場合には、ＣＰＵ２５は、ステップ５０２〜ステップ５０６の各工程を繰り返し、ステップ５０８では、ｐ＜Ｉの場合には、ステップ５０２〜ステップ５０７の各工程を繰り返し、ステップ３０２で生成された最適平滑化スプライン曲面x(s，t)を、s座標をp(p=1,2,・・・，I)で固定して生成された複数の曲線x(p，t)のすべてについてエッジを検出する。

（デジタル画像のエッジ検出例）
図１７及び図１８は、本発明の第３実施形態の最適平滑化スプラインによる極値検出方法に用いたデジタル画像のエッジ検出方法を用いて、実際にエッジ検出を行った例を示す図である。図１７は、原画像（ａ）にノイズが含まれない例を示し、図１８は、原画像（ａ）にノイズが含まれる例を示す。

図１７（ｂ）と図１８（ｂ）は、本実施形態のエッジ検出法を適用して生成した画像である。また、比較例として、cannyフィルタ及びsobelフィルタを用いてエッジ検出した画像を、図１７（ｃ）と図１８（ｃ）（cannyフィルタ）及び図１７（ｄ）と図１８（ｄ）（sobelフィルタ）に示した。

図からもわかるように、特に、原画像（ａ）にノイズが含まれる場合において、他のフィルタと比較して、精度良くエッジの検出ができることがわかる。よって、本発明の最適平滑化スプラインによる極値検出方法に用いたデジタル画像のエッジ検出方法の有効性が証明された。

以上説明したように、本実施形態の最適平滑化スプラインによる極値検出方法によれば、スプライン曲面の設計が画像情報全体の情報に基づいた曲面近似となるため、デジタル画像のエッジ検出は大域的な方法となる。従来の局所的な（画像の一部の情報に基づいた）方法とはまったく異なる方法であるといえ、本実施形態の最適平滑化スプラインによる極値検出方法は、結果的に画像のノイズの影響を受けにくい。

なお、上記各実施形態においては、スプライン曲線（１次元）の場合、スプライン曲面（２次元）の場合について説明をしたが、本発明はこれに限られず、さらなる高次元の場合に拡張して適用することが可能である。

最適平滑化スプライン曲線を概略的に示す図である。 本発明の第１実施形態の最適平滑化スプラインによる極値検出方法に用いられる処理装置を示すブロック図である。 本発明の第１実施形態の最適平滑化スプラインによる極値検出方法の処理の流れを示すフローチャートである。 本発明の第１実施形態の最適平滑化スプラインによる極値検出方法の処理の流れにおいて、最適平滑化スプラインが３次である場合を示すフローチャートである。 本発明の第１実施形態の最適平滑化スプライン曲線とその極値を示す図である。 空間における離散データから求められた最適な制御点（重み係数）の分布を示す図である。 最適平滑化スプライン曲面を概略的に示す図である。 本発明の第２実施形態の最適平滑化スプラインによる極値検出方法に用いたデジタル画像のエッジ検出処理装置を示すブロック図である。 本発明の第２実施形態の最適平滑化スプラインによる極値検出方法及び極値検出プログラムの処理の概要を示すフローチャートである。 図８のフローチャートのステップ２０３における処理の詳細を示すフローチャートである。 図８のフローチャートのステップ２０６における処理の詳細を示すフローチャートである。 図８のフローチャートのステップ２１１における処理の詳細を示すフローチャートである。 本発明の第２実施形態の最適平滑化スプラインによる極値検出方法及び極値検出プログラムを用いて、実際に極値検出を行った例を示す図である。 本発明の第３実施形態の最適平滑化スプラインによる極値検出方法及び極値検出プログラムを用いたデジタル画像のエッジ検出処理装置を示すブロック図である。 本発明の第３実施形態の最適平滑化スプラインによるデジタル画像のエッジ検出方法の処理の概要を示すフローチャートである。 図１４のフローチャートのステップ３０３における処理の詳細を示すフローチャートである。 図１４のフローチャートのステップ３０４における処理の詳細を示すフローチャートである。 本発明の第３実施形態の最適平滑化スプラインによる極値検出方法を用いたデジタル画像のエッジ検出方法を用いて、実際にエッジ検出を行った例を示す図である。 本発明の第３実施形態の最適平滑化スプラインによる極値検出方法を用いたデジタル画像のエッジ検出方法を用いて、実際にエッジ検出を行った例を示す図である。

符号の説明

１：最適平滑化スプラインによる極値検出処理装置
３：バス
４：メモリ装置
５：ＣＰＵ
６：入力装置
７：ディスプレイ装置
８：プリンタ装置
９：ネットワーク装置

Claims (20)

1. データが入力される入力手段と、プログラムを格納する記憶手段と、前記プログラムを実行し前記データを処理する演算手段と、前記演算手段の処理結果を出力する出力手段とを備える処理装置で実行される極値検出方法であって、
   前記入力手段から入力された平面における離散データ群又は関数から、前記演算手段が最適平滑化スプライン曲線を生成するステップと、
   前記演算手段が、前記最適平滑化スプライン曲線又はその導関数の各節点区間について各々極値を求めることにより、前記最適平滑化スプライン曲線又はその導関数の極値のすべてを求め、該極値を前記出力手段へ出力するステップと、
   を有することを特徴とする最適平滑化スプライン曲線による極値検出方法。
2. 前記最適平滑化スプライン曲線又はその導関数の極値のすべてを求め、該極値を前記出力手段へ出力するステップは、
   前記演算手段が、前記最適平滑化スプライン曲線又はその導関数の各節点区間について各々正規化して、単位区間における正規化された多項式関数又はその導関数を生成するステップと、
   前記演算手段が、前記単位区間における正規化された多項式関数又はその導関数について極値を求めることにより、前記最適平滑化スプライン曲線又はその導関数の極値のすべてを求めるステップと、
   を有することを特徴とする請求項１に記載の最適平滑化スプライン曲線による極値検出方法。
3. 前記演算手段が前記最適平滑化スプライン曲線を生成するステップは、
   前記演算手段が前記平面における離散データ群又は関数から節点ごとの最適重み係数を求めるステップと、
   前記演算手段が前記節点ごとの正規化された一様なＢスプライン関数を、前記節点ごとの最適重み係数により重みづけして加え合わせて構成することにより、最適平滑化スプライン曲線を生成するステップと、
   を有することを特徴とする請求項１または２に記載の最適平滑化スプライン曲線による極値検出方法。
4. 前記スプライン関数は、３次のスプライン関数であることを特徴とする請求項３に記載の最適平滑化スプライン曲線による極値検出方法。
5. データが入力される入力手段と、プログラムを格納する記憶手段と、前記プログラムを実行し前記データを処理する演算手段と、前記演算手段の処理結果を出力する出力手段とを備える処理装置で実行される極値検出方法であって、
   前記入力手段から入力された空間における離散データ群又は関数から前記演算手段が最適平滑化スプライン曲面を生成するステップと、
   前記演算手段が前記最適平滑化スプライン曲面又はその偏導関数の各節点領域について各々極値を求めることにより、前記最適平滑化スプライン曲面又はその偏導関数の極値のすべてを求め、該極値を前記出力手段へ出力するステップと、
   を有することを特徴とする最適平滑化スプライン曲面による極値検出方法。
6. 前記最適平滑化スプライン曲面又はその偏導関数の極値のすべてを求め、該極値を前記出力手段へ出力するステップは、
   前記演算手段が前記最適平滑化スプライン曲面又はその偏導関数の各節点領域について各々正規化して、単位領域における正規化された多項式関数又はその偏導関数を生成するステップと、
   前記演算手段が前記単位領域における正規化された多項式関数又はその偏導関数について極値を求めることにより、前記最適平滑化スプライン曲面又はその偏導関数の極値のすべてを求めるステップと、
   を有することを特徴とする請求項５に記載の最適平滑化スプライン曲面による極値検出方法。
7. 前記演算手段が前記最適平滑化スプライン曲面を生成するステップは、
   前記演算手段が前記空間における離散データ群又は関数から節点ごとの最適重み係数を求めるステップと、
   前記演算手段が前記節点ごとの正規化された一様なＢスプライン関数を、前記節点ごとの最適重み係数により重みづけして加え合わせて構成することにより、最適平滑化スプライン曲面を生成するステップと、
   を有することを特徴とする請求項５または６に記載の最適平滑化スプライン曲面による極値検出方法。
8. ピクセルデータが入力される入力手段と、プログラムを格納する記憶手段と、前記プログラムを実行し前記ピクセルデータを処理する演算手段と、前記演算手段の処理結果を出力する出力手段とを備える画像処理装置で実行されるデジタル画像のエッジ検出方法であって、
   前記入力手段から入力されたピクセルの座標及び輝度値からなるピクセルデータから、前記演算手段が最適平滑化スプライン曲面を生成するステップと、
   前記演算手段が前記最適平滑化スプライン曲面を第１の座標軸上の複数の座標で固定した複数の曲線からなる第１の曲線群を生成し、該第１の曲線群を構成する複数の曲線の１次導関数の各節点区間について各々極値を求めることにより、前記第１の曲線群の１次導関数の極値のすべてを求めるステップと、
   前記演算手段が前記最適平滑化スプライン曲面を第２の座標軸上の複数の座標で固定した複数の曲線からなる第２の曲線群を生成し、該第２の曲線群を構成する複数の曲線の１次導関数の各節点区間について各々極値を求めることにより、前記第２の曲線群の１次導関数の極値のすべてを求めるステップと、
   前記演算手段が前記第１の曲線群の１次導関数の極値又は前記第２の曲線群の１次導関数の極値の絶対値が所定のしきい値より大きい箇所を、デジタル画像のエッジ部分と判定し、該判定結果を前記出力手段へ出力するステップと、
   を有することを特徴とするデジタル画像のエッジ検出方法。
9. 前記第１の曲線群の１次導関数の極値のすべてを求めるステップは、
   前記演算手段が前記第１の曲線群を構成する複数の曲線の１次導関数の各節点区間について各々正規化して、単位区間における正規化された第１の曲線群を構成する複数の曲線の１次導関数の極値を求めることにより、前記第１の曲線群の１次導関数の極値のすべてを求めるステップと、を有し、
   前記第２の曲線群の１次導関数の極値のすべてを求めるステップは、
   前記演算手段が前記第２の曲線群を構成する複数の曲線の１次導関数の各節点区間について各々正規化して、単位区間における正規化された第２の曲線群を構成する複数の曲線の１次導関数の極値を求めることにより、前記第２の曲線群の１次導関数の極値のすべてを求めるステップと、
   を有することを特徴とする請求項８に記載のデジタル画像のエッジ検出方法。
10. 前記演算手段が前記最適平滑化スプライン曲面を生成するステップは、
    前記演算手段が前記ピクセルの座標及び輝度値からなるピクセルデータから節点ごとの最適重み係数を求めるステップと、
    前記演算手段が前記節点ごとの正規化された一様なＢスプライン関数を、前記節点ごとの最適重み係数により重みづけして加え合わせて構成することにより、最適平滑化スプライン曲面を生成するステップと、
    を有することを特徴とする請求項８または９に記載のデジタル画像のエッジ検出方法。
11. データが入力される入力手段と、プログラムを格納する記憶手段と、前記プログラムを実行し前記データを処理する演算手段と、前記演算手段の処理結果を出力する出力手段とを備える処理装置に、
    前記入力手段から入力された平面における離散データ群又は関数から、前記演算手段が最適平滑化スプライン曲線を生成するステップと、
    前記演算手段が前記最適平滑化スプライン曲線又はその導関数の各節点区間について各々極値を求めることにより、前記最適平滑化スプライン曲線又はその導関数の極値のすべてを求め、該極値を前記出力手段に出力するステップと、
    を実行させるための最適平滑化スプライン曲線による極値検出プログラム。
12. 前記最適平滑化スプライン曲線又はその導関数の極値のすべてを求め、該極値を前記出力手段に出力するステップは、
    前記演算手段が前記最適平滑化スプライン曲線又はその導関数の各節点区間について各々正規化して、単位区間における正規化された多項式関数又はその導関数を生成するステップと、
    前記演算手段が前記単位区間における正規化された多項式関数又はその導関数について極値を求めることにより、前記最適平滑化スプライン曲線又はその導関数の極値のすべてを求めるステップと、
    を有することを特徴とする請求項１１に記載の最適平滑化スプライン曲線による極値検出プログラム。
13. 前記演算手段が前記最適平滑化スプライン曲線を生成するステップは、
    前記演算手段が前記平面における離散データ群又は関数から節点ごとの最適重み係数を求めるステップと、
    前記演算手段が前記節点ごとの正規化された一様なＢスプライン関数を、前記節点ごとの最適重み係数により重みづけして加え合わせて構成することにより、最適平滑化スプライン曲線を生成するステップと、
    を有することを特徴とする請求項１１または１２に記載の最適平滑化スプライン曲線による極値検出プログラム。
14. 前記スプライン関数は、３次のスプラインであることを特徴とする請求項１３に記載の最適平滑化スプライン曲線による極値検出プログラム。
15. データが入力される入力手段と、プログラムを格納する記憶手段と、前記プログラムを実行し前記データを処理する演算手段と、前記演算手段の処理結果を出力する出力手段とを備える処理装置に、
    前記入力手段から入力された空間における離散データ群又は関数から、前記演算手段が最適平滑化スプライン曲面を生成するステップと、
    前記演算手段が前記最適平滑化スプライン曲面又はその偏導関数の各節点領域について各々極値を求めることにより、前記最適平滑化スプライン曲面又はその偏導関数の極値のすべてを求め、該極値を前記出力手段に出力するステップと、
    を実行させるための最適平滑化スプライン曲面による極値検出プログラム。
16. 前記最適平滑化スプライン曲面又はその偏導関数の極値のすべてを求め、該極値を前記出力手段に出力するステップは、
    前記演算装置が前記最適平滑化スプライン曲面又はその偏導関数の各節点領域について各々正規化して、単位領域における正規化された多項式関数又はその導関数を生成するステップと、
    前記演算装置が前記単位領域における正規化された多項式関数又はその偏導関数について極値を求めることにより、前記最適平滑化スプライン曲面又はその偏導関数の極値のすべてを求めるステップと、
    を有することを特徴とする請求項１５に記載の最適平滑化スプライン曲面による極値検出プログラム。
17. 前記演算手段が前記最適平滑化スプライン曲面を生成するステップは、
    前記演算手段が空間における離散データ群又は関数から節点ごとの最適重み係数を求めるステップと、
    前記演算手段が前記節点ごとの正規化された一様なＢスプライン関数を、前記節点ごとの最適重み係数により重みづけして加え合わせて構成することにより、最適平滑化スプライン曲面を生成するステップと、
    を有することを特徴とする請求項１５または１６に記載の最適平滑化スプライン曲面による極値検出プログラム。
18. ピクセルデータが入力される入力手段と、プログラムを格納する記憶手段と、前記プログラムを実行し前記ピクセルデータを処理する演算手段と、前記演算手段の処理結果を出力する出力手段とを備える画像処理装置に、
    前記入力手段から入力されたピクセルの座標及び輝度値からなるピクセルデータから前記演算手段が最適平滑化スプライン曲面を生成するステップと、
    前記演算手段が前記最適平滑化スプライン曲面を第１の座標軸上の複数の座標で固定した複数の曲線からなる第１の曲線群を生成し、該第１の曲線群を構成する複数の曲線の１次導関数の各節点区間について各々極値を求めることにより、前記第１の曲線群の１次導関数の極値のすべてを求めるステップと、
    前記演算手段が前記最適平滑化スプライン曲面を第２の座標軸上の複数の座標で固定した複数の曲線からなる第２の曲線群を生成し、該第２の曲線群を構成する複数の曲線の１次導関数の各節点区間について各々極値を求めることにより、前記第２の曲線群の１次導関数の極値のすべてを求めるステップと、
    前記演算手段が前記第１の曲線群の１次導関数の極値又は前記第２の曲線群の１次導関数の極値の絶対値が所定のしきい値より大きい箇所を、デジタル画像のエッジ部分と判定し、該判定結果を前記出力手段へ出力するステップと、
    を実行させるためのデジタル画像のエッジ検出プログラム。
19. 前記第１の曲線群の１次導関数の極値のすべてを求めるステップは、
    前記演算手段が前記第１の曲線群を構成する複数の曲線の１次導関数の各節点区間について各々正規化して、単位区間における正規化された前記第１の曲線群を構成する複数の曲線の１次導関数の極値を求めることにより、前記第１の曲線群の１次導関数の極値のすべてを求めるステップと、を有し、
    前記第２の曲線群の１次導関数の極値のすべてを求めるステップは、
    前記演算手段が前記第２の曲線群を構成する複数の曲線の１次導関数の各節点区間について各々正規化して、単位区間における正規化された前記第２の曲線群を構成する複数の曲線の１次導関数の極値を求めることにより、前記第２の曲線群の１次導関数の極値のすべてを求めるステップと、
    を有することを特徴とする請求項１８に記載のデジタル画像のエッジ検出プログラム。
20. 前記演算手段が前記最適平滑化スプライン曲面を生成するステップは、
    前記演算手段が前記ピクセルの座標及び輝度値からなるピクセルデータから節点ごとの最適重み係数を求めるステップと、
    前記演算手段が前記節点ごとの正規化された一様なＢスプライン関数を、前記節点ごとの最適重み係数により重みづけして加え合わせて構成することにより、最適平滑化スプライン曲面を生成するステップと、
    を有することを特徴とする請求項１８または１９に記載のデジタル画像のエッジ検出プログラム。