

pythonを使って画像の相関を取るプログラムを書いています。  
計算を始めると最初は動くのですが、しばらくすると何のerrorも示さないまま計算が中断されていしまいます。  
自分で調べてみると、PCのメモリが不足している可能性があるとのことでしたが、それ以外の原因は考えられないのでしょうか？  
どなたか回答のほどよろしくお願いいたします。

コードは以下の通りです。  
129行目「for r in
range(int(subset_interpolation.shape[0]/2)):」以下の処理の最中に中断してしまっている状況です。

**ソースコード**

[code]

    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)
    '''
    
[/code]


CC BY-SA 4.0
2020-02-05T06:34:55.210
62798
2020-02-05T06:34:55.210
37096
question
python3 numpy
pythonの計算が突然中断されるが、errorが出ないので原因がわかりません。