JP4218824B2 - 木材の年輪箇所検出方法および年輪幅計測方法 - Google Patents
木材の年輪箇所検出方法および年輪幅計測方法 Download PDFInfo
- Publication number
- JP4218824B2 JP4218824B2 JP2002330131A JP2002330131A JP4218824B2 JP 4218824 B2 JP4218824 B2 JP 4218824B2 JP 2002330131 A JP2002330131 A JP 2002330131A JP 2002330131 A JP2002330131 A JP 2002330131A JP 4218824 B2 JP4218824 B2 JP 4218824B2
- Authority
- JP
- Japan
- Prior art keywords
- annual ring
- measurement line
- point
- signal
- wood
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
- G01N23/02—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
- G01N23/04—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
- G06F2218/10—Feature extraction by analysing the shape of a waveform, e.g. extracting parameters relating to peaks
Landscapes
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Biochemistry (AREA)
- General Health & Medical Sciences (AREA)
- General Physics & Mathematics (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Complex Calculations (AREA)
- Image Analysis (AREA)
- Length Measuring Devices By Optical Means (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
- Image Processing (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、波形上の特徴点の値または特徴点間距離が不規則な波形信号の特徴点を精度良く検出する方法に関し、更に詳しくは、例えば木材年輪数、年輪幅等の計測に好適に用いることができる木材の年輪箇所検出方法および年輪幅計測方法に関する。
【０００２】
【従来の技術】
年輪年代学では、年代確定の基準となる標準的な年輪幅変動のデータベースと照合することにより、各年輪層が何年に形成されたものかを一年単位で確定することができる。このデータベースは、奈良文化財研究所による精力的な調査研究によって、現在のところ日本では檜材で紀元前912年、杉材で紀元前1313年まで遡ることができる。ちなみに、年輪年代学の先進地域であるドイツでは、約一万年にわたる標準データベースが作成されている。年輪年代学では、主に以下のような事例が取り扱われている。
【０００３】
（ａ）木材の伐採年の推定
（ｂ）木造文化財（建築物、仏像、美術工芸品等）の制作年代や修理経緯の推定、真贋鑑定等
（ｃ）過去の長期にわたる気候変動の研究、地球環境温暖化の研究等
年輪年代学に使用される各年輪幅の時系列データは、誤検出（年輪以外のものを年輪として認識する誤り）も、非検出（年輪を認識しない誤り）も共に０という究極の検出性能が必要とされる。
【０００４】
このため、年輪幅の計測は、計測顕微鏡を用いた専用システムによって目視で行われており、その作業には熟練した技術と膨大な時間（300層クラスの試料で約１時間）を要していた。また、システム構成が大掛かりなものになるという問題もかかえていた。
【０００５】
計測作業を自動化するためには、デジタルカメラやスキャナ等の画像取得装置を用いて取得した年輪画像から、パソコンレベルのコンピュータを用いて各年輪幅の計測を行う方法が考えられており、これまでにもいくつかの試みがなされているが、前述の検出性能やシステム規模、価格などの点で問題があり、研究手段として一般的に広く用いられるには至っていないのが現状である。
【０００６】
とりわけ檜材については、文化財等に広く利用されていて年輪年代学上重要な樹種であるにもかかわらず、年輪幅が狭小であること、杉材に比べて年輪が不明瞭なことなどの条件から、自動計測の実用化が極めて困難であった。
【０００７】
本特許案件に類似した公知の技術の主要なもの以下に示す。
【０００８】
（１）「Win DENDRO」 カナダ Regent Instruments社 1988年
（http://www.regent.qc.ca.products/dendro/DENDRO.html参照）
ケベック大学のDr. Rejean Gagnon と Dr. Hubert Morin が設計し、Regent Instruments社が商品化した年輪年代学用ソフトである。このソフトは、年輪画像の輝度変化情報をもとに、年輪計測や木材組織の解析を行うことができる。このソフトのアルゴリズムの詳細は不明であるが、カタログ文面から察するところ、ウェーブレット処理および複数計測線情報の統合技術は用いていないと思われる。
【０００９】
（２）「画像処理システムを用いた年輪幅計測」 野田真人 1990年 樹木年輪研究会にて発表
この発表では、杉年輪の計測は可能であったが、檜年輪は計測不可能と結論づけている。この発表に係る計測方法においては、ウェーブレット処理と複数計測線情報の統合技術は用いられていない。
【００１０】
（３）特開平11−232427号公報
年輪数測定のために画像輝度情報を用いる点が記載されているが、同技術は上記（２）によりすでに公知である。この公報記載の公知技術においても、ウェーブレット処理と複数計測線情報の統合技術は何ら用いられていない。
【００１１】
【発明が解決しようとする課題】
年輪画像から画像上の計測線に沿って画素情報を取得することにより、輝度変化情報の波形信号または濃度変化情報の波形信号を得ることができる。濃度波形の極大箇所（輝度波形の場合は極小箇所）は各年輪層の各年輪層のうちの最も濃い部分（晩材最高濃度部）に対応する。したがって、濃度信号波形の極大箇所または輝度信号波形の極小箇所を認識していくことにより、各年輪層を識別することができる。
【００１２】
一方、濃度波形信号をさらに微分処理すると、暗部から明部へと移行する点（微分波形の極小箇所）が各年輪層の終端点（晩材終端）に対応することなる。したがって、計測線上でこの微分波形信号の極小箇所間距離を計測することにより、年輪幅をより正確に計測することができる。
【００１３】
年輪画像から画像上の計測線に沿って画素情報を取得することにより、輝度変化情報の波形信号または濃度変化情報の波形信号等を得ることは、上記した公知文献等により公知である。
【００１４】
しかしながら、木材の年輪幅は一般に不規則であり、最大年輪幅が最小年輪幅の数倍から約１００倍程度まで異なることも少なくない、したがって、年輪画像から得た波形信号に対する特徴点（極大箇所であるピーク点または極小箇所であるトラフ点）の検出精度を小さい特徴点間距離のレベルに合わせると、特徴点間距離の大きい領域では年輪に無関係なノイズを拾いやすくなる。一方、波形信号に対する特徴点（ピーク点またはトラフ点）の検出精度を大きい特徴点間距離のレベルに合わせると、特徴点間距離の小さい領域では特徴点の検出漏れが生じやすくなる。したがって、最小年輪幅と最大年輪幅との差が大きい計測対象物に対し年輪数や年輪幅を正確に計測することが困難となっている。
【００１５】
さらに、年輪画像においては、濃度のレベルが一定ではないため、年輪画像から得られる上記輝度波形信号、濃度波形信号、微分波形信号等はいずれも計測区間全体にわたりうねったような変化特性を有する場合が多い。このため、一定の閾値を用いて波形信号における特徴点（極大箇所であるピーク点または極小箇所であるトラフ点）を検出しようとすると、検出漏れが生じる原因となり、年輪数や年輪幅を正確に計測できない原因となる。
【００１６】
したがって、年輪年代学にとって最も基礎的なデータとなる各年輪幅の時系列データを迅速かつ高精度に得ることができる年輪箇所検出方法および年輪幅計測方法が要望されている。
【００１７】
本発明は、上記の要望に応え得る木材の年輪箇所検出方法および年輪幅計測方法を提供することを目的とするものであり、年輪画像から得た波形信号の特徴点をなす年輪箇所の値または該年輪箇所間の距離が不規則であっても、年輪画像から得た波形信号から波形の特徴点をなす年輪箇所を迅速かつ高精度に検出することができる木材の年輪箇所検出方法および該年輪箇所検出方法を用いて行う年輪幅計測方法を提供しようとするものである。
【００１８】
【課題を解決するための手段】
上記課題を解決するために、本発明は、第１発明として、撮像装置によって取得した木材の表面又は内部の年輪画像をコンピュータ装置により処理するコンピュータ処理工程を有する木材の年輪箇所検出方法であって、前記コンピュータ処理工程が、前記年輪画像上に開始点と終端点とで区間を指定して設定した計測線上に沿って前記年輪画像から画素情報を取得する画素情報取得工程と、前記画素情報を輝度情報信号または濃度情報信号である波形信号に変換する波形信号取得工程と、所定のマザーウェーブレットと、ゼロまたは負の整数値からなる複数段階のスケールレベルとを用いて前記波形信号をウェーブレット変換することにより各スケールレベル毎のウェーブレット変換信号を生成するウェーブレット変換工程と、前記計測線上の各点において前記ウェーブレット変換信号から前記各スケールレベルに応じて予め設定した区間幅における区間二乗平均値を算出する区間二乗平均値算出工程と、前記計測線上の各点において前記区間二乗平均値が最大となるウェーブレット変換信号のスケールレベルをその点における支配レベルとして決定する支配レベル決定工程と、前記計測線上の各点において前記各支配レベルに応じて予め設定した区間幅における前記波形信号が最大値となる点を前記計測線上の年輪箇所として決定する年輪箇 所決定工程と、を有することを特徴とする木材の年輪箇所検出方法を提供する。
【００１９】
さらに、第２発明として、撮像装置によって取得した木材の表面又は内部の年輪画像をコンピュータ装置により処理するコンピュータ処理工程を有する木材の年輪箇所検出方法であって、前記コンピュータ処理工程が、前記年輪画像上に開始点と終端点とで区間を指定して設定した計測線上に沿って前記年輪画像から画素情報を取得する波形信号取得工程と、前記画素情報を輝度情報信号または濃度情報信号に変換した後更に微分処理して波形信号を得る波形信号取得工程と、所定のマザーウェーブレットと、ゼロまたは負の整数値からなる複数段階のスケールレベルとを用いて前記波形信号をウェーブレット変換することにより各スケールレベル毎のウェーブレット変換信号を生成するウェーブレット変換工程と、前記計測線上の各点において前記ウェーブレット変換信号から前記各スケールレベルに応じて予め設定した区間幅における区間二乗平均値を算出する区間二乗平均値算出工程と、前記計測線上の各点において前記区間二乗平均値が最小となるウェーブレット変換信号のスケールレベルをその点における支配レベルとして決定する支配レベル決定工程と、前記計測線上の各点において前記各支配レベルに応じて予め設定した区間幅における前記波形信号が最小値となる点を前記計測線上の年輪箇所として決定する年輪箇所決定工程と、を有することを特徴とする木材の年輪箇所検出方法を提供する。
【００２０】
さらに第３発明として、上記第１発明又は第２発明の構成を有する木材の年輪箇所検出方法において、前記ウェーブレット変換に次式（６）、すなわち、
【数３】
ここで、ｆ(x)は波形信号、ψ(x)はマザーウェーブレット、ｂjはスケールパラメータであり、ｂは定数（ただし、ｂ＞１）、ｊは0または負の整数からなるスケールレベル、ｋはトランスレートパラメータである
を用い、さらに、前記区間二乗平均の演算に、次式（９）、すなわち、
【数４】
ここで、ｊは式（６）で用いたスケールレベル、ｋはトランスレートパラメータ、ｐjはスケールレベルｊが低レベルになるほど大きくなるようにスケールレベルｊに応じて設定される定数である
を用いることを特徴とする木材の年輪箇所検出方法を提供する。
【００２１】
上記構成において、好ましくは、前記区間二乗平均の演算式（９）におけるｐjが次式（１０）、すなわち、
【００２２】
【数５】
ここで、ａは式（６）で用いたマザーウエーブレットψ(x)のサポートによって決定される定数、ｂは式（６）で用いた定数、ｊは式（６）で用いたスケールレベルである
で定義されるものである。
【００２３】
さらに、第４発明として、上記第２発明の構成を有する木材の年輪箇所検出方法において、前記微分処理は、前記計測線上における複数画素分の間隔を隔てた複数の画素同士における輝度または濃度情報信号の差分演算処理であることを特徴とする請求項２記載の木材の年輪箇所検出方法を提供する。
【００２４】
さらに、第５発明として、上記第１発明又は第２発明の構成を有する木材の年輪箇所検出方法において、前記計測線が主計測線と、該主計測線の両側に各々等間隔を隔てて並列に設定される複数の副計測線とからなり、前記年輪箇所決定工程は、該主計測線および副計測線の始端から各々略同一距離の点において年輪候補点が検出されたとき、少なくとも該年輪候補点の数が主計測線および副計測線の本数に対し過半数であることを条件の１つとして、その年輪候補点を前記年輪箇所として決定することを特徴とする木材の年輪箇所検出方法を提供する。
【００２５】
さらに、第６発明として、上記第５発明の構成を有する木材の年輪箇所検出方法において、前記年輪箇決定工程は、前記主計測線の両側に各々２本の副計測線を設定し、該主計測線および副計測線の始端から各々略同一距離の点において年輪候補点が検出されたとき、少なくとも主計測線上と該主計測線に近い位置にある２本の副計測線のうちのいずれか一方上に年輪候補点があり且つ他のいずれかの副計測線上に年輪候補点がある場合、または、主計測線に近い位置にある２本の副計測線上と他のいずれかの副計測線上に年輪候補点がある場合に、その年輪候補点を前記年輪箇所として決定するものであることを特徴とする木材の年輪箇所検出方法を提供する。
【００２６】
さらに、第７発明として、上記第１発明又は第２発明の構成を有する木材の年輪箇所検出方法において、前記波形信号取得工程は、前記計測線に沿って取得した前記波形信号を前記計測線の周辺の画素の情報を用いて平滑化する工程を含むことを特徴とする年輪箇所検出方法を提供する。
【００２７】
さらに、第８発明として、上記第１発明又は第２発明の構成を有する木材の年輪箇所検出方法を用いて行う木材の年輪幅計測方法であって、前記コンピュータ処理工程が、前記年輪箇所決定工程により決定された前記計測線上の前記各年輪箇所間の画素数と１画素の大きさとから前記各年輪箇所間の年輪幅を算出する年輪幅算出工程を有することを特徴とする木材の年輪幅計測方法を提供する。
【００２８】
【発明の実施の形態】
以下、図面を参照して本発明の好ましい実施形態を説明する。
【００２９】
図１および図２は、本発明に係る木材の年輪箇所検出方法および年輪幅計測方法の一実施形態を示すフローチャートである。図１および図２を参照すると、木材の年輪箇所検出および年輪幅計測の手順は次のプロセスからなる。
【００３０】
１．年輪画像の取得（ステップＳ１）
２．計測箇所の指定（ステップＳ２）
３．画素情報の取得（ステップＳ３）
４．濃度変換（ステップＳ４、Ｓ５）
５．平滑化処理（ステップＳ６、Ｓ７）
６．ピーク信号またはエッジ信号の取得（ステップＳ８〜Ｓ１２）
７．マザーウェーブレットの指定（ステップＳ１３）
８．ウェーブレット変換（ステップＳ１４）
９．区間二乗平均の演算（ステップＳ１５、Ｓ１６）
10．支配レベルの決定（ステップＳ１７）
11．区間内最大値（エッジ信号の場合は区間内最小値）の決定（ステップＳ１８〜Ｓ２０）
12．主計測線と副計測線の相互参照（ステップＳ２１）
13．副計測線間の相互参照（ステップＳ２２）
14．年輪箇所の決定（ステップＳ２３）
15．年輪幅の計測および結果出力（ステップＳ２４、Ｓ２５）
また、本発明の上記実施形態は、上記ステップＳ2からステップＳ２５までのプロセスを実行することができるプログラムを、基本システムがインストールされたパソコンなどのコンピュータ装置にインストールすることにより、実行可能となる。
【００３１】
この実施形態におけるプログラムは、計測箇所の指定手段、画素情報取得手段、濃度変換手段、平滑化処理手段、ピーク信号またはエッジ信号取得手段、マザーウェーブレット指定手段、ウェーブレット変換手段、区間二乗平均の演算手段、支配レベル決定手段、区間内最大値（エッジ信号の場合は区間内最小値）の決定手段、主計測線と副計測線の相互参照手段、副計測線間の相互参照手段、年輪箇所の決定手段、年輪幅の計測手段および結果出力手段を有する。
【００３２】
次に、上記各プロセスについて詳述する。
【００３３】
１．年輪画像の取得（ステップＳ１）
デジタルカメラやスキャナ等の撮像装置を用いて年輪画像を取得する。例えば小型の移動可能な試料の場合にはスキャナを使用し、大型の試料には高解像度タイプのデジタルカメラを使用する。この際、最も注意すべき点は、所定の画像解像度を得ることである。檜材では、最も年輪間隔の狭小な部分で1層あたり約0.1mmという箇所を生ずることがある。このような部分を識別するためには、ナイキストの標本化定理により、１画素あたりの大きさをその半分の0.05mm以下にすべきである。
【００３４】
さらに安全と計測誤差を考慮すると、１画素あたりの大きさ約0.02mmに相当する1200dpiの解像度で画像取得を行うことが好ましい。大面積の試料でこの解像度を計測面全域にわたって保持するためには、一回あたりの撮影を小領域に限定し、画像を継ぎ合わせるという手法を用いることが好ましい。
【００３５】
本発明で扱う年輪画像としては、上記のようなデジタル画像取得装置を用いて試料表面から直接取得した画像のほかに、銀塩写真（通常のフィルムを使用する写真）等をデジタル変換した画像も同様に扱うことができる。また、X線写真、X線CT画像（X線断層写真）、MRI画像のような試料内部の年輪画像も扱うことができる。
【００３６】
２．計測箇所の指定（ステップＳ２）
年輪幅計測は、図３に示すように、少なくとも１本の計測線１に沿って行われる。木口（横断面）円盤上では、通常、計測線は樹木中心（図３中左方向）から樹皮方向（図３中右方向）に向かって設定される。また、柾目（縦断面）では、木理と直交するように樹木中心方向から樹皮方向に向かって設定される。本発明においては、好ましくは、認識性能の向上を図るため、１本の主計測線１に対し、それと平行に適切な間隔を置いて主計測線１の両側に各々複数の副計測線２、３、４、５を設定する（図３参照）。
【００３７】
副計測線の数は、木口で2本、柾目面で4本が適当であることが我々の実験によって確かめられている。ただし、試料による個体差などもあり、必ずしもこの限りではない。実際のプログラム上では、副計測線を設定する間隔と本数、および年輪画像上で主計測線１の開始点１ａと終端点１ｂの2点を指定することにより、計測箇所を指定可能とすることが好ましい。
【００３８】
３．画素情報の取得（ステップＳ３）
上記したように、主計測線１および副計測線２〜５に沿って各画像の画素情報（通常のパソコンの場合には、BGR3チャネルの8ビットデジタル信号）を取得する。画像取得に関しては、計測線上の画素情報のみでなく、後述の平滑化処理のために計測線周辺画素（具体的な周辺画素数は、平滑化の程度によって決定される）の情報も同時に取得することが好ましい。
【００３９】
なお、本発明方法による年輪認識では色情報を使用しないため、BGR信号の平均信号や混合信号、あるいは特定チャネルの信号のみ（たとえばG信号のみ）の輝度に基づく信号を使用することができる。通常は、B：G：R=1：1：1の割合で平均混合した信号を使用すればよい。
【００４０】
４．濃度変換（ステップＳ５）
上記方法で取得された画素情報は、各画素の輝度に起因しているため、値が大きくなるほど明るくなるという特性を有する。さらに、パソコンシステム上のRGB信号には、通常、CRTモニタの特性を補正するためにガンマ補正という非線形変換が施されている。一方、年輪の濃淡変化は年輪の細胞密度に起因しているため、密度あるいは密度に関連している画像濃度をもとに濃淡を記述するのが、年輪年代学の慣例である。
【００４１】
濃度記述では、値が大きくなるほど暗くなる。本発明による実施形態では、こうした年輪年代学の慣例に従い、RGB平均信号を濃度信号に変換している（図４参照）。変換方法としては、画像中の濃度既知部分（通常はステップタブレットチャート）から導出される変換式を用いることができる。変換式の一例をあげると、
【００４２】
【数６】
また、非線形の変換ではなく、下記のような線形変換で代用してもよい。
【００４３】
【数７】
なお、２値化等の画像処理を施す場合や、従来の方法では年輪識別のための閾値設定の必要性から、年輪信号波形の振幅特性をどのように記述するかが重要であった。
【００４４】
しかし、本発明では後述するとおり、年輪信号波形の値そのものとは無関係に年輪識別を行うことができるため、RGB、輝度あるいは濃度いずれの信号を用いても、年輪識別性能に大差はない。したがって、処理ステップや時間短縮を図るためには、このプロセスは省略可能である。ただし、このプロセスを省略すると年輪信号波形は反転することに留意する必要がある。
【００４５】
５．平滑化処理（ステップＳ６、Ｓ７）
上記プロセスによって取得した主計測線１および副計測線２〜５上の各画素値に対し、周辺画素の情報を用いて平滑化を行う。平滑化を施すことによって、試料自体や測定機器に起因するノイズによる誤認識を軽減することができる（図５参照）。
【００４６】
平滑化方法には、移動平均処理または、移動メジアン処理を用いることができる。周辺画素数としては、例えば1200dpi解像度の場合、我々の実験では木口面で、計測線直交方向に5画素、計測線方向に3画素程度、柾目面で計測線直交方向に15画素、計測線方向に5画素程度が適切であったが、これらの画素数に限定されない。なお、平滑化処理を行わない場合には、平滑化処理範囲を計測線直交方向および計測線方向に各々１画素として設定する。
【００４７】
６．ピーク信号またはエッジ信号の取得（ステップＳ８〜Ｓ１２）
上記プロセスによって取得した信号は、波形の極大箇所が各年輪層の最も濃い部分（晩材最高濃度部）を表している。この濃度情報信号を、以下、ピーク信号と称することにする。このピーク信号の特徴点であるピーク点（波形の極大箇所）を認識していくことにより、各年輪層を識別することができる。
【００４８】
一方、このピーク信号の微分信号（実際の演算は差分信号）における特徴点であるトラフ点（波形の極小箇所）は、年輪画像の暗部（高濃度部）から明部（低濃度部）へと移行する点に対応し、各年輪層の終端点（晩材終端）に対応する。この微分信号を、以下、エッジ信号と称することにする。
【００４９】
通常の微分信号は、隣接画素どうしの差分演算によって得られるが、隣接差分演算ではノイズの影響を受けやすくなる。そこで、本発明の好ましい実施形態では、複数画素の間隔を隔てた複数の画素同士の差分演算を行うことにより、ノイズの軽減を図っている（図６中央の波形説明参照）。
【００５０】
ピーク信号（図６上部の波形）とエッジ信号（図６下部の波形）を比較すると、波形ノイズの影響などにより、各年輪の識別にはピーク信号の方が検出性能は良好であるが、年輪幅を正確に計測するという観点では晩材終端に対応しているエッジ信号の方が精度は高い。そこで、本発明の好ましい実施形態では、このような特徴を活かして、ピーク信号による認識、エッジ信号による認識、両者を併用した認識をメニュー上で選択可能にしている。
【００５１】
後述の「８．ウェーブレット変換」の箇所で述べるウェーブレット変換（畳み込み積分演算）では、ある１点の画素に対して周辺画素を含む演算を行う。したがって、計測線両端部（図３の開始点１ａおよび終端点１ｂ）では、画素不足による演算不能を防止するため、信号自身を折り返してダミーの信号（図７の波形の両端の点線部分参照）を作成しておくことが望ましい。ダミー信号として必要な長さは、後述の「７．マザーウェーブレットの指定」で指定するマザーウェーブレットのサポートと「８．ウェーブレット変換」で指定するレベル数に依存するが、最も長い端部所理の長さは最低位のレベルｊの時であるから、これに合わせてダミー信号を延長しておけば十分である。
【００５２】
７．マザーウェーブレットの指定（ステップＳ１３）
マザーウェーブレットとして用いることができる関数ψ(x)は、
（ａ）端部が０または０に収束する関数であること
（ｂ）全区間での総和（積分値）が０であること
の２つの条件を満たす関数である。
【００５３】
また次の条件は、必ずしも満たす必要はないが、可能であれば満たしていることが望ましい。
【００５４】
（ｃ）サポート（関数の値が０ではない区間）がコンパクトに纏まっていること
マザーウェーブレットの具体例として、次のフレンチハット型の関数（式(１３)）とメキシカンハット型の関数（式(１４)）をあげることができる（図８（ａ）、（ｂ）参照）。
【００５５】
【数８】
【００５６】
【数９】
本発明の好ましい実施形態においては、生成が容易でサポートもコンパクトにまとまっているため、式(１３)のフレンチハット型関数をデフォルトとして用いている。もちろん、指定すればメキシカンハット型関数やその他のマザーウェーブレットも使用可能である。
【００５７】
８．ウェーブレット変換（ステップＳ１４）
１次元の画像信号ｆ(x)に対するウェーブレット変換を次式(１５)で定義する。
【００５８】
【数１０】
式（１４）において、ｆ(x)は計測線上のピーク信号またはエッジ信号である。また、上記「６．ピーク信号またはエッジ信号の取得」で述べたように、信号両端部は、信号自身を折り返すことによって適切に処理されている。ψ(x)は「７．マザーウェーブレットの指定」で述べたマザーウェーブレットである。２jはスケールパラメータで、レベル（スケールレベル）ｊは０または負の整数である。このスケールパラメータによって、マザーウエーブレットの幅および高さを同一面積に保持したまま２のべき乗則にしたがって変化させることができる（図９参照）。ｋはトランスレートパラメータを表している。なお、式(１５)ではスケーリングパラメータを２jとしているが、これに限られず、ｂj（ただし、ｂ＞１）としておいて次式(１６)のように表すこともできる。
【００５９】
【数１１】
式(１６)では、１＜ｂ＜２のとき、スケーリングパラメータのステップ間隔は式(１４)の時よりも細かくなり、２＜ｂのとき、スケーリングパラメータのステップ間隔は式(１５)の時よりも粗くなる。
【００６０】
我々の実験で確認されているように、通常は式(１５)で表される２のべき乗ステップで良好な検出性能を得ている。
【００６１】
式(１５)や式(１６)の意味するところは、大きさ（スケーリングパラメータによって規定される）と位置（トランスレートパラメータによって規定される）を変化させたマザーウエーブレットψ(x)と画像から形成される信号ｆ(x)の畳み込み積分演算である。この演算によって生成されるｄj(x)は、ｆ(x)の局所的な周期とマザーウェーブレットψ(x)のサポート（ψが０でない区間）が概ね一致したときに激しく波立つという特徴がある（図１０参照）。
【００６２】
図９や図１０では、説明および理解を容易にするためにレベルｊを３段階にしか変化させていないが、理論的にはレベルｊを８段階に変化させると、
【００６３】
【数１２】
となり、100倍以上の年輪幅に対応することができる。我々の実験でも、８段階に設定したレベルｊによって、年輪幅0.1mmから１cmを超える箇所まで対応できることが確認された。したがって、レベルｊの段階数は、計測対象物の特性に応じて任意に設定することができる。
【００６４】
９．区間二乗平均の演算（ステップＳ１６）
「８．ウェーブレット変換」で生成されたｄj(x)に対して、次式(１８)に示すように、区間二乗平均を演算する（図１１参照）。
【００６５】
【数１３】
ただし、ａはマザーウエーブレットψ(x)のサポートによって決定される定数である。
【００６６】
式(１８)では、二乗平均を求める区間は、ｘ−２-jａからｘ+２-jａまでになっているが、これに限定されず、次式(１９)のように定義してもよい。
【００６７】
【数１４】
すなわち、ｐjはスケールレベルｊが低レベルになるほど大きくなるようにスケールレベルｊに応じて設定されるので、区間二乗平均の演算は、レベルｊの高い高周波のｄj(x)に対しては狭い区間[ｘ−ｐj, ｘ＋ｐj]で、レベルｊの低い低周波のｄj(x)に対しては広い区間[ｘ−ｐj, ｘ＋ｐj]でそれぞれ求めればよい。式(１８)に示すｘ−２-jａからｘ+２-jａまでというのは、あくまでも標準となる積分区間である。積分区間[ｘ−ｐj, ｘ＋ｐj]をこれよりも狭くすると、局所的な濃度変化に対して鋭敏になりやすいがノイズの影響も受けやすくなる。一方、積分区間[ｘ−ｐj, ｘ＋ｐj]をこれよりも広くすると、局所的な濃度変化に対しては鋭敏ではなくなるものの、ノイズに対して強い安定志向になる傾向がある。したがって、積分区間[ｘ−ｐj, ｘ＋ｐj]を上記の標準的な状態から加減してやることによって、検出の感度調節を図ることもできる。
【００６８】
１０．支配レベルの決定（ステップＳ１７）
ある点ｘにおいて、式(１８)または式(１９)で定義されたｇj(ｘ)が最大となるレベルｊを点ｘにおける支配レベルｊdとする。これにより、計測線上のすべての点が、ある一つの支配レベルを持つことになる（図１２参照）。支配レベルｊdは、年輪幅が広い部分では低くなる傾向があるのに対し、年輪幅が狭い部分や晩材周辺の信号が著しく変化する部分では高くなる傾向がある。
【００６９】
１１．区間内最大値（エッジ信号の場合は区間内最小値）の決定（ステップＳ１８〜Ｓ２０）
図１３に示すように、
（ａ）ｆ(x)がピーク信号の場合には、ｆ(x)の値が区間[ｘ−２-jdａ, ｘ＋２-jdａ]における ｆ(x)の最大値ｆmax(x)と等しければ、点ｘを年輪層内の最大濃度箇所と決定する。
【００７０】
（ｂ）ｆ(x)がエッジ信号の場合には、ｆ(x)の値が区間[ｘ−２-jdａ, ｘ＋２-jdａ]におけるｆ(x)の最小値ｆmin(x)と等しければ、点ｘ年輪層内の晩材終端と決定する。
【００７１】
これを全範囲にわたって行うことで、計測線上のすべての年輪を検出することができる。
【００７２】
なお、上記（ａ）、（ｂ）は上記したスケールパラメータｂjにおけるｂが２の場合であるが、下記のように表すこともできる。
【００７３】
（ｃ）ｆ(x)がピーク信号の場合には、ｆ(x)の値が区間[ｘ−ｂ-jd, ｘ＋ｂ-jd]におけるｆ(x)の最大値ｆmax(x)と等しければ、点ｘを年輪層内の最大濃度箇所と決定する。
【００７４】
（ｄ）ｆ(x)がエッジ信号の場合には、ｆ(x)の値が区間[ｘ−ｂ-jd, ｘ＋ｂ-jd]におけるｆ(x)の最小値ｆmin(x)と等しければ、点ｘを年輪層内の晩材終端と決定する。
【００７５】
これをさらに一般化すると、支配レベルｊdに応じた定数ｑj ｄを用いて下記のように表すことができる。
【００７６】
（ｅ）ｆ(x)がピーク信号の場合には、ｆ(x)の値が区間[ｘ−ｑj ｄ, ｘ＋ｑj ｄ]におけるｆ(x)の最大値ｆmax(x)と等しければ、点ｘを年輪層内の最大濃度箇所と決定する。
【００７７】
（ｆ）ｆ(x)がエッジ信号の場合には、ｆ(x)の値が区間[ｘ−ｑj ｄ, ｘ＋ｑj ｄ]におけるｆ(x)の最小値ｆmin(x)と等しければ、点ｘを年輪層内の晩材終端と決定する。
【００７８】
１２．主計測線と副計測線の相互参照（ステップＳ２１）
上記「２．計測箇所の指定」で定めた主計側線１上の年輪が不明瞭であった場合には、上記「３．画素情報の取得」から「１１．区間内最大値（エッジ信号の場合は区間内最小値）の決定」の手順で各層の最大濃度箇所や晩材終端を検出しようとしても検出されないことがある。また、年輪層の部分的な不明瞭は、ランダムに発生する場合が多い。この状態を非検出と呼ぶことにする。
【００７９】
また、主計側線１上に木材組織の不均一な箇所があったり、割れや研磨痕などの本来の年輪以外濃度変化（以下ノイズと称する）があったりした場合には、最大濃度箇所や晩材終端以外の箇所を検出してしまうことがある。この状態を誤検出と呼ぶことにする。
【００８０】
通常、年輪は木口では同心円状に、柾目では平行状に生じるのに対し、ノイズの発生パターンは不確定的である。これらの違いを利用して、主計測線上に非検出や誤検出を生じても、正しい検出結果を出力できるように、「２．計測箇所の指定」で設定した副計測線上の検出結果（「３．画素情報の取得」から「１１． 区間内最大値（エッジ信号の場合は区間内最小値）の決定」のプロセスをそれぞれの副計測線２〜５上に対して行ったもの）と相互参照しながら、年輪として確定してもよいか否かを検証する。その具体的な手順を以下に記す。
【００８１】
実際のプログラム上では相互参照を効率よく進めるため、あらかじめ、すべての主および副計測線上の検出点の位置と計測線番号を確定候補点としてメモリーに記憶しておく。
【００８２】
木口面では、通常年輪は同心円状であるが、微小な幅の中で考えた場合には、ほぼ平行とみなすことができる。また、柾目では平行状である。ところが、いずれも計測線方向に直交する完全な平行状態ではないため、主計測線上の検出位置と副計測線上の検出位置とで若干のズレを生じることがある。本来同一年輪層であるべき検出点が、主または副計測線設定位置の違いによって生じる誤差を吸収するために、すなわち、各計測線の始端から略同一距離の範囲内にある検出点か否かを決定するために、同一とみなしうる幅ｈを設定する。
【００８３】
次に、複数存在する主および計測線上のうち、何箇所以上に検出点を生じた場合に、その検出点を正しいものとして確定してよいかを定義する。我々の実験では、国産檜材の場合、木口では主および副計測線合計３本中２箇所以上に検出点を生じたもの、柾目では主および副計測線合計５本中３箇所以上に検出点を生じたものとした場合、最も良好な結果が得られている。より一般的には、同一とみなしうる幅±ｈ以内で主および副計測線合計ｎ本中ｍ箇所以上に検出点のあるものと、定義することができる。
【００８４】
上記規則に従って、主計測線上の全検出点に対し、幅±ｈ以内で主および副計測線合計ｎ本中ｍ以上に検出点のあるものを、正しい年輪（各層内最大濃度箇所または各層内晩材終端）として確定していく(図１５参照)。これにより、
（１）主計測線上の検出点が正しい年輪であった場合には、相互参照により正しいものとして確定される。例えば、図１５の例では、検出点ａ1〜ａ3においては主計測線上に検出点があり且つ５本の計測線中３つ以上の検出点があるので、検出点が確定する。
【００８５】
（２）主計測線上の検出点が誤検出であった場合には、相互参照により排除される。例えば、図１５の例では、検出点ａ5は５本の計測線中３箇所未満であるため、確定しないで次の「１３．副計測線間の相互参照」を行う。
【００８６】
（３）主計測線上に非検出があった場合（図１５における検出点ａ4参照）には、この過程では考慮されないので、次の「１３．副計測線間の相互参照」を行う。
【００８７】
「１３．副計測線間の相互参照」においては、図１６に示すように、検出点ａ4においては３箇所以上に検出点があるので、優先順位１を検索中に確定する。しかし、検出点ａ5においては５本の計測線中３箇所未満であるので確定しない。
【００８８】
図１７に示す段階においては、検索点ａ1〜ａ4は確定済みであり、一方、検出点ａ5は検出点として確定せずに相互参照が終了する。
【００８９】
なお、次の「１３．副計測線間の相互参照」を効率よく行うために「１２．主計測線と副計測線の相互参照」の過程で正しい年輪として確定されたものは、確定候補点から順次消去する。
【００９０】
１３．副計測線間の相互参照（ステップＳ２２）
複数存在する副計測線２〜５に、あらかじめ優先順位をつける。通常は、主計測線１に近いほど、高い優先順位とする。このとき、同順位のないようにしておく。
【００９１】
まず、確定候補点として残っている検出点のなかで最も優先順位の高い副計測線上の全検出点を、「１２．主計測線と副計測線の相互参照」と同様の規則に従って正しい年輪として確定していく（図１６における検出点ａ4参照）。なお、この過程で正しい年輪として確定されたものは、確定候補点から順次消去する。
【００９２】
以下、優先順位の高いものから順に同様のことを行う。これらの一連のプロセスは、確定候補点がなくなった時点、または最も優先順位の低い副計測線に対して行った時点で終了する（図１７参照）。なお、実際には、ｎ本中ｍ以上の規則の場合には、最低順位を含む順位下位のｍ−１本に対しては無意味であるから、この順位で打切りにしてよい。
【００９３】
上記したように、「１２．主計測線と副計測線の相互参照」のプロセスで主計測線上に非検出があった場合にも、「１３．副計測線間の相互参照」を行うことで正しい年輪として確定される(図1５、図１６における検出点ａ4参照)。
【００９４】
１４．年輪箇所の決定（ステップＳ２３）
「１２．主計測線と副計測線の相互参照」と「１３．副計測線間の相互参照」の年輪確定点を統合して、年輪箇所（各層内最大濃度箇所または各層内晩材終端）とする（図１４参照）。なお、「１３．副計測線間の相互参照」の過程で確定された副計測線上の確定点は、統合の際に主計測線に射影しておく。
【００９５】
１５．年輪幅の計測および結果の出力（ステップＳ２４、Ｓ２５）
すべての年輪確定点の情報は主計測線１上に集められているため、主計測線上の各年輪確定点の座標値から各層の画素数をカウントしていく。「１．年輪画像の取得」で設定した画像解像度によって画素１個あたりの大きさが決められるので、
【００９６】
【数１５】
によって各層の年輪幅を求めることができる。この結果を、出力することにより、一連の全プロセスが完了する。
【００９７】
ところで、年輪年代学に有用なデータは、年輪幅だけではない。層内最大濃度（すなわち層内最大密度）、層内最小濃度（すなわち層内最小密度）、早材晩材比なども、重要なデータとなる。
【００９８】
上記したプロセスによる木材の年輪箇所検出方法の効果を図１８および図１９に示す。図１８は、波形信号に上記濃度信号を用いた場合の最大濃度点検出結果を従来方法による検出結果と比較し表にして示したものである。また、図１９は、波形信号に上記微分信号を用いた場合の最大晩材終端点検出結果を従来方法による検出結果と比較し表にして示したものである。図１８および図１９において、本発明による新しい方法は、ウェーブレット変換のみ（複数の計測線の相互参照は省略）の場合と、複数の計測線の相互参照のみ（ウェーブレット変換は無し）の場合と、ウェーブレット変換および複数の計測線の相互参照を併用した場合とに分けている。
【００９９】
また、図１８および図１９に示す検出率、誤検出率、非検出率は、図２０に示す、Ｓｓ、Ｓｎ、Ｎｓ、Ｎｎの関係から、次式（２１）、（２２）、（２３）のように定義したものである。
【０１００】
【数１６】
図１８および図１９からわかるように、本発明の検出方法を用いれば、いずれの場合も従来の方法に比べて検出率の増加、誤検出率および非検出率の低下が確認され、検出性能が全般的に向上することが実証された。
【０１０１】
以上、木材の年輪箇所を検出し年輪幅を計測する場合について詳述したが、「１４．年輪箇所の決定」までを正確に行うことができれば、ピーク信号やエッジ信号からこれら年輪幅以外の特徴量も容易に求めることができる。
【０１０２】
また、本発明方法の応用例として、例えば個人認証のための指紋、音紋、網膜等の計測、電子部品における配線パターンの計測、脳波信号等の生体信号の計測等にも適用することができる。
【０１０３】
【発明の効果】
以上の説明から明らかなように、本発明による木材の年輪箇所検出方法および年輪幅計測方法によれば、年輪画像から得られる輝度情報信号又は濃度情報信号が特徴点間すなわち年輪間の距離が不規則な波形信号であっても、それら輝度情報信号又は濃度情報信号である波形信号（ピーク信号）又はそれらを微分処理してなる波形信号（エッジ信号）の特徴点すなわち年輪箇所を正確に検出することができるので、木材の年輪箇所の検出性能および年輪幅計測性能を著しく向上することができる。
【０１０４】
すなわち、計測法にウェーブレット変換を導入する前は、ｆ(x)上にある閾値を定めて閾値以上か閾値以下かで区別する単純二値化で年輪検出を行っていたため、ノイズや不鮮明な年輪に対して非常に弱く、また閾値が個体差に左右されやすい欠点を持っていたが、本発明によるウェーブレット変換を用いた計測法を導入することにより、その場その場その場の年輪の細かさ（あるいは粗さ）や年輪の鮮明さ（あるいは不鮮明さ）に自動的に適応しながら年輪検出を行えるようになり、これが検出性能の向上に寄与するのである。
【０１０５】
また、この検出性能は、閾値設定不要なため、個体差に左右されにくい頑健な面も持ちあわせており、この点においてもきわめて画期的かつ効果的である。
【図面の簡単な説明】
【図１】 本発明の一実施形態に係る木材の年輪箇所検出方法および年輪幅計測方法のプロセス前半部分を示すフローチャートである。
【図２】 図１のフローチャートに続くプロセス後半部分を示すフローチャートである。
【図３】 主計測線及び副計測線の設定例を模式的に示す説明図である。
【図４】 計測線に沿って取得したＲＧＢ信号から濃度信号へ変換する様子を示す説明図であり、（ａ）はＲＢＧ信号の波形説明図、（ｂ）はＲＢＧ信号から作成した濃度信号の波形説明図、（ｃ）はＲＧＢ信号から濃度信号への変換に用いられるＲＧＢ−濃度変換曲線の説明図である。
【図５】 （ａ）は濃度信号の変化を示す波形説明図であり、（ｂ）は平均化処理した濃度信号を示す波形説明図である。
【図６】 （ａ）は濃度信号の波形説明図であり、（ｂ）は濃度信号の微分処理により作成された微分信号の波形説明図、（ｃ）は微分処理（差分処理）方法を説明する説明図である。
【図７】 濃度信号の両端部に延長ダミー信号を付加した例を示す波形説明図である。
【図８】 （ａ）はウェーブレット変換に用いられるフレンチハット型のマザーウェーブレットの代表例を示す説明図であり、（ｂ）はメキシカンハット型マザーウェーブレットの代表例を示す説明図である。
【図９】 フレンチハット型のマザーウェーブレット波形がレベル値に応じて段階的に変化する様子を示す説明図である。
【図１０】 濃度信号を３段階のレベル値を用いてウェーブレット変換することにより得られる信号波形を示す説明図である。
【図１１】 図１０のウェーブレット変換により得られる信号波形をレベル値に応じて区間二乗平均処理した結果を示す説明図である。
【図１２】 図１１の区間二乗平均処理結果から支配レベルを決定演算した結果を示す説明図である。
【図１３】 図１２の支配レベル決定演算結果を用いて濃度信号波形上の区間内最大値および微分信号波形上の区間内最小値を検索し決定するプロセスを模式的に示す説明図であり、（ａ）は支配レベル説明図、（ｂ）は濃度信号の波形説明図、（ｃ）は微分信号の波形説明図である。
【図１４】 （ａ）は濃度信号からの最大濃度点（特徴点）検出結果を示す説明図であり、（ｂ）は微分信号からの晩材終端点（特徴点）検出結果を示す説明図である。
【図１５】 主計測線および副計測線の相互参照による検出点（特徴点）確定プロセスを模式的に示す説明図である。
【図１６】 図１５のプロセスに続けて副計測線の相互参照による検出点（特徴点）確定プロセスを模式的に示す説明図である。
【図１７】 図１６のプロセスに続く相互参照による検出点（特徴点）確定プロセスの最終結果を模式的に示す説明図である。
【図１８】 濃度信号を用いた場合の本発明方法による最大濃度点（特徴点）検出結果を従来方法による検出結果と比較して示す説明図である。
【図１９】 濃度信号を用いた場合の本発明方法による晩材終端点（特徴点）検出結果を従来方法による検出結果と比較して示す説明図である。
【図２０】 検出率と誤検出率と非検出率の定義方法を説明するための説明図である。
【符号の説明】
１ 主計測線
２〜５ 副計測線
【発明の属する技術分野】
本発明は、波形上の特徴点の値または特徴点間距離が不規則な波形信号の特徴点を精度良く検出する方法に関し、更に詳しくは、例えば木材年輪数、年輪幅等の計測に好適に用いることができる木材の年輪箇所検出方法および年輪幅計測方法に関する。
【０００２】
【従来の技術】
年輪年代学では、年代確定の基準となる標準的な年輪幅変動のデータベースと照合することにより、各年輪層が何年に形成されたものかを一年単位で確定することができる。このデータベースは、奈良文化財研究所による精力的な調査研究によって、現在のところ日本では檜材で紀元前912年、杉材で紀元前1313年まで遡ることができる。ちなみに、年輪年代学の先進地域であるドイツでは、約一万年にわたる標準データベースが作成されている。年輪年代学では、主に以下のような事例が取り扱われている。
【０００３】
（ａ）木材の伐採年の推定
（ｂ）木造文化財（建築物、仏像、美術工芸品等）の制作年代や修理経緯の推定、真贋鑑定等
（ｃ）過去の長期にわたる気候変動の研究、地球環境温暖化の研究等
年輪年代学に使用される各年輪幅の時系列データは、誤検出（年輪以外のものを年輪として認識する誤り）も、非検出（年輪を認識しない誤り）も共に０という究極の検出性能が必要とされる。
【０００４】
このため、年輪幅の計測は、計測顕微鏡を用いた専用システムによって目視で行われており、その作業には熟練した技術と膨大な時間（300層クラスの試料で約１時間）を要していた。また、システム構成が大掛かりなものになるという問題もかかえていた。
【０００５】
計測作業を自動化するためには、デジタルカメラやスキャナ等の画像取得装置を用いて取得した年輪画像から、パソコンレベルのコンピュータを用いて各年輪幅の計測を行う方法が考えられており、これまでにもいくつかの試みがなされているが、前述の検出性能やシステム規模、価格などの点で問題があり、研究手段として一般的に広く用いられるには至っていないのが現状である。
【０００６】
とりわけ檜材については、文化財等に広く利用されていて年輪年代学上重要な樹種であるにもかかわらず、年輪幅が狭小であること、杉材に比べて年輪が不明瞭なことなどの条件から、自動計測の実用化が極めて困難であった。
【０００７】
本特許案件に類似した公知の技術の主要なもの以下に示す。
【０００８】
（１）「Win DENDRO」 カナダ Regent Instruments社 1988年
（http://www.regent.qc.ca.products/dendro/DENDRO.html参照）
ケベック大学のDr. Rejean Gagnon と Dr. Hubert Morin が設計し、Regent Instruments社が商品化した年輪年代学用ソフトである。このソフトは、年輪画像の輝度変化情報をもとに、年輪計測や木材組織の解析を行うことができる。このソフトのアルゴリズムの詳細は不明であるが、カタログ文面から察するところ、ウェーブレット処理および複数計測線情報の統合技術は用いていないと思われる。
【０００９】
（２）「画像処理システムを用いた年輪幅計測」 野田真人 1990年 樹木年輪研究会にて発表
この発表では、杉年輪の計測は可能であったが、檜年輪は計測不可能と結論づけている。この発表に係る計測方法においては、ウェーブレット処理と複数計測線情報の統合技術は用いられていない。
【００１０】
（３）特開平11−232427号公報
年輪数測定のために画像輝度情報を用いる点が記載されているが、同技術は上記（２）によりすでに公知である。この公報記載の公知技術においても、ウェーブレット処理と複数計測線情報の統合技術は何ら用いられていない。
【００１１】
【発明が解決しようとする課題】
年輪画像から画像上の計測線に沿って画素情報を取得することにより、輝度変化情報の波形信号または濃度変化情報の波形信号を得ることができる。濃度波形の極大箇所（輝度波形の場合は極小箇所）は各年輪層の各年輪層のうちの最も濃い部分（晩材最高濃度部）に対応する。したがって、濃度信号波形の極大箇所または輝度信号波形の極小箇所を認識していくことにより、各年輪層を識別することができる。
【００１２】
一方、濃度波形信号をさらに微分処理すると、暗部から明部へと移行する点（微分波形の極小箇所）が各年輪層の終端点（晩材終端）に対応することなる。したがって、計測線上でこの微分波形信号の極小箇所間距離を計測することにより、年輪幅をより正確に計測することができる。
【００１３】
年輪画像から画像上の計測線に沿って画素情報を取得することにより、輝度変化情報の波形信号または濃度変化情報の波形信号等を得ることは、上記した公知文献等により公知である。
【００１４】
しかしながら、木材の年輪幅は一般に不規則であり、最大年輪幅が最小年輪幅の数倍から約１００倍程度まで異なることも少なくない、したがって、年輪画像から得た波形信号に対する特徴点（極大箇所であるピーク点または極小箇所であるトラフ点）の検出精度を小さい特徴点間距離のレベルに合わせると、特徴点間距離の大きい領域では年輪に無関係なノイズを拾いやすくなる。一方、波形信号に対する特徴点（ピーク点またはトラフ点）の検出精度を大きい特徴点間距離のレベルに合わせると、特徴点間距離の小さい領域では特徴点の検出漏れが生じやすくなる。したがって、最小年輪幅と最大年輪幅との差が大きい計測対象物に対し年輪数や年輪幅を正確に計測することが困難となっている。
【００１５】
さらに、年輪画像においては、濃度のレベルが一定ではないため、年輪画像から得られる上記輝度波形信号、濃度波形信号、微分波形信号等はいずれも計測区間全体にわたりうねったような変化特性を有する場合が多い。このため、一定の閾値を用いて波形信号における特徴点（極大箇所であるピーク点または極小箇所であるトラフ点）を検出しようとすると、検出漏れが生じる原因となり、年輪数や年輪幅を正確に計測できない原因となる。
【００１６】
したがって、年輪年代学にとって最も基礎的なデータとなる各年輪幅の時系列データを迅速かつ高精度に得ることができる年輪箇所検出方法および年輪幅計測方法が要望されている。
【００１７】
本発明は、上記の要望に応え得る木材の年輪箇所検出方法および年輪幅計測方法を提供することを目的とするものであり、年輪画像から得た波形信号の特徴点をなす年輪箇所の値または該年輪箇所間の距離が不規則であっても、年輪画像から得た波形信号から波形の特徴点をなす年輪箇所を迅速かつ高精度に検出することができる木材の年輪箇所検出方法および該年輪箇所検出方法を用いて行う年輪幅計測方法を提供しようとするものである。
【００１８】
【課題を解決するための手段】
上記課題を解決するために、本発明は、第１発明として、撮像装置によって取得した木材の表面又は内部の年輪画像をコンピュータ装置により処理するコンピュータ処理工程を有する木材の年輪箇所検出方法であって、前記コンピュータ処理工程が、前記年輪画像上に開始点と終端点とで区間を指定して設定した計測線上に沿って前記年輪画像から画素情報を取得する画素情報取得工程と、前記画素情報を輝度情報信号または濃度情報信号である波形信号に変換する波形信号取得工程と、所定のマザーウェーブレットと、ゼロまたは負の整数値からなる複数段階のスケールレベルとを用いて前記波形信号をウェーブレット変換することにより各スケールレベル毎のウェーブレット変換信号を生成するウェーブレット変換工程と、前記計測線上の各点において前記ウェーブレット変換信号から前記各スケールレベルに応じて予め設定した区間幅における区間二乗平均値を算出する区間二乗平均値算出工程と、前記計測線上の各点において前記区間二乗平均値が最大となるウェーブレット変換信号のスケールレベルをその点における支配レベルとして決定する支配レベル決定工程と、前記計測線上の各点において前記各支配レベルに応じて予め設定した区間幅における前記波形信号が最大値となる点を前記計測線上の年輪箇所として決定する年輪箇 所決定工程と、を有することを特徴とする木材の年輪箇所検出方法を提供する。
【００１９】
さらに、第２発明として、撮像装置によって取得した木材の表面又は内部の年輪画像をコンピュータ装置により処理するコンピュータ処理工程を有する木材の年輪箇所検出方法であって、前記コンピュータ処理工程が、前記年輪画像上に開始点と終端点とで区間を指定して設定した計測線上に沿って前記年輪画像から画素情報を取得する波形信号取得工程と、前記画素情報を輝度情報信号または濃度情報信号に変換した後更に微分処理して波形信号を得る波形信号取得工程と、所定のマザーウェーブレットと、ゼロまたは負の整数値からなる複数段階のスケールレベルとを用いて前記波形信号をウェーブレット変換することにより各スケールレベル毎のウェーブレット変換信号を生成するウェーブレット変換工程と、前記計測線上の各点において前記ウェーブレット変換信号から前記各スケールレベルに応じて予め設定した区間幅における区間二乗平均値を算出する区間二乗平均値算出工程と、前記計測線上の各点において前記区間二乗平均値が最小となるウェーブレット変換信号のスケールレベルをその点における支配レベルとして決定する支配レベル決定工程と、前記計測線上の各点において前記各支配レベルに応じて予め設定した区間幅における前記波形信号が最小値となる点を前記計測線上の年輪箇所として決定する年輪箇所決定工程と、を有することを特徴とする木材の年輪箇所検出方法を提供する。
【００２０】
さらに第３発明として、上記第１発明又は第２発明の構成を有する木材の年輪箇所検出方法において、前記ウェーブレット変換に次式（６）、すなわち、
【数３】
を用い、さらに、前記区間二乗平均の演算に、次式（９）、すなわち、
【数４】
を用いることを特徴とする木材の年輪箇所検出方法を提供する。
【００２１】
上記構成において、好ましくは、前記区間二乗平均の演算式（９）におけるｐjが次式（１０）、すなわち、
【００２２】
【数５】
で定義されるものである。
【００２３】
さらに、第４発明として、上記第２発明の構成を有する木材の年輪箇所検出方法において、前記微分処理は、前記計測線上における複数画素分の間隔を隔てた複数の画素同士における輝度または濃度情報信号の差分演算処理であることを特徴とする請求項２記載の木材の年輪箇所検出方法を提供する。
【００２４】
さらに、第５発明として、上記第１発明又は第２発明の構成を有する木材の年輪箇所検出方法において、前記計測線が主計測線と、該主計測線の両側に各々等間隔を隔てて並列に設定される複数の副計測線とからなり、前記年輪箇所決定工程は、該主計測線および副計測線の始端から各々略同一距離の点において年輪候補点が検出されたとき、少なくとも該年輪候補点の数が主計測線および副計測線の本数に対し過半数であることを条件の１つとして、その年輪候補点を前記年輪箇所として決定することを特徴とする木材の年輪箇所検出方法を提供する。
【００２５】
さらに、第６発明として、上記第５発明の構成を有する木材の年輪箇所検出方法において、前記年輪箇決定工程は、前記主計測線の両側に各々２本の副計測線を設定し、該主計測線および副計測線の始端から各々略同一距離の点において年輪候補点が検出されたとき、少なくとも主計測線上と該主計測線に近い位置にある２本の副計測線のうちのいずれか一方上に年輪候補点があり且つ他のいずれかの副計測線上に年輪候補点がある場合、または、主計測線に近い位置にある２本の副計測線上と他のいずれかの副計測線上に年輪候補点がある場合に、その年輪候補点を前記年輪箇所として決定するものであることを特徴とする木材の年輪箇所検出方法を提供する。
【００２６】
さらに、第７発明として、上記第１発明又は第２発明の構成を有する木材の年輪箇所検出方法において、前記波形信号取得工程は、前記計測線に沿って取得した前記波形信号を前記計測線の周辺の画素の情報を用いて平滑化する工程を含むことを特徴とする年輪箇所検出方法を提供する。
【００２７】
さらに、第８発明として、上記第１発明又は第２発明の構成を有する木材の年輪箇所検出方法を用いて行う木材の年輪幅計測方法であって、前記コンピュータ処理工程が、前記年輪箇所決定工程により決定された前記計測線上の前記各年輪箇所間の画素数と１画素の大きさとから前記各年輪箇所間の年輪幅を算出する年輪幅算出工程を有することを特徴とする木材の年輪幅計測方法を提供する。
【００２８】
【発明の実施の形態】
以下、図面を参照して本発明の好ましい実施形態を説明する。
【００２９】
図１および図２は、本発明に係る木材の年輪箇所検出方法および年輪幅計測方法の一実施形態を示すフローチャートである。図１および図２を参照すると、木材の年輪箇所検出および年輪幅計測の手順は次のプロセスからなる。
【００３０】
１．年輪画像の取得（ステップＳ１）
２．計測箇所の指定（ステップＳ２）
３．画素情報の取得（ステップＳ３）
４．濃度変換（ステップＳ４、Ｓ５）
５．平滑化処理（ステップＳ６、Ｓ７）
６．ピーク信号またはエッジ信号の取得（ステップＳ８〜Ｓ１２）
７．マザーウェーブレットの指定（ステップＳ１３）
８．ウェーブレット変換（ステップＳ１４）
９．区間二乗平均の演算（ステップＳ１５、Ｓ１６）
10．支配レベルの決定（ステップＳ１７）
11．区間内最大値（エッジ信号の場合は区間内最小値）の決定（ステップＳ１８〜Ｓ２０）
12．主計測線と副計測線の相互参照（ステップＳ２１）
13．副計測線間の相互参照（ステップＳ２２）
14．年輪箇所の決定（ステップＳ２３）
15．年輪幅の計測および結果出力（ステップＳ２４、Ｓ２５）
また、本発明の上記実施形態は、上記ステップＳ2からステップＳ２５までのプロセスを実行することができるプログラムを、基本システムがインストールされたパソコンなどのコンピュータ装置にインストールすることにより、実行可能となる。
【００３１】
この実施形態におけるプログラムは、計測箇所の指定手段、画素情報取得手段、濃度変換手段、平滑化処理手段、ピーク信号またはエッジ信号取得手段、マザーウェーブレット指定手段、ウェーブレット変換手段、区間二乗平均の演算手段、支配レベル決定手段、区間内最大値（エッジ信号の場合は区間内最小値）の決定手段、主計測線と副計測線の相互参照手段、副計測線間の相互参照手段、年輪箇所の決定手段、年輪幅の計測手段および結果出力手段を有する。
【００３２】
次に、上記各プロセスについて詳述する。
【００３３】
１．年輪画像の取得（ステップＳ１）
デジタルカメラやスキャナ等の撮像装置を用いて年輪画像を取得する。例えば小型の移動可能な試料の場合にはスキャナを使用し、大型の試料には高解像度タイプのデジタルカメラを使用する。この際、最も注意すべき点は、所定の画像解像度を得ることである。檜材では、最も年輪間隔の狭小な部分で1層あたり約0.1mmという箇所を生ずることがある。このような部分を識別するためには、ナイキストの標本化定理により、１画素あたりの大きさをその半分の0.05mm以下にすべきである。
【００３４】
さらに安全と計測誤差を考慮すると、１画素あたりの大きさ約0.02mmに相当する1200dpiの解像度で画像取得を行うことが好ましい。大面積の試料でこの解像度を計測面全域にわたって保持するためには、一回あたりの撮影を小領域に限定し、画像を継ぎ合わせるという手法を用いることが好ましい。
【００３５】
本発明で扱う年輪画像としては、上記のようなデジタル画像取得装置を用いて試料表面から直接取得した画像のほかに、銀塩写真（通常のフィルムを使用する写真）等をデジタル変換した画像も同様に扱うことができる。また、X線写真、X線CT画像（X線断層写真）、MRI画像のような試料内部の年輪画像も扱うことができる。
【００３６】
２．計測箇所の指定（ステップＳ２）
年輪幅計測は、図３に示すように、少なくとも１本の計測線１に沿って行われる。木口（横断面）円盤上では、通常、計測線は樹木中心（図３中左方向）から樹皮方向（図３中右方向）に向かって設定される。また、柾目（縦断面）では、木理と直交するように樹木中心方向から樹皮方向に向かって設定される。本発明においては、好ましくは、認識性能の向上を図るため、１本の主計測線１に対し、それと平行に適切な間隔を置いて主計測線１の両側に各々複数の副計測線２、３、４、５を設定する（図３参照）。
【００３７】
副計測線の数は、木口で2本、柾目面で4本が適当であることが我々の実験によって確かめられている。ただし、試料による個体差などもあり、必ずしもこの限りではない。実際のプログラム上では、副計測線を設定する間隔と本数、および年輪画像上で主計測線１の開始点１ａと終端点１ｂの2点を指定することにより、計測箇所を指定可能とすることが好ましい。
【００３８】
３．画素情報の取得（ステップＳ３）
上記したように、主計測線１および副計測線２〜５に沿って各画像の画素情報（通常のパソコンの場合には、BGR3チャネルの8ビットデジタル信号）を取得する。画像取得に関しては、計測線上の画素情報のみでなく、後述の平滑化処理のために計測線周辺画素（具体的な周辺画素数は、平滑化の程度によって決定される）の情報も同時に取得することが好ましい。
【００３９】
なお、本発明方法による年輪認識では色情報を使用しないため、BGR信号の平均信号や混合信号、あるいは特定チャネルの信号のみ（たとえばG信号のみ）の輝度に基づく信号を使用することができる。通常は、B：G：R=1：1：1の割合で平均混合した信号を使用すればよい。
【００４０】
４．濃度変換（ステップＳ５）
上記方法で取得された画素情報は、各画素の輝度に起因しているため、値が大きくなるほど明るくなるという特性を有する。さらに、パソコンシステム上のRGB信号には、通常、CRTモニタの特性を補正するためにガンマ補正という非線形変換が施されている。一方、年輪の濃淡変化は年輪の細胞密度に起因しているため、密度あるいは密度に関連している画像濃度をもとに濃淡を記述するのが、年輪年代学の慣例である。
【００４１】
濃度記述では、値が大きくなるほど暗くなる。本発明による実施形態では、こうした年輪年代学の慣例に従い、RGB平均信号を濃度信号に変換している（図４参照）。変換方法としては、画像中の濃度既知部分（通常はステップタブレットチャート）から導出される変換式を用いることができる。変換式の一例をあげると、
【００４２】
【数６】
【００４３】
【数７】
【００４４】
しかし、本発明では後述するとおり、年輪信号波形の値そのものとは無関係に年輪識別を行うことができるため、RGB、輝度あるいは濃度いずれの信号を用いても、年輪識別性能に大差はない。したがって、処理ステップや時間短縮を図るためには、このプロセスは省略可能である。ただし、このプロセスを省略すると年輪信号波形は反転することに留意する必要がある。
【００４５】
５．平滑化処理（ステップＳ６、Ｓ７）
上記プロセスによって取得した主計測線１および副計測線２〜５上の各画素値に対し、周辺画素の情報を用いて平滑化を行う。平滑化を施すことによって、試料自体や測定機器に起因するノイズによる誤認識を軽減することができる（図５参照）。
【００４６】
平滑化方法には、移動平均処理または、移動メジアン処理を用いることができる。周辺画素数としては、例えば1200dpi解像度の場合、我々の実験では木口面で、計測線直交方向に5画素、計測線方向に3画素程度、柾目面で計測線直交方向に15画素、計測線方向に5画素程度が適切であったが、これらの画素数に限定されない。なお、平滑化処理を行わない場合には、平滑化処理範囲を計測線直交方向および計測線方向に各々１画素として設定する。
【００４７】
６．ピーク信号またはエッジ信号の取得（ステップＳ８〜Ｓ１２）
上記プロセスによって取得した信号は、波形の極大箇所が各年輪層の最も濃い部分（晩材最高濃度部）を表している。この濃度情報信号を、以下、ピーク信号と称することにする。このピーク信号の特徴点であるピーク点（波形の極大箇所）を認識していくことにより、各年輪層を識別することができる。
【００４８】
一方、このピーク信号の微分信号（実際の演算は差分信号）における特徴点であるトラフ点（波形の極小箇所）は、年輪画像の暗部（高濃度部）から明部（低濃度部）へと移行する点に対応し、各年輪層の終端点（晩材終端）に対応する。この微分信号を、以下、エッジ信号と称することにする。
【００４９】
通常の微分信号は、隣接画素どうしの差分演算によって得られるが、隣接差分演算ではノイズの影響を受けやすくなる。そこで、本発明の好ましい実施形態では、複数画素の間隔を隔てた複数の画素同士の差分演算を行うことにより、ノイズの軽減を図っている（図６中央の波形説明参照）。
【００５０】
ピーク信号（図６上部の波形）とエッジ信号（図６下部の波形）を比較すると、波形ノイズの影響などにより、各年輪の識別にはピーク信号の方が検出性能は良好であるが、年輪幅を正確に計測するという観点では晩材終端に対応しているエッジ信号の方が精度は高い。そこで、本発明の好ましい実施形態では、このような特徴を活かして、ピーク信号による認識、エッジ信号による認識、両者を併用した認識をメニュー上で選択可能にしている。
【００５１】
後述の「８．ウェーブレット変換」の箇所で述べるウェーブレット変換（畳み込み積分演算）では、ある１点の画素に対して周辺画素を含む演算を行う。したがって、計測線両端部（図３の開始点１ａおよび終端点１ｂ）では、画素不足による演算不能を防止するため、信号自身を折り返してダミーの信号（図７の波形の両端の点線部分参照）を作成しておくことが望ましい。ダミー信号として必要な長さは、後述の「７．マザーウェーブレットの指定」で指定するマザーウェーブレットのサポートと「８．ウェーブレット変換」で指定するレベル数に依存するが、最も長い端部所理の長さは最低位のレベルｊの時であるから、これに合わせてダミー信号を延長しておけば十分である。
【００５２】
７．マザーウェーブレットの指定（ステップＳ１３）
マザーウェーブレットとして用いることができる関数ψ(x)は、
（ａ）端部が０または０に収束する関数であること
（ｂ）全区間での総和（積分値）が０であること
の２つの条件を満たす関数である。
【００５３】
また次の条件は、必ずしも満たす必要はないが、可能であれば満たしていることが望ましい。
【００５４】
（ｃ）サポート（関数の値が０ではない区間）がコンパクトに纏まっていること
マザーウェーブレットの具体例として、次のフレンチハット型の関数（式(１３)）とメキシカンハット型の関数（式(１４)）をあげることができる（図８（ａ）、（ｂ）参照）。
【００５５】
【数８】
【数９】
【００５７】
８．ウェーブレット変換（ステップＳ１４）
１次元の画像信号ｆ(x)に対するウェーブレット変換を次式(１５)で定義する。
【００５８】
【数１０】
【００５９】
【数１１】
【００６０】
我々の実験で確認されているように、通常は式(１５)で表される２のべき乗ステップで良好な検出性能を得ている。
【００６１】
式(１５)や式(１６)の意味するところは、大きさ（スケーリングパラメータによって規定される）と位置（トランスレートパラメータによって規定される）を変化させたマザーウエーブレットψ(x)と画像から形成される信号ｆ(x)の畳み込み積分演算である。この演算によって生成されるｄj(x)は、ｆ(x)の局所的な周期とマザーウェーブレットψ(x)のサポート（ψが０でない区間）が概ね一致したときに激しく波立つという特徴がある（図１０参照）。
【００６２】
図９や図１０では、説明および理解を容易にするためにレベルｊを３段階にしか変化させていないが、理論的にはレベルｊを８段階に変化させると、
【００６３】
【数１２】
【００６４】
９．区間二乗平均の演算（ステップＳ１６）
「８．ウェーブレット変換」で生成されたｄj(x)に対して、次式(１８)に示すように、区間二乗平均を演算する（図１１参照）。
【００６５】
【数１３】
【００６６】
式(１８)では、二乗平均を求める区間は、ｘ−２-jａからｘ+２-jａまでになっているが、これに限定されず、次式(１９)のように定義してもよい。
【００６７】
【数１４】
【００６８】
１０．支配レベルの決定（ステップＳ１７）
ある点ｘにおいて、式(１８)または式(１９)で定義されたｇj(ｘ)が最大となるレベルｊを点ｘにおける支配レベルｊdとする。これにより、計測線上のすべての点が、ある一つの支配レベルを持つことになる（図１２参照）。支配レベルｊdは、年輪幅が広い部分では低くなる傾向があるのに対し、年輪幅が狭い部分や晩材周辺の信号が著しく変化する部分では高くなる傾向がある。
【００６９】
１１．区間内最大値（エッジ信号の場合は区間内最小値）の決定（ステップＳ１８〜Ｓ２０）
図１３に示すように、
（ａ）ｆ(x)がピーク信号の場合には、ｆ(x)の値が区間[ｘ−２-jdａ, ｘ＋２-jdａ]における ｆ(x)の最大値ｆmax(x)と等しければ、点ｘを年輪層内の最大濃度箇所と決定する。
【００７０】
（ｂ）ｆ(x)がエッジ信号の場合には、ｆ(x)の値が区間[ｘ−２-jdａ, ｘ＋２-jdａ]におけるｆ(x)の最小値ｆmin(x)と等しければ、点ｘ年輪層内の晩材終端と決定する。
【００７１】
これを全範囲にわたって行うことで、計測線上のすべての年輪を検出することができる。
【００７２】
なお、上記（ａ）、（ｂ）は上記したスケールパラメータｂjにおけるｂが２の場合であるが、下記のように表すこともできる。
【００７３】
（ｃ）ｆ(x)がピーク信号の場合には、ｆ(x)の値が区間[ｘ−ｂ-jd, ｘ＋ｂ-jd]におけるｆ(x)の最大値ｆmax(x)と等しければ、点ｘを年輪層内の最大濃度箇所と決定する。
【００７４】
（ｄ）ｆ(x)がエッジ信号の場合には、ｆ(x)の値が区間[ｘ−ｂ-jd, ｘ＋ｂ-jd]におけるｆ(x)の最小値ｆmin(x)と等しければ、点ｘを年輪層内の晩材終端と決定する。
【００７５】
これをさらに一般化すると、支配レベルｊdに応じた定数ｑj ｄを用いて下記のように表すことができる。
【００７６】
（ｅ）ｆ(x)がピーク信号の場合には、ｆ(x)の値が区間[ｘ−ｑj ｄ, ｘ＋ｑj ｄ]におけるｆ(x)の最大値ｆmax(x)と等しければ、点ｘを年輪層内の最大濃度箇所と決定する。
【００７７】
（ｆ）ｆ(x)がエッジ信号の場合には、ｆ(x)の値が区間[ｘ−ｑj ｄ, ｘ＋ｑj ｄ]におけるｆ(x)の最小値ｆmin(x)と等しければ、点ｘを年輪層内の晩材終端と決定する。
【００７８】
１２．主計測線と副計測線の相互参照（ステップＳ２１）
上記「２．計測箇所の指定」で定めた主計側線１上の年輪が不明瞭であった場合には、上記「３．画素情報の取得」から「１１．区間内最大値（エッジ信号の場合は区間内最小値）の決定」の手順で各層の最大濃度箇所や晩材終端を検出しようとしても検出されないことがある。また、年輪層の部分的な不明瞭は、ランダムに発生する場合が多い。この状態を非検出と呼ぶことにする。
【００７９】
また、主計側線１上に木材組織の不均一な箇所があったり、割れや研磨痕などの本来の年輪以外濃度変化（以下ノイズと称する）があったりした場合には、最大濃度箇所や晩材終端以外の箇所を検出してしまうことがある。この状態を誤検出と呼ぶことにする。
【００８０】
通常、年輪は木口では同心円状に、柾目では平行状に生じるのに対し、ノイズの発生パターンは不確定的である。これらの違いを利用して、主計測線上に非検出や誤検出を生じても、正しい検出結果を出力できるように、「２．計測箇所の指定」で設定した副計測線上の検出結果（「３．画素情報の取得」から「１１． 区間内最大値（エッジ信号の場合は区間内最小値）の決定」のプロセスをそれぞれの副計測線２〜５上に対して行ったもの）と相互参照しながら、年輪として確定してもよいか否かを検証する。その具体的な手順を以下に記す。
【００８１】
実際のプログラム上では相互参照を効率よく進めるため、あらかじめ、すべての主および副計測線上の検出点の位置と計測線番号を確定候補点としてメモリーに記憶しておく。
【００８２】
木口面では、通常年輪は同心円状であるが、微小な幅の中で考えた場合には、ほぼ平行とみなすことができる。また、柾目では平行状である。ところが、いずれも計測線方向に直交する完全な平行状態ではないため、主計測線上の検出位置と副計測線上の検出位置とで若干のズレを生じることがある。本来同一年輪層であるべき検出点が、主または副計測線設定位置の違いによって生じる誤差を吸収するために、すなわち、各計測線の始端から略同一距離の範囲内にある検出点か否かを決定するために、同一とみなしうる幅ｈを設定する。
【００８３】
次に、複数存在する主および計測線上のうち、何箇所以上に検出点を生じた場合に、その検出点を正しいものとして確定してよいかを定義する。我々の実験では、国産檜材の場合、木口では主および副計測線合計３本中２箇所以上に検出点を生じたもの、柾目では主および副計測線合計５本中３箇所以上に検出点を生じたものとした場合、最も良好な結果が得られている。より一般的には、同一とみなしうる幅±ｈ以内で主および副計測線合計ｎ本中ｍ箇所以上に検出点のあるものと、定義することができる。
【００８４】
上記規則に従って、主計測線上の全検出点に対し、幅±ｈ以内で主および副計測線合計ｎ本中ｍ以上に検出点のあるものを、正しい年輪（各層内最大濃度箇所または各層内晩材終端）として確定していく(図１５参照)。これにより、
（１）主計測線上の検出点が正しい年輪であった場合には、相互参照により正しいものとして確定される。例えば、図１５の例では、検出点ａ1〜ａ3においては主計測線上に検出点があり且つ５本の計測線中３つ以上の検出点があるので、検出点が確定する。
【００８５】
（２）主計測線上の検出点が誤検出であった場合には、相互参照により排除される。例えば、図１５の例では、検出点ａ5は５本の計測線中３箇所未満であるため、確定しないで次の「１３．副計測線間の相互参照」を行う。
【００８６】
（３）主計測線上に非検出があった場合（図１５における検出点ａ4参照）には、この過程では考慮されないので、次の「１３．副計測線間の相互参照」を行う。
【００８７】
「１３．副計測線間の相互参照」においては、図１６に示すように、検出点ａ4においては３箇所以上に検出点があるので、優先順位１を検索中に確定する。しかし、検出点ａ5においては５本の計測線中３箇所未満であるので確定しない。
【００８８】
図１７に示す段階においては、検索点ａ1〜ａ4は確定済みであり、一方、検出点ａ5は検出点として確定せずに相互参照が終了する。
【００８９】
なお、次の「１３．副計測線間の相互参照」を効率よく行うために「１２．主計測線と副計測線の相互参照」の過程で正しい年輪として確定されたものは、確定候補点から順次消去する。
【００９０】
１３．副計測線間の相互参照（ステップＳ２２）
複数存在する副計測線２〜５に、あらかじめ優先順位をつける。通常は、主計測線１に近いほど、高い優先順位とする。このとき、同順位のないようにしておく。
【００９１】
まず、確定候補点として残っている検出点のなかで最も優先順位の高い副計測線上の全検出点を、「１２．主計測線と副計測線の相互参照」と同様の規則に従って正しい年輪として確定していく（図１６における検出点ａ4参照）。なお、この過程で正しい年輪として確定されたものは、確定候補点から順次消去する。
【００９２】
以下、優先順位の高いものから順に同様のことを行う。これらの一連のプロセスは、確定候補点がなくなった時点、または最も優先順位の低い副計測線に対して行った時点で終了する（図１７参照）。なお、実際には、ｎ本中ｍ以上の規則の場合には、最低順位を含む順位下位のｍ−１本に対しては無意味であるから、この順位で打切りにしてよい。
【００９３】
上記したように、「１２．主計測線と副計測線の相互参照」のプロセスで主計測線上に非検出があった場合にも、「１３．副計測線間の相互参照」を行うことで正しい年輪として確定される(図1５、図１６における検出点ａ4参照)。
【００９４】
１４．年輪箇所の決定（ステップＳ２３）
「１２．主計測線と副計測線の相互参照」と「１３．副計測線間の相互参照」の年輪確定点を統合して、年輪箇所（各層内最大濃度箇所または各層内晩材終端）とする（図１４参照）。なお、「１３．副計測線間の相互参照」の過程で確定された副計測線上の確定点は、統合の際に主計測線に射影しておく。
【００９５】
１５．年輪幅の計測および結果の出力（ステップＳ２４、Ｓ２５）
すべての年輪確定点の情報は主計測線１上に集められているため、主計測線上の各年輪確定点の座標値から各層の画素数をカウントしていく。「１．年輪画像の取得」で設定した画像解像度によって画素１個あたりの大きさが決められるので、
【００９６】
【数１５】
【００９７】
ところで、年輪年代学に有用なデータは、年輪幅だけではない。層内最大濃度（すなわち層内最大密度）、層内最小濃度（すなわち層内最小密度）、早材晩材比なども、重要なデータとなる。
【００９８】
上記したプロセスによる木材の年輪箇所検出方法の効果を図１８および図１９に示す。図１８は、波形信号に上記濃度信号を用いた場合の最大濃度点検出結果を従来方法による検出結果と比較し表にして示したものである。また、図１９は、波形信号に上記微分信号を用いた場合の最大晩材終端点検出結果を従来方法による検出結果と比較し表にして示したものである。図１８および図１９において、本発明による新しい方法は、ウェーブレット変換のみ（複数の計測線の相互参照は省略）の場合と、複数の計測線の相互参照のみ（ウェーブレット変換は無し）の場合と、ウェーブレット変換および複数の計測線の相互参照を併用した場合とに分けている。
【００９９】
また、図１８および図１９に示す検出率、誤検出率、非検出率は、図２０に示す、Ｓｓ、Ｓｎ、Ｎｓ、Ｎｎの関係から、次式（２１）、（２２）、（２３）のように定義したものである。
【０１００】
【数１６】
【０１０１】
以上、木材の年輪箇所を検出し年輪幅を計測する場合について詳述したが、「１４．年輪箇所の決定」までを正確に行うことができれば、ピーク信号やエッジ信号からこれら年輪幅以外の特徴量も容易に求めることができる。
【０１０２】
また、本発明方法の応用例として、例えば個人認証のための指紋、音紋、網膜等の計測、電子部品における配線パターンの計測、脳波信号等の生体信号の計測等にも適用することができる。
【０１０３】
【発明の効果】
以上の説明から明らかなように、本発明による木材の年輪箇所検出方法および年輪幅計測方法によれば、年輪画像から得られる輝度情報信号又は濃度情報信号が特徴点間すなわち年輪間の距離が不規則な波形信号であっても、それら輝度情報信号又は濃度情報信号である波形信号（ピーク信号）又はそれらを微分処理してなる波形信号（エッジ信号）の特徴点すなわち年輪箇所を正確に検出することができるので、木材の年輪箇所の検出性能および年輪幅計測性能を著しく向上することができる。
【０１０４】
すなわち、計測法にウェーブレット変換を導入する前は、ｆ(x)上にある閾値を定めて閾値以上か閾値以下かで区別する単純二値化で年輪検出を行っていたため、ノイズや不鮮明な年輪に対して非常に弱く、また閾値が個体差に左右されやすい欠点を持っていたが、本発明によるウェーブレット変換を用いた計測法を導入することにより、その場その場その場の年輪の細かさ（あるいは粗さ）や年輪の鮮明さ（あるいは不鮮明さ）に自動的に適応しながら年輪検出を行えるようになり、これが検出性能の向上に寄与するのである。
【０１０５】
また、この検出性能は、閾値設定不要なため、個体差に左右されにくい頑健な面も持ちあわせており、この点においてもきわめて画期的かつ効果的である。
【図面の簡単な説明】
【図１】 本発明の一実施形態に係る木材の年輪箇所検出方法および年輪幅計測方法のプロセス前半部分を示すフローチャートである。
【図２】 図１のフローチャートに続くプロセス後半部分を示すフローチャートである。
【図３】 主計測線及び副計測線の設定例を模式的に示す説明図である。
【図４】 計測線に沿って取得したＲＧＢ信号から濃度信号へ変換する様子を示す説明図であり、（ａ）はＲＢＧ信号の波形説明図、（ｂ）はＲＢＧ信号から作成した濃度信号の波形説明図、（ｃ）はＲＧＢ信号から濃度信号への変換に用いられるＲＧＢ−濃度変換曲線の説明図である。
【図５】 （ａ）は濃度信号の変化を示す波形説明図であり、（ｂ）は平均化処理した濃度信号を示す波形説明図である。
【図６】 （ａ）は濃度信号の波形説明図であり、（ｂ）は濃度信号の微分処理により作成された微分信号の波形説明図、（ｃ）は微分処理（差分処理）方法を説明する説明図である。
【図７】 濃度信号の両端部に延長ダミー信号を付加した例を示す波形説明図である。
【図８】 （ａ）はウェーブレット変換に用いられるフレンチハット型のマザーウェーブレットの代表例を示す説明図であり、（ｂ）はメキシカンハット型マザーウェーブレットの代表例を示す説明図である。
【図９】 フレンチハット型のマザーウェーブレット波形がレベル値に応じて段階的に変化する様子を示す説明図である。
【図１０】 濃度信号を３段階のレベル値を用いてウェーブレット変換することにより得られる信号波形を示す説明図である。
【図１１】 図１０のウェーブレット変換により得られる信号波形をレベル値に応じて区間二乗平均処理した結果を示す説明図である。
【図１２】 図１１の区間二乗平均処理結果から支配レベルを決定演算した結果を示す説明図である。
【図１３】 図１２の支配レベル決定演算結果を用いて濃度信号波形上の区間内最大値および微分信号波形上の区間内最小値を検索し決定するプロセスを模式的に示す説明図であり、（ａ）は支配レベル説明図、（ｂ）は濃度信号の波形説明図、（ｃ）は微分信号の波形説明図である。
【図１４】 （ａ）は濃度信号からの最大濃度点（特徴点）検出結果を示す説明図であり、（ｂ）は微分信号からの晩材終端点（特徴点）検出結果を示す説明図である。
【図１５】 主計測線および副計測線の相互参照による検出点（特徴点）確定プロセスを模式的に示す説明図である。
【図１６】 図１５のプロセスに続けて副計測線の相互参照による検出点（特徴点）確定プロセスを模式的に示す説明図である。
【図１７】 図１６のプロセスに続く相互参照による検出点（特徴点）確定プロセスの最終結果を模式的に示す説明図である。
【図１８】 濃度信号を用いた場合の本発明方法による最大濃度点（特徴点）検出結果を従来方法による検出結果と比較して示す説明図である。
【図１９】 濃度信号を用いた場合の本発明方法による晩材終端点（特徴点）検出結果を従来方法による検出結果と比較して示す説明図である。
【図２０】 検出率と誤検出率と非検出率の定義方法を説明するための説明図である。
【符号の説明】
１ 主計測線
２〜５ 副計測線
Claims (8)
- 撮像装置によって取得した木材の表面又は内部の年輪画像をコンピュータ装置により処理するコンピュータ処理工程を有する木材の年輪箇所検出方法であって、
前記コンピュータ処理工程が、
前記年輪画像上に開始点と終端点とで区間を指定して設定した計測線上に沿って前記年輪画像から画素情報を取得する画素情報取得工程と、
前記画素情報を輝度情報信号または濃度情報信号である波形信号に変換する波形信号取得工程と、
所定のマザーウェーブレットと、ゼロまたは負の整数値からなる複数段階のスケールレベルとを用いて前記波形信号をウェーブレット変換することにより各スケールレベル毎のウェーブレット変換信号を生成するウェーブレット変換工程と、
前記計測線上の各点において前記ウェーブレット変換信号から前記各スケールレベルに応じて予め設定した区間幅における区間二乗平均値を算出する区間二乗平均値算出工程と、
前記計測線上の各点において前記区間二乗平均値が最大となるウェーブレット変換信号のスケールレベルをその点における支配レベルとして決定する支配レベル決定工程と、
前記計測線上の各点において前記各支配レベルに応じて予め設定した区間幅における前記波形信号が最大値となる点を前記計測線上の年輪箇所として決定する年輪箇所決定工程と、
を有することを特徴とする木材の年輪箇所検出方法。 - 撮像装置によって取得した木材の表面又は内部の年輪画像をコンピュータ装置により処理するコンピュータ処理工程を有する木材の年輪箇所検出方法であって、
前記コンピュータ処理工程が、
前記年輪画像上に開始点と終端点とで区間を指定して設定した計測線上に沿って前記年輪画像から画素情報を取得する波形信号取得工程と、
前記画素情報を輝度情報信号または濃度情報信号に変換した後更に微分処理して波形信号を得る波形信号取得工程と、
所定のマザーウェーブレットと、ゼロまたは負の整数値からなる複数段階のスケールレベルとを用いて前記波形信号をウェーブレット変換することにより各スケールレベル毎のウェーブレット変換信号を生成するウェーブレット変換工程と、
前記計測線上の各点において前記ウェーブレット変換信号から前記各スケールレベルに応じて予め設定した区間幅における区間二乗平均値を算出する区間二乗平均値算出工程と、
前記計測線上の各点において前記区間二乗平均値が最小となるウェーブレット変換信号のスケールレベルをその点における支配レベルとして決定する支配レベル決定工程と、
前記計測線上の各点において前記各支配レベルに応じて予め設定した区間幅における前記波形信号が最小値となる点を前記計測線上の年輪箇所として決定する年輪箇所決定工程と、
を有することを特徴とする木材の年輪箇所検出方法。 - 前記ウェーブレット変換に次式（１）、すなわち、
を用い、さらに、前記区間二乗平均の演算に、次式（４）、すなわち、
を用いることを特徴とする請求項１又は２記載の木材の年輪箇所検出方法。 - 前記微分処理は、前記計測線上における複数画素分の間隔を隔てた複数の画素同士における輝度または濃度情報信号の差分演算処理であることを特徴とする請求項２記載の木材の年輪箇所検出方法。
- 前記計測線が主計測線と、該主計測線の両側に各々等間隔を隔てて並列に設定される複数の副計測線とからなり、前記年輪箇所決定工程は、該主計測線および副計測線の始端から各々略同一距離の点において年輪候補点が検出されたとき、少なくとも該年輪候補点の数が主計測線および副計測線の本数に対し過半数であることを条件の１つとして、その年輪候補点を前記年輪箇所として決定することを特徴とする請求項１又は２記載の木材の年輪箇所検出方法。
- 前記年輪箇決定工程は、前記主計測線の両側に各々２本の副計測線を設定し、該主計測線および副計測線の始端から各々略同一距離の点において年輪候補点が検出されたとき、少なくとも主計測線上と該主計測線に近い位置にある２本の副計測線のうちのいずれか一方上に年輪候補点があり且つ他のいずれかの副計測線上に年輪候補点がある場合、または、主計測線に近い位置にある２本の副計測線上と他のいずれかの副計測線上に年輪候補点がある場合に、その年輪候補点を前記年輪箇所として決定するものであることを特徴とする請求項５記載の木材の年輪箇所検出方法。
- 前記波形信号取得工程は、前記計測線に沿って取得した前記波形信号を前記計測線の周辺の画素の情報を用いて平滑化する工程を含むことを特徴とする請求項１または２記載の木材の年輪箇所検出方法。
- 請求項１または２に記載した木材の年輪箇所検出方法を用いて行う木材の年輪幅計測方法であって、前記コンピュータ処理工程が、前記年輪箇所決定工程により決定された前記計測線上の前記各年輪箇所間の画素数と１画素の大きさとから前記各年輪箇所間の年輪幅を算出する年輪幅算出工程を有することを特徴とする木材の年輪幅計測方法。
Priority Applications (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002330131A JP4218824B2 (ja) | 2002-11-13 | 2002-11-13 | 木材の年輪箇所検出方法および年輪幅計測方法 |
| US10/535,068 US20060140484A1 (en) | 2002-11-13 | 2003-05-13 | Method and equipment for measuring feature points of wave signal |
| US12/380,358 US20100067804A1 (en) | 2002-11-13 | 2009-02-26 | Method and equipment for measuring feature points of wave signal |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002330131A JP4218824B2 (ja) | 2002-11-13 | 2002-11-13 | 木材の年輪箇所検出方法および年輪幅計測方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2004164346A JP2004164346A (ja) | 2004-06-10 |
| JP4218824B2 true JP4218824B2 (ja) | 2009-02-04 |
Family
ID=32807893
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2002330131A Expired - Fee Related JP4218824B2 (ja) | 2002-11-13 | 2002-11-13 | 木材の年輪箇所検出方法および年輪幅計測方法 |
Country Status (2)
| Country | Link |
|---|---|
| US (2) | US20060140484A1 (ja) |
| JP (1) | JP4218824B2 (ja) |
Families Citing this family (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2007139699A (ja) * | 2005-11-22 | 2007-06-07 | Tokyo Electric Power Co Inc:The | 周波数解析方法 |
| US7643953B2 (en) * | 2007-03-07 | 2010-01-05 | Microsoft Corporation | Determination of incremental value in server processed data |
| US7945098B2 (en) * | 2007-05-17 | 2011-05-17 | The United States Of America As Represented By The Secretary Of Agriculture | Method for characterizing the density and cross-section morphology of trees |
| US8065529B2 (en) * | 2007-05-21 | 2011-11-22 | Ut-Battelle, Llc | Methods for using a biometric parameter in the identification of persons |
| JP5384316B2 (ja) * | 2009-12-24 | 2014-01-08 | 株式会社パスコ | 変位計測装置、変位計測方法、及び変位計測プログラム |
| CN102394064B (zh) * | 2011-08-30 | 2013-06-05 | 浙江大学 | 一种组合式植物声频调控方法 |
| US8755578B1 (en) | 2011-10-17 | 2014-06-17 | Igor Vladimir Smolyar | System and method for quantification of size and anisotropic structure of layered patterns |
| WO2016008529A1 (de) * | 2014-07-17 | 2016-01-21 | Hamberger Wolfgang | Identifizierung von gegenständen aus holz mittels der vermessung von jahresringen |
| US10819881B1 (en) | 2015-03-12 | 2020-10-27 | Igor Vladimir Smolyar | System and method for encryption/decryption of 2-D and 3-D arbitrary images |
| CN105160688B (zh) * | 2015-10-22 | 2018-05-25 | 上海斐讯数据通信技术有限公司 | 应用于移动终端的方向判断系统及方法 |
| US20200378933A1 (en) * | 2016-05-30 | 2020-12-03 | Shimadzu Corporation | Peak detection method and data processing device |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3624462B2 (ja) * | 1995-03-27 | 2005-03-02 | アイシン精機株式会社 | 車両の乗員保護装置 |
| JPH09126932A (ja) * | 1995-10-31 | 1997-05-16 | Aisin Seiki Co Ltd | タイヤの空気圧検出方法及び空気圧検出装置 |
| JP3648023B2 (ja) * | 1997-07-31 | 2005-05-18 | 株式会社山武 | データ変換装置およびプログラム記憶媒体 |
| EP1024435A1 (en) * | 1999-01-28 | 2000-08-02 | Atr Human Information Processing Research Laboratories | A mellin-transform information extractor for vibration sources |
-
2002
- 2002-11-13 JP JP2002330131A patent/JP4218824B2/ja not_active Expired - Fee Related
-
2003
- 2003-05-13 US US10/535,068 patent/US20060140484A1/en not_active Abandoned
-
2009
- 2009-02-26 US US12/380,358 patent/US20100067804A1/en not_active Abandoned
Also Published As
| Publication number | Publication date |
|---|---|
| US20100067804A1 (en) | 2010-03-18 |
| JP2004164346A (ja) | 2004-06-10 |
| US20060140484A1 (en) | 2006-06-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US20100067804A1 (en) | Method and equipment for measuring feature points of wave signal | |
| JP5117982B2 (ja) | 情報抽出装置、情報抽出方法、プログラムおよび記録媒体 | |
| EP1624412B1 (en) | Biological information measuring device | |
| EP1958158B1 (en) | Method for detecting streaks in digital images | |
| JP2017053819A (ja) | コンクリートのひび割れ検出方法及び検出プログラム | |
| US20100277774A1 (en) | Image quality indicator responsive to image processing | |
| EP2104059A1 (en) | Apparatus and method for processing image | |
| JP6823486B2 (ja) | ひび割れ検出方法 | |
| JP4870016B2 (ja) | ひび割れ検出方法 | |
| KR910017328A (ko) | 위치인식장치 | |
| JP5680458B2 (ja) | 原稿傾き判定装置、原稿傾き判定方法、プログラムおよびその記録媒体 | |
| US7046837B2 (en) | System and method for locating irregular edges in image data | |
| US11493453B2 (en) | Belt inspection system, belt inspection method, and recording medium for belt inspection program | |
| JP3696212B2 (ja) | パターン認識における照合に用いられる画像の生成ならびに同画像を用いたパターン認識のための方法、装置、およびプログラム | |
| JP2008116206A (ja) | パターン寸法測定装置、方法及びプログラム | |
| JP5413283B2 (ja) | 焦点ズレ検出装置、焦点ズレ検出方法およびプログラム | |
| CN114140620A (zh) | 一种物体直线轮廓检测方法 | |
| US6597805B1 (en) | Visual inspection method for electronic device, visual inspecting apparatus for electronic device, and record medium for recording program which causes computer to perform visual inspecting method for electronic device | |
| CN104318580B (zh) | 一种模式搜索方法及装置 | |
| KR101334029B1 (ko) | 두경부 근육 추출 방법 및 ｒｏｉ분석 방법 | |
| JP2638921B2 (ja) | 画像評価装置 | |
| JP2021189238A (ja) | ベルト検査システムおよびベルト検査プログラム | |
| JP2017534949A (ja) | 走査動作中にターゲットオブジェクトが存在したかどうかの決定 | |
| CN117237617B (zh) | 一种架空输电线图像识别方法、装置、设备及存储介质 | |
| CN111899237B (zh) | 标尺精度测量方法、装置、计算机设备和存储介质 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20051104 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20080516 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20080714 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20081105 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20081106 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111121 Year of fee payment: 3 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121121 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20151121 Year of fee payment: 7 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| LAPS | Cancellation because of no payment of annual fees |