«前の日(08-31) 最新 次の日(09-02)» 追記

ema log


2007年09月01日 この日を編集

_ [最近] 微妙に不調

なんか、たとえるならばバイオリズムの下の方にいる感じ。まぁ、いつも通りといえばいつも通りなのですが。

なんとかレポートとかも仕上げますよ。AS3 いじくるぐらいの元気はありますから。Squeak はわけわかめ。情報の探し方が皆目わからない。ようやく、クラスの作り方とメソッドの追加方法っぽいのがわかったぐらい。

AS3 にしても、(ありがちな入門的なものを除いては)情報はなくて実質的にはリファレンスと格闘しながらの手探りなんですが、とりあえず リファレンス と 文法 と trace(いわゆる printf 的なもの)の使い方さえ身につければ、手の探りようがあります。しかも、わかってきたら、流石に良くできてます。型付けが面倒なぐらい。

Squeak は リファレンス も 文法 も デバッグ出力 もまだわからないw。なんかよくわからないけど、動いたよレベルから進めない。独自の世界過ぎる。わかれば良くできた環境なんだろうという手応えがあるだけに、力不足が残念。

_ [Linux] いい加減、再度 coLinux に挑戦するかなぁ。

しかし、coLinux なり導入して Linux 環境で作業した方が効率が良いのではないだろうか。vim (not emacs) をマスターした方が良いのではないだろうか。

そう思い続けてもはや一年近い。慣れというのは最大の敵。

_ [最近] ウィーナフィルタとか

頑張って数式追っかけて、最後の結論は、「理論的には最善だけど、パラメータは環境依存でサンプリング以外で決定できない。だから試行錯誤してね」って結論だったのが悲しい。流石、工学。


2008年09月01日 この日を編集

_ [最近]Welcome back! Ubuntu!!

CentOS → Ubuntu に戻した。NISまわりの設定にてこずるも二時間ほどでほぼ復旧完了。 テキストインストーラじゃないと入れれなかったので、alternate の方をダウンロード。

tex やら dvi やら、勝手しったる環境と違いすぎて困ったからなのですが。 それにしても、yum より aptitude のほうが反応早いなぁ。

_ [Programming]画像の拡大縮小(サンプリングに基づく補間アルゴリズム比較)

ディジタル画像処理

以前は抵抗があった「重み付き加算」だとか「サンプリングに基づく*1」だとか、そういうアレな単語を使うのに抵抗を覚えない。ちょっと残念だ。もっとまろやかに書けたらいいのに。でも、まとまったページがなかったんだ。堅い書き方でごめんなさい。背景理論はさておき、ざっくりと統一感のあるサンプル付きで紹介が目的。自分用まとめ。

絵がないとわかりにくいですね。ごめんなさい。絵については「ディジタル画像処理」か「WhiteForest/白ノ森」さんかを参照のこと*2

サンプリングによる方法には、ニアレストネイバー、バイリニア、バイキュービック、Lanczos2、Lanczos3 などがある。難しそうな名前がついているが、それぞれ 1点、4点、16点、16点、36点を元に変換後の画素値を決定する。その際の重み付けの違いが各手法。「Lanczos」は「ディジタル画像処理」に載ってなかった。理論的な話が知りたい人は、「ディジタル画像処理」や「まるも製作所 - 縮小アルゴリズム」辺りをヒントにあさると良いと思います。縮小限定だと「平均画素法」というのもあるみたいです。

Nearest NeighborBiLinearBiCubicLanczos2
点数141616
重みなし距離Bicubic_spline(距離)lanczos(距離)

バイリニア以上だと、補間材料が画像の範囲外になることがあるので注意。といっても、その辺の対処方法について書いてる本/ページがないんだよなぁ。使える点で補間するのが良いんだろうけど。

後、基本的に「無い/無くなった」情報を再現しようというのだから、それぞれの方法の善し悪しはあくまでもケースバイケースになる。困ったもんだ。

基本

変換先から、変換元の画素座標を計算する。拡大縮小以外にも応用可能だが、拡大縮小に限定して話を進める。

変換元画像srcの幅と高さをsw,sh、変換先画像dstの幅と高さをdw,dhとし、ある dx,dy の点の画素値を計算することを考える。このとき、dx, dy を変換元画像の座標 sx,syに変換すると、以下の変換になる。

float sx = dx * sw / dw;
float sy = dy * sh / dh;

諸事情により、サンプルは擬似コードとCで、Cの場合OpenCVで読み書きしてるのを前提に BGR の順に画素が連続したメモリ配置を仮定している。以下、これらの変数を使って説明。limitRGB は次のような関数。

unsigned char limitRGB(int color){
  if (color < 0x00) return 0x00;
  if (color > 0xff) return 0xff;
  return color;
}

ちなみに、ElectricFence などを使うと、範囲外アクセスを簡単にチェックできます。

ニアレストネイバー

重み付けせず、最も近い画素値を採用する。補間がないので、はっきり言って汚い。

dst(dx,dy) = src( round(sx), round(sy) ); // round は四捨五入。 floor(sx+0.5) と思いなせえ。

ただ、ドット絵の編集とか処理結果を見るときとかはこれじゃないと困ります。Windows の画像ビュワーや Eye of GNOME なんかだと、補間されちゃって困ることが。

使わないからCのソースはなし。以下、Cのソース付き。BGRの順に並んだ1画素ずつループ内で計算しているので、dst[hoge] = src[hoge]... の式が 3 つ並んでます。

バイリニア

近辺4画素を元に、線形補間する。線形補間と書くと難しそうだが、距離に応じて4点を重み付けして加算するだけ。縦横2方向の線形補間だからバイリニア。

x1 = floor(sx); y1 = floor(sx);
x2 = x1 + 1; y2 = y1 + 1;

をの 2 x 2 の組み合わせで 4 点を使う。点の取り方は自由度があるはず。

sx==x1, sy==y1 の点は、ニアレストネイバーと同じになる。中途半端なところで、「間を補う」わけです。

欠点としては、拡大だと平滑化*3がかかり、縮小だとモアレがでたりする。例えば拡大だと引き延ばすだけだから、伸びちゃうイメージ。画像だとピンとこないかもしれないけど、音声だとわかりやすい。線形補間で半分の速度で再生すると低くなるのと一緒。逆に単純に縮めると、高くなっちゃうのね。

というわけで、周波数領域でみるとよろしくない(再現性が低い)わけです。縮小したときに、モアレと呼ばれる縞が出たりします。cf.画像を縮小するときにモアレを防ぐには - 日経トレンディネット

以下、点数を増やすと同時に重み付けについても改善が図られています(といっても、2倍に拡大とかだと以下の手法とほとんど見分けがつかないのも事実。後、縮小の方が差がわかりやすい)。

void BilinearInterpolation( unsigned char * src, unsigned char * dst, int sw, int sh, int dw, int dh )
{
  float sx,sy;
  int dx,dy; // source xy, destination xy
  int bx,by;
  float x1,x2,y1,y2;
  unsigned int idx, step;
  unsigned int base1, base2;
  idx = 0;
  step = sw*3;
  for(dy=0; dy<dh; ++dy) {
    for(dx=0; dx<dw; ++dx, idx+=3) {
      sx = (float)dx*sw/dw;
      sy = (float)dy*sh/dh;
      bx = (int)floor((double)sx); by = (int)floor((double)sy);
      x1 = bx-sx+1.0; y1 = by-sy+1.0;
      x2 = sx-bx;     y2 = sy-by;
      base1 = by*step + bx*3;
      base2 = base1 + step;

      dst[idx+0] = limitRGB((int)(
        y1 * (x1*(float)src[base1+0] + x2*(float)src[base1+3+0]) +
        y2 * (x1*(float)src[base2+0] + x2*(float)src[base2+3+0])
      ));
      dst[idx+1] = limitRGB((int)(
        y1 * (x1*(float)src[base1+1] + x2*(float)src[base1+3+1]) +
        y2 * (x1*(float)src[base2+1] + x2*(float)src[base2+3+1])
      ));
      dst[idx+2] = limitRGB((int)(
        y1 * (x1*(float)src[base1+2] + x2*(float)src[base1+3+2]) +
        y2 * (x1*(float)src[base2+2] + x2*(float)src[base2+3+2])
      ));
    }
  }
}

バイキュービック

x1 = x1 - 1; y1 = y1 - 1;
x2 = floor(sx); y2 = floor(sx);
x3 = x1 + 1; y3 = y1 + 1;
x4 = x1 + 2; y4 = y1 + 2;

上記、近辺16画素を元に補間する。サンプリング定理に基づいた画像に対して、理論上最適な変換らしい。よく分からない。「拡大縮小アルゴリズム(バイキュービックとLanczos) (きっちん)」のグラフを見ると雰囲気が分かった気になれて良い感じ。

Bicubic splineはパラメータとしてαをもっていて

という式を使う様子(Bicubic interpolation - Wikipedia, the free encyclopediaより引用)。「ディジタル画像処理」の式はα=-1の式かな。Wikipediaには-0.5〜-0.75ぐらいで使うとあるんだけどなぁ。

シャープな画像が得られるけど、シャープすぎたりすることもある様子。

float bicubic_spline(float t)
{
  float tmp = fabs(t);
  if ( tmp <= 1.0 ) return ( tmp*tmp*tmp - 2.0 *tmp*tmp          + 1.0);
  if ( tmp <= 2.0 ) return (-tmp*tmp*tmp + 5.0 *tmp*tmp -8.0*tmp + 4.0);
  return 0.0;
}

void BicubicInterpolation( unsigned char * src, unsigned char * dst, int sw, int sh, int dw, int dh )
{
  float sx,sy;
  int dx,dy; // source xy, destination xy
  int bx,by;
  float x1,x2,x3,x4,y1,y2,y3,y4;
  float hx1, hx2, hx3, hx4, hy1, hy2, hy3, hy4;
  unsigned int idx, step;
  unsigned int base1,base2,base3,base4;
  idx = 0;
  step = sw*3;
  for(dy=0; dy<dh; ++dy) {
    for(dx=0; dx<dw; ++dx, idx+=3) {
      sx = (float)dx*sw/dw;
      sy = (float)dy*sh/dh;
      bx = (int)floor((double)sx); by = (int)floor((double)sy);

      x1 = 1.0 + sx - bx; y1 = 1.0 + sy - by;
      x2 =       sx - bx; y2 =       sy - by;
      x3 = bx + 1.0 - sx; y3 = by + 1.0 - sy;
      x4 = bx + 2.0 - sx; y4 = by + 2.0 - sy;

      hx1 = bicubic_spline(x1); hy1 = bicubic_spline(y1);
      hx2 = bicubic_spline(x2); hy2 = bicubic_spline(y2);
      hx3 = bicubic_spline(x3); hy3 = bicubic_spline(y3);
      hx4 = bicubic_spline(x4); hy4 = bicubic_spline(y4);

      base2 = by*step + bx*3;
      base1 = base2 - step;
      base3 = base2 + step;
      base4 = base3 + step;

      dst[idx+0] = limitRGB((int)(
        hy1 * ( hx1*src[base1-3+0] + hx2*src[base1+0] + hx3*src[base1+3+0] + hx4*src[base1+6+0] ) +
        hy2 * ( hx1*src[base2-3+0] + hx2*src[base2+0] + hx3*src[base2+3+0] + hx4*src[base2+6+0] ) +
        hy3 * ( hx1*src[base3-3+0] + hx2*src[base3+0] + hx3*src[base3+3+0] + hx4*src[base3+6+0] ) +
        hy4 * ( hx1*src[base4-3+0] + hx2*src[base4+0] + hx3*src[base4+3+0] + hx4*src[base4+6+0] )
      ));
      dst[idx+1] = limitRGB((int)(
        hy1 * ( hx1*src[base1-3+1] + hx2*src[base1+1] + hx3*src[base1+3+1] + hx4*src[base1+6+1] ) +
        hy2 * ( hx1*src[base2-3+1] + hx2*src[base2+1] + hx3*src[base2+3+1] + hx4*src[base2+6+1] ) +
        hy3 * ( hx1*src[base3-3+1] + hx2*src[base3+1] + hx3*src[base3+3+1] + hx4*src[base3+6+1] ) +
        hy4 * ( hx1*src[base4-3+1] + hx2*src[base4+1] + hx3*src[base4+3+1] + hx4*src[base4+6+1] )
      ));
      dst[idx+2] = limitRGB((int)(
        hy1 * ( hx1*src[base1-3+2] + hx2*src[base1+2] + hx3*src[base1+3+2] + hx4*src[base1+6+2] ) +
        hy2 * ( hx1*src[base2-3+2] + hx2*src[base2+2] + hx3*src[base2+3+2] + hx4*src[base2+6+2] ) +
        hy3 * ( hx1*src[base3-3+2] + hx2*src[base3+2] + hx3*src[base3+3+2] + hx4*src[base3+6+2] ) +
        hy4 * ( hx1*src[base4-3+2] + hx2*src[base4+2] + hx3*src[base4+3+2] + hx4*src[base4+6+2] )
      ));
    }
  }
}

Lanczos2

読めない。こちらも近辺16画素を元に補間する。Lanczos 窓関数で、画素の周囲 ±2(つまり4x4で16画素)を使うので、Lanczos2 という名前のようだ。Lanczos3 は 6x6 の 36 点になる。「Lanczos resampling - Wikipedia, the free encyclopedia」の式でいうと、a というのに、この 2 や 3 が入る。

L(x) = sinc(pi*x) * sinc(pi*x/a) = sin(pi*x)*sin(pi*x/a) / (pi^2 * x^2)

が基本の式になる(ただし、x=0 のとき L(x)=1、|x| > a のとき L(x) = 1)。ひとまず、近似*4せず実装した。

バイキュービックの sinc を lanczos に置き換えればいい。

#define PI 3.14159265358979
#define LIMIT 0.0000001
// @see http://en.wikipedia.org/wiki/Lanczos_resampling
float lanczos(float x)
{
  float pi_x;
  if ( -LIMIT<x && x<LIMIT ) return 1.0;
  if ( x < -2.0 || x > 2.0 ) return 0.0;

  pi_x = PI*x;
  return 2.0*sin(pi_x) * sin(pi_x/2.0) / (pi_x*pi_x);
}

void BicubicInterpolation( unsigned char * src, unsigned char * dst, int sw, int sh, int dw, int dh )
{
  float sx,sy;
  int dx,dy; // source xy, destination xy
  int bx,by;
  float x1,x2,x3,x4,y1,y2,y3,y4;
  float hx1, hx2, hx3, hx4, hy1, hy2, hy3, hy4;
  unsigned int idx, step;
  unsigned int base1,base2,base3,base4;
  idx = 0;
  step = sw*3;
  for(dy=0; dy<dh; ++dy) {
    for(dx=0; dx<dw; ++dx, idx+=3) {
      sx = (float)dx*sw/dw;
      sy = (float)dy*sh/dh;
      bx = (int)floor((double)sx); by = (int)floor((double)sy);

      x1 = 1.0 + sx - bx; y1 = 1.0 + sy - by;
      x2 =       sx - bx; y2 =       sy - by;
      x3 = bx + 1.0 - sx; y3 = by + 1.0 - sy;
      x4 = bx + 2.0 - sx; y4 = by + 2.0 - sy;

      hx1 = lanczos(x1); hy1 = lanczos(y1);
      hx2 = lanczos(x2); hy2 = lanczos(y2);
      hx3 = lanczos(x3); hy3 = lanczos(y3);
      hx4 = lanczos(x4); hy4 = lanczos(y4);

      base2 = by*step + bx*3;
      base1 = base2 - step;
      base3 = base2 + step;
      base4 = base3 + step;

      dst[idx+0] = limitRGB((int)(
        hy1 * ( hx1*src[base1-3+0] + hx2*src[base1+0] + hx3*src[base1+3+0] + hx4*src[base1+6+0] ) +
        hy2 * ( hx1*src[base2-3+0] + hx2*src[base2+0] + hx3*src[base2+3+0] + hx4*src[base2+6+0] ) +
        hy3 * ( hx1*src[base3-3+0] + hx2*src[base3+0] + hx3*src[base3+3+0] + hx4*src[base3+6+0] ) +
        hy4 * ( hx1*src[base4-3+0] + hx2*src[base4+0] + hx3*src[base4+3+0] + hx4*src[base4+6+0] )
      ));
      dst[idx+1] = limitRGB((int)(
        hy1 * ( hx1*src[base1-3+1] + hx2*src[base1+1] + hx3*src[base1+3+1] + hx4*src[base1+6+1] ) +
        hy2 * ( hx1*src[base2-3+1] + hx2*src[base2+1] + hx3*src[base2+3+1] + hx4*src[base2+6+1] ) +
        hy3 * ( hx1*src[base3-3+1] + hx2*src[base3+1] + hx3*src[base3+3+1] + hx4*src[base3+6+1] ) +
        hy4 * ( hx1*src[base4-3+1] + hx2*src[base4+1] + hx3*src[base4+3+1] + hx4*src[base4+6+1] )
      ));
      dst[idx+2] = limitRGB((int)(
        hy1 * ( hx1*src[base1-3+2] + hx2*src[base1+2] + hx3*src[base1+3+2] + hx4*src[base1+6+2] ) +
        hy2 * ( hx1*src[base2-3+2] + hx2*src[base2+2] + hx3*src[base2+3+2] + hx4*src[base2+6+2] ) +
        hy3 * ( hx1*src[base3-3+2] + hx2*src[base3+2] + hx3*src[base3+3+2] + hx4*src[base3+6+2] ) +
        hy4 * ( hx1*src[base4-3+2] + hx2*src[base4+2] + hx3*src[base4+3+2] + hx4*src[base4+6+2] )
      ));
    }
  }
}
References
まるも製作所 - 縮小アルゴリズム - Diary 2000-12
まるも製作所 - 縮小アルゴリズム - Diary 2001-7
陶見の窓 画像縮小アルゴリズム(?)
拡大縮小アルゴリズム(バイキュービックとLanczos) (きっちん)
Lanczos resampling - Wikipedia, the free encyclopedia
Bicubic interpolation - Wikipedia, the free encyclopedia
ディジタル画像処理
画像を縮小するときにモアレを防ぐには - 日経トレンディネット
WhiteForest/白ノ森
Lanczos関数による画像の拡大縮小

*1 根本的に方向性の異なる方法としては、周波数領域に変換してから、戻すってのもある。けど、ここでは扱わない。周波数領域でやるのが多分、一番綺麗なんだと思う。JPEG なんかで圧縮された画像の拡大縮小なんかだと、デコード時に拡大縮小できないのかな?その辺よく分からない。

*2 後者のページ、参考文献があげられてないけど、たぶん「ディジタル画像処理」が元ネタ

*3 ニアレストネイバー + 単純平均による平滑化と一緒ですから

*4 参考にしたページには書いてなかったけど、丸ごとテイラー展開するなり、バイキュービックのときの近似式を元に近似式作ればいいかな?後者はバイキュービックの式の誤差の項がないから、どれぐらいの精度か分からないけど。Lanczos2 とかならテーブル作っちゃうのも有りだろうなぁ。

_ [Programming]補間の補完

ついでに蛇足かもしれないけど補足。

拡大縮小って、飛び飛びの点を一回連続的なデータに近似して、再度離散値に戻すわけで。だもんで、たぶん、周波数領域でやるのが一番いい。ただ、どう繋ぐべきかは、何点考えるかや、曲がり具合の調整が色々あると。そういう問題。

ドット絵は連続的な領域に持って行きたくないから、ニアレストネイバーがいいわけです。

後、高速化についてはメモリ周りのケアが重要かと。prefetch 入れるだけでも変わるんじゃないかなぁ。SSE2 の整数演算は MMX の命令を二回実行するだけとか、ショボーンな実装だったりするみたいですし。

個人的な趣味ですが、冗長でも繰り返しが見えたり、並列性のある部分が目に見えるように書くようにしています。

本日のツッコミ(全2件) [ツッコミを入れる]

_ ねこ [始めまして、ソースを参考にさせて頂きました。 バイリニアは問題なく縮小できたのですが、 バイキュービックの dst..]

_ ema [諸般の事情で、動作確認しているソースそのままではありません。今ざっとみたら、確かにループの回し方がまずくて、確実に範..]