Cプログラミングによる画像処理(1)

これまでは、ImageMagickを使用した画像処理を練習してきました。ImageMagick自体には使いこなせないくらいの多彩な機能がある一方、どうしても自分が望むような処理を実現するコマンドがない場合もあります。その場合には、自分自身でプログラムを作成するしか方法がありません。

リモートセンシング画像の処理においては、古くはfortran、現在はC言語でのプログラミングが多く、インターネット上で入手可能なプログラムソースもC言語で記述されたものをよく見かけます。私自身、C言語での経験が長いため、ここでもC言語でのプログラミングに限定して話を進めていきます。

このページの内容:

画像におけるデータ構成

手始めに、入力画像の輝度値を2倍にして出力する簡単なプログラムを作成しましょう。この演習では、汎用フォーマット(rawフォーマット)の画像を使用するという前提で話を進めていきます。

そもそも、画像ファイルにはどのような仕組みでデータが記録されているのでしょうか?画像は1つ1つの単位である画素(ピクセル)の集合体で、下記に示しますように、画像そのものは行と列から構成されているように見えます。リモートセンシングの分野では、画像上の行はライン(line)、列はピクセル(pixel)と呼ばれることもあります。「ピクセル」という表現は、画素と列のどちらを意味するか、文脈から判断するようにしましょう。

画像のライン、ピクセル

しかし実際のファイルには、ライン数、ピクセル数に関係なく、単に一列のデータが並んでいるだけです。こちらがわが画像の幅であるピクセル数を指定することで、画像処理ソフトはピクセル数ごとに折り返してくれるので、行と列から構成される画像として表示されるというわけです。

ということで、大事な点としては、

画素単位での処理の場合には、「ファイルの終わりまで処理する」という命令を指示すればいいですから、特にライン数、ピクセル数は要りません。ライン数、ピクセル数を明示し、処理しても構いませんが、必ずしも要求されるわけではない、ということです。

一方、「近傍での処理」と書きましたが、ある画素を中心とした3×3や5×5の領域において、平均や差分などを計算することもよくあります。その場合には、画像における縦と横がわかっていないといけませんので、ライン数、ピクセル数は必要不可欠になります。

このように状況によって、求められる変数(パラメータ)が異なるという点を理解しましょう。慣れれば、そんなに大したことではなくなります。

輝度値を2倍へ変換する処理

それでは実際の処理に移ります。まずは下記の画像を保存して下さい。

演習用画像データ(rawフォーマット):
[バンド3] [バンド4] (右クリックして、「保存」を選択して、各々の画像を保存してください)

次にプログラムソースの作成に移ります。テキストエディタ(「秀丸」や「メモ帳」など)を起動し、下記のプログラムソースを入力して、twice.cというファイル名で保存しましょう。左側に記載してある「1:」などの数字は、解説の必要上記述した行番号ですので、実際には入力する必要はありませんので、注意してください。

 1: #include <stdio.h>
 2: #include <stdlib.h>
 3:
 4: int main(int argc, char *argv[])
 5: {
 6:   FILE  *fpi, *fpo;
 7:   unsigned char  idat;
 8:
 9:   /* 引数のチェック */
10:   if (argc != 3) {
11:     fprintf(stderr, "Usage: %s [input] [output]\n", argv[0]);
12:     exit(1);
13:   }
14:
15:   /* 入力画像のオープン */
16:   if((fpi=fopen(argv[1], "rb")) == NULL){
17:     fprintf(stderr, "input file open error\n");
18:     exit(1);
19:   }
20:
21:   /* 出力画像のオープン */
22:   if((fpo=fopen(argv[2], "wb")) == NULL){
23:     fprintf(stderr, "output file open error\n");
24:     exit(1);
25:   }
26:
27:   /* 入力画像の読込み */
28:   while (fread(&idat, sizeof(unsigned char), 1, fpi) == 1){
29:
30:     /* 2倍の変換 */
31:     if (idat * 2 > 255) {
32:       idat = 255;
33:     } else {
34:       idat = idat * 2;
35:     } 
36:   
37:     /* 変換データの書出し */
38:     if(fwrite(&idat, sizeof(unsigned char), 1, fpo) != 1){
39:       fprintf(stderr, "data write error\n");
40:       exit(1);
41:     }
42:
43:   }
44:
45:   fclose(fpi);
46:   fclose(fpo);
47:
48:   return (0);
49: }

上記のプログラム処理の大まかな流れを説明します。

まず、入出力ファイルを開く処理ですが、読み込みの場合には「r」、書き込みの場合には「w」を指定します。画像ファイルは、テキストファイルで開いて中身を確認できないバイナリと呼ばれる形式で保存されています。したがって、画像ファイルの読み書きの際には、バイナリ形式用の「b」をつけることになっています。

次に、1画素ごとのデータの読み取りですが、バイナリ形式のデータ読み込みに 使用される関数として、ここではfreadを使用しています。 freadの書式(使い方)は下記の通りです。freadは正常に読み取りれれば、 読み取ったデータの数を返してくれます。

fread(読む込んだデータの格納先, データの大きさ, データの数, 入力ファイルのポインタ)

データを読み取った後ですが、単純に2倍して出力しているわけではありません。 入力画像も出力画像も1バイトという前提ですが、1バイト とは8ビットのことで、28 = 256の階調を持っています。実際には、 各画素は0から255までの値を取り得ます。2倍した結果、 255より大きくなった場合には、255という値を与える処理を、 if文を使って実現しています。if文については、 説明のページ演習問題のページが用意されていますので、 興味がある方は参照してください。

最後になりますが、freadという関数は、データの読み取り に成功すれば、読み取ったデータの数を返してくれます。 ピクセルごとに1つずつデータを読み取りますので、1を返してくれるはずです。 その返してくれる値を元に読み取りを継続するかどうか、while文の中に記述し ています。全てのピクセルでの処理が終わり、ファイルの終わりに達すると、 while文での繰り返し処理は終わります。

それでは、上記の「twice.c」というプログラムを入力できたら、次はコンパイルする必要があります。コンピュータが実行可能な形式のファイルへ変換する作業で、下記のように行います。コンパイルを含むC言語の基本的な処理に関しては、こちらのページに説明がありますので、参照してください。

コンパイル

(方法1)$ gcc twice.c
(方法2)$ gcc -o twice twice.c

実行

$ (実行ファイル) (入力画像ファイル) (出力画像ファイル)
(方法1でコンパイルした場合)
$ ./a.out tm931028b3.raw tm931028b3_twice.raw
(方法2でコンパイルした場合)
$ ./twice tm931028b3.raw tm931028b3_twice.raw

JPEG画像への変換
画像処理ソフトであるImageMagickのうちの、convertというコマンドを 使用。convertの他の使い方は
こちらのページへ。

$ convert   -depth 8   -size 2000x1500   gray:tm931028b3_twice.raw   tm931028b3_twice.jpg
(演習問題)上記のプログラムでは、255の値を超えるかどうかで条件分岐を行っていたが、条件分岐をしない(つまりif文をなくす)プログラムをtwice2.cとして作成し、出力画像(JPEG画像)を確認せよ。2倍するくらいでは差がわからない場合もあり、その場合には3倍、4倍など倍率を変えて試してみよ。

バンド3、4データからのNDVIの生成

上記では1枚の画像を読み込み処理しましたが、今度は複数枚の画像を使った 処理を練習してみます。Landsat TMのバンド3は可視光線の赤の波長帯 (630 nm - 690 nm)、バンド4は近赤外線の波長帯(750 nm - 900 nm)に 設定されています。植生(植物の集合体)の 電磁波の波長帯ごとの反射率スペクトルカーブ(分光曲線)と呼ばれますが、 スペクトルカーブは植生の種類、また同一の植生でも季節によって変動すると 知られています。

植生においては赤から近赤外におけるカーブの上昇具合が、植生の光合成の状態(光合成が盛んかどうか)と相関が高いと言われています。ですから、単にバンド3、バンド4のデータを各々単独で用いるよりは、スペクトルカーブのおおよその傾きに相当するバンド4とバンド3の差が使われたり、さらにバンド4とバンド3の差を、バンド4とバンド3の和で割り標準化されたNDVI (normalized difference vegetation index: 正規化植生指標)と呼ばれる指標が、植生分布の解析によく使用されます。

NDVI = (band4 - band3) / (band4 + band3)

このNDVIが取り得る値の範囲は、-1(band4=0のとき)から1(band3=0のとき)までと言われています。

これから演習に入りますが、先の演習で作成したtwice.cを基本にして、NDVIを計算して出力するndvi.cを作成していきます。twice.cをコピーしてもいいですし、別名保存ででも結構ですので、先にndvi.cを作成してください。その後、NDVIの計算を実現するように、プログラムソースを書き直していきます。

もう一度、ここでの状況を整理しますと、

ということです。ここには正解のソースを書きませんので、自分自身で考えて、要求を満足するndvi.cを作成してください。

作成し終えたら、コンパイルが成功するか確認してください。C言語の文法上の誤りがあれば、この段階でエラーメッセージが出ます。また、コンパイルに成功しても処理が正しかったかどうか、画像を確認してください。上記の演習と同様、rawフォーマット画像をJPEGフォーマットへ変換した後に、画像を開いて確認してください。画像輝度値の範囲が狭いと暗く見えたりしますが、その場合には下記の「画像強調」を試してください。

画像強調
onvertのオプションとして、-equalizeを使用。

$ convert   -equalize   tm931028b3_ndvi.jpg   tm931028b3_ndvi_eql.jpg

どうしてもわからない場合には、こちらのプログラム ソースを参照してください。NDVI画像作成用のプログラムソースの例が 記載されています。


Cygwin、ImageMagickによる画像処理のページに戻る
須崎純一トップページに戻る
須崎純一 京都大学大学 工学研究科都市環境工学専攻 環境情報学講座