distortionカーブをpython上で実装する方法

Python

はじめに

distortionとはザイデルの5収差の一つで、一般的にカメラで撮った画像がゆがむのはこの特性が影響しています。
(最近カメラレンズのdistortionはほとんど0またはcalibration技術の発達によってそこまで歪んでいるようには見えませんが…)
こちらについて原理から詳しく知りたい場合はボルンの光学の原理あたりを読むのをおすすめします。
また簡単な理解ではエドモンドのDitortionのページに書かれているものを読むとある程度の理解はできます。

しかし、これらの説明を受けても具体的に性能欄にあるレンズのDistortion値で数10%と言われても直感的に理解できません。
また、もっと詳しい資料になると入射角(またはsensorの像高:image height)ごとのDisotion値、いわゆるDistortion curveを受け取ることもしばしば…

これらについて、手元でいじくることで直感的に理解できるようにしたいと思いPythonで実装してみました。

数学的な理解

数学的な理解といっても難しいことはありません。
原理的な話となると5収差はh = image heihtのhの関数として表せて…というザイデルの5収差そのものとそれに関連する数学の話になりますが、今回は実装上の話だけをします。
Distortion前のimage heightを\(AD\)、Distortion後のimage heightを\(PD\)とするとDistortionの割合\(D\)は以下の関係で定義されます。

$$D(AD) (\%) \equiv \frac{AD-PD}{PD}　\times 100 (\%)$$

ここで\(D\)はレンズの設計データから一般的に得ることができて、ADに依存した値になっています。
今回使う数式はこれだけになります。
注意すべきはデータ処理上この数式を使うと計算がしやすいだけで、物理的なものに基づくものではありません。
(そもそもDistorionはimage heightの3乗に比例しますし…)

Pythonでの実装

まず、今回はサンプルのために適当にDistortion関数を作りました。
本来はデータからcurve_fit関数なんかを使ってDistorion関数を推定してあげてください。

import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit
import numpy as np

def distortion_function(theta,a,b,c,d,e):
    return a*theta**4+b*theta**3+c*theta**2+d*theta+e

theta = np.arange(0,45,1)
#適当にDistortion関数を定義
popt = [0.000000045,0.000003,0.000003,0.004,0]
plt.plot(theta,distortion_function(theta,*popt))
plt.xlabel("angle (deg)")
plt.ylabel("distortion")
plt.xlim(0)
plt.ylim(0)

これで次のようなDistortionカーブとなります。

また、続いて格子をわかりやすく表記するためにbold関数を定義してあげます。
その後格子を描写してみます。

def bold(array,x_center,y_center,width,intensity):
    point = [-int(width//2)+i for i in range(width)]
    for i in point:
        for j in point:
            try:
                array[x_center+i,y_center+j] = intensity
            except:
                continue
    return array

# 1mm×1mm(1um間隔、-0.5~0.5mm)の座標を作る, unit:um
axis = np.arange(-500,501,1)
grid = np.meshgrid(axis,axis)[0]

### distorionのもととなるグリッド線を作る
#0のデータを格納するarrayを作る
grid_line = np.zeros_like(grid)
#grid間隔
interval = 100

for i,data in enumerate(grid):
    for j,_ in enumerate(data):
        x,y = i-500,j-500
        if (x%interval==0 or y%interval==0):
            grid_line = bold(grid_line,i,j,5,255)

#grid plot
plt.imshow(grid_line,cmap="gray")

#横軸縦軸のラベルを整える
ran = int(1000/interval)+1
ticks = [i * 100 for i in range(ran)]
tick_labels = [-500 + i *100 for i in range(ran)]
plt.xticks(ticks, tick_labels)
plt.yticks(ticks, tick_labels)
#ラベリング
plt.xlabel("horizontal (um)")
plt.ylabel("vertical (um)")


この格子にdistortionをかけてあげると次のようになります。
ここで端にdistortion関数をかけるとグリッド外にでてしまうのでif文で弾いてあります。

# 1mm×1mm(1um間隔、-0.5~0.5mm)の座標を作る, unit:um
axis = np.arange(-500,501,1)
grid = np.meshgrid(axis,axis)[0]

### distorionデータを格納する箱を作る
#0のデータを格納するarrayを作る
distortion_line = np.zeros_like(grid)
#grid間隔
interval = 100

for i,data in enumerate(grid):
    for j,_ in enumerate(data):
        #sensor上でのx,y座標に変換
        x,y = i-500,j-500
        #print(x,y)
        if ((i >= 50 and i <= 950) and (j >= 50 and j <= 950)):
            if (x%interval==0 or y%interval==0):
                #ゆがむ前のラインを保存する、ここで薄くするために100を代入
                #distortion_line = bold(distortion_line,i,j,5,50)

                #EFL = 3mmとするとsensor上での座標はr(x,y)=3mm*tan(theta)で表せる
                #よってtheta = arctan(r um/3000 um)
                theta = np.arctan(np.sqrt(x**2+y**2)/3000)
                theta = np.rad2deg(theta)
                m = distortion_function(theta,*popt)
                #ゆがんだ座標=m*ゆがむ前の座標+ゆがむ前の座標
                x_dis = int(m*x+x)
                y_dis = int(m*y+y)
                #データに格納するためにデータ座標に変換
                i_dis,j_dis = x_dis+500,y_dis+500
                distortion_line = bold(distortion_line,i_dis,j_dis,5,250)


#grid plot
plt.imshow(distortion_line,cmap="gray")

#横軸縦軸のラベルを整える
ran = int(1000/interval)+1
ticks = [i * 100 for i in range(ran)]
tick_labels = [-500 + i *100 for i in range(ran)]
plt.xticks(ticks, tick_labels)
plt.yticks(ticks, tick_labels)
#ラベリング
plt.xlabel("horizontal (um)")
plt.ylabel("vertical (um)")

 

全コード

import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit
import numpy as np

def distortion_function(theta,a,b,c,d,e):
    return a*theta**4+b*theta**3+c*theta**2+d*theta+e

def bold(array,x_center,y_center,width,intensity):
    point = [-int(width//2)+i for i in range(width)]
    for i in point:
        for j in point:
            try:
                array[x_center+i,y_center+j] = intensity
            except:
                continue
    return array


theta = np.arange(0,45,1)
popt = [0.000000045,0.000003,0.000003,0.004,0]
plt.plot(theta,distortion_function(theta,*popt))
plt.xlabel("angle (deg)")
plt.ylabel("distortion")
plt.xlim(0)
plt.ylim(0)
plt.savefig("distortion_cuve.png")
plt.show()

# 1mm×1mm(1um間隔、-0.5~0.5mm)の座標を作る, unit:um
axis = np.arange(-500,501,1)
grid = np.meshgrid(axis,axis)[0]

### distorionのもととなるグリッド線を作る
#0のデータを格納するarrayを作る
grid_line = np.zeros_like(grid)
#grid間隔
interval = 100

for i,data in enumerate(grid):
    for j,_ in enumerate(data):
        x,y = i-500,j-500
        if (x%interval==0 or y%interval==0):
            grid_line = bold(grid_line,i,j,5,255)

#grid plot
plt.imshow(grid_line,cmap="gray")

#横軸縦軸のラベルを整える
ran = int(1000/interval)+1
ticks = [i * 100 for i in range(ran)]
tick_labels = [-500 + i *100 for i in range(ran)]
plt.xticks(ticks, tick_labels)
plt.yticks(ticks, tick_labels)
#ラベリング
plt.xlabel("horizontal (um)")
plt.ylabel("vertical (um)")
plt.show()

# 1mm×1mm(1um間隔、-0.5~0.5mm)の座標を作る, unit:um
axis = np.arange(-500,501,1)
grid = np.meshgrid(axis,axis)[0]

### distorionデータを格納する箱を作る
#0のデータを格納するarrayを作る
distortion_line = np.zeros_like(grid)
#grid間隔
interval = 100

for i,data in enumerate(grid):
    for j,_ in enumerate(data):
        #sensor上でのx,y座標に変換
        x,y = i-500,j-500
        #print(x,y)
        if ((i >= 50 and i <= 950) and (j >= 50 and j <= 950)):
            if (x%interval==0 or y%interval==0):
                #ゆがむ前のラインを保存する、ここで薄くするために100を代入
                #distortion_line = bold(distortion_line,i,j,5,50)

                #EFL = 3mmとするとsensor上での座標はr(x,y)=3mm*tan(theta)で表せる
                #よってtheta = arctan(r um/3000 um)
                theta = np.arctan(np.sqrt(x**2+y**2)/3000)
                theta = np.rad2deg(theta)
                m = distortion_function(theta,*popt)
                #ゆがんだ座標=m*ゆがむ前の座標+ゆがむ前の座標
                x_dis = int(m*x+x)
                y_dis = int(m*y+y)
                #データに格納するためにデータ座標に変換
                i_dis,j_dis = x_dis+500,y_dis+500
                distortion_line = bold(distortion_line,i_dis,j_dis,5,250)


#grid plot
plt.imshow(distortion_line,cmap="gray")

#横軸縦軸のラベルを整える
ran = int(1000/interval)+1
ticks = [i * 100 for i in range(ran)]
tick_labels = [-500 + i *100 for i in range(ran)]
plt.xticks(ticks, tick_labels)
plt.yticks(ticks, tick_labels)
#ラベリング
plt.xlabel("horizontal (um)")
plt.ylabel("vertical (um)")
plt.show()

参考、その他グラフ用
https://florisryo.com/2021/07/10/soturon-python/

コメント

タイトルとURLをコピーしました