LoginSignup
1
1

More than 5 years have passed since last update.

WACULAdvent Calendar 2015Day 23
@mmaioe

Scalaで始める行列計算 - 第一回

Last updated at Posted at 2015-12-19

概要

機械学習、画像処理、コンピューターアニメーション等、いろいろな場面で使用される行列計算。その行列計算をあえて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")
  }
}
1
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
1