(農業)情報工学・画像3　(2019/12/19)

今回の講義では、まず、カラー画像の基本的な特徴について理解する。
次に、Image-Jのマクロ処理を利用する。
その後、リモートセンシングの入門演習を行う。
(なお、リモートセンシングの作業では、作業スペースとして500MB〜1GB程度の余裕が好ましい。)
自分の持っているディスクスペースを確認して、不要なファイルを削除しておく必要がある。

画像処理入門(3)

カラー画像の表現方法
これまでは、主に白黒の画像（グレイスケール画像）を取り扱ってきた。今回はカラー画像を取り扱う。
色情報はスペクトル分布を持っているが、任意の色を再現する場合には互いに独立な３つの色があれば良いことが知られて（グラスマンの法則）いる。その中で、コンピュータで広く利用されているカラー表現は、ＲＧＢ表色系である。これは、赤（Ｒ）緑（Ｇ）青（Ｂ）の３つの色を光の３原色として混合して色空間を構成するもので、各色の階調（色強度）を８ビット＝２５６階調とすることによって、２４ビット＝16,777,216色の色が再現できる。
この２４ビットカラーのことを、通常、フルカラーと称している。人間の識別可能な色数は、７〜８００万色程度といわれているため、通常は２４ビットの色再現で十分である。（ただし、一部、金色などの色はコンピュータディスプレイでは再現が難しい。）
一方、現在ではこの２４ビットにソフトウェアで制御可能な情報８ビットを加えた３２ビットで「トゥルーカラー」と称する場合（Windows系）や、各色10ビットずつの約10億色を再現させて「トゥルーカラー」や「リアルカラー」などとして区別する例などがあり、表現が混乱している場合もある。
ＲＧＢ以外には、Cyan,Mazenta,Yelowの３色（実際の印刷などでは黒(K)も利用した４色）を用いるＣＭＹ表色系や、色相(H) 明度(V)彩度(C)で色を表現するマンセル表色系（HV/C）などが使われているが、本講義では取り扱わない。
今回、最初に取り扱うカラー画像は、夏ミカンのカラー画像である。
フルカラー画像とRGB成分
Summer Orange01画像1 は、夏ミカンの画像である。この画像を元に、ミカンの個数を数える処理を行ってみよう。
この画像は、単純に白黒化した場合、葉の白っぽい部分とミカンの黄色の部分との明るさの差が小さく、うまくミカンの画像のみを抽出することが出来ない。（画像の2値化を行って確認してみる事。(Process->Binary->Make Binary)）
そのため、果実部分のみを抽出するために、色情報を利用する。ここでは黄色の画像を強調するために、R画像からB画像成分を引くことが有効である。（単色情報だけでは、うまくミカン画像が抽出できないことを確認してみよう。）
続いて、ミカン画像を数える手順を以下に示す。
  1. ミカンの画像を適当な場所に保存する。
  2. ImageJで画像データを読み込む
  3. RGB成分に分解する。(Image->Color->Spilit Channels)
  4. 画像計算機を呼び出し、(Process->Image Calculater)
    黄色成分を計算する。即ち、[R画像]から[B画像]を引く（[R画像]　−　[B画像]）。
    (Image1->S01.jpg (red), Operation->Subtract, Image2->S01.jpg (blue))
    この際に、create new window, 32-bit を選択しておく。
    結果を確認する。(Result of S01.jpg)
  5. 8Bitイメージに変換する。(Image->Type>8 bit)
  6. 画像を2値化する。(Process->Binary->Make Binary)
  7. ミカンの個数を数える。
    　(Analyze->Analyze Particles)
    　Size (pixel^2): 3000 (ミカンと見なす最小のピクセル数。これでノイズを除去できる。)
    　Circularity: 0.00-1.00 (真円度。デフォルトで良い) Show:Elcipsesを指定。
    　Display Results,Clear Results, Summarize, Include Holesをチェックして「OK」する。
  8. 結果(Summary)ウィンドウは、最後の処理画面の下に隠れてるので、クリックして表示すること。ここで、ミカンの数(3個)が、[Count]として、表示されていれば成功である。適切に抽出されているかをチェックする。
上記のステップ２からステップ７を、[ミカンcount手順]と定義しておく。
ImageJのマクロ処理
上の手順では２．から８．の7段階の作業を行って、ミカンの個数を数えることができた。
この作業を画像2に適用する場面を考えてみよう。
ImageJでは、同じ処理を続けて実行する場合、マクロ処理という機能がある。これは一連のコマンドを実行するためのプログラムの一種であるが、最も簡便なマクロの作成手法としては、操作手順の記録である。ここでは、まず、その手順を行う。
  1. ミカンの画像を作業用の場所に保管する。
  2. ImageJで、マクロの記録モードにする。
    　Plugins->Macros->Records
    　マクロ記録用ウィンドウが表示される。
  3. 一連の手順（[ミカンcount手順]）を実行する。
  4. マクロを作成・保存する。(Createのボタンを押す。デフォルトではMacro.jimが作成される。)
  5. マクロには、適当な名前（例：O-count.jim）をつけて、保存すると良い。File->Save
この手順で、おそらく、次のようなマクロが作成されているはずである。
  1. open("/Users/nori/Desktop/Image03/S02.jpg");
  2. run("Split Channels");
  3. imageCalculator("Subtract create 32-bit", "S02.jpg (red)","S02.jpg (blue)");
  4. selectWindow("Result of S02.jpg (red)");
  5. run("8-bit");
  6. setOption("BlackBackground", false);
  7. run("Make Binary");
  8. run("Analyze Particles...", "size=3000 show=Ellipses display clear include summarize");
  9. selectWindow("Summary");
  1. ファイルのオープン(※)
  2. ＲＧＢ成分に分解
  3. イメージ画像の演算(red - blue)
  4. 結果画面の選択
  5. 8ビット化
  6. 8ビット化のオプション
  7. ２値化
  8. 粒子個数のカウント
  9. Summaryウィンドウの表示
(※)　ファイルを保存した場所の違いによって、ディレクトリ名（１．のOpenのディレクトリ名 /Users/nori/Desktop/image/S02.jpg　の部分）が異なっている。
このマクロ（プログラム）は、(1)ファイルのオープン,(2)ＲＧＢ成分に分解,(3)イメージ画像の演算,(4)結果画面の選択,(5)8ビット化,(6)8ビット化のオプション,(7)２値化,(8)粒子個数のカウント,(9)結果の表示を順に実行するものである。
保存したマクロを実行するには、

ここまでで、ミカンの個数カウントマクロが完成した。しかし、このままでは、いつも同じファイル（S02.jpg）に対してだけ＜個数カウント＞を実行するマクロとなってしまい、ほとんど実用性は無い。
そこで、このマクロに修正を加えて、 としたのが、次のマクロプログラムである。(O-count2.jim)
  1. // - Macro for Count the Number of Oranges -
  2. dir="/Users/nori/Desktop/Image03/";
  3. file=getString("File Name->", "S01.jpg");
  4. min=getNumber("Minimum Particle Size->", 1000);
  5. //
  6. filename=dir+file;
  7. Img1=file+" (red)";
  8. Img2=file+" (blue)";
  9. Win1="Result of "+file+" (red)"
  10. Cmd1="size="+min+" circularity=0.00-1.00 show=Ellipses display clear include summarize"
  11. //
  12. open(filename);
  13. run("Split Channels");
  14. imageCalculator("Subtract create 32-bit", Img1, Img2);
  15. selectWindow(Win1);
  16. run("8-bit");
  17. setOption("BlackBackground", false);
  18. run("Make Binary");
  19. run("Analyze Particles...", Cmd1);
  20. selectWindow("Summary");
  1. コメント行
  2. 変数"dir"の指定(自分の環境に合わせて変更)
  3. ファイル名の入力(デフォルトはS01.jpg)
  4. 最小画素数を入力(デフォルトは1000)
  5. //
  6. ファイル名をフルパス名で指定。(2+3を合成)
  7. イメージ演算のファイル１((red) の前に半角スペース)
  8. イメージ演算のファイル２((blue)の前に〃)
  9. 演算結果ウィンドウの文字列
  10. ミカン数コマンドを文字列(Cmd1)として合成
  11. //
  12. ミカン画像ファイルオープン
  13. ＲＧＢ成分に分解
  14. イメージ画像の演算（Img1 - Img2）
  15. 結果画面の選択
  16. 8ビット化
  17. 8ビット化のオプション
  18. ２値化
  19. 粒子個数のカウント
  20. Summaryウィンドウの表示
マクロの注意点。
  • ６行目でファイル名をフルパス名で指定。(2行目で指定したディレクトリ名と3行目でキー入力したファイル名を合成して、ファイル名としている。（FQDNみたいですね(^_^)）
  • 7,8行目は、image calculater用の入力ファイル名１と２を指定している（ここでも、文字列の合成である。）
    ここで、(red) (blue)の前に、半角スペースが空いていることに注意
  • 9行目は、演算結果ウィンドウを選択するための文字列
  • 10行目は、ミカン数カウントコマンドに指定するコマンドを文字列(Cmd1)として合成している。
    • ここで、このようにコマンドを作成するという手順を、理解してほしい。
    すなわち、最小粒子サイズ変数minが、[ "size="+min ]という文字列として埋め込まれている。
    実際に、ImageJの実行を行っているのは、12行目以降である。
  • 12〜20行目が、O-count.ijmでの１〜9行目に相当する。 このプログラムを作成し、O-count2.ijmとして保存して、実行してみる。
    上の例（画像１）でミカンが３個、画像2でミカンが4個と表示される事を確認したら画像3を使って、ミカンの数を数えてみよう。
    最小の画素数をいくつにすれば良いであろうか。
    他にも、黄色いミカンの木の画像などをネットで探して、実行してみるのも良い。
    ・ミカンの木画像の検索←Googleの画像検索
    リンゴの画像の場合は、どのような工夫をすれば良いであろうか。
    ・リンゴの赤い果実は、G-Y画像で抽出できる。なお、Gは緑画像。Yは輝度画像であるが、通常のカラー画像を単純に8ビット化することによって、近似画像が得られる。簡略化した手順は、次の通り。
    1. 元画像を色分解して、G画像を準備する。
    2. 元画像を再度呼び出し、8ビット化する。(Y画像）
    3. G-Yを演算し、ミカン画像と同じ手順で2値化、個数カウントする。
    リモートセンシング画像と擬似カラー表現
    リモートセンシング(Remote Sensing)とは遠隔から測定するという意味で、隔測と呼ばれることもある。地表の状況を大規模に測定する事を目的として、人工衛星から得られた画像を利用する。
    主な人工衛星の名称とその特徴については、次の通りである。
    衛星名(国・呼称) 主なセンサと解像度 衛星高度
    ALOS-2(エイロス)　(日：だいち2号)
    2014年    		 PALSAR−2　(1-3m) 628km
    14日周期
    ALOS(エイロス)　(日：だいち)
    2006年    		 AVNIR-2　(10m) 692km
    46日周期
    PRISM　(2.5m)
    PALSAR　(10m)
    LANDSAT4,5　(米：ランドサット)
    LANDSAT7　(米:ランドサット)
    		 MSS (80m) 705km
    TM (30m)
    SPOT　(仏：スポット） HRV-XS　(20m) 832km
    HRV-P　(10m)
    MOS　(日：モス) MESSER　(50m) 909km
    JERS　(日：ジェイ・エルス) OPS　(18) 568km
    SAR
    ADEOS　(日：アデオス) AVNIR　(8m,16m) 800km
    OCT　(LACデータ：700m)
    OCT　(GACデータ：4km)
    EOS Terra　(米：イオス・テラ) ASTER　(15m,30m,90m) 705km
    CERES　(25km)
    他
    NOAA　(米：ノア) AVHRR　(1.1km) 850km
    IKONOS-1　(米民間：イコノス) 白黒(0.82m) 680km
    マルチスペクトル(3.3m)
    Quick bird-2　(米民間：クイックバード) 白黒(0.61m) 450km
    マルチスペクトル(2.44m)
    RADARSAT-1,2　カナダ:レーダーサット) SAR-10　(25m×28m)
    (SyntheticApertureRadar)    		 793km〜821km
    ここでは、NOAA(National Oceanic and Atmospheric Administration)のデータを元にして植生指標を求める計算を行う。
    具体的には、この手順に従って正規化植生指標（NDVI:Normalized Difference Vegetation Index）を作成する。
    植生指標とは、簡易な計算式で植生の状況を把握することを目的として考案されたもので、植物の量や活力を表す。植物は通常緑色に見えている。すなわち太陽から届く光の中で、緑色の光を特に反射しやすく（反射率が高い）、逆に赤色の光は反射しにくい（反射率が低い）という特徴があるためである。この結果、人間の目には緑色の光がより多く届き、植物が緑色に見えている。
    また、植物体は、一部の赤外線（近赤外域(800nm〜1300nm)）で高い反射率を示すという特徴もある。この特徴を元に計算されたものが植生指標である。
    ここで用いるNOAAのデータは、改良型高分解能放射計(AVHRR:Advanced Very High Resolution Radiometer)のチャンネル１(Ch1.)とチャンネル２(Ch2.)画像である。具体的な計算式を下に示す。

    (1)植生指標(r_ndvi)の計算
    r_ndvi = (ch.2 - ch.1)/(ch.2 + ch.1)
    ch.1: AVHRRのチャンネル１(可視域(赤色))
    ch.2: AVHRRのチャンネル２(近赤外域)
    r_ndvi: -1〜+1の実数による植生指標
    ※ここでr_ndviは実数値であるため、画像データとして利用しやすくしたのが下記のNDVI値であるする。
    (2)正規化植生指標(NDVI)
    　NDVI = (r_ndvi + 1.0) * 100
    　NDVI: 0〜200の整数(8bitに収まる)にストレッチされた植生指標

    具体的には、上に見るように近赤外域と可視域との差を元にした計算式で、数値が大きいほど植生が活発な点であるとみなしている。
    今回は、この指標を計算してみる。
    1997年7月のデータ:CH1, CH2
    2004年7月のデータ:CH1, CH2
    データは圧縮してあるが、解凍すると、フォーマットの無いRAW(生)データが得られる。一つのファイルサイズが12メガとなり、かなり大きいので注意が必要である。画像形式は、3600x3600ドットで、各点8ビット(256階調)になっている。
    これを、ImageJで読み込み、上の演算を行うと、NDVIデータが得られる。
    得られたNDVIの結果は、上記に示すように0〜２００の整数化された数値となっている。この数値は標準ではグレースケールであるが、白黒の濃淡だけでは細かな数値の違いが見にくいため、数値ごとに適当な色をあてはめてディスプレイに表示させることがある。天気予報の際に、気温が高いほど赤く低いほど青く色を付ける作業と同じ原理である。これを擬似カラー（Pnuedo Color)と呼ぶ。どの数値に何色を割り当てるかの対応表の事を、ルックアップテーブル(LUT)と呼んでいる。ImageJでは、標準でいくつかのLUTが存在するが、残念ながら植生を表す場合にはあまり適当な物が無い。今回は、赤から連続的に緑色に移行するRed/GreenのNDVI2.lutを適用した。ただし、単純にRed/GreenのLUTを適用すると、数値が０となっている水域が真っ赤に表現されるため、数値０の場合には白を割り当てている。

    処理手順は次の通りである。
    1. リモートセンシングデータ(上記のCH1とCH2)をダウンロードし、適当なディレクトリに保存する。その後、解凍しておく。
    2. 同様に、ルックアップテーブルデータ(NDVI2.lut)も保存しておく。
    3. ImageJでCH1のデータを読み込む。
      File->Import->Raw->(ファイルの指定)->8bits,3600,3600と指定。
    4. 同様にCH2のデータを読み込む。
      File->Import->Raw->(ファイルの指定)->8bits,3600,3600と指定。
    5. CH1とCH2の数値の合計を計算する。その際に、結果を３２ビット値とするのを忘れないようにする。（数秒の時間がかかる）
      Process->Image Calculater->add(CH1とCH2を指定)
    6. 結果を、add.rawとして一時的に保存する。このファイルは４９メガあるので不要になったら、なるべく削除すること。
      File->Save as->Raw Data..
    7. 同様に、CH2 - CH1 を計算し、結果をsub.rawとして保存する。
    8. 元のCH1とCH2の画像を、一端、閉じる。（メモリー消費対策のため）
    9. sub.raw÷add.rawを実行する。 Process->Image Calculater->devide(subとaddを指定)
    10. 結果のデータ(r_ndvi値：-1.0〜1.0）に1を加え、その後、１００倍する。 Process->Math->add(1)
      Process->Math->Multiply(100)
    11. 結果を、８ビット化する。
      Image->Type->8bits
    12. ＬＵＴを適用する。
      File->Import->LUT(NDVI2.lutを指定)
    13. このような画像が得られる筈である。
    少し複雑な手順であると思う。しかし、内容を理解しながら、確実に行って欲しい。
    この作業で利用可能なNOAAのデータは、先ほども述べた植生指標データのダウンロードリンクから入手することが出来る。(今回は1kmメッシュとする)ここのCH1とCH2のデータを元に、NDVIデータを作成することができる。
    （2019年12月現在、CH1とCH2のデータはダウンロード不能で、計算済みのNDVIのデータのみがダウンロードできる。）
    【参考】ここまでに紹介した手法を用いて、LANDSATのデータに対してNDVIを求めてみる。 手順は、次のリンクに従った。
    Landsat 8画像でNDVI算出
    NOAAのデータは1ｋｍメッシュ（現在は250ｍメッシュ）サイズであるが、LANDSATのデータは30ｍメッシュである。LANDSATのデータは、産業総合研究所で公開されている。
    NDVIを計算するためには、Landsat8のBand4(赤:R)とBand5(近赤外:IR)のデータが必要である。また、雲の少ない画像を入手することが好ましい。

    以下に、入手方法を示す。
    サンプル：Band4, Band5
    入手した二つの画像(Band5とBand4)を用いて、NDVIを求める
    但しLandsat8のデータは16ビットで符号化されているので、そのままではImageJで開けることが出来ない。そのため、一旦、PreViewでファイルを開けてから、8ビット化して保存する。（ファイルフォーマットはTIFFで良い。Jpegにはしないこと。Jpeg形式はデータ圧縮の際にデータの非可逆的変換を行うので、元の数値が変化してしまう。） NDVI=(band5-Band4)/(Band5+Band4)　
    LookUpTableはNDVI3.lut(0と255を白にしている)を用いると良い。 なお上の式で得られるNDVI値は、32ビット化したときに実数値となるが、今回は、そのまま8ビット化する。今回のデータは海洋のデータも緑色になるので、数値を見て、海洋には色が付かないように調整しても良い。(Prosess-Mathで調整する。)
    今回の課題
    以上を、ITC-LMSのレポート機能を用いて提出する。健闘を祈る。
    リンク
    色々雑学（コニカミノルタ）：色に関する説明のページ。この中でも「表色系」の説明は秀逸である。
    国土環境モニタリング：ここには、今回の講義内容に相当するデータが保管されている。