ハフ変換

今回はハフ変換について。
ハフ変換の説明はウィキペディアとかまあそういうので調べると良いんじゃないかと思います。
私は研究室にあった本を使いましたがね。

参考にしたサイト
http://codezine.jp/a/article.aspx?aid=153

実行結果はのせるのが面倒なので、ソースコードだけ貼ります。

import java.util.Date;

/**
 * 2次元配列に対するHought変換を行うための手段を提供するクラスです。
 * @author MineAP
 */
public class HoughTransform {

    //直線検出用
    private final double PI = Math.PI;
    private final int MAX_THETA = 360;
    private final int MAX_RHO = 720;
    private final int MAX_BRIGHTNESS = 255; 
    
    //円検出用
    private final int XMAX;
    private final int YMAX;
    private final int RMAX;
    private final int THETA_MAX = 1024;
    private final float PIK = (float) Math.PI / THETA_MAX;

    //三角関数テーブル(サイン)
    float[] sn;
    //三角関数テーブル(コサイン)
    float[] cs;
    //半径計算用斜線長テーブル
    short[][] diagonal;
    
    
    private int mode = 0;
    
    private int x_size1;
    private int y_size1;
    private int x_size2;
    private int y_size2;
    private int x_size3, y_size3, z_size3;
    
    private int[][] image1;
    private int[][] image2d;    
    private int[][][] image3d;
    
    /**
     * コンストラクタ
     * 引数で得た2次元配列を使って初期化します。
     * @param image1 
     */
    public HoughTransform(int[][] image1){
        this.image1 = image1;
        this.x_size1 = image1.length;
        this.y_size1 = image1[0].length;
        
        //円検出で使用
        this.YMAX = this.y_size1;
        this.XMAX = this.x_size1;
        this.RMAX = this.YMAX/4;
        
    }
    
    /**
     * モードをセットします。
     * モードはデフォルトでは0です。
     * mode = 0 直線の検出
     * mode = 1 円の検出
     * @param mode 
     */
    public void setMode(int mode){
        this.mode = mode;
    }
    
    //////////////////// 直線の検出用 ////////////////////////////
    
    /**
     * image2[x][y]上の2点(x1,y1)、(x2,y2)間の直線上の配列要素の値をインクリメントする。
     * 
     */
    private void drawStraightLine(int x1, int y1, int x2, int y2) {
        
        double distance,step,t;
        int x,y;
        
        //2点間の距離
        distance = Math.sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
        //発生させる内分点の数
        step = 1.0/(1.5*distance);
        
        t = 0.0;
        while(t<1.0){
            x = (int)(t*x2+(1.0-t)*x1);
            y = (int)(t*y2+(1.0-t)*y1);
            if(x >= 0 && x < x_size2 && y >= 0 && y < y_size2 && image2d[x][y]<MAX_BRIGHTNESS){
                image2d[x][y]++;
                //System.out.println("inc"+x+","+y);
            }
            t = t+step;
        }
    }
    
    /**
     * ρを計算するメソッド
     */
    private int calcRho(double px, double py, double r_mg, int theta) {
        return (int)(r_mg*(px*Math.cos(theta*PI/MAX_THETA)+py*Math.sin(theta*PI/MAX_THETA)) + MAX_RHO/2.0);
    }
    
    /**
     * ρ=x cosθ + y sinθ の曲線を配列image2[rho][theta]に描く
     */
    private void drawCurve(double px, double py){
        //ρ方向の最大値[画素]
        double rho_max = 1.0/1.5*x_size1 + 1.0/1.5*y_size1;
        //ρ方向の係数
        double rho_c = MAX_RHO/2.0/rho_max;
        
        for(int theta = 0; theta < MAX_THETA; theta++){
            //System.out.println("drawStraightLine()");
            drawStraightLine(theta, calcRho(px, py, rho_c, theta), theta+1, calcRho(px, py, rho_c, theta+1));
        }
    }
    
    //////////////////// 直線の検出用 ここまで ////////////////////////////
    
    //////////////////// 円の検出用 ////////////////////////////
    
    /**
     * 円を検出するためのハフ変換を行う
     */
    private void houghTransformForCircle(){
        
        sn = new float[THETA_MAX];
        cs = new float[THETA_MAX];
        diagonal = new short[YMAX][XMAX];
        
        //三角関数テーブルを作成
        for (int i = 0; i < THETA_MAX; i++) {
            sn[i] = (float) Math.sin(PIK*i);
            cs[i] = (float) Math.cos(PIK*i);
        }
        
        //斜線長テーブルを作成
        for (int y = 0; y < YMAX; y++) {
            for (int x = 0; x < XMAX; x++) {
                diagonal[y][x] = (short) (Math.sqrt(y*y+x*x)+0.5);
            }
        }
        
        int centerX,centerY,distX,distY,radius;
        //円検出用頻度カウンタ
        image3d = new int[YMAX][XMAX][RMAX];
        int count = 0;
        //y座標を捜索
        for (int y = 0; y < YMAX; y++) {
            //x座標を捜索
            for (int x = 0; x < XMAX; x++) {
                //元画像の明度が最大の点(0xff?)を見つけた場合の処理
                if (image1[y][x] == MAX_BRIGHTNESS) {
                    //明度が最大の点を通る円を使ってimage3dに写像する
                    //y座標について捜索する
                    for (centerY = 0; centerY < YMAX; centerY++) {
                        distY = Math.abs(y - centerY);
                        //半径がimage3dの空間に収まらなかったら中断してやり直し
                        if (distY > RMAX) {
                            continue;
                        }
                        //x座標について捜索する
                        for (centerX = 0; centerX < XMAX; centerX++) {
                            distX = Math.abs(x - centerX);
                            radius = diagonal[distY][distX];
                            //対角線の長さがimage3dの空間に収まらなかったら中断してやり直し
                            if (radius >= RMAX) {
                                continue;
                            }
                            if (image3d[centerY][centerX][radius] < MAX_BRIGHTNESS) {
                                image3d[centerY][centerX][radius]++;
                                count++;
                            }
                        }
                    }
                }
            }
        }
        
        System.out.println(count + "個の座標を写像しました。");
    }
    
    //////////////////// 円の検出用 ここまで ////////////////////////////
    
    
    /**
     * コンストラクタで初期化されたimage1[x][y]をHough変換し、
     * 結果をθρ平面image2[theta][rho]に代入する。
     * ただし、θρ平面の横軸・縦軸の大きさはMAX_THETA,MAX_RHOで予め決められている。
     * また、MAX_BRIGHTNESS以上の累積は考慮しない。
     * 
     * @return image1[][]をHough変換した値
     */
    public int[][] getHoughTransformedArray() {
        
        if(mode != 0){
            System.out.println("直線検出モードではないため、ハフ変換できません。: mode = " + this.mode);
            return null;
        }
        
        x_size2 = MAX_THETA;
        y_size2 = MAX_RHO;
        image2d = new int[x_size2][y_size2];
        int count = 0;
        
        //配列の初期化
        for(int y=0; y < y_size2; y++){
            for(int x=0; x < x_size2; x++){
                image2d[x][y] = 0;
            }
        }
        
        //処理の実施
        System.out.println("ハフ変換を実施します。");
        for(int y=0; y<y_size1; y++){
            for(int x=0; x<x_size1; x++){
                if(image1[x][y] == MAX_BRIGHTNESS){
                    //System.out.println("("+x+"," +y+") drawCurve()");
                    drawCurve(x, y);
                    count++;
                }
            }
            //System.out.println("");
        }
        System.out.println(count + "個の座標をθρ平面に写像しました。");
        
        return image2d;
    }
    
    /**
     * コンストラクタで初期化されたimage1[x][y]をHough変換し、
     * 結果をxxyyzz平面image3d[xx][yy][r]に代入する。
     * MAX_BRIGHTNESS以上の累積は考慮しない。
     * 
     * @return image1[][]をHough変換した値
     */
    public int[][][] getHoughTransformedArray3D() {
        
        if(mode != 1){
            System.out.println("円検出モードではないため、ハフ変換できません。: mode = " + this.mode);
            return null;
        }
        
        x_size3 = 0;
        y_size3 = 0;
        z_size3 = 0;
        image3d = new int[x_size3][y_size3][z_size3];
        System.out.println("ハフ変換開始...");
        long time = (new Date()).getTime();
        //処理の実施
        this.houghTransformForCircle();
        
        System.out.println("処理にかかった時間[ms]:" + ((new Date()).getTime() - time));
        
        return image3d;
        
    }
    
}


以上、ハフ変換をして2次元配列に格納して返すクラス。

次、逆ハフ変換のためのクラス。

import java.awt.Color;
import java.awt.Graphics;
import java.awt.image.BufferedImage;
import java.util.Date;

/**
 *
 * @author MineAP
 */
public class InvHoughTransform {

    //直線検出用
    private final double PI = Math.PI;
    private final int MAX_THETA = 360;
    private final int MAX_RHO = 720;
    private final int MAX_BRIGHTNESS = 255; 
    
    //円検出用
    private int XMAX;
    private int YMAX;
    private int RMAX;
    private int THETA_MAX = 1024;
    private float PIK = (float) Math.PI / THETA_MAX;
    
    
    private int x_size1;
    private int y_size1;
    private int x_size2;
    private int y_size2;
    private int x_size3,y_size3,r_size;
    
    private int[][] image1;
    private int[][] image2;  
    private int[][][] image3d;
    
    /**
     * コンストラクタ。
     * 渡された引数を使って画像を初期化します。
     * @param image1 
     */
    public InvHoughTransform(int[][] image1){
        this.image1 = image1;
        this.x_size1 = image1.length;
        this.y_size1 = image1[0].length;
    }
    
    /**
     * コンストラクタ。
     * 渡された引数(円をハフ変換した3次元配列)をつかって配列を初期化します。
     * @param image3d 
     */
    public InvHoughTransform(int[][][] image3d){
        this.image3d = image3d;
        this.x_size3 = image3d.length;
        this.y_size3 = image3d[0].length;
        this.r_size = image3d[0][0].length;
    }
    
    /**
     * θρ平面(image1)上の点(xs,ys)を画像上の直線に逆変化し、その結果をimage2[x][y]に描きます。
     */
    private void extractLine(int xs, int ys) {
        double theta,rho,rho_max,sin,cos;
        
        theta = 180.0*xs/MAX_THETA;
        sin = Math.sin(theta/180.0*PI);
        cos = Math.cos(theta/180.0*PI);
        rho_max = 1.0/1.5*x_size2 + 1.0/1.5*y_size2;
        rho = (ys - MAX_RHO/2.0) * rho_max/MAX_RHO * 2.0;
        /* rho = xcos+ysin -> y=rho/sin-xcos/sin */
        if(45.0<=theta && theta <= 135.0) {
            //System.out.println("45°<=θ<=135°");
            for(int x=0; x<x_size2;x++){
                int y=(int)((rho-x*cos)/sin);
                if(0<=y && y<y_size2) {
                    //System.out.println("set");
                    image2[x][y] = MAX_BRIGHTNESS;
                }
            }
        }else{
            //System.out.println("0°<=θ<45°,135°<θ<360°");
            for(int y=0;y<y_size2;y++){
                int x=(int)((rho-y*sin)/cos);
                if(0 <= x && x<x_size2) {
                    //System.out.println("set");
                    image2[x][y] = MAX_BRIGHTNESS;
                }
            }
        }
    }
    
    /**
     * image1をθρ平面と見なし、閾値より上の点の画像空間に逆変換して画像中に直線を描きます。
     * 出力はimage2に行われます。
     */
    private void inverseHoughTransform(int threshold){
        
        int counter;
        int x,y;
        //String str;
        
        System.out.println("しきい値:" + threshold);
        counter = 0;
        for (y = 0; y < MAX_RHO; y++) {
            for (x = 0; x < MAX_THETA; x++) {
                if (image1[x][y] >= threshold) {
                    counter++;
                    extractLine(x, y);
                }
            }
        }
        System.out.println("検出された直線数:" + counter);
        
    }
    
    /**
     * データが正しいかどうかをチェックします。
     * 正しくない場合false,正しい場合trueを返します。
     */
    private boolean checkData() {
       
        if (x_size1 == MAX_THETA || y_size1 == MAX_RHO) {
            return true;
        }

        return false;
    }
    
    /**
     * 指定された引数(閾値、出力先画像のサイズ)を使って予め初期化された画像を逆ハフ変換し、
     * その画像を返します。
     * @param threshold 閾値
     * @param x_size 出力画像の幅
     * @param y_size 出力画像の縦
     * @return 逆ハフ変換後の画像(配列)
     */
    public int[][] getInverseHoughTransformedArray(int threshold, int x_size, int y_size) {
        if(this.checkData()){
            this.image2 = new int[x_size][y_size];
            this.x_size2 = x_size;
            this.y_size2 = y_size;
            for(int x=0; x<x_size; x++) {
                for(int y=0; y<y_size; y++){
                    image2[x][y] = 0;
                }
            }
            System.out.println("逆ハフ変換を実行します。");
            
            this.inverseHoughTransform(threshold);
            
            return this.image2;
        }else{
            System.out.println("初期化されている画像はサポートされていません。");
            return null;
        }
    }

    /**
     * 円を検出するために行ったハフ変換の結果を格納した3次元配列に対して逆ハフ変換を行い、その結果を
     * BufferedImageに格納して返します。
     * 検出された図形は白い背景に黒い線で描かれます。
     * @param threshold 
     * @param x_size 
     * @param y_size 
     * @return 
     */
    public BufferedImage getInverseHoughTransformedImageForCircle(int threshold, int x_size, int y_size) {
        
        int X0=0,Y0=0,X1=x_size,Y1=y_size;
        
        System.out.println("逆ハフ変換開始...");
        long time = (new Date()).getTime();
        
        this.XMAX = x_size;
        this.YMAX = y_size;
        this.RMAX = this.YMAX / 4;
        
        BufferedImage img = new BufferedImage(XMAX, YMAX, BufferedImage.TYPE_3BYTE_BGR);
        
        for(int x=0; x<X1; x++){
            for(int y=0; y<Y1; y++){
                img.setRGB(x, y, 0xffffffff);
            }
        }
        
        Graphics g = img.getGraphics();
        
        g.setColor(Color.BLACK);
        
        //繰り返しを終了させるフラグ
        int end_flag;
        //検出された直線または円の個数カウンタ
        int count;
        //
        int counter1_max;   

        int centerX_max = 0;
        int centerY_max = 0;
        int radius_max = 0;
        
        end_flag = 0;
        count = 0;

        /*
         * 大きい円から順に検出していくアルゴリズムになっている。 
         * 果たして指先を検出したいと思った時に正しく働くか。
         */
        do {
            //r^2 = (x - centerX)^2 + (y - centerY)^2
            counter1_max = 0;
            //counter1が最大になるcenterX_max、centerY_maxとradius_maxを求める
            for (int centerY = 0; centerY < YMAX; centerY++) {
                for (int centerX = 0; centerX < XMAX; centerX++) {
                    for (int radius = 0; radius < RMAX; radius++) {
                        if (image3d[centerX][centerY][radius] > counter1_max) {
                            counter1_max = image3d[centerY][centerX][radius];
                            //引数で指定された大きさになったら検出を終了
                            if (counter1_max <= threshold) {
                                end_flag = 1;
                                
                            } else {
                                end_flag = 0;
                            }
                            centerY_max = centerY;
                            centerX_max = centerX;
                            radius_max = radius;
                        }
                    }
                }
            }
            System.out.println("counter1_max = " + counter1_max);
            System.out.println("\tcenterY:"+centerY_max);
            System.out.println("\tcenterX:"+centerX_max);
            System.out.println("\tr:"+radius_max);
            
            //検出された円の近くにある円を消す
            for (int k = -5; k <= 5; k++) {
                if (centerY_max + k >= YMAX || centerY_max + k < 0) {
                    continue;
                }
                for (int j = -5; j <= 5; j++) {
                    if (centerX_max + j >= XMAX || centerX_max + j < 0) {
                        continue;
                    }
                    for (int i = -5; i <= 5; i++) {
                        if(radius_max+i>=RMAX || radius_max+i<0) continue;
                        if (radius_max + i < 0) {
                            continue;
                        }
                        try{
                            image3d[centerX_max + k][centerY_max + j][radius_max + i] = 0;
                        }catch(ArrayIndexOutOfBoundsException e){
                            /*
                            System.out.println("x = " + centerX_max);
                            System.out.println("y = " + centerY_max);
                            System.out.println("radius = " + radius_max);
                            e.printStackTrace();
                            System.exit(0);
                             */
                        }
                    }
                }
            }

            //検出した円の描画
            //if(counter1_max > threshold){
                g.drawOval(centerX_max - radius_max, centerY_max - radius_max, radius_max * 2, radius_max * 2);
                count ++;
            //}
            //大きさが100ピクセル以下か、円が5個検出されたら終了
        } while (end_flag == 0 && count < 10);
        
        g.dispose();
        
        System.out.println("検出された円:"+count);
        
        System.out.println("処理にかかった時間[ms]:" + ((new Date()).getTime() - time));
        
        return img;
    }
    
}

以上、逆ハフ変換してBufferedImageを返すクラス。
円検出用のハフ変換を半円に応用して指先を検出しようと思ったんだけど、全然リアルタイム性が出せなくて断念しました。
おかげでこっち方面にあまり興味がなくなってきた。

Processingの本を読んでグラフィカルな方向に戻ろうかな、というところ。