pythonを使って画像の相関を取るプログラムを書いています。
計算を始めると最初は動くのですが、しばらくすると何のerrorも示さないまま計算が中断されていしまいます。
自分で調べてみると、PCのメモリが不足している可能性があるとのことでしたが、それ以外の原因は考えられないのでしょうか?
どなたか回答のほどよろしくお願いいたします。
コードは以下の通りです。
129行目「for r in range(int(subset_interpolation.shape[0]/2)):」以下の処理の最中に中断してしまっている状況です。
ソースコード
import cv2
import numpy as np
from numba.decorators import jit
from numba import float64,int64
import time
def numerical_diff(f,x,i,subset,d_subset,split,H):
h = 1e-4
h_vec = np.zeros_like(x)
h_vec[i] = h
return (f(x + h_vec,subset,d_subset,split,H) - f(x - h_vec,subset,d_subset,split,H)) / (2*h)
def hessian(f,x,subset,d_subset,split,H,i,j):
h = 1e-4
h_i = np.zeros_like(x)
h_j = np.zeros_like(x)
h_i[i] = h
h_j[j] = h
if i == j:
fpp = (f(x + 2*h_i,subset,d_subset,split,H) - f(x,subset,d_subset,split,H))/(2*h)
fmm = (f(x,subset,d_subset,split,H) - f(x - 2*h_i,subset,d_subset,split,H))/(2*h)
return (fpp - fmm)/(2*h)
else:
fpp = (f(x + h_i + h_j,subset,d_subset,split,H) - f(x- h_i + h_j,subset,d_subset,split,H))/(2*h)
fmm = (f(x + h_i - h_j,subset,d_subset,split,H) - f(x - h_i - h_j,subset,d_subset,split,H))/(2*h)
return (fpp - fmm)/(2*h)
@jit(float64[:,:](float64[:,:],int64,float64),nopython = True)
def Heaviside(subset,r,theta):
subset_center = (subset.shape[0]/2,subset.shape[1]/2)
split = np.empty_like(subset)
for x in range(subset.shape[0]):
for y in range(subset.shape[1]):
R = (x-subset_center[1])*np.cos(theta*np.pi/180.) + (y-subset_center[0])*np.sin(theta*np.pi/180.) - r
if R < 0:
split[y][x] = 0
if R >=0:
split[y][x] = 1
return split
@jit(float64(float64[:],float64[:,:],float64[:,:],float64[:,:],int64),nopython = True)
def cross_correlation_coficient(u,subset,d_subset,split,H):
fg = 0
f2 = 0
g2 = 0
F = subset
G = np.empty_like(subset)
for id_h in list(zip(*np.where(split == H))):
x = id_h[1]
y = id_h[0]
xp = (x + u[1] + u[2]*np.abs(d_subset.shape[0]/2 - x) + u[3]*np.abs(d_subset.shape[1]/2 - y) + (d_subset.shape[0]/2 - subset.shape[0]/2))
yp = (y + u[0] + u[4]*np.abs(d_subset.shape[0]/2 - x) + u[5]*np.abs(d_subset.shape[1]/2 - y) + (d_subset.shape[1]/2 - subset.shape[1]/2))
xp_int = int(xp)
yp_int = int(yp)
dx = xp - xp_int
dy = yp - yp_int
G[y][x] = (1-dy)*(1-dx)*d_subset[yp_int][xp_int] + dy*(1-dx)*d_subset[(yp_int)+1][xp_int] + (1-dy)*dx*d_subset[yp_int][(xp_int)+1] + dy*dx*d_subset[(yp_int)+1][(xp_int)+1]
fg = fg + F[y][x]*G[y][x]
f2 = f2 + (F[y][x])**2
g2 = g2 + (G[y][x])**2
#print(F[y][x])
if fg == 0:
c = 0
else:
c = fg/np.sqrt(f2*g2)
c = float(c)
return c
t0=time.time()
reference = cv2.imread("C:\\Users\\nakada\\Desktop\\academic\\sample\\sample0.jpg",0)
deformed = cv2.imread("C:\\Users\\nakada\\Desktop\\academic\\sample\\sample1.jpg",0)
subset_center = (300,300)#(y,x)
subset_size = (25,25)
#seek =13
#print(reference)
subset = reference[subset_center[0] - int(subset_size[0]/2):subset_center[0] + int(subset_size[0]/2),subset_center[1] - int(subset_size[1]/2):subset_center[1] + int(subset_size[1]/2)]
result = cv2.matchTemplate(deformed, subset, cv2.TM_CCORR_NORMED)#(探査範囲,探査部位,探査方法)
min_val, max_val, min_loc, max_loc = cv2.minMaxLoc(result)
top_left = max_loc #(x,y)
print('左上',top_left)
deformed_center = (top_left[0]+int(subset.shape[0]/2) , top_left[1]+int(subset.shape[1]/2))
print('変形前中心',(subset_center[1],subset_center[0]))
print('変形後中心',deformed_center)
U0 = (deformed_center[0]-subset_center[1],deformed_center[1]-subset_center[0])
d_subset = deformed[deformed_center[0] - int(subset_size[0]):deformed_center[0] + int(subset_size[0]),deformed_center[1] - int(subset_size[1]):deformed_center[1] + int(subset_size[1])]
print('raw(Ux,Uy)',U0)#(x,y)
scale = 5
subset_interpolation = cv2.resize(subset,(int(scale*subset.shape[1]),int(scale*subset.shape[0])))
subset_interpolation = subset_interpolation.astype(np.float64)
subset_d_interpolation = cv2.resize(d_subset,(int(scale*d_subset.shape[0]),int(scale*d_subset.shape[1])))
subset_d_interpolation = subset_d_interpolation.astype(np.float64)
up = np.array([U0[1]*scale,U0[0]*scale,0,0,0,0],dtype = np.float64)
upD = np.array([U0[1]*scale,U0[0]*scale,0,0,0,0],dtype = np.float64)
history =[]
c_max = 0
for r in range(int(subset_interpolation.shape[0]/2)):
for theta in np.arange(0,180,10):
hist = np.empty(10)
split = Heaviside(subset_interpolation,r,theta)
H = 0
Hesse = np.empty((6,6))
d_up = up
for i in range(20):
for x in range(6):
for y in range(6):
Hesse[x][y] = hessian(cross_correlation_coficient,d_up,subset_interpolation,subset_d_interpolation,split,H,x,y)
#print(Hesse)
H_inv = np.linalg.pinv(Hesse)
grad = np.array([numerical_diff(cross_correlation_coficient,d_up,i,subset_interpolation,subset_d_interpolation,split,H) for i in range(6)])
d_up = d_up - np.dot(H_inv, grad)
grad_norm = np.linalg.norm(grad)
print(grad_norm)
if grad_norm <= 1.0:
c0 = cross_correlation_coficient(d_up,subset_interpolation,subset_d_interpolation,split,H)
hist[0] = c0 #H=0側のc
hist[3] = d_up[0] #H=0側のu_x
hist[4] = d_up[1] #H=0側のu_y
hist[9] = 0 #収束した
break
if grad_norm > 1.0:
c0 = cross_correlation_coficient(d_up,subset_interpolation,subset_d_interpolation,split,H)
hist[0] = c0
hist[3] = d_up[0]
hist[4] = d_up[1]
hist[9] = 1 #収束してない
H = 1
Hesse = np.empty((6,6))
d_upD = upD
for i in range(20):
for x in range(6):
for y in range(6):
Hesse[x][y] = hessian(cross_correlation_coficient,d_up,subset_interpolation,subset_d_interpolation,split,H,x,y)
H_inv = np.linalg.pinv(Hesse)
gradD = np.array([numerical_diff(cross_correlation_coficient,d_up,i,subset_interpolation,subset_d_interpolation,split,H) for i in range(6)])
d_upD = d_upD - np.dot(H_inv, grad)
grad_normD = np.linalg.norm(gradD)
if grad_normD < 1.0:
c1 = cross_correlation_coficient(d_upD,subset_interpolation,subset_d_interpolation,split,H)
hist[1] = c1 #H=1側のc
hist[5] = d_upD[0] #H=1側のu_x
hist[6] = d_upD[1] #H=1側のu_y
hist[9] = hist[9] + 0 #収束した
break
if grad_normD > 1.0:
c1 = cross_correlation_coficient(d_upD,subset_interpolation,subset_d_interpolation,split,H)
hist[3] = c1
hist[5] = d_upD[0]
hist[6] = d_upD[1]
hist[9] = hist[9] + 1 #収束してない
if hist[9] == 0:
hist[2] = hist[0] + hist[1]
if hist[2] > c_max:
c_max = hist[2]
hist =np.array(hist)
history.append(hist)
print(history)
t1 = time.time()
print('計算時間',t1-t0)
'''