赤いグラフは周波数220Hzのsinwave、緑のグラフは直流を示している。 青い点は、sample周波数220Hzでsampleをとったときに とられる値を示している。 この図を見れば分かるように、ちょうどsinwaveが1になったときに sample値がとられていることがわかる。そのため 得られるサンプル値はすべて1である。 直流の方はもちろん1である。
ようやく手をつけることにした。以後、レポートのもととなる 研究やら解説やらをしてゆくことにする。 最終的にTeXソースにまとめてレポートとするつもりである(その レポートも公開する予定)
まずはレポートの周辺の状況などを。 学校のimagekouというユーザーのHOMEに行くと、テスト用の画像と レポートの課題が入っているのでそれを見る。 また、(一応)portableな(つまりHP-UX以外でも使える)画像ローダも 入っているので、必要ならコンパイルしなおして使うとよい。
なお、~imagekou/README.txtには、
注1: ここのimagesの下にある画像を扱うときには、リンクを張って コピーせずに使ってください。自分のディレクトリーにコピーすると、 計算機の中にメモリーをたくさん取ることになります。などと書かれているが(適当に改行した)、リンクを張るにはいくらかの 方法がある。
% cd ~ % ln -s ~imagekou ./imagekouまあ普通な方法。でもこれだと ~/imagekouの中で作業ができないので (実体は~imagekouなので書き込み権限がない)、 おれならlndirを使うことをすすめる。 lndirはXの配布ファイルの中に含まれているスクリプトだったと思う (HP上にもちゃんと置いてあった)。
% cd ~ % lndir ~imagekou ./imagekou違いは使い分けて実感すればよい :-)
さて、ここでの目的は、まず「自宅のFreeBSD上で研究する環境を 整える」である。いちいち学校の計算機室の使いづらいシステムを使ってまで 研究などしたくない :-)
上にもちょっと書いた通り、~imagekou/bin にあるxgpwmというプログラムを 使えば ~imagekou/images にある画像を見ることができる。 こいつは ~imagekou/src にソースがあるので、それをもらってきて コンパイルすればFreeBSD上でも問題なく使える。
ただし、カラーマップを全て奪って色を確保する(xboingみたいな)やりかたを しているので、TrueColorではエラーが出て使えない(と思う)。 うちの16bppの環境では使えなかった。まあ研究のときだけ8bppの環境で やればいいんだけどそれではメンドウだ。 TrueColorの環境で使うためのローダを作ろう。...と、できあがり。 xgpwm-true.19981026.tar.gz を利用可能である。
build方法はかんたん。上記ソースを展開した後
% xmkmf % makeでおっけーだ。あとはできあがったxgpwm-trueをpathの通ったディレクトリ なんかに置けばいつでも使えるはずだ。 以前作ったHGローダのソースを少々いじってでっちあげたものなので かなり恥ずかしい作りをしているが、まあ実用には問題ないだろう。 色は最初に256コ確保するというみっともない方法を使っているので 8bppだとうまく動かないと思う。8bppでもなんとか表示させる方法を 実装するのはそんなに難しくないけど需要ないので無視している。 8bppの環境ならオリジナルのものを使えばいいと思う。
しかしまあ、オリジナルのソースは「monocrome」だし 「X Windows」だし、ほんとにもう...。
画像はいわゆる「raw形式」などと呼ばれるものだ。 矩形画像の左上から右下までの濃度分布を 8ビットで格納してある。つまり0が一番暗くて(ふつうは黒) 255が一番明るい(ふつうは白)ということである。 それ以外のヘッダなどは一切ついていないので、 ファイルの中身から画像のサイズやらの情報を得ることはできない。 仕方ないのでオリジナルのソースでもやっているように 「ファイルのサイズを得てそれの平方根を縦と横のサイズとする」という やりかたをするしかない。 ~imagekou/images以下の画像はすべて正方形を仮定しているようなので、 これでうまく動く。おれとしては こんなコードは美観と激しく対立するので書きたくないんだけど。
いわゆるraw形式ということで、rawtopbm(netpbm)やconvert(Image)などを 使えば他の形式(GIFや(E)PSなど)に変換することができる。 例えば、~imagekou/images/child(256x256)をpsにするには、
% rawtopgm 256 256 child | convert PGM:- child.psなどとすればよい。いちいちサイズ指定をしなければならないのは 面倒なので、簡単なスクリプトを作ってみた。
% rawtops.pl childなどとすると child.psを得ることができる。 psファイルになったら、そいつをpsutilsなどで編集するとかすれば 印刷にもちゃんと使える。形式かえたかったら中身を適当に 書換えてね。EPSの方がTeXでレポート作るぶんには便利そうだけど。
以上のような話、 ドキュメントやレポート問題の中にはまったくとりあげられていない。 やる気あんのか、あんたたち。 まあローダのソースを読めばすぐわかるが(読めない奴の方が絶対多いと 思うけど。どうやってレポートでっちあげてくるのかねえ :-))、 まさか説明なしでも対処できるかどうかで実力をはかる、とか ビンボーくさいこと考えてるんじゃねえだろうな。
ようやく課題にとりくむ準備ができた。課題の内容は ~imagekou/report/imrepo98A1にある。
平成10年度 画像工学A レポート問題１ 10/13/'98出題 以下の1から３までを、10月末日までに提出のこと。 １．細かい構造は要らないので、細かい構造を持った画像から必要な 細かさの構造を持った画像に変換したい。単純に間引いてサンプリン グをすることでは、本来なかった構造が見えてくるエリアシングが 起きることを、１０２４x１０２４画素の画像（ディレクトリImagesに おいてあるsms0など）を粗くサンプルした画像で確かめよ。 ２．エリアシングが起きないようにするにはどうすればよいかを考えよ。 考えに従って、元の画像を前処理してからのサンプリングを実際に 行って示せ。普段からFT(フーリエ変換)、DCT(離散コサイン変換)、Wavelet変換などに 親しんでいる人間にとっては、何を今更...人をバカにしたような問題だが :-)
おおっと、そろそろ学校に行かねば。続きは帰ってきてからだな。
まずはダウンサンプリングをしませう。
こんなのは...さっさと作ってしまおう
(sqrt()使ってるからコンパイルには-lmおぷしょんをわすれずに)。
元の画像とそれをダウンサンプルした画像をのっけようと思うんだが。
rawだとでかいなあ。
というわけでgifにしてみる。
→元の画像(1024x1024)
→ダウンサンプルした画像(256x256)
元の画像にすでにAliasingが生じている。こまったもんだ。
さて問題にある通り、Aliasingが起きないように処理するのは どうすればよいかを考えるために、なぜAliasingが起きるのかということを 考えてみよう。といっても数学的に厳密な議論は 「FFTの使い方」(安居院猛・中嶋正之共著 秋葉出版)あたりを 見てもらうとして、多少ルーズでも感覚的にわかりやすい説明を めざそう(別に厳密な議論ができないわけではないのだが :-) 厳密にするにはそれなりのパワーが必要なので) 以後Sampling(標本化)はすべて一定の周期をもつ(いわゆるlinear)ものとする。
下に示すのはごくふつうのsinwaveである。
赤いグラフは周波数220Hzのsinwave、緑のグラフは直流を示している。
青い点は、sample周波数220Hzでsampleをとったときに
とられる値を示している。
この図を見れば分かるように、ちょうどsinwaveが1になったときに
sample値がとられていることがわかる。そのため
得られるサンプル値はすべて1である。
直流の方はもちろん1である。
sinwaveと直流、まったく違う原波形をsamplingしているのに、 得られるsampleデータが同じになる...逆に言うと、 この場合、sample値を見ただけでは 直流なのかsinwaveなのかの区別がつかないということになる。 これがAliasingの大雑把な原理である。 Aliasingが起こると、高い周波数成分が低い周波数成分に化ける (上の波形で言えば220Hzのsinwaveが直流に化ける)ことになる。 上のdownsampleの例で言うと、downsampleした画像になんか変な 模様が出たのがそのAliasingでできた模様ということになる。
このような現象がおきるのは「sample間隔が長いせいで、もっと
短くすればそんなことはなくなる」ということが
すぐに想像できると思う。実際のところ、
sample間隔をもっと狭く、例えば4.4kHzくらいにしてみると...
この通り、できあがったsample値がsinwaveか直流かは一目瞭然である。
sample間隔を狭くすればAliasingが起きないことがわかったが、 Aliasingが起きないようにするためには どのくらいまで狭くすればよいのだろう。逆に言うと、 どのくらいまでsample間隔を広げるとAliasingが起こるのだろう。 そこらへんをきっちりと定式化したのが 「Shannon-染谷のSampling定理」と呼ばれるものである (単に「ShannonのSampling定理」と呼ばれることの方が多い)。 簡単にまとめると、「元信号のもつ周波数成分の最大値fmaxの 2倍以上のsample周波数fsをもってsamplingすれば、 元の信号を完全に復元できる」ということである。 つまり、fs ≧ 2 * fmax を満たしていれば Aliasingを起こさずにsamplingができるということになる。
逆に言うと、fsでsamplingした信号には fmax = fs / 2 より大きい周波数を持つ 成分は存在しないことになる。 元の信号に fmax(今後nyquist周波数と呼ぶ) より大きい周波数成分が含まれていた場合、 nyquist周波数以下の周波数成分に折り返されることになる(Aliasing)
では1問目の問題でなぜAliasingが起きたのかということを 考えてみると... 例えば、4点おきにdownsampleするということは、sample周波数を 1/4にするということとおなじである(詳しくは「空間周波数」という 話をしなきゃならないけどまあ感覚的にわかると思う)。 そのため、元の信号に含まれた高い周波数成分が、downsampleで sample周波数が1/4になったことでsample周波数の2倍を越えてしまい Aliasingが起きたというふうに考えることができる。
ここまでくると、2問目の「Aliasingが起きないようにするには どうすればよいか」という答えもわかってくるだろう。 downsampleの前にAliasingが起きないように高い周波数成分を あらかじめ除去するようなフィルタ(lpf -- low pass filter)を 通しておけばよいのである。 このlpfのアルゴリズムは...
フーリエ変換の係数を...めんどうになったからやめた(^^;
やる気も戻ってそれなりに成果も出たので続きを書くかのぉ。
いきなりだけど...Octaveを 使うことにした。 OctaveはMATLABに よく似ている数学的処理をするツールである。 MATLABのソースは大抵Octaveでも通るのでとてもありがたい。 MATLABは商品だけど、OctaveはGPLにしたがって配布されている。
情報源はとりあえず以下のようなところがある。
数学的処理に色気を出した処理系といえば、 MathematicaとかU-BASICあたりがよく知られているだろうが、 行列演算を中心とした分野(信号処理など)では MATLABの方が向いているし、利用者数も上だと思う。 画像/音声処理の論文などでは、C(++)ソースと MATLABのソースがほとんど同じくらいの頻度で出現する(おれが 手に入れた限りでは :-))
上に書いたようなページをざっと見ればわかるだろうけど、 ここでも一応軽くまとめておくかあ。
めんどいので例だけ。
> A = [ 1 2 3 ;
4 5 6 ;
7 8 9 ];
といった感じである。なお、最後のセミコロンを
つけないと、最後に行列の内容を表示する。
行列の定義に限らず、最後のセミコロンは結果の表示をしないという
意味があるらしい。
from:step:to という形式で等差数列をつくることができる。例えば、 1:4:32 などとすると、「1 5 9 13 17 21 25 29」という 数列を得ることができる。 stepを省略することもできる(1:4 なら「1 2 3 4」となる)
A(i,j)で行列Aのi行j列の要素を得ることができる。 また、上に書いたような数列を引数にすることもできるので、 A(1,1:20)とすると、1行の1〜20列の要素を得ることができる。 また、「:」単体で用いるとすべての行(列)を得ることができる。 例えばA(1,:)は行列Aの1行目の要素だけを得る。
fft, ifftという組み込み関数を使えばよい。 2次元行列のときはfft2, ifft2 を使う。
スカラー倍するときは 2 * A などとする。 画像の表示は...めんどうだからあとで。
% raw2img.pl child > child.imgなお逆変換だが、こいつはOctaveのsaveimage関数を使えば ppmやps形式で出力できるし、imageやimshow関数で 起動されるxvからsave asしてもよいし...というわけなので 用意しない。
ここまでの説明を読めばわかると思うが、Octaveを使うといっこめの課題、 つまりダウンサンプリングの処理は...
1> [A, map] = loadimage("sms0.img");
2> B = A(1:4:1024, 1:4:1024);
これでおしまいである。1行目で画像を行列Aに読み込む。
mapはカラーマップだ。
sms0.imgはsms0をOctaveの形式になおしたもの。
2行目でデータの間引き。できた画像を見るなら
> colormap(map); > image(B, 1);でよろしい。
High Cut Filterともいう。時間(または空間領域)で フィルタを設計することもできるが (ここらへんは既に出ている次の課題で扱うようだ)、ここでは 感覚的にもわかりやすい周波数領域でのフィルタについて。
> [A, map] = loadimage("sms0.img");
> B = fft2(A);
Fourier変換すると行列Bに変換係数が入る。いわゆる複素Fourier変換なので
行列の要素Bは複素数。
さて。この変換係数によって周波数成分というやつがわかるのだが、
具体的にわ...この画像は1024x1024なので、1行(1列)が直流成分、
真中がNyquist周波数の成分、1024行(列)がサンプリング周波数より1サンプル間隔だけ
短かいもの。
画像処理系の本では、行列の真中が直流成分で外側に行くほど
高い周波数成分をあらわす...としていることが多い(光学的配置とかいうらしい)。
で、成分はNyquist成分を中心に対称になっとるので、...うーん
説明はめんどうだあ。
> B(256:767,:) = 0 + 0i; > B(:,256:767) = 0 + 0i; > C = real(ifft2(B));これで高域成分は除去できましたあ(;_;)
> image(A(1:4:1024,1:4:1024),1); > image(C(1:4:1024,1:4:1024),1);上が除去前、 下が除去後... Aliasingはかなりの部分おさえられているが、 もともとAliasingがおきていた部分はボケてしまって少々みっともないなあ。
この作業はめちゃメモリを食うので(うちで試したときは120MB食った) メモリの少ないマシンでは動かないかもしれない。 また、仮想メモリがいっぱいあるのに動かないという場合は shellのLimitを外そう。tcshなどなら
% unlimit datasizeでおっけーである。swapがきびしいけど、ある程度は動くと思う。 うちの仮想メモリサイズは200MB弱なので、これからさらに フクザツな処理をさせようと思うとちょっときびしいかもしれない。 ここらへん、グラフィックの処理を外部プログラムに任せっきりにする Octaveの弱点と言えるかな。
とりあえず今日じゅうにレポートをある程度作っておかないと 不幸なことになりそうなので、ちょっとがんばってレポート書くとしよう。 といっても絵を印刷して説明加えるだけなんだけど。
...と、いうわけでなんとかレポートをでっちあげた。 あんまり納得できないけど。基本がなってないんだよな :-) このドキュメントももうちょっと気合い入ったものにするつもりだったんだけど なあ...今後がんばるとしよう。