複天一流:どんな手を使ってでも問題解決を図るブログ

宮本武蔵の五輪書の教えに従い、どんな手を使ってでも問題解決を図るブログです(特に、科学、数学、工学の問題についてですが)

画像処理にチャレンジ(part 1):彗星観測の失敗と無謀な挑戦

Pons-Brook彗星の観測失敗

今年は大きめの彗星が2つほどやってくると噂に聞いた。実は一つ目はすでに観測可能であり、ネットを検索すると綺麗な緑青色のコマをもった彗星の画像が溢れている。

www.theguardian.com

(Arundel城の頭上に光るこの写真はおそらく合成された偽写真だと思うが)こんな感じで見えたら最高である!ということで、遅まきながら私も観測に参加することにした。ちなみに、この彗星は19世紀の初めに発見されたが、最初がフランス、その70年後が英国の天文学者による観測だったという。つまり70年ほどの周期で太陽の周りを周回する楕円軌道の彗星である。

彗星の観測といえば、吉田氏のHPである。光度グラフをみると5月1日あたりが最大光度のようだが、関東周辺の場合、現在は夕方になると地平線の上方にわずかに見える位置にある。おそらく4月に入ると夜間は地平線に沈み、昼間の空に彗星の位置が移動するだろうから、(北半球からは)観測できなくなるであろう。ということで、「今しかない」のである。

さっそく観測してみたが、現在5等星(最大光度の場合でも4.5等星)なので肉眼でみるのはきついが、観測写真なら十分捕捉可能な明るさである。しかし、現在の彗星存在地点が馴染みの薄い場所であり探しにくいという問題点が(私の場合は)ある。

ということで、広角レンズで、少し広めに夜空を切り取ってみた。 今、もっとも目立つのが木星である(写真の中央からやや下めに位置する明るい点)。その少し上にプレアデス星団(すばる)が見える。すばるの左斜め上にあるV字型の星座が牡牛座でアルデバランの赤い点が比較的よく目立つ。昨年は、アルデバランのすぐ横にZTF2023E3彗星が位置した時期があり、撮影に成功した。アルデバランの赤(朱色)と彗星の緑青のコントラストが見事であった。そのさらに右上にずーっと視線を動かすとオリオン座が右下を覗き込むように傾いて光っている。

しかし、今回のPons-Brooks彗星は「逆側」の領域にある。逆という意味は、木星とスバルをむすぶ垂直線に対して右側という意味である。右側に目を写し、上方を仰ぎ見ると比較的明るめの恒星カペラがみえる(ぎょしゃ座)。写真の中央右上にぎょしゃ座の5角形は丸ごと写っている。カペラの下にあるのがペルセウス座で、星がぐちゃぐちゃっと写っている領域である。その手にメデューサの首ならぬアルゴル変光星が握られている。カペラとアルゴルを結ぶ直線の先に2つの星が縦に並んでいるのが見える。これが牡羊座のハマル(上)とシェラタン(下)である。写真中央部の地平線近くにのっぽの針葉樹が見えるが、そのすぐ左上にある、縦に並んだ2つの恒星である。この星の右上に、鋭角三角形が下向きに位置している。三角座である。彗星は、牡羊座と三角座の中間に位置しているはずである。...が、残念ながらこの写真には写っていない。

拡大写真

彗星は動き回る天体なので、毎回毎回その位置ぎめが大変な作業である。星座によく習熟しておかないと、さっぱりの場合もある(今回がそれである)。いつも思うのが、星のパターンを自動認識してくれるソフトウェアがないかな、という願望である。それよりなにより、この写真に写っている星はいったい幾つあるのだろうか?明るい星はどのくらいで、暗い星は何個あるだろうか?そういう分析こそ、画像処理が得意ではないのか?

....と思ったのがきっかけである。デジタル画像処理について興味が湧いたので、さっそく自分でできる範囲までやってみようと思い立ったのである(無謀?)

画像処理の教科書を納屋から引っ張り出す

遥か以前、画像処理をやってみようと思い立った時期があっった。あの頃は、変光星に興味があり、定期的に観測して、その光度がどのように時間変化するのか興味があったのだ。フリーウェアがいくつか見つかったのだが、そのほとんどがWindows向けのものであり、LinuxからmacOSに軸足を移しつつあった私は、その試みを諦めたのである。吉田氏の彗星光度分析にあるような画像解析をやりたかったのだった。

そこで、自分でなんとかならないかと考え、いろいろな書籍を購入してみたのだが、そのほとんどが(やはり)windows向けの書籍であり、内容を読んでもさっぱりわからず、再び挫折してしまったのである。

この時購入したいくつかの書籍が埃をかぶって納屋に眠っているのを思い出し、もしかすると今なら多少はわかるようになっているかも、と根拠のない楽観に基づき、再再度挑戦することにしたのである。

doomとnanosaursの経験が生きているかも

この間チャレンジしたdoomnanosaursというゲームのコンパイルをやっている最中に、アプリケーションのコードがどのように書かれているのか、多少見聞きする機会があった。特にlinuxxdoomでは、C言語のプログラムのバグを手で修正したので、その構造や書き方について(少しではあるが)その慣習を知ることができた。

昔買った教科書のひとつが「画像処理プログラミング」というものである。以前はこのプログラムの意味がまったくわからず、コンパイルも通らず途方に暮れたのだが、今見ると、どうしてうまくいかなかったのかすぐにわかった。関数のコードリストはそのままコンパイルできないこともあるし、includeファイルの設定や、main()の記述などもどこかでやらねばならないのである。

doomの場合には、mainはほとんど形式的な利用になっていて、コードを見ると、すぐにD_DoomMain()を呼び出して実行的な処理はそちらでやるようになっていた(main関数はi_main.cにある)。

int main( int argc,  char **argv ) 
{ 
    myargc = argc; 
    myargv = argv; 
 
    D_DoomMain (); 

    return 0;
} 

さて、教科書に掲載されたプログラムであるが、まずはmain関数が載っているプログラムを探さないといけない。ずっとページを追いかけていくと、第1章の終わりに近いところにあった。ただ、ヘッダーファイルのリストが教科書のどこにも載っていないことがわかった。これではコンパイルが通らないわけである。

よくよく調べると、この教科書のコードはすべてソフトバンクのサーバーで公開されているそうで、まずはそちらをダウンロードしてみることにした。すると、そこにヘッダーファイルが含まれていたのであった。しかし、これだけではまだ不十分である。ダウンロードパッケージにはmakefileが含まれていなかったので、自分で書く必要があった。ファイルの依存関係を見出すのがなかなか骨折りであったが、mainを含むコードを読んで、必要な関数が含まれているコードを探し、それを個別にオブジェクトファイルにコンパイルできるかどうかチェックし、最後にリンクしてみたら、うまくいった!

できあがったコードを走らせると、小さなBMPファイルが生成されていた。教科書で解説していた市松模様の画像であった!

灰色の四角形に見えるが、拡大すると市松模様になっているBMP画像ファイル

これで、やっとこの教科書の読み方がわかった。画像処理の技術はひとつひとつが関数として分割して書かれており、見栄えが悪いヘッダーファイルは教科書では省略されてしまっている。コードの書き方だけを学んでね、という書き方である。

実際にコンパイルしてコンピュータで動かすには、ヘッダーファイルを書き足し、必要なメイン関数を書いて、関数ファイルとリンクさせる必要がある。主な作業はダウンロードパッケージに含まれたファイルを借用すれば省けるが、ファイルの依存性だけは自分でやらないとだめである。コードの中身を読み、教科書を読みながら、自分で分析してmakefileを書く。これらを突破し最後までは辿り着けば、実行ファイルが生成されて、BMP画像が作られるのである。

教科書をまねる

どっちみち自分でコードの構造を理解し、分析しないと動かないのであるから、思い切って同じ内容のプログラムを自分で書き直してみることにした。

まずは、画像データの先頭に組み込まれているデータ情報の読み込みをやってみることにした。画像ファイルをプログラムで読み込めなければ、なにも始まらない。

この教科書はC言語で画像処理を実装している。画像形式としては圧縮しないタイプのBMPが選ばれている。jpegとかpngとかをやりたいところであるが、これら圧縮された画像データの処理はかなり複雑であり、ライブラリを介してやらないと大変骨が折れる(ライブラリを使っても大変だと思うが)。画像処理プログラミングの最初の一歩としては(現在はあまり利用されなくなったとはいえ)BMPが適しているのかもしれない。BMP利用が減少した最大の理由がその大きなファイルサイズであったが、近年の技術革新により、BMPサイズの大きさがあまり気にならなくなってきたというのもある。Windows3.1やOS/2が利用されていた「大昔」においては4MBのフロッピーディスクに画像を何枚しまえるかどうかが問題だったが、今では4TBの外付けHDDや256GBのSDカードが1万円から2万円程度で買えるようになった。もはや画像はBMPのままでも問題ないのかもしれない(ネット利用でなければjpegは必要ないかも)。

BMPファイルを用意する

まずはBMPファイルを用意する必要がある。教科書のプログラムを使って生成した(上の)市松模様のファイルがその一つである。もう一つは、macOS Montereyのバージョンを示す画像をpngスクリーンショットした後、imagemagickをつかってBMPに変換したものである。

まずは、imagemagickをつかって、用意した画像のプロパティを確認しておこう。

% identify out.bmp

out.bmp BMP3 100x100 100x100+0+0 8-bit sRGB 30054B 0.010u 0:00.000

いろいろな情報が表示されたが、わかるのは「形式がBMP3、大きさが100x100ピクセル、ファイルサイズは30054byte」ということであろう。

pngからbmpに変換したものについても同じようにやってみると、

%identify monterey.bmp

monterey.bmp BMP 654x400 654x400+0+0 8-bit sRGB 789826B 0.010u 0:00.003

となった。ファイル形式は(BMP3ではなく)BMPとある。いろいろなバージョンがBMPにもあるのであろう。「大きさは654x400ピクセル、そしてファイルサイズが789826byte」であった。

その他の情報については、BMPについて詳しくなれば次第にわかるようになるであろう。

ファイルを読み込むコード

BMPはバイナリファイルであるから、まずはC言語でバイナリファイルを開いてから読み込むことになる。

K&Rを久しぶりに読み直してみると、「ファイルはポインタで定義される」とある。バッファの位置や、その現在の文字位置、加えて読み書きどちらなのかといった情報、ファイルの終端(EOF)が来たのか来てないのか、エラーのあるなしなど、ファイルに関する様々なプロパティが構造体の形でまとめられているのだという。そしてその情報の先頭が置かれているアドレスがファイルポインタである。したがって、まずはファイルポインタを宣言し、そのインスタンスを生成する。ちなみに、このファイル構造体の定義は、stdio.hに書かれているという(見たことはないが)。

ということで、C言語でファイルの読み出しを行う第一歩は次のようになる。

#include <stdio.h>

int main(){
    FILE *fp;
    return 0;
}

この状態ですでにコンパイルでき、実行できる(なにも起きないが)。

% gcc openFile.c

次に、BMPデータファイルの情報を、このファイルポインタ構造体に結びつける。このために必要なのが、stdioライブラリに収められているfopenである。ファイル名とファイルアクセスモード(読み込みか書き込みかなど)を指定して実行する。

FILE *fp;
fp = fopen("out.bmp","rb");

rは読み込み(read)、bはバイナリ(binary)形式、を意味する。BMPはバイナリファイルなのでbが必要である。

「ファイルオープン」つまりファイルとファイルポインタの結びつきが失敗すると、ポインタfpにはNULLが結びつく。それ以外の場合は成功と見做される。K&Rでは、ファイルオープンの成功失敗の検査も込みした例文が書き下されているので、それを採用してみた。

#include <stdio.h>

int main(){
    FILE *fp;
    if( (fp = fopen("out.dat","rb")) == NULL){
        printf("File not found. Failsure.\n");
        return 1;
    }else{
        printf("File found. Success.\n");
   }
    return 0;
}

最後にファイルの中身にアクセスしてみる。といっても最初2バイトのみとする。これは、BMPデータファイルの最初の2バイトには、識別のための"B"と"M"が必ず書き込まれているからである(そういう仕様になっていて「マジックナンバー(魔法数)」と呼ばれる)。

ASCII文字コード

C言語は米国で開発されたプログラミング言語であるため、扱う文字は米国系のアルファベットと関連する記号である。ドイツ語だったり、北欧系、東欧系の諸国で利用されている文字や記号は採用されてない。したがって、128通りの区別があれば、米国系の英文に限ってならコンピュータで表現できる(らしい)。ということで、文字を表す表現としては27=128の情報をもつASCIIコードが採用された。たとえば、Aは整数65、Bは整数66....によって識別する。ちなみに、「1」は49, 「2」は50、...といった具合に「数字」についても「文字」として認識され、ASCIIコードが割り振られている。

BMPの最初の2バイトに書き込まれているBMという2つの文字は、ASCIIでは66と77に相当する(16進数では0x42と0x4dである)。

現代のコンピューターアーキテクチャーでは、1byte=8bitであるから、28=256文字を処理することができる。ASCIIコードの最大値は127=27-1、つまり7ビットであるから、大幅に「余り」が出てしまう。この余の部分に各国の言語特有の記号を埋め込んだ新しいASCIIコードというものが制定されてはいるようだが、今考えるべきはオリジナルのASCIIコード7ビット分、つまり0.5バイトの情報量である。しかし、0.5バイトというのは中途半端な数字である。ファイルに文字を書き込む時、実際はどのようにやっているのか疑問に思うのである。

putc()関数

ファイルポインタfpで定義されたファイルの中へ一文字つまり1バイトのASCIIデータを書き込む関数はputc()である。

似たような関数にfputc()があるため、どちらを使ってよいか初心者は非常に困惑する。色々読んで調べてみたが、結局はどちらでもよいということになった(少なくとも初心者目線では)。教科書ではfputcが採用されているが、あえてputcを使って書き換えてみることにした。どちらでも同じ、というならば、ちゃんと動くだろう。

次のようなコードを書いて、ファイルに1文字書き込んでみた。果たして0.5byteになるのか、それとも1byteになるのか?

#include <stdio.h>

int main(){
  int c;
  FILE *fp;
  fp = fopen("sample.dat","wb");
  c = putc('B',fp);
  fclose(fp);
  return 0;
}

wbというのはバイナリ形式での(新規ファイルへの)書き込みモード、という意味である。書き込む先のファイルは「sample.dat」と名付けた。1文字を表すときは一重引用符で囲む。つまり'B'とする。二重引用符で囲む、つまり"B"とすると、C言語では文字列扱いになってしまうので、終端記号'\0'が自動的に付加されて2バイトになってしまう。putcが書き込めるのは、あくまで「一文字」だけであり、文字列は別の関数を使わねばならない。

さて出来上がったファイルを見ると、

% ls -l sample.dat

-rw-r--r--  1 user  group  1  3 26 12:50 sample.dat

% xxd sample.dat

00000000: 42                                       B

となっていて、1バイト、文字Bひとつだけが書き込まれたファイルになっている。では2文字書き込んだらどうなるか?

#include <stdio.h>

int main(){
  int c;
  FILE *fp;
  fp = fopen("sample.dat","wb");
  c = putc('B',fp);
  c = putc('M',fp);
  fclose(fp);
  return 0;
}

結果は、

% xxd sample.dat

00000000: 424d                                     BM

% ls -l sample.dat

-rw-r--r--  1 user  group  2  3 26 12:53 sample.dat

たしかに2バイトである。

この実験により明らかになったのは、ASCII文字コードは7bit仕様(128文字)であるが、現代のコンピュータアーキテクチャでは、8bit(256文字)=1byteとして扱っているということである。かなり無駄をしているように見えるが、無駄が無駄にならないくらい技術が進歩したのであろう。

ちなみに、「無駄」な情報といえば、生体アミノ酸を思い出す。生体アミノ酸は20種類あるそうだが、塩基配列RNAの遺伝子情報)では4進数表現(G,C,A,U)の3ビット(例えばGCCはアラニンを表す)で表現されているそうである。43=64なので、20種類のアミノ酸を表すには「多すぎる」のだが、この冗長さが進化論的に有利に働いたのであろう。これは、コロナウイルスが変異しやすい理由などと若干関係しているのかもしれないが、よくしらないのでこれ以上の深入はやめておこう.....

さて、putc()はint型なので、返値は整数で与えられる。前回確認したようにC言語の整数表現は4byte、つまり32ビット表現がデフォルトである。putcが返す整数値は、書き込んだASCIIコードになっているらしいことが実験をしてみるとわかる。しかし、ASCIIコードは0.5byteで十分だから、4byteの整数表現では「ゆるゆる」ではないだろうか? 

まずは、文字がASCIIコード(つまり0-127の整数)と等価であることを確認しておこう。putc('B',fp)の代わりに、BのASCIIコードである66を使っても同じ結果になる。つまり、

c = putc(66, fp);

c = putc('B',fp);

は同じ結果となるのである。K&Rにも

int putc( int, FILE *)

と書いてある。返値cを印字させると66と出力される。予想通り、書き込んだASCIIコードが返値になっているのである。

ここで実験である。ASCIIコードをはるかに超える大きな数、たとえば578を使ってファイルに書き込んでみたらどうなるだろうか?つまり、

c = putc(578, fp);

である。コンパイルはうまく通るので、これは間違いではないらしい。書き込んだファイルを覗いてみると、Bが書き込まれており、さらにcを出力すると66となっているではないか。つまり578は66なのである(?)どういうことかというと、現代のマシンではASCIIコードが8ビット(つまり256種)に拡張されており、256を超える大きな数は「周期的に」再定義されているらしいのである。つまり、

578 % 256 = 66

として元に戻ってくるように見えるのだ。ちなみに、putc(834,fp)もputc(21314,fp)も、Bを書き込む命令putc(66,fp)と等価であることを確認した。どれも余りが66となる整数である。

BMPの最初の2byteのデータを読み込む

ということで、bmp画像ファイルの最初の2バイトに書き込んであるという「魔法数」を読み取り、それがBMになっていることを確認してみよう。プログラムする前に、先ほど使ったxxdでBMPファイルをダンプしてみた。

% xxd out.bmp

00000000: 424d 6675 0000 0000 0000 3600 0000 2800  BMfu......6...(.
00000010: 0000 6400 0000 6400 0000 0100 1800 0000  ..d...d.........
00000020: 0000 0000 0000 2c01 0000 2c01 0000 0000  ......,...,.....
00000030: 0000 0000 0000 0000 00ff ffff 0000 00ff  ................
00000040: ffff 0000 00ff ffff 0000 00ff ffff 0000  ................
00000050: 00ff ffff 0000 00ff ffff 0000 00ff ffff  ................
00000060: 0000 00ff ffff 0000 00ff ffff 0000 00ff  ................
00000070: ffff 0000 00ff ffff 0000 00ff ffff 0000  ................
00000080: 00ff ffff 0000 00ff ffff 0000 00ff ffff  ................
00000090: 0000 00ff ffff 0000 00ff ffff 0000 00ff  ................
000000a0: ffff 0000 00ff ffff 0000 00ff ffff 0000  ................
000000b0: 00ff ffff 0000 00ff ffff 0000 00ff ffff  ................
000000c0: 0000 00ff ffff 0000 00ff ffff 0000 00ff  ................
000000d0: ffff 0000 00ff ffff 0000 00ff ffff 0000  ................
000000e0: 00ff ffff 0000 00ff ffff 0000 00ff ffff  ................
000000f0: 0000 00ff ffff 0000 00ff ffff 0000 00ff  ................
00000100: ffff 0000 00ff ffff 0000 00ff ffff 0000  ................
(以下省略)

確かに最初の2byteがBMになっている。しばらくはヘッダー情報が続くが、やがて0000と00ffの繰り返しになっていて、これが市松模様に対応していることがわかる。

では読み取ってみよう。

#include <stdio.h>

int main(){
  int c;
  FILE *fp;
  if( (fp = fopen("out.bmp","rb")) == NULL) {
     printf("File not found.\n");
     return 1;
  }
  printf("File found.\n");
  c = getc(fp);
  printf("%c %d %x\n", c, c, c);
  c = getc(fp);
  printf("%c %d %x\n", c, c, c);
  fclose(fp);
  return 0;
}

コンパイルして実行すると、次の結果を得る。

% a.out 

File found.
B 66 42
M 77 4d

out.bmpBMP画像データであることが確認できた!

BMPファイルの初めの14byteには、このような識別データがいくつか書き込まれているとBMPの仕様書には書いてある。2つ目のデータは4byteで構成されるファイルサイズのデータである。3つ目と4つ目は2byteずつ確保されているそうだが、零値が充てあられる規約となっているという。5つ目が4byteの情報で、画像データがどこから始まるか書いてあるそうである。

これらのデータは、BMという2バイトの識別データを読み込んだ要領と同様に読み取っていけばよい。

まずは、画像ファイルのヘッダー部分の読み取りができるようになった。この要領で、画像データ本体も読み取ってみたいと思うが、それは次回のテーマである。