読者です 読者をやめる 読者になる 読者になる

PILのImage.transformの使い方 (PERSPECTIVE)

Python

PILのImage.transformの使い方 (QUAD, MESH) - umejanのブログのつづき。

Image.PERSPECTIVEの変形

遠近感のある形に変形する。

Image.PERSPECTIVEのdata指定方法

data : 8 tupples (a, b, c, d, e, f, g, h)

(x0, y0) = ((ax+by+c)/(gx+hy+1), (dx+ey+f)/(gx+hy+1))

で元の画像のピクセルで見るとのこと。


自力でa,b,c,d,e,f,hを計算するのは難しい。

変換元座標から変換先座標に変換する a,b,c,d,e,f を求める処理がネットにあったのでそれを利用する。 ( python - How does perspective transformation work in PIL? - Stack Overflow ) ※中身わからず。方程式を解いている?


遠近感のある画像から正面から見たように変換する場合は、変換前と変換後の4座標は決められる。

しかし、正面の画像を遠近感のある画像に変換する場合は、本来であれば変換後の4点を計算する必要があるだろう。

そこで、サンプルでは、視野変換、透視投影変換で計算した座標を求める処理も用意した。 詳しくはサンプル参照のこと。

視野変換は、 http://marina.sys.wakayama-u.ac.jp/~tokoi/?date=20090902 を参考にした。

透視投影変換は、 http://marupeke296.com/DXG_No70_perspective.html を参考にした。

Image.PERSPECTIVEのサンプル

直接座標を指定して変形したサンプルと
視野変換、透視投影変換で計算して変換後の座標を求めて変形したサンプル。

# -*- coding: utf-8 -*-
from __future__ import print_function
from PIL import Image
from PIL import ImageDraw
import numpy as np
import math
import argparse

parser = argparse.ArgumentParser(description="Image.transformのサンプル")

parser.add_argument("--output","-o", required=True, type=str
                   , help="出力ファイル")

args = parser.parse_args()

# 変換元の画像をロード
def load_image2() :
  #im = Image.open("xxxx.png")

  # サンプルとして、82x66の画像を作成。グレースケール256階調
  # 16ドットごとに格子を描画
  im = Image.new("L", (82, 65), color=0)
  draw = ImageDraw.Draw(im)
  for i in range(6) :
    draw.line([(i*16, 0), (i*16, 65)], fill=255, width=2)
  for i in range(5) :
    draw.line([(0, i*16), (82, i*16)], fill=255, width=2)
  
  # 四角を左上の格子の中に描画
  draw.rectangle([6, 6, 12, 12], fill=255)
  
  del draw
  #im.save("moto2.png")

  return im


im = load_image2()

#
# Image.PERSPECTIVE
# 
# data : 8 tupples (a, b, c, d, e, f, g, h)
# 
# (x0, y0) =
#   ((ax+by+c)/(gx+hy+1), (dx+ey+f)/(gx+hy+1))
# 
# (1) 変換元から変換先に変換するパラメータ取得処理
# 変換元座標から変換先座標に変換する a,b,c,d,e,f を求める処理が
# ネットにあったのでそれを利用する。
#
# ただ、長方形を遠近感のある四角形に変形するには、
# 事前に変形後の座標を求める必要がある。
# 四角形の4点はでたらめに指定しても、
# 正しく遠近感のある四角形にならない気がする。
#
# (2) 変換先座標を求める処理
# 3次元(x,y,z)でのカメラの位置、透視投影スクリーンの位置、注目点等を与えて、
# カメラの位置を移動させることで変換先座標を求めることにする。
#
# 変換前のカメラの位置  :
#     (0, 0, 0)
# 変換前の画像の位置 :
#     カメラの位置から真正面にあるものとし、距離はnZ。z座標は、-nZ。
#     左上、左下、右下、右上のそれぞれの (x, y)は、下記とする。
#       (-width/2, -height/2)
#       ( width/2, -height/2)
#       ( width/2,  height/2)
#       (-width/2,  height/2)
# 変換前の注目点(視線とスクリーンの交差する座標) :
#     画像中央の位置とする。(0, 0, -nZ)
# 変換前の透視投影スクリーン :
#     変換前の画像と一致しているものとする。
#     べたっとくっついてる状態で枠も画像サイズと同じ。
# 
# 上記の設定で、以下のカメラの移動を行い、座標変換(x1,y1,z1系)を行う。
# - カメラ位置の移動。(0, 0, 0) -> (ex, ey, ez)。
# - カメラの上方向のベクトルを指定。(ux, uy, uz)。通常、(0, 1, 0)。
# - 目標点の移動。(tx, ty, tz)。画像中央にする場合は、(0, 0, -nZ)。
# (補足) この変換を視野変換というらしい。
# (参考) : http://marina.sys.wakayama-u.ac.jp/~tokoi/?date=20090902
# 
# 視野変換後、カメラ位置から目標点(z軸)方向に透視投影スクリーンがあるものとして、
# その透視投影スクリーンに投影された画像の4隅の座標を取得する。(x2,y2,z2系)
# なお、描画空間の奥行が必要だが適当に以下にしておく。
#   nZ + max(width,height) * 2 + カメラ位置のz移動量としておく。
#
# (補足) この変換を透視投影変換というらしい。
# (参考) : http://marupeke296.com/DXG_No70_perspective.html
#
# この透視投影変換で得られた座標を(1)の変換後の座標として使用する。
#


#######################################################
# 関数
#######################################################

##############################################
# a b c d e f g の係数を取得する関数
# pa : 変換先の座標系
# pb : 現在の座標系
# コピー元:
# http://stackoverflow.com/questions/14177744/how-does-perspective-transformation-work-in-pil
##############################################
def find_coeffs(pa, pb):
    matrix = []
    for p1, p2 in zip(pa, pb):
        matrix.append([p1[0], p1[1], 1, 0, 0, 0, -p2[0]*p1[0], -p2[0]*p1[1]])
        matrix.append([0, 0, 0, p1[0], p1[1], 1, -p2[1]*p1[0], -p2[1]*p1[1]])

    A = np.matrix(matrix, dtype=np.float)
    B = np.array(pb).reshape(8)

    res = np.dot(np.linalg.inv(A.T * A) * A.T, B)
    ret = np.array(res).reshape(8)
    wk = []
    for v in ret :
      wk.append(round(v, 3))
    #print("coeffs", wk)
    (alpha,beta) = (wk[6]+1, wk[7]+1)
    (x0,y0) = (wk[2], wk[5])
    (x1,y1) = (round((wk[0]+x0)/alpha,3), round((wk[3]+y0)/alpha,3))
    (x2,y2) = (round((wk[1]+x0)/beta,3), round((wk[4]+y0)/beta,3))
    #print("alpha,beta", wk[6]+1,wk[7]+1, "x0,y0", wk[2],wk[5])
    #print("x1,y1", x1, y1, "x2,y2", x2,y2)
    #print("pa", pa)
    #print("pb", pb)
    
    return ret

# 視野変換
def getSiyaHenkanMatrix(ev, tv, uv) :
  k1 = np.sqrt(np.power(ev[0]-tv[0],2) +
               np.power(ev[1]-tv[1],2) +
               np.power(ev[2]-tv[2],2))
  
  zv = ((ev[0]-tv[0])/k1, (ev[1]-tv[1])/k1, (ev[2]-tv[2])/k1)
  #print("zv", zv)
  
  wv = np.cross(uv, zv)
  k2 = np.sqrt(np.power(wv[0], 2) +
               np.power(wv[1], 2) +
               np.power(wv[2], 2))
  xv = wv / k2
  yv = np.cross(zv, xv)
  
  A = np.matrix([
    [xv[0], xv[1], xv[2], -ev[0]*xv[0]-ev[1]*xv[1]-ev[2]*xv[2]],
    [yv[0], yv[1], yv[2], -ev[0]*yv[0]-ev[1]*yv[1]-ev[2]*yv[2]],
    [zv[0], zv[1], zv[2], -ev[0]*zv[0]-ev[1]*zv[1]-ev[2]*zv[2]],
    [0.0  , 0.0  , 0.0  , 1.0                                 ]
  ], dtype=float)
  return A.T

# 透視投影変換
#  (x y z 1)*M = (x' y' z' 1)
#  P = (x'/z y'/z z' 1)
def getTosiHenkanMatrix(size, nZ, ev) :
  (w, h) = (float(size[0]), float(size[1]))
  fZ = float(nZ) + max(w, h) * 2.0 + float(ev[2])
  cot_h_theta = 2.0 * float(nZ) / w / (h/w)
  print("cot_h_theta:", cot_h_theta, "nZ:", nZ, "fZ:", fZ)
  B = np.matrix([
    [h/w * cot_h_theta, 0.0        , 0.0        , 0.0],
    [0.0              , cot_h_theta, 0.0        , 0.0],
    [0.0              , 0.0        , -1/(fZ-nZ) , 0.0],
    [0.0              , 0.0        , -nZ/(fZ-nZ), 1.0]
  ])
  return B


#######################################################
# サンプルのメイン
#######################################################

#----------------------------------
# 1. 直接指定する場合
#----------------------------------
# 右端を起点に左に視線を向ける
(w, h) = im.size

# 魔法の関数を呼ぶ。
coeffs = find_coeffs(
        [(0, 0), (w/2, w/4), (w/2, h - w/4), (0, h)],
        [(0, 0), (w  ,  0 ), (w  , h)      , (0, h)])

print("coeffs:", coeffs)

subnm = "directly"
im2 = im.transform(im.size, Image.PERSPECTIVE, coeffs, Image.BILINEAR)
filenm = args.output + "_" + subnm + ".png"
print("filenm:", filenm)
im2.save(filenm)


#----------------------------------
# 2. 視野変換、透視投影変換を行う場合
#----------------------------------
#------------------------------------
# 視野変換、透視投影変換のパラメータ
# nZ : 透視投影スクリーンまでの距離
# ev : カメラの移動ベクトル
# tv : 注目点の移動ベクトル
# uv : カメラの上ベクトル
#------------------------------------
nZ = im.size[0]/2
ev = (-im.size[0]/2, 0,  im.size[0]/4)

tv = (im.size[0]/5, 0, -nZ)
uv = (0, 1, 0)

#---------------------

params = (nZ, ev[0],ev[1],ev[2],tv[0],tv[1],tv[2],uv[0],uv[1],uv[2])

A = getSiyaHenkanMatrix(ev, tv, uv)
print("A:",A)

B = getTosiHenkanMatrix(im.size, nZ, ev)
print("B:",B)

(w, h) = (float(im.size[0]), float(im.size[1]))
P = np.matrix([
  [-w/2, -h/2, -nZ, 1],
  [ w/2, -h/2, -nZ, 1],
  [ w/2,  h/2, -nZ, 1],
  [-w/2,  h/2, -nZ, 1]
], dtype=float)

Pv = P * A
Pwk = P * A * B
#print("P:",P)
#print("Pwk:",Pwk)
#print("P * A", P * A)

Pn = []
for i in range(Pwk.shape[0]) :
  Pn.append((int(round(Pwk[i, 0]/(-Pv[i,2])*w/2 + w/2)),
             int(round(Pwk[i, 1]/(-Pv[i,2])*h/2 + h/2))))
  
  #print("p1:", Pwk[i, :])
  #print(" p1x/z:", Pwk[i, 0]/(-Pv[i,2]), "p1y/z:", Pwk[i,1]/(-Pv[i,2]))
  #print(" p1x' :", Pwk[i, 0]/(-Pv[i,2])*w/2,
  #       "p1y' :", Pwk[i, 1]/(-Pv[i,2])*h/2)


print("Pn:", Pn)

# 魔法の関数を呼ぶ。
coeffs = find_coeffs(
        [Pn[0], Pn[1], Pn[2], Pn[3]],
        [(0, 0), (w, 0), (w, h), (0, h)]
)

#print("coeffs:", coeffs)

subnm = "bysitenido"
im2 = im.transform(im.size, Image.PERSPECTIVE, coeffs, Image.BICUBIC)
filenm = args.output + "_" + subnm + "_" + "_".join(map(str, params)) + ".png"
print("filenm:", filenm)
im2.save(filenm)

変換元の画像

f:id:umejan:20160427223126p:plain ※ふちが白枠。

変換後の画像

  • 直接変形
    f:id:umejan:20160427223204p:plain

  • 視野変換、透視投影変換して変形
    f:id:umejan:20160427223222p:plain

コメント

このImage.PERSPECTIVEにはまる。

最初、変換後の座標を適当に決めてa,b,c,d,e,f,g,hを求めようとしたが、 何度も解を間違え、あげくには検証するための計算も間違っていたという情けないことに。

行列の計算もあやしく、行列の計算で大きな勘違いもあった。

あらためて自分が数学が苦手であることを認識させられた。勉強すれば何とかなる部分があるにしても、集中力の足りなさを実感する。

それにしても、このa,b,c,d,e,f,g,hを計算する処理をネットにあげた人に感謝である。 これがなければあきらめていた。

関連

PILのImage.transformの使い方 (EXTENT) - umejanのブログ
PILのImage.transformの使い方 (AFFINE) - umejanのブログ
PILのImage.transformの使い方 (QUAD, MESH) - umejanのブログ