Java画像処理をはじめよう。初めての縮小(・∀・)

今までこのブログでも何回か拡大縮小の話が出てきたことがありますが、もう一度やっておきましょう。
今回は縮小についてです。

最も簡単な縮小処理　間引き法（最近傍法）

　画像を小さくするとピクセルも一緒に小さくなってくれたら嬉しいのですが、そんなわけにはいきません………。
　画像を小さくするとピクセル数が減るので表現できる情報が少なくなってしまいます。したがって、いったいどの情報を残して、どの情報を削除するのか、が一番重要なことになります。
　最初に説明するのは間引き法という物です。これは単純に、容量オーバーしているピクセル分まるっと削除するやり方です。
　例えば、3×3ピクセルの画像を2×2ピクセルの画像に縮小する場合、各方向1ピクセル多いので下の図の様に削除して2×2の画像を得ます。

　あ（ー△ー；）　今、この方法馬鹿にしたでしょ。「ふふん、こんなの使えねーよ」って言ったでしょ(-_-)
　あまいな〜。この方法、実は最近の手法に繋がるかもしれないんだよ？
　いったいどこをどのように間引くのか。これをどういう関数にするかっていうだけで最近の技術になるんだよ。

　ま、今回はただ単純に、縮小後の座標を元の画像の座標に戻したとき、もっともそこに近い画素を残して、他を削除する最近傍法でやります。この方法は全く持って新しくないよ(´Д⊂ヽ
　何？どーでも良いから最新技術を教えろ？
　嫌だ！
Seam carvingで検索

/**
 * 間引き法による縮小
 * @param img 元画像
 * @param zoom 縮小率（0<zoom<1）
 * @return 縮小画像
 */
public static BufferedImage thinningout(BufferedImage img,double zoom){
    if(zoom <=0 || zoom >=1)
        throw new IllegalArgumentException("zoom:"+zoom);
    int w = img.getWidth();
    int h = img.getHeight();
    //縮小後の大きさ
    int zoomw=(int)ceil(zoom*w);
    int zoomh=(int)ceil(zoom*h);

    BufferedImage ret = new BufferedImage(zoomw,zoomh,TYPE_INT_RGB);
    for(int y=0;y<zoomh;y++){
        //元の画像での座標yを求めます
        int imgy=(int)(y/zoom);
        for(int x=0;x<zoomw;x++){
            //元の画像での座標xを求めます
            int imgx=(int)(x/zoom);
            int color = img.getRGB(imgx, imgy);
            ret.setRGB(x, y, color);
        }
    }
    return ret;
}

縮小前

縮小後（50％）

面積平均法で縮小する

　さっきの結果はけっこう汚いです。間引いたデータは無かった物として扱っているから縮小結果ががたがたになってしまいます。
　これを改善する為に、「ピクセルの色の平均値を使う」という発想は自然に出てくると思う。いわゆる面積平均法という奴だ。下の図の赤い枠は、縮小後のピクセルが元の画像のピクセルに対応する範囲を表す。
　

　この赤枠の大きさは、2分の一に縮小するなら2×2ピクセル、3分の一に縮小するなら3×3ピクセルだ。

　これをそのまま実装してみると下のようになります。
もう二度と実装したくねー。思った以上に面倒くさい。(´Д⊂ヽ一発で動かなかせなかったｗ

public static BufferedImage areaAvarage(BufferedImage img,double zoom){
    if(zoom <=0 || zoom >=1)
        throw new IllegalArgumentException("zoom:"+zoom);
    int w = img.getWidth(),h = img.getHeight();
    int zoomw=(int)ceil(zoom*w),zoomh=(int)ceil(zoom*h);
    BufferedImage ret = new BufferedImage(zoomw,zoomh,TYPE_INT_RGB);
    
    final double length = 1/zoom;//赤い四角の長さ
    //ここの部分は後で修正するよ。
    for(int i=0;i<zoomw*zoomh;i++){
        int x= i%zoomw;
        int y= i/zoomw;
        
        double imgx = x/zoom;
        double imgy = y/zoom;
        double r=0,g=0,b=0,
                s=0;//面積
        for(double dy=imgy;dy<imgy+length;){
            int Y = (int)floor(dy);
            if(Y>=h)break;
            double nextdy=(int)floor(dy+1);
            if(nextdy>imgy+length)nextdy=imgy+length;
            double ry = nextdy-dy;
            
            for(double dx=imgx;dx<imgx+length;){
                if((int)floor(dx) >=w)break;
                double nextdx=(int)floor(dx+1);
                if(nextdx>imgx+length)nextdx=imgx+length;
                double rx = nextdx-dx;
                //このピクセルが赤い枠に入っている面積
                double S = rx*ry;
                int c = img.getRGB((int)floor(dx), Y);
                r += r(c)*S;//面積の重みをかけて和をとる。
                g += g(c)*S;
                b += b(c)*S;
                s +=S;
                dx = nextdx;
            }
            dy=nextdy;
        }
        r/=s;//平均を出すために面積で割る。
        g/=s;
        b/=s;
        int rgb = rgb((int)r,(int)g,(int)b);
        ret.setRGB(x, y, rgb);
    }
    return ret;
}

　え？読みにくい？私も読みにくいと思うぜ。

　先のようにギザギザした感じが消えましたね。

分解を利用して面積平均法を簡略、高速化

　さて、ここら辺から少し難しい話に入っていきましょう。
　面積平均法というのは、元の画像のそれぞれのピクセルの値を面積重みの平均で再度定義した後、間引き法で画素を間引く事と同等です。（正確には「間引く」ではなく「リサンプリング」です。）
　間引く前の、面積重みの平均で定義するということを式で表すと
$I_1(X,Y)=\frac{\sum_y\sum_xS(x,y)I_0(x,y)}{len^2}$
　となります。（式中のlenはソースコード内のlength（=1/zoom）の略です。）
　ここで、I(x,y)は点x,yにおける画像Iの色を表します。（今後In(x,y)など、大文字Iの後に番号nがきて、(x,y)などと書いてある式は画像Inの点x,yにおける色を表すことにします。）

　I₀は元の画像を表し、I₁は面積重みの平均で再度定義された画像になります。（なお、間引く前の話ですので画像の大きさはどちらも同じです。）
　S(x,y)は先の説明で用いた赤枠の面積です。

　Σのx,yの範囲は書くと非常に読みにくくなったので省略していますが、X-length/2からX+length/2までです。yも同様です。

　ピクセルは1という幅を持っており、座標(1.2,2.5)などといった小数点を含む座標のピクセルは存在しません。（離散化されていると言います）
　ですが、ここでは仮に小数点の座標でも画像の色が存在すると仮定しましょう。（連続と言います）
　すると先ほどの式は

　となります。
　ここで矩形関数R(x)を導入してみましょう。
$R(x)=\left\{\begin{array}{ll}1&if |x|\leq len/2\\0&else\end{array}\right.$
　これを用いると
$I _{1} \left(X,Y\right) =\frac{1}{len^2} \int_{- \infty }^{ \infty } R \left(y\right) \int_{- \infty }^{ \infty } R \left(x\right) I _{0} \left(x+X,y+Y\right) dxdy$
$I _{1} \left(X,Y\right) =\frac{1}{len^2} \int_{- \infty }^{ \infty } R \left(y\right) \left\{\int_{- \infty }^{ \infty } R \left(x\right) I _{0} \left(X-x,Y-y\right) dx\right\}dy$
　この式の形は畳み込みと言います。
　面積平均法とは、結局の所画像に対して矩形関数をたたみ込むという処理なのです。

　まぁ、今は畳み込みがどうのこうのは放っておきましょう。
　上の式を見ていただくと気がつかれると思いますが、xとyの積分はそれぞれ綺麗に分離されています。xの積分をするときはyは定数と見なして積分した後、yについて積分すればいいと言うことです。
　先ほどのプログラムソースコードではx軸とy軸を同時に計算していましたが、それぞれの方向に分解して実行してやればよいと言うことです。書くのが楽になるだけでなく、計算量も減るので高速になります。(元の計算量はlength^2に比例するのに対して、分離すると2*lengthに比例します。普通二乗より2倍の方が明らかに小さいですよね。)

　この様にたたみ込む関数（フィルタ）を
$A(x)\times B(y)$
　と分離できるフィルタは常にx軸とy軸のそれぞれの畳み込みに分離できます。

　さて、長くなりましたが面積平均法をx,y軸に分離してみましょう。
　先ほどに比べて非常に長くなっていますが、それぞれのループは中身が簡単になっています。

　簡単になっていますが、同じピクセルを2度読み込んだりしないように効率よく書いていますため、結構読みにくいかもしれません。簡単に説明すると「縮小画像の座標から元画像の座標を求める」のではなく、「元画像の座標から縮小画像の座標を求める」という処理になっています。これは、データの読み込み回数を減らした高速化のテクニックの一つです。

　Javaは配列読み込みが遅いので、読み込み回数を減らすと劇的に速度上がりますよ。

public static BufferedImage areaAvarage(BufferedImage img,double zoom){
    if(zoom <=0 || zoom >=1)
        throw new IllegalArgumentException("zoom:"+zoom);
    int w = img.getWidth(),h = img.getHeight();
    int zoomw=(int)ceil(zoom*w),zoomh=(int)ceil(zoom*h);
    //中間画像が必要になります
    BufferedImage imp = new BufferedImage(zoomw,h,TYPE_INT_RGB);
    BufferedImage ret = new BufferedImage(zoomw,zoomh,TYPE_INT_RGB);
    
    double length = 1/zoom;
    
    //まずはx方向に縮小します
    for(int y=0;y<h;y++){
        double l=0;
        double r=0,g=0,b=0;
        int impx = 0;
        for(int x=0;x<w;x++){
            int c = img.getRGB(x, y);
            if(l+1 >= length){
                r=(r+r(c)*(length-l))/length;
                g=(g+g(c)*(length-l))/length;
                b=(b+b(c)*(length-l))/length;
                int rgb = rgb((int)r,(int)g,(int)b);
                imp.setRGB(impx++, y, rgb);
                l = l+1-length;
                r = r(c)*l;
                g = g(c)*l;
                b = b(c)*l;
            }else{
                r = r(c);
                g = g(c);
                b = b(c);
                l+=1;
            }
        }
        //残っちゃった部分を処理
        if(impx<zoomw){
            r/=length;
            g/=length;
            b/=length;
            int rgb = rgb((int)r,(int)g,(int)b);
            imp.setRGB(impx, y, rgb);
        }
    }
    
    //同様にy方向にimpを縮小します
    img = imp;
    imp = ret;//コピペで済ますために入れ替えました (-_-)
    for(int x=0;x<zoomw;x++){
        double l=0;
        double r=0,g=0,b=0;
        int impy = 0;
        for(int y=0;y<h;y++){
            int c = img.getRGB(x, y);
            if(l+1 >= length){
                r=(r+r(c)*(length-l))/length;
                g=(g+g(c)*(length-l))/length;
                b=(b+b(c)*(length-l))/length;
                int rgb = rgb((int)r,(int)g,(int)b);
                imp.setRGB(x,impy++,rgb);
                l = l+1-length;
                r = r(c)*l;
                g = g(c)*l;
                b = b(c)*l;
            }else{
                r = r(c);
                g = g(c);
                b = b(c);
                l+=1;
            }
        }
        //残っちゃった部分を処理
        if(impy<zoomh){
            r/=length;
            g/=length;
            b/=length;
            int rgb = rgb((int)r,(int)g,(int)b);
            imp.setRGB(x, impy, rgb);
        }
    }
    return ret;
}