概要

機械学習、画像処理、コンピューターアニメーション等、いろいろな場面で使用される行列計算。その行列計算をあえてScala(一部、Javaです)でやりたというというのが今回の趣旨です。今回は、数ある行列計算のオープンソースの中からLAPACKを使ってみました。

LAPACK

ちなみに、LAPACKは、大規模な疎行列の計算には向いてません（＊私個人的に思っただけなので、もし間違えてたら、ご指摘お願いします。）

LAPACKの依存関係

SBTを使っていれば、以下をbuild.sbtに追加するだけで、使用できます(*注意: 適宜適切なバージョンをしてしてください)

"com.googlecode.netlib-java" % "netlib-java" % "1.0"

活用例

今回は、画像処理で基本の射影変換をLAPACKを使って、求めてみたいと思います。

以上のGitにサンプルプログラムを乗せました。IntelliJを使用している方は、cloneして、Sbtプロジェクトとして、Importすれば使えます。

ホモグラフィ行列

ホモグラフィ行列に関する記事は、ネット上に多数存在し、多くの応用例があります。以下のブログなどでも、その応用例がみれます。

http://blog.livedoor.jp/ryo_ogawa/archives/3109950.html
http://www.cuspy.org/diary/2012-06-29

4点からホモグラフィ行列計算

ホモグラフィ行列の計算は、４つの点があれば計算できます。このことは、多くのネット上の文献でもみれるので、是非検索してみてください。

Homography.scala

package mmaioe.com.cameracalibration

import mmaioe.com.linearsystem.LinearSolver

/**
 * Created by maoito on 12/19/15.
 */
object Homography {
  /**
   *
   * p1_* belongs to the samge image
   * p2_* belongs to the same image
   *
   * p1_i and p2_i should be a pair of corresponding points
   * @param p1_1
   * @param p1_2
   * @param p1_3
   * @param p1_4
   * @param p2_1
   * @param p2_2
   * @param p2_3
   * @param p2_4
   */
  def solve(
             p1_1: {val x:Double; val y:Double; },
             p1_2: {val x:Double; val y:Double; },
             p1_3: {val x:Double; val y:Double; },
             p1_4: {val x:Double; val y:Double; },
             p2_1: {val x:Double; val y:Double; },
             p2_2: {val x:Double; val y:Double; },
             p2_3: {val x:Double; val y:Double; },
             p2_4: {val x:Double; val y:Double; }
             ):List[List[Double]]={
    val homography = LinearSolver.solveWithSquareMatrix(
      Array(
        Array(p1_1.x, p1_1.y, 1, 0,      0,      0, -p1_1.x*p2_1.x, -p1_1.y*p2_1.x),
        Array(0,      0,      0, p1_1.x, p1_1.y, 1, -p1_1.x*p2_1.y, -p1_1.y*p2_1.y),
        Array(p1_2.x, p1_2.y, 1, 0,      0,      0, -p1_2.x*p2_2.x, -p1_2.y*p2_2.x),
        Array(0,      0,      0, p1_2.x, p1_2.y, 1, -p1_2.x*p2_2.y, -p1_2.y*p2_2.y),
        Array(p1_3.x, p1_3.y, 1, 0,      0,      0, -p1_3.x*p2_3.x, -p1_3.y*p2_3.x),
        Array(0,      0,      0, p1_3.x, p1_3.y, 1, -p1_3.x*p2_3.y, -p1_3.y*p2_3.y),
        Array(p1_4.x, p1_4.y, 1, 0,      0,      0, -p1_4.x*p2_4.x, -p1_4.y*p2_4.x),
        Array(0,      0,      0, p1_4.x, p1_4.y, 1, -p1_4.x*p2_4.y, -p1_4.y*p2_4.y)
      )
      ,
      Array(p2_1.x,p2_1.y,p2_2.x,p2_2.y,p2_3.x,p2_3.y,p2_4.x,p2_4.y)
    )

    return List[List[Double]](
      List(homography(0),homography(1),homography(2)),
      List(homography(3),homography(4),homography(5)),
      List(homography(6),homography(7),1)
    )
  }

  def multiply(homography:List[List[Double]], point: List[Double]): List[Int] ={
    var threePoint:collection.mutable.ListBuffer[Double] = collection.mutable.ListBuffer[Double]()++point
    threePoint += 1

    val new_point:collection.mutable.ListBuffer[Int] = collection.mutable.ListBuffer[Int]()
    (0 until 2).foreach{
      i =>
        var sum = 0.0
        (0 until 3).foreach{
          j =>
            sum += homography(i)(j) * threePoint(j);
        }
        if(sum - sum.toInt > 0.5) sum = sum.toInt+1
        else                      sum = sum.toInt

        new_point += sum.toInt
    }

    return new_point.toList
  }
}

Image.java

package mmaioe.com.image;

import javax.imageio.ImageIO;
import java.awt.*;
import java.awt.image.BufferedImage;
import java.io.File;
import java.io.FileInputStream;
import java.io.IOException;
import java.io.InputStream;

/**
 * Created by maoito on 12/19/15.
 */
public class Image {
    public BufferedImage image;
    public String fileName;


    public Image(String fileName){
        InputStream is = null;
        this.fileName = fileName;
        try {
            is = new FileInputStream(fileName);
            image = ImageIO.read(is);
        } catch (IOException e) {
            throw new RuntimeException("failed to read "+fileName);
        } finally {
            if (is != null) try { is.close(); } catch (IOException e) {
                throw  new RuntimeException("failed to close FileInputStream");
            }
        }
    }
    public Image(int width, int height){
        this.image = new BufferedImage(width,height,BufferedImage.TYPE_INT_RGB);
    }

    public Color getRGBColor(int x,int y){
        return new Color(image.getRGB(x, y));
    }

    public void setColor(int x,int y, Color c){
        image.setRGB(x, y, c.getRGB());
    }

    public int width(){
        return image.getWidth();
    }
    public int height(){
        return image.getHeight();
    }

    public void writeImage(String fileName){
        try{
            ImageIO.write(image, "jpg", new File(fileName));
        }catch(Exception e){
            throw new RuntimeException("failed to write image to "+fileName);
        }
    }
}

LinearSolver.java

package mmaioe.com.linearsystem;

import org.netlib.lapack.Dgesv;
import org.netlib.util.intW;

/**
 * Created by maoito on 12/19/15.
 */
public class LinearSolver {
    /**
     *
     * @param _A
     * @param _B
     * @return
     */
    public static double[] solveWithSquareMatrix(double[][] _A, double[] _B) {
        int N = _A.length;
        int NRHS = 1;
        int LDA = _A.length;
        int LDB = _B.length;

        double[] A = new double[_A.length*_A.length];
        double[] B = _B;

        intW info = new intW(0);

        int[] IPIV = new int[_A.length];
        for(int i=0;i<IPIV.length;i++) IPIV[i] = 0;

        for(int i=0;i<_A.length;i++){
            for(int j=0;j<_A[i].length;j++){
                A[i+j*_A.length] = _A[i][j];
            }
        }

        //Solve Ah = B
        Dgesv.dgesv(N, NRHS, A, 0, LDA, IPIV, 0, B, 0, LDB, info);

        return B;
    }

}

適用例

変換した画像を見せられないのは残念ですが、以下に使用例のプログラムを貼ります。今回のプログラムでは、単純に射影変換したのみで、補間などはおこなっていません。

Test.scala


import mmaioe.com.cameracalibration.Homography
import mmaioe.com.image.Image

/**
 * Created by maoito on 12/19/15.
 */
object Test {
  def main(args: Array[String]) {
    val image = new Image("/Users/XXXX/DSC01018.JPG")

    println(image.height()+","+image.width())


    val homography = Homography.solve(
      new {val x:Double = 0; val y:Double=0; },
      new {val x:Double = 0; val y:Double=5; },
      new {val x:Double = 5; val y:Double=0; },
      new {val x:Double = 5; val y:Double=5; },
      new {val x:Double = 0; val y:Double=0; },
      new {val x:Double = 2; val y:Double=6; },
      new {val x:Double = 4; val y:Double= -1; },
      new {val x:Double = 6; val y:Double=4; }
    )

    val newImage = new Image(image.width()*2,image.height()*2)

    val base_x=900
    val base_y=900
    (0 until image.width()).foreach{
      x =>
        (0 until image.height()).foreach{
          y =>
            val new_point = Homography.multiply(homography,List(x,y))

            if(0 <= new_point(0)+base_x && new_point(0)+base_x < newImage.width() && 0 <= new_point(1)+base_y && new_point(1)+base_y < newImage.height()){
              newImage.setColor(new_point(0)+base_x,new_point(1)+base_y,image.getRGBColor(x,y))
            }
        }
    }

    newImage.writeImage("homographied.png")
  }
}