0

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)
'''
3
  • 難問ですね。sample0.jpgは、縦横、どれぐらいの画素数の画像ですか?小さい画像でも発生しますか。(興味深い問題ではあります。)
    – boundary1
    2020年2月5日 9:05
  • 自己解決しました。 処理が途中で止まってしまう原因はやはり処理落ちでした。おそらく逆行列の計算がうまくいかずに、処理落ちをおこしていたみたいです。スペックの高い解析用PCを使ったところ計算が終了した時があったので、そちらの結果を使用します。
    – NNL
    2020年2月6日 5:27
  • 解決に至った情報をもしよければ「自己回答」として投稿してみてください。分かる範囲で構いませんので、NGだった環境とスペックの高い環境とでどの程度違いがあるのかなどの情報もあると、より参考になると思います。
    – cubick
    2020年2月10日 1:12

0

この質問に回答するには、ログインする必要があります。

のタグが付いた他の質問を参照する。