広島工業大学紀要研究編
第 44 巻（2010）pp.291-296

論文

稲の生育診断のための圃場写真の変換と合成
田中健太郎 *・山下英生 **（平成21年10月29日受理）

Transformation and Composition of Field Pictures
for Estimation of Plant Growth in Paddy Field
Kentaro TANAKA and Hideo YAMASHITA(Received Oct. 29, 2009)

Abstract
It is important to agricultural production to estimate a growth state of rice in paddy field.
The system that diagnoses the growth of rice by the image analysis is proposed. This systemtransforms bird’s-eye view images of paddy fields to top views of a parallel projection and analyzes them. Image analysis is performed for each image. Therefore, there is a possibility that
we cannot diagnose all farms within a limited period in case a farm to diagnose is wide. Analytical time is shortened by composing plural plane views to one piece of image.
In this paper, transformation and composition methods of field pictures are proposed. The
angle of depression of the camera is automatically calculated from the field picture and themap, and the field images are transformed to a parallel projected image. Parallel projected field
image of target area is made by transforming them to be suitable for the map, and composingit.
Key Words: image transformation, image composition, parallel projected image, growth diagnosis

１を要するという問題があった。一方，衛星画像を利用する

はじめに方法では，撮影された画像がユーザの手元に届くまで時間

稲の生育診断を行うことは農業において重要なことである。稲の生育状況を把握し，適切な時期に適切な施肥を行がかかるため，施肥時期を逃してしまうことから生育診断には適していなかった。

うことにより，高品質な米を生産することができる。また，

そこで現在，圃場をその脇からデジタルカメラで俯瞰撮収穫量の増加も見込める。この施肥を行う適切な時期とは，影した圃場写真を用いて，稲の生育診断を行う手法が提案

減数分裂期と呼ばれる出穂前約 10 ～ 15 日の生育時期であされている１）。これは，圃場写真を画像解析することにより，とても短い期間である。

り，稲穂の生育状態を診断するものである。圃場の窒素含

これまで，稲の生育状況を把握するためにいろいろな方有率分布をわかりやすく把握するために，圃場写真を真上法が用いられてきた。ハンディタイプの計測器を用いる方から見た平面画像に変換して画像解析する。現在，画像解法や，衛星画像を利用する方法などである。しかし，計測析は写真１枚ずつについて行っている。しかし，生育診断器を用いる場合，圃場全体の生育状況の把握に多大な労力を行いたい圃場が広範囲である場合，施肥を行う期間内に＊＊＊広島工業大学大学院工学系研究科情報システム科学専攻＊＊＊

広島工業大学情報学部知的情報システム学科― 291 ―

田中健太郎・山下英生

全圃場の生育診断を行えない可能性があるため，解析時間して次式が成り立つ。の短縮が要求されている。

（2.1）

そこで複数の圃場写真を合成することにより，１枚の画像にして解析を行う。これにより，解析時間を短縮することができると考える。本論文では，複数の圃場写真と地図（2.1）式のパラメータ a，b を最小二乗法により求める。

を用いて広範囲の圃場平面画像を作成する手法を提案す

歪曲収差はカメラやレンズによって異なるため，使用する

２）る 。圃場写真を平面に変換した後，地図に合わせて合成カメラ，レンズごとにパラメータを算出する必要がある。

することにより対象圃場の平面画像を作成する。平面画像

補正パラメータ a，b を最小二乗法により求めるためには，への変換には撮影カメラの俯角が必要となる。文献１では（u，v），（u'，v'）の複数の組が必要である。本研究では，レーザ距離計を用いて平面座標を取得し，変換パラメータ画像座標が正確に取得でき，補正後の座標が推定できると算出を行っているが，本研究では，平面図である地図を利いう理由から方眼紙を真上から撮影し，その画像座標を使用して算出を行う。

用した。図２⒜の撮影には焦点距離 9.2 ㎜のレンズを用い２た。図２⒜に示す 20 点（黄色点）を用いてパラメータを圃場写真の変換・合成算出し補正を行った結果を図２⒝に示す。図中の赤点は画広範囲の圃場画像を作成するために，本研究では撮影カ像中心を示す。

メラの向きを変えながら撮影した複数の圃場写真のそれぞれについて平面変換を行う。撮影された写真には歪曲収差という歪みが存在する。これは平面変換の結果に影響を与えるため，補正を行う必要がある。また，稲の生育診断には圃場の平面画像を用いるが，この変換には撮影カメラの俯角が必要である。この俯角を写真撮影時に計測するには特殊な装置が必要となる。本論文では，写真と地図を用いて俯角の算出を行う。次に，平面画像への変換を行った後，地図と合成を行う。このとき，地図に合わせて画像に拡大・縮小，回転の処理を行う。

⒜

撮影した方眼紙画像⒝

歪曲収差補正画像

以下に，歪曲収差補正，カメラの俯角算出，平面変換，地図と平面変換画像の合成の４つの処理について説明する。
2.1

歪曲収差補正カメラで撮影した写真には，通常，図１のような糸巻形や樽形の歪みが存在する（歪曲収差）。
図２

歪曲収差補正補正後画像を見ると歪曲収差が補正されていることがわかる。算出されたパラメータの値は，a＝－8.50×10－9，b＝1.00 であった。
⒜

糸巻形収差
図１⒝

樽形収差

2.2

歪曲収差の例

俯角算出

撮影した写真を平面画像へ変換するには，撮影時の俯角
この収差は，平面変換後の結果に影響を与えるので，撮が必要となる。そこで，地図と撮影画像を用いて俯角を求影画像に歪曲収差除去処理の補正を行った後，平面変換をめる。撮影画像と地図に共通して存在する点（特徴点）を行う。今，画像の座標系の原点を画像中心とし，補正前の３点選び，これらの点間の距離情報からカメラの俯角を算画像座標（u，v）と原点との距離を r，補正後の画像座標（u'，出する。以下で，図３を用いて俯角算出手順を示す。レンv'）と原点との距離を r' とすると，図１のような収差に対ダリングソフト（Softimage）を用いて作成した圃場モデ― 292 ―稲の生育診断のための圃場写真の変換と合成

ルをコンピュータ上で撮影した圃場写真が図３⒜である。
撮影時の俯角は 20 度に設定した。また，圃場モデルの平（2.2）

面図を地図の代わりとして用いた（図３⒝）。

この変換式に θ ＝ 20.01 を代入し，図４⒜に示す赤枠部分を平面変換した結果を図４⒝に示す。平面図である地図のように平面画像になっていることがわかる。

⒜

圃場写真⒝
図３

地図

圃場モデル⒜圃場写真⒝図４

以下に俯角算出のアルゴリズムを示す。

⒜の平面変換結果

平面変換

Step １：撮影画像と地図上の特徴点を３点（a，b，c）選2.4

ぶ。
Step ２：地図上における ab 間，bc 間の距離を求め，距離比 ab ／ bc を求める（これを α0 とする）。

地図と平面変換画像の合成

地図と平面画像を合成するために，2.3 節で作成した平
面画像を地図画像と大きさ，方向が一致するように拡大・Step ３：俯角は未知であるので，適当に２つの角度（θ1，縮小，回転する。平面画像に行う拡大・縮小，回転変換はθ2）を初期の俯角として定め，それぞれ 90°－ θ1，90°－（2.3）式で表される。（x''，y''）は，拡大・縮小，回転後のθ2 平面方向に撮影画像を回転変換する。それぞれの画画像座標，k は拡大・縮小率，β は回転角度，dx，dy は座像における ab ／ bc を求める。それぞれの距離比を α1，標の平行移動量を示す。

α2 とする。
Step ４：俯角と距離比 ab ／ bc で構成される２次元直交（2.3）

座標系において，（θ1，α1），（θ2，α2）の２点を通る直線の方程式を求める。これに α0 を代入し，近似解 θ3 を求（2.3）式の変換パラメータ k，βを，平面画像と地図のめる。

特徴点２点を用いて算出する。パラメータの算出手順を以Step ５：角度 θ3 のときの距離比 α3 を求め，（θ1，α1）と（θ2，下に示す。なお，算出には以下の評価関数を用いる。

α2）のどちらか一方を（θ3，α3）と置き換える。

（2.4）

Step ６：Step ４，５の処理を繰り返し，¦α0 － α3¦ が一定値以下になったときの解を最終的な俯角として決定す（mxi，myi は地図画像における特徴点の座標を表す。）

る。

Step １：パラメータの初期値 β0，k0 を適当に定める。
Step ２：β ＝ β0，k ＝ k0 のときの評価関数 E の値 E0 と，上記の手順により，図３に示す特徴点を用いて撮影カメ

β0 を 0.1 増減させたときの評価関数 E+，E －を求める。

ラの俯角を求めたところ 20.01 度と算出され，真値に近い

Step ３：β と E で定義される２次元直交座標系において，値が求まった。

３点（β －，E － ），（β0，E0），（β+，E+）を通る２次式は次式で表される。

2.3

平面変換

（2.5）

生育診断には，圃場写真を真上から見た平面画像に変換して使用する。ここでは，2.2 節で求めた俯角 θ を用いて画像の平面変換を行う。（x，y）を写真画像座標，（x'，y'）を平面変換後の画像座標とすると，平面変換は次式で表さ（2.5）式の係数 a，b，c を以下の式より求める。

a＝れる３）。ここで，P（θ） は画像を平面方向へ θ 度回転を行う４×３の透視投影行列である。

b＝― 293 ―田中健太郎・山下英生

般に使用する点の数が多いほど精度は良くなる。そこで，

c＝使用する点の数によってパラメータ算出結果がどのように

Step ４：β' ＝と，β ＝ β' のときの E ＝を求める。

変化するかを調べた。使用する点は，図中の黄色点４点から始め，赤点，青点，緑点の順に各エリアに２点ずつ増や次に，β' を β0 と置き換える。
Step ５：k についても同様に２次式近似を用いて k' を求し，算出結果を調べた（①のエリアに示す点と同様に他エリアにも点を増やす）。使用点数とパラメータ算出結果をめ k0 と置き換える。
Step ６：β，k について交互に Step ２～４を繰り返し，１表１に示す。

ステップ前の E と現ステップの E の値が変わらなくなったとき，β'，k' を最終的なパラメータとして決定する。
画像変換後，地図と合成を行う。合成は，地図の画素を平面画像の画素に置き換えることにより行う。上述の変換で地図と平面画像の大きさを一致させているので，両画像の座標は１：１の対応になっている。前処理として地図中の田んぼの部分を赤色に変更する。特徴点２点のうち１点を基準として決め，そこから上下左右に地図の画素を調べ図６ていく。画素が赤であるときには，その座標の画素を対応する平面画像の画素に置き換える。例として，図５⒜に図表１３⒝の右手前の田んぼを赤色にした地図の一部を示す。また，図５⒝に，2.3 節の平面画像を拡大・縮小，回転処理

方眼紙画像と４つのエリア

使用点数とパラメータ算出結果使用した点の数

ab

した画像と図５⒜を合成した結果を示す。パラメータ算出には図５に示す特徴点を使用し，基準点として１番の点を使用した。

４点（黄点）－ 4.18 Ｅ－ 08

1.07

12 点（赤点）－ 8.57 Ｅ－ 09

1.00

20 点（青点）－ 8.50 Ｅ－ 09

1.00

28 点（緑点）－ 8.24 Ｅ－ 09

1.00

使用点数とパラメータ a のグラフを以下に示す。

⒜

合成用地図（一部）
図５⒝平面画像合成結果

圃場モデル写真の合成

図５⒝から正しく合成できていることがわかり，地図に合わせる変換も正常に行われている。写真１枚では，図５⒝のように田んぼ一面の画像を作成することができないため，他の圃場写真を同様の手法で変換・合成することによ

図７

パラメータ a と使用点数のグラフ

り，田んぼ一面の画像を作成する。

３グラフから，特徴点の数が 10 数点以上あれば，ほぼ精写真画像の変換についての精度検証度が確保できているといえる。

パラメータ算出時には誤差が生じる。算出に用いる特徴点の数やその位置により，算出精度がどのように影響を受

3.2

特徴点の選び方による俯角算出精度けるかを検証する。使用する点の数による歪曲収差量の算

撮影カメラの俯角は，圃場写真と地図の特徴点を３点使出精度と，特徴点の選び方による俯角算出精度の検証を用し，２点間の距離比を用いて算出を行う。本節では特徴行った。

点の選び方が，算出精度にどのように影響するか検証する。
図８のように特徴点の候補を設定する。選択した特徴点の3.1

歪曲収差量の算出精度組と算出された俯角を表２に示す。なお，図８は俯角 20

本研究では歪曲収差量を最小二乗法により算出する。一度で撮影した画像である。

― 294 ―

稲の生育診断のための圃場写真の変換と合成㎜レンズのときの俯角は 6.76 度と求まった。
4.4

平面変換

前節で求めた俯角を用いて，圃場写真をそれぞれ平面画像に変換する。4.2 節で例に挙げた画像（図 10）の赤枠部分を平面変換した結果を以下に示す（図 11）。また比較用として歪曲収差補正前の平面変換結果も示す（図 12）。図
図８

11，12 から，歪んでいた田んぼのあぜ部分が真っすぐに俯角算出を行う特徴点の候補補正されており，歪曲収差補正の効果を確認することがで表２きた。

選んだ点の組み合わせと俯角算出結果

点の組み合わせ（距離比）

算出角度（度）a，c，e（ac ／ ce）

20.01

a，d，b（ad ／ db）

19.0

a，d，f（ad ／ df）

19.33

a，c，b（ac ／ cb）

19.23

平面変換時，視線方向に垂直な線分は距離の変化が緩やかであり，視線方向に平行な線分は変化が激しい。表２か⒜圃場写真Ａ⒝圃場写真Ｂ⒞圃場写真Ｃ⒟圃場写真Ｄ⒠圃場写真Ｅ⒡圃場写真Ｆら，ac と ce のように距離の変化が緩やかな線分と変化が激しい線分を距離比算出に選ぶと，算出精度が良くなることがわかる。

４

圃場写真への適用例

実際の圃場写真を用いて変換・合成を行った結果を以下に示す。
4.1

圃場写真・地図

今回使用した圃場写真と地図を図９に示す。地図に示す
田んぼの番号と写真に示す番号は対応しており，同じ田んぼであることを表す。図９⒜，⒝は焦点距離 27 ㎜のレンズ，
図９⒞，⒟，⒠，⒡は焦点距離 105 ㎜のレンズで撮影されたものである。撮影は，カメラを水平に回転させながら行っているため，図９⒜，⒝のカメラの俯角は等しく，図９⒞，⒟，⒠，⒡のカメラの俯角は等しい。図９⒢の地図の縮尺は１／ 2500 である。
4.2

歪曲収差補正各圃場写真に歪曲収差補正を行う。例として，図９⒜に歪曲収差補正を行った画像を図 10 に示す。目視では，補正の前後で違いはほとんどわからないが，4.4 節の平面変⒢換の結果に影響が出ている。
4.3

図９俯角算出図９⒝，⒡に示す特徴点から 27 ㎜レンズで撮影時と，105 ㎜レンズで撮影時のカメラの俯角をそれぞれ算出した。その結果，27 ㎜レンズのときの俯角は 19.36 度，105― 295 ―地図

圃場写真と地図

田中健太郎・山下英生

図 10

歪曲収差補正画像（図９⒜）図 13５

地図と平面画像の合成結果

おわりに

本論文では，稲の生育診断を圃場写真を用いて行うための画像変換と合成について述べた。すなわち，複数の圃場
図 11

図 10 の平面変換結果

写真に歪曲収差補正や平面変換を行い，それらを地図と合成することにより広範囲の圃場平面画像を作成する方法について述べた。使用する点の数による歪曲収差パラメータの算出精度，特徴点の選び方による俯角算出精度の検証を行い，最適な算出法を確かめた。また，実際の圃場写真へ変換・合成手法を適用し，良い合成結果が得られることを確かめた。今後の課題として，段差のある田んぼへの適用，明るさの補正などが挙げられる。

図 12

文

歪曲収差補正前画像の平面変換結果

献１）松田真典・大澤昭一・保坂幸男・金田和文・山下英生，
4.5

画像合成

「近接リモートセンシングによる水稲の生育診断‐デ

4.4 節で平面変換された各圃場写真について回転などのジタルカメラによる葉身窒素含有率の測定‐」，日本変換を行い，地図と合成を行った結果を以下に示す。図９

リ モ ー ト セ ン シ ン グ 学 会 誌，Vol.23，No.5，pp.506-⒢の青枠で囲まれた部分は図９⒜，⒝の平面画像を，赤枠515，（2003）で囲まれた部分は図９⒞，⒟，⒠の平面画像を合成したも

２）田中健太郎・山下英生，「稲の生育診断のための圃場写真の平面画像への変換」，平成 21 年度のである。

通信学会中国支部

平面に変換された圃場画像と地図の形が合致しており，

第 60 回連合大会

電気・情報
講演論文集，p.182,（2009）

良い合成結果が得られている。

３）ディジタル画像処理編集委員会，「ディジタル画像処理」，CG-ARTS 協会，pp.252-256，（2004）― 296 ―