Ads by Google
新しい記事を書く事で広告が消せます。
自然画像に対する球面調和関数の展開
3Dで使われることがある球面調和関数の勉強として、
球面調和関数を使用し、単純な画像変換を行ってみました。
球面調和関数とは？：
直感的に分かりやすい紹介をしてくれている
【コラム】3Dグラフィックス・マニアックス (68) 事前計算放射輝度伝搬(PRT)〜PRTの基本。静的PRT(3) | パソコン | マイコミジャーナル
によると（引用）
『球面調和関数とは、かなり極端に簡略化して説明すると、球っぽいモノから任意の長さの針が突き出た「いびつなウニ」みたいな物体の形状を数学的な関数で表すためのもの……といえる(ウニの針が逆に凹んでいる場合もある)。』
『MPEGやJPEGの圧縮に用いられる離散コサイン変換の「球体バージョン」ということができるかもしれない。』
ということである。（ウニという表現が非常にGoodだと思います。）
この球面調和関数を使うことで、球面上の位置を使って定義される関数を球面調和関数展開し、
近似関数をつくることができる。
高次の係数（ほぼ0に近くなる）を無視すれば、情報の圧縮が可能になります。
情報の圧縮が目的です。
本当に綺麗に近似できるのかと思ったので、
今回、次のような自然画像を球面上の二つの半球面に
それぞれ正射影で貼り付けられているモノとして考え、
球面調和関数展開を行って係数を計算し、
その係数を作った近似関数から画像を再構築してみました。
計算式は、
t-pot『テイラー、フーリエ、球面調和関数』で
紹介されていた式を離散的なものに修正して使用しました。
[結果]
元画像:
l=0の時(係数は1個) ｌ=5の時(係数は36個) ｌ=20の時(係数は441個)
l=50の時(係数は2601個)
「う〜ん・・」な結果となりました。
記事では、解像度の低い32×32テクセルの6面体のキューブ環境マップの例を使い、
「実用レベルで係数は16個程度でも十分」と書かれていますが、
今回使ったのは、512*512の画像です。
結果を見る限り、少ない係数では悲惨な状況になっており、
解像度が高い場合は、多くの係数を必要とすることが
分かりました。
*記事で係数は16個程度でも十分と書かれている理由は、
32*32の解像度の低い自然画像を用いており、解像度の低い画像は、
詳細部分が失われ周波数の低い画像となっているためだと考えられます。
そもそも、解像度の高い球面全体の自然画像に対して、球面調和関数展開しても、
球面全体の自然画像では、周波数が高くなり、高次の係数が必要になってきます。
ですので、解像度の高い球面全体の自然画像に対して、
球面調和関数展開するのは、よくないと思いました。
今回の実験を行ってみて、フーリエ級数の「球体バージョン」ということを
感じ取ることができました。
もっと調べてみるとおもしろそうです。
仕事忙しいけど、こつこつやってみよう！
参考：
【コラム】3Dグラフィックス・マニアックス (68) 事前計算放射輝度伝搬(PRT)〜PRTの基本。静的PRT(3) | パソコン | マイコミジャーナル
http://journal.mycom.co.jp/column/graphics/068/index.html
t-pot『テイラー、フーリエ、球面調和関数』
http://www.t-pot.com/program/88_SH/index.html
球面調和関数を使用し、単純な画像変換を行ってみました。
球面調和関数とは？：
直感的に分かりやすい紹介をしてくれている
【コラム】3Dグラフィックス・マニアックス (68) 事前計算放射輝度伝搬(PRT)〜PRTの基本。静的PRT(3) | パソコン | マイコミジャーナル
によると（引用）
『球面調和関数とは、かなり極端に簡略化して説明すると、球っぽいモノから任意の長さの針が突き出た「いびつなウニ」みたいな物体の形状を数学的な関数で表すためのもの……といえる(ウニの針が逆に凹んでいる場合もある)。』
『MPEGやJPEGの圧縮に用いられる離散コサイン変換の「球体バージョン」ということができるかもしれない。』
ということである。（ウニという表現が非常にGoodだと思います。）
この球面調和関数を使うことで、球面上の位置を使って定義される関数を球面調和関数展開し、
近似関数をつくることができる。
高次の係数（ほぼ0に近くなる）を無視すれば、情報の圧縮が可能になります。
情報の圧縮が目的です。
本当に綺麗に近似できるのかと思ったので、
今回、次のような自然画像を球面上の二つの半球面に
それぞれ正射影で貼り付けられているモノとして考え、
球面調和関数展開を行って係数を計算し、
その係数を作った近似関数から画像を再構築してみました。
計算式は、
t-pot『テイラー、フーリエ、球面調和関数』で
紹介されていた式を離散的なものに修正して使用しました。
[結果]
元画像:
l=0の時(係数は1個) ｌ=5の時(係数は36個) ｌ=20の時(係数は441個)
l=50の時(係数は2601個)
「う〜ん・・」な結果となりました。
記事では、解像度の低い32×32テクセルの6面体のキューブ環境マップの例を使い、
「実用レベルで係数は16個程度でも十分」と書かれていますが、
今回使ったのは、512*512の画像です。
結果を見る限り、少ない係数では悲惨な状況になっており、
解像度が高い場合は、多くの係数を必要とすることが
分かりました。
*記事で係数は16個程度でも十分と書かれている理由は、
32*32の解像度の低い自然画像を用いており、解像度の低い画像は、
詳細部分が失われ周波数の低い画像となっているためだと考えられます。
そもそも、解像度の高い球面全体の自然画像に対して、球面調和関数展開しても、
球面全体の自然画像では、周波数が高くなり、高次の係数が必要になってきます。
ですので、解像度の高い球面全体の自然画像に対して、
球面調和関数展開するのは、よくないと思いました。
今回の実験を行ってみて、フーリエ級数の「球体バージョン」ということを
感じ取ることができました。
もっと調べてみるとおもしろそうです。
仕事忙しいけど、こつこつやってみよう！
参考：
【コラム】3Dグラフィックス・マニアックス (68) 事前計算放射輝度伝搬(PRT)〜PRTの基本。静的PRT(3) | パソコン | マイコミジャーナル
http://journal.mycom.co.jp/column/graphics/068/index.html
t-pot『テイラー、フーリエ、球面調和関数』
http://www.t-pot.com/program/88_SH/index.html
3次元点集合を使ったドロネー分割
ドロネー分割は、頂点群を与えるだけで、三角ポリゴンのメッシュを
作ることができるというもの。
ドロネー分割を使えば、モデリングとか自動化できて楽になるかも？と思ったのと、
勉強のために3次元のドロネー分割を行うプログラムを作ってみた。
3次元空間に頂点群を与えると、それに応じたメッシュを作ってくれる。
こんなデータを与えると、
D3DXVECTOR3 v[10];
v[0] = D3DXVECTOR3(0.0f, 0.0f, 0.0f);
v[1] = D3DXVECTOR3(0.0f, 0.0f, 5.0f);
v[2] = D3DXVECTOR3(5.0f, 0.0f, 5.0f);
v[3] = D3DXVECTOR3(5.0f, 0.0f, 0.0f);
v[4] = D3DXVECTOR3(0.0f, 5.0f, 0.0f);
v[5] = D3DXVECTOR3(0.0f, 5.0f, 5.0f);
v[6] = D3DXVECTOR3(5.0f, 5.0f, 5.0f);
v[7] = D3DXVECTOR3(5.0f, 5.0f, 0.0f);
v[8] = D3DXVECTOR3(8.0f, 2.5f, 2.5f);
v[9] = D3DXVECTOR3(10.0f, 2.5f, 5.5f);
こんな感じになった。
形状を三角錐の集合で、つくることができるみたい。
頂点群を与えるだけで、メッシュを作ってくれるのは、おもしろい。
しかし、次のことも気ががついた。
「意図する形状モデルを作るためには、意図する形状モデルを作るための
頂点群を自分で考える必要がある。」
「テクスチャ貼るときは、結局、自分で座標を設定してやる必要がある」
「凹の部分を表現できない。」
いろいろ大変そう。
モデル作るなら、素直にモデリングソフト使った方がよさそうだと思った。
作ることができるというもの。
ドロネー分割を使えば、モデリングとか自動化できて楽になるかも？と思ったのと、
勉強のために3次元のドロネー分割を行うプログラムを作ってみた。
3次元空間に頂点群を与えると、それに応じたメッシュを作ってくれる。
こんなデータを与えると、
D3DXVECTOR3 v[10];
v[0] = D3DXVECTOR3(0.0f, 0.0f, 0.0f);
v[1] = D3DXVECTOR3(0.0f, 0.0f, 5.0f);
v[2] = D3DXVECTOR3(5.0f, 0.0f, 5.0f);
v[3] = D3DXVECTOR3(5.0f, 0.0f, 0.0f);
v[4] = D3DXVECTOR3(0.0f, 5.0f, 0.0f);
v[5] = D3DXVECTOR3(0.0f, 5.0f, 5.0f);
v[6] = D3DXVECTOR3(5.0f, 5.0f, 5.0f);
v[7] = D3DXVECTOR3(5.0f, 5.0f, 0.0f);
v[8] = D3DXVECTOR3(8.0f, 2.5f, 2.5f);
v[9] = D3DXVECTOR3(10.0f, 2.5f, 5.5f);
こんな感じになった。
形状を三角錐の集合で、つくることができるみたい。
頂点群を与えるだけで、メッシュを作ってくれるのは、おもしろい。
しかし、次のことも気ががついた。
「意図する形状モデルを作るためには、意図する形状モデルを作るための
頂点群を自分で考える必要がある。」
「テクスチャ貼るときは、結局、自分で座標を設定してやる必要がある」
「凹の部分を表現できない。」
いろいろ大変そう。
モデル作るなら、素直にモデリングソフト使った方がよさそうだと思った。
ハールのウェーブレット変換
ハールの基底を使った画像のウェーブレット変換の
プログラムを作ってみた。DCTより単純で楽。
記録を残しておく。
順変換のアルゴリズム：
(1)
画像を用意する。
(2)
画像の一行分のピクセルに対して、平均と差分を計算し、
左側に平均をおき、右側に差分を置いていく。
この処理を上の行から、下へ1行ずつ順番に行っていく。
このような処理を行い、左側に平均、右側に差分が集まった画像を新しく作る。
[ a, b, c, d ] →[ ( a+b)/2, (c+d)/2, (a-b)/2, (c-d)/2]
(3)
(2)の結果の画像における、一列分のピクセルに対して、平均と差分を計算し、
上側に平均をおき、下側に差分を置いていく。
この処理を左の列から、右の列へ1列ずつ順番に行っていく。
このような処理を行い、上側に平均、下側に差分が集まった画像を新しく作る。
(4)
結果として、次のような画像が作られる。
・左上には、元の画像の行方向の平均をし、列方向に平均をした画像
(元の画像を1/2にした画像)。
・左下には、元の画像の行方向の平均をし、列方向の差分をした画像。
・右上には、元の画像の行方向の差分をし、列方向の平均をした画像。
・右下には、元の画像の行方向の差分をし、列方向の差分をした画像。
これで、一つの階層の処理の終了となる。
さらに、階層を深くすることができる。
さらに階層を深くする場合は、左上の元の画像の行方向の平均をし、
列方向に平均をした画像に対して、(2)と(3)の処理を行っていけばよい。
before after
逆変換のアルゴリズム：
順変換と逆の工程を行えば、元の画像に戻すことができる。
列の走査をして、行の走査をすればよい。
[ ( a+b)/2, (c+d)/2, (a-b,)/2 (c-d)/2] == [A, B, C, D] とすると、
a = A+C;
b = A-C;
c = B+D;
d = B-D;
を行うことで、もとのデータ[ a, b, c, d ] に戻せる。
圧縮について：
順変換の結果画像に対して、適切な量子化を行い、
各種符号化すればよい。適切な量子化については、まだ調査中。
例：順変換の結果画像の高周波成分(左下、右上、右下の画像の部分)を0にするとか。
思ったこと：
・ハール以外の基底があるので、他の基底も試して見る必要がある。
・ハールは、整数で演算できるから良いと思う。
・SURFやいろんな所でも使われてるみたいだから、勉強の価値あり。
・適切な量子化の方法を勉強しなければいけない。
参考：
『画像圧縮アルゴリズム (8) ウェーブレット変換 -1-』
http://www2.starcat.ne.jp/~fussy/algo/algo8-8.htm
『画像圧縮アルゴリズム (9) ウェーブレット変換 -2-』
http://www2.starcat.ne.jp/~fussy/algo/algo8-9.htm
『直交ウェーブレット変換について』
http://laputa.cs.shinshu-u.ac.jp/~yizawa/InfSys1/advanced/daubechies/
『ウェーブレット変換の基礎と応用』
http://wwwmp.jwu.ac.jp/mine/SemPPT/2005sem/2005FoundationOfWaveletTransformRevised.pdf
わかりやすいスライドでした。
『USC-SIPI Image Database』
http://sipi.usc.edu/database/
レナとかの画像処理で使われる画像がおいてあるサイト。
プログラムを作ってみた。DCTより単純で楽。
記録を残しておく。
順変換のアルゴリズム：
(1)
画像を用意する。
(2)
画像の一行分のピクセルに対して、平均と差分を計算し、
左側に平均をおき、右側に差分を置いていく。
この処理を上の行から、下へ1行ずつ順番に行っていく。
このような処理を行い、左側に平均、右側に差分が集まった画像を新しく作る。
[ a, b, c, d ] →[ ( a+b)/2, (c+d)/2, (a-b)/2, (c-d)/2]
(3)
(2)の結果の画像における、一列分のピクセルに対して、平均と差分を計算し、
上側に平均をおき、下側に差分を置いていく。
この処理を左の列から、右の列へ1列ずつ順番に行っていく。
このような処理を行い、上側に平均、下側に差分が集まった画像を新しく作る。
(4)
結果として、次のような画像が作られる。
・左上には、元の画像の行方向の平均をし、列方向に平均をした画像
(元の画像を1/2にした画像)。
・左下には、元の画像の行方向の平均をし、列方向の差分をした画像。
・右上には、元の画像の行方向の差分をし、列方向の平均をした画像。
・右下には、元の画像の行方向の差分をし、列方向の差分をした画像。
これで、一つの階層の処理の終了となる。
さらに、階層を深くすることができる。
さらに階層を深くする場合は、左上の元の画像の行方向の平均をし、
列方向に平均をした画像に対して、(2)と(3)の処理を行っていけばよい。
before after
逆変換のアルゴリズム：
順変換と逆の工程を行えば、元の画像に戻すことができる。
列の走査をして、行の走査をすればよい。
[ ( a+b)/2, (c+d)/2, (a-b,)/2 (c-d)/2] == [A, B, C, D] とすると、
a = A+C;
b = A-C;
c = B+D;
d = B-D;
を行うことで、もとのデータ[ a, b, c, d ] に戻せる。
圧縮について：
順変換の結果画像に対して、適切な量子化を行い、
各種符号化すればよい。適切な量子化については、まだ調査中。
例：順変換の結果画像の高周波成分(左下、右上、右下の画像の部分)を0にするとか。
思ったこと：
・ハール以外の基底があるので、他の基底も試して見る必要がある。
・ハールは、整数で演算できるから良いと思う。
・SURFやいろんな所でも使われてるみたいだから、勉強の価値あり。
・適切な量子化の方法を勉強しなければいけない。
参考：
『画像圧縮アルゴリズム (8) ウェーブレット変換 -1-』
http://www2.starcat.ne.jp/~fussy/algo/algo8-8.htm
『画像圧縮アルゴリズム (9) ウェーブレット変換 -2-』
http://www2.starcat.ne.jp/~fussy/algo/algo8-9.htm
『直交ウェーブレット変換について』
http://laputa.cs.shinshu-u.ac.jp/~yizawa/InfSys1/advanced/daubechies/
『ウェーブレット変換の基礎と応用』
http://wwwmp.jwu.ac.jp/mine/SemPPT/2005sem/2005FoundationOfWaveletTransformRevised.pdf
わかりやすいスライドでした。
『USC-SIPI Image Database』
http://sipi.usc.edu/database/
レナとかの画像処理で使われる画像がおいてあるサイト。
DirectXのシェーダ(HLSL)における勾配関数ddx ddyの挙動調査
DirectXのシェーダ(HLSL)で使う勾配関数ddx ddyが行う、
具体的な演算方法がよく分からないので、調査した。（DirectX9）
DirectXのヘルプを読んでも、どのように傾きを求めているのかがわからない。
調査の方法：
浮動小数テクスチャに調べやすそうな任意の値を書き込んで、
デバッグを使って、地道にデータを見て、調べた。
結果：
次のような結果となった。
ddxの場合：
入力データ[0, 1, 4, 9, 16, 25, 36, 49] → 出力データ[1, 1, 5, 5, 9, 9, 13, 13]
この結果より、二つのピクセル値を使って、
傾き(ただの差分である。右 - 左)を求め、
二つのピクセルにその傾きを格納している。
(上の結果の 0, 1 は、差分を取ると1-0 == 1となる。
その1を対応する二つの出力ピクセルに入れている。)
(上の結果の 4, 9 は、差分を取ると9-1 == 5となる。
その5を対応する二つの出力ピクセルに入れている。)
ddyの場合：
ddxと同様な処理をy軸方向に行っていた。
まとめ：
・ddxとddy関数は、隣り合う二つのピクセルの値を使い、
その二つのピクセルの差分(右 - 左)を取り、
その結果を対応する二つの出力ピクセルに格納する関数である。
・ddx,ddyは、ベクトルデータ(float2,float3,float4とか)を与えることもできる。
ベクトルを与えた場合は、各要素で独立して差分を求める。
具体的な演算方法がよく分からないので、調査した。（DirectX9）
DirectXのヘルプを読んでも、どのように傾きを求めているのかがわからない。
調査の方法：
浮動小数テクスチャに調べやすそうな任意の値を書き込んで、
デバッグを使って、地道にデータを見て、調べた。
結果：
次のような結果となった。
ddxの場合：
入力データ[0, 1, 4, 9, 16, 25, 36, 49] → 出力データ[1, 1, 5, 5, 9, 9, 13, 13]
この結果より、二つのピクセル値を使って、
傾き(ただの差分である。右 - 左)を求め、
二つのピクセルにその傾きを格納している。
(上の結果の 0, 1 は、差分を取ると1-0 == 1となる。
その1を対応する二つの出力ピクセルに入れている。)
(上の結果の 4, 9 は、差分を取ると9-1 == 5となる。
その5を対応する二つの出力ピクセルに入れている。)
ddyの場合：
ddxと同様な処理をy軸方向に行っていた。
まとめ：
・ddxとddy関数は、隣り合う二つのピクセルの値を使い、
その二つのピクセルの差分(右 - 左)を取り、
その結果を対応する二つの出力ピクセルに格納する関数である。
・ddx,ddyは、ベクトルデータ(float2,float3,float4とか)を与えることもできる。
ベクトルを与えた場合は、各要素で独立して差分を求める。
