【C#】画像処理 ~ 幾何補正 / 擬似2次元アフィン変換編 ~

■ はじめに

https://blogs.yahoo.co.jp/dk521123/38069294.html
の続き。
今回は、「擬似2次元アフィン変換(共一次変換)」について、扱う。

■ サンプル

 * 以下で実装してみる
  + 補正方法 : 擬似2次元アフィン変換
  + 補間方法 : 最近傍補間(ニアレストネイバー)
 * ベースは、以下の関連記事の2次元アフィン変換と変わらない
  => 変わるのは、係数の数と等式のみ。
https://blogs.yahoo.co.jp/dk521123/38069294.html
 * 等式については、以下のサイトを参照
http://www.civil.kumamoto-u.ac.jp/matsu/geome.pdf

Form1.cs

using System;
using System.Collections.Generic;
using System.Drawing;
using System.Windows.Forms;

namespace SampleForm
{
  public partial class Form1 : Form
  {
    private const int DataNumberForPseudoAffineCoefficients = 8;

    private Bitmap originalBitmap = null;

    public Form1()
    {
      InitializeComponent();
    }

    private void Form1_Load(object sender, EventArgs e)
    {
      // 画像ファイルのImageオブジェクトを作成する
      this.originalBitmap = new Bitmap(@"C:\temp\20161215052204.gif");
      this.pictureBox1.Image = this.originalBitmap;
    }

    private void button1_Click(object sender, EventArgs e)
    {
      var inputs = new List<PointF>()
      {
        new PointF(0, 0),
        new PointF(this.originalBitmap.Width, this.originalBitmap.Height),
        new PointF(0, this.originalBitmap.Height),
        new PointF(this.originalBitmap.Width, 0),
      };
      var outputs = new List<PointF>()
      {
        new PointF(0, 0),
        new PointF((float)(this.originalBitmap.Width * 1.5), (float)(this.originalBitmap.Height * 1.5)),
        new PointF(0, (float)(this.originalBitmap.Height * 1.5)),
        new PointF((float)(this.originalBitmap.Width * 1.5), 0),
      };

      var pseudoAffineCoefficients = CalculatePseudoAffineCoefficients(inputs, outputs);
      if (pseudoAffineCoefficients == null)
      {
        MessageBox.Show("Error...", "Error", MessageBoxButtons.OK, MessageBoxIcon.Error);
        return;
      }

      // 逆変換
      var transferedBitmap = new Bitmap(this.originalBitmap.Width, this.originalBitmap.Height);
      for (int x = 0; x < this.originalBitmap.Width; ++x)
      {
        for (int y = 0; y < this.originalBitmap.Height; ++y)
        {
          Point correctedPoint = GetCorrectPointByPseudoAffineTransfer(pseudoAffineCoefficients, x, y);

          Color color;
          if (correctedPoint.X >= 0 && correctedPoint.X < this.originalBitmap.Width &&
            correctedPoint.Y >= 0 && correctedPoint.Y < this.originalBitmap.Height)
          {
            color = this.originalBitmap.GetPixel(correctedPoint.X, correctedPoint.Y);
          }
          else
          {
            color = Color.Black;
          }

          transferedBitmap.SetPixel(x, y, color);
        }
      }

      this.pictureBox2.Image = transferedBitmap;
    }

    /// <summary>
    /// 擬似アフィン変換係数を求める
    /// </summary>
    /// <param name="observedPoints">観測点</param>
    /// <param name="transferedPoints">変換点</param>
    /// <returns></returns>
    private double[] CalculatePseudoAffineCoefficients(
      List<PointF> observedPoints, List<PointF> transferedPoints)
    {
      if (observedPoints.Count < 4 || observedPoints.Count != transferedPoints.Count)
      {
        // 4点未満の際は、計算不可なので、nullを返す
        return null;
      }

      var matrixValues = new double[DataNumberForPseudoAffineCoefficients / 2][]
      {
        new double[] { 0.0, 0.0, 0.0, 0.0 },
        new double[] { 0.0, 0.0, 0.0, 0.0 },
        new double[] { 0.0, 0.0, 0.0, 0.0 },
        new double[] { 0.0, 0.0, 0.0, 0.0 }
      };
      var vectors = new double[DataNumberForPseudoAffineCoefficients] 
      {
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
      };

      int index = 0;
      foreach (var observedPoint in observedPoints)
      {
        var observedX = (double)observedPoint.X;
        var observedY = (double)observedPoint.Y;
        var transferedX = (double)transferedPoints[index].X;
        var transferedY = (double)transferedPoints[index].Y;
        var xy = observedX * observedY;
        var x2 = observedX * observedX;
        var y2 = observedY * observedY;
        var x2y = x2 * observedY;
        var xy2 = observedX * y2;
        var x2y2 = x2 * y2;

        matrixValues[0][0] += 1;
        matrixValues[0][1] += observedX;
        matrixValues[0][2] += observedY;
        matrixValues[0][3] += xy;
        matrixValues[1][0] += observedX;
        matrixValues[1][1] += x2;
        matrixValues[1][2] += xy;
        matrixValues[1][3] += x2y;
        matrixValues[2][0] += observedY;
        matrixValues[2][1] += xy;
        matrixValues[2][2] += y2;
        matrixValues[2][3] += xy2;
        matrixValues[3][0] += xy;
        matrixValues[3][1] += x2y;
        matrixValues[3][2] += xy2;
        matrixValues[3][3] += x2y2;

        vectors[0] += transferedX;
        vectors[1] += transferedX * observedX;
        vectors[2] += transferedX * observedY;
        vectors[3] += transferedX * observedX * observedY;
        vectors[4] += transferedY;
        vectors[5] += transferedY * observedX;
        vectors[6] += transferedY * observedY;
        vectors[7] += transferedY * observedX * observedY;

        index++;
      }

      // 逆行列を計算する
      if (!CalculateInverseMatrix(matrixValues, out double[][] inverseMatrixValues))
      {
        // 逆行列算出に失敗したため、null
        return null;
      }

      var pseudoAffineCoefficients = new double[DataNumberForPseudoAffineCoefficients];
      pseudoAffineCoefficients[0] =
        inverseMatrixValues[0][0] * vectors[0] +
        inverseMatrixValues[0][1] * vectors[1] +
        inverseMatrixValues[0][2] * vectors[2] +
        inverseMatrixValues[0][3] * vectors[3];
      pseudoAffineCoefficients[1] =
        inverseMatrixValues[1][0] * vectors[0] +
        inverseMatrixValues[1][1] * vectors[1] +
        inverseMatrixValues[1][2] * vectors[2] +
        inverseMatrixValues[1][3] * vectors[3];
      pseudoAffineCoefficients[2] =
        inverseMatrixValues[2][0] * vectors[0] +
        inverseMatrixValues[2][1] * vectors[1] +
        inverseMatrixValues[2][2] * vectors[2] +
        inverseMatrixValues[2][3] * vectors[3];
      pseudoAffineCoefficients[3] =
        inverseMatrixValues[3][0] * vectors[0] +
        inverseMatrixValues[3][1] * vectors[1] +
        inverseMatrixValues[3][2] * vectors[2] +
        inverseMatrixValues[3][3] * vectors[3];
      pseudoAffineCoefficients[4] =
        inverseMatrixValues[0][0] * vectors[4] +
        inverseMatrixValues[0][1] * vectors[5] +
        inverseMatrixValues[0][2] * vectors[6] +
        inverseMatrixValues[0][3] * vectors[7];
      pseudoAffineCoefficients[5] =
        inverseMatrixValues[1][0] * vectors[4] +
        inverseMatrixValues[1][1] * vectors[5] +
        inverseMatrixValues[1][2] * vectors[6] +
        inverseMatrixValues[1][3] * vectors[7];
      pseudoAffineCoefficients[6] =
        inverseMatrixValues[2][0] * vectors[4] +
        inverseMatrixValues[2][1] * vectors[5] +
        inverseMatrixValues[2][2] * vectors[6] +
        inverseMatrixValues[2][3] * vectors[7];
      pseudoAffineCoefficients[7] =
        inverseMatrixValues[3][0] * vectors[4] +
        inverseMatrixValues[3][1] * vectors[5] +
        inverseMatrixValues[3][2] * vectors[6] +
        inverseMatrixValues[3][3] * vectors[7];

      return pseudoAffineCoefficients;
    }

    // 逆行列算出。詳細は、以下の関連記事を参照
    // https://blogs.yahoo.co.jp/dk521123/38068366.html
    private static bool CalculateInverseMatrix(
      double[][] targetMatrixes, out double[][] outInverseMatrixes)
    {
      var dataNumber = targetMatrixes.GetLength(0);
      outInverseMatrixes = new double[dataNumber][];
      //単位行列を作る
      for (int i = 0; i < dataNumber; i++)
      {
        outInverseMatrixes[i] = new double[dataNumber];
        for (int j = 0; j < dataNumber; j++)
        {
          outInverseMatrixes[i][j] = (i == j) ? 1.0f : 0.0f;
        }
      }

      // 掃き出し法
      for (int i = 0; i < dataNumber; i++)
      {
        if (targetMatrixes[i][i] == 0)
        {
          return false;
        }
        var tempValue = 1.0 / targetMatrixes[i][i];
        for (int j = 0; j < dataNumber; j++)
        {
          targetMatrixes[i][j] *= tempValue;
          outInverseMatrixes[i][j] *= tempValue;
        }
        for (int j = 0; j < dataNumber; j++)
        {
          if (i != j)
          {
            tempValue = targetMatrixes[j][i];
            for (int k = 0; k < dataNumber; k++)
            {
              targetMatrixes[j][k] -= targetMatrixes[i][k] * tempValue;
              outInverseMatrixes[j][k] -= outInverseMatrixes[i][k] * tempValue;
            }
          }
        }
      }
      return true;
    }


    /// <summary>
    /// 擬似アフィン変換による補正座標を返す
    /// </summary>
    /// <param name="pseudoAffineCoefficients">擬似アフィン係数</param>
    /// <param name="targetPoint">対象座標</param>
    /// <returns>擬似アフィン変換による補正座標</returns>
    private Point GetCorrectPointByPseudoAffineTransfer(
      double[] pseudoAffineCoefficients, int x, int y)
    {
      // u = a0xy + a1x + a2y + a4
      var correctedX = 
        pseudoAffineCoefficients[0] * x * y + 
        pseudoAffineCoefficients[1] * x + 
        pseudoAffineCoefficients[2] * y + 
        pseudoAffineCoefficients[3];
      // v = b0xy + b1x + b2y + b4
      var correctedY = 
        pseudoAffineCoefficients[4] * x * y + 
        pseudoAffineCoefficients[5] * x + 
        pseudoAffineCoefficients[6] * y + 
        pseudoAffineCoefficients[7];

      // 座標にするので四捨五入
      return new Point(
        (int)Math.Round(correctedX, MidpointRounding.AwayFromZero),
        (int)Math.Round(correctedY, MidpointRounding.AwayFromZero));
    }
  }
}

関連記事

C#】画像処理 ~ 幾何補正 / 2次元アフィン変換編 ~

https://blogs.yahoo.co.jp/dk521123/38069294.html

画像処理 ~ 幾何補正について ~

https://blogs.yahoo.co.jp/dk521123/38073336.html