気泡の輪郭の移動量を求める

前提・実現したいこと

ここに質問の内容を詳しく書いてください。

図のような2つの気泡の輪郭データ（frame1オレンジ、frame2緑）があります。
この気泡の界面の速度を求めるために、移動量を求めようとしています。
方法として、
１：frame1(オレンジ)の輪郭座標の隣接する2点を結ぶ直線（傾き、切片）と、その2点間の中点を求め、これを気泡の界面を示す直線とします。
２：その2点間の中点を通り、界面を示す直線の法線（傾き、切片）を求めます。
３：手順1と同様にしてframe2(緑)の界面を示す直線を求めます。
４：手順2で求めた法線と、手順３で求めたframe2の界面を示す直線の交点を求めます。
５：手順1の中点と、手順4で求めた交点の距離が気泡の移動量となる。

発生している問題・エラーメッセージ

上記のような手順で気泡の界面の移動量を求めようとして下記のようなプログラムを組みました。
しかし、手順4で求めた交点は図中の赤色でプロットされてたような箇所にまとまってしまいました。（本来なら、ほぼframe2の輪郭と重なるようにプロットされるはず）

どのようにプログラムを修正してよいか分からず困っています。少しでも修正箇所が見つかるだけでも非常に助かります。
助言のほどよろしくお願いいたします。

該当のソースコード

import cv2
import numpy as np
import csv
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib.patches import Polygon
import glob
import pylab
import math

#画像読み込み
img1 = cv2.imread('D:/github/sample/python/opencv/laplacian/bubble1.tif')
img2 = cv2.imread('D:/github/sample/python/opencv/laplacian/bubble2.tif')

#輪郭検出
def contours(img):
    #グレイスケール変換
    gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
    #grayを閾値処理
    thresh = cv2.threshold(gray, 210, 255, cv2.THRESH_BINARY)[1]
    #ラプラシアンフィルター
    dst = cv2.Laplacian(thresh, cv2.CV_64F, ksize=3)
    #dstをFloat64からuint8に変換(findContoursにFloat64入らない)
    dst_uint8 = dst.astype(np.uint8)
    #dstから輪郭検出
    contours,hierarchy = cv2.findContours(dst_uint8.copy(), cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)
    
    return contours, hierarchy

#輪郭の座標データをx，yに分けてlistに保存
def xylist(contours):
    x_list = []
    y_list = []
    for i in range(0, len(contours)):
            buf_np = contours[i].flatten()# numpyの多重配列になっているため、一旦展開する。
            #print(buf_np)
            for j, elem in enumerate(buf_np):
                if j%2==0:
                    x_list.append(elem)
                else:
                    y_list.append(elem*(-1))
                                 
    return x_list, y_list

#輪郭関数呼び出し
contours1, hierarchy1 = contours(img1)
contours2, hierarchy2 = contours(img2)



x_list1, y_list1 = xylist(contours1)
x_list2, y_list2 = xylist(contours2)
#print(x_list1)



#輪郭座標を0列目x、1列目yのndarrayに格納
def xy_arr(x_list, y_list):
    x_arr_pre = np.array(x_list)
    x_arr = np.reshape(x_arr_pre, [len(x_arr_pre), 1])
    xy_arr = np.insert(x_arr, 1, y_list, axis=1)
    
    return xy_arr

xy_arr1 = xy_arr(x_list1, y_list1)
xy_arr2 = xy_arr(x_list2, y_list2)


'''
#重心検出関数
def center_of_gravity(contours):
    x_CoG = []#重心Center of Gravity
    y_CoG = []
    for i in range(0, len(contours)):
        cnt = len(contours[i])#contours1[i]内の座標数をカウント((x,y)で一つとカウント)
        CoG_sum = sum(contours[i]).flatten()#一旦展開する
        CoG = list(map(lambda x: x/cnt, CoG_sum))#重心=要素総和/要素数
        for i, elem in enumerate(CoG):#重心座標をx,yごとリストに保存
            if i%2==0:
                x_CoG.append(elem)
            else:
                y_CoG.append(elem*(-1)) 
                
    return x_CoG, y_CoG

x_CoG1, y_CoG1 = center_of_gravity(contours1)
x_CoG2, y_CoG2 = center_of_gravity(contours2)
'''


#バブル界面ベクトル，隣接輪郭座標の中点導出関数
def InterfaceVelocity_and_ContousCenter(xy_arr):
    #バブルの界面を表すベクトルと隣接する輪郭座標の中点を求める
    #vector_arr=隣接する輪郭座標どうしのベクトル（バブル界面を表すベクトル）
    vector_arr = np.zeros((len(xy_arr), 2))#sort_theta_x_y_npの要素数　× 2の0埋めndarray作成
    #contours_center_arr=隣接する輪郭座標の中点[Xc, Yc]
    contours_center_arr = np.zeros((len(xy_arr), 2))
    for i in range(0,len(xy_arr)):
        if i+1 < len(xy_arr):#i+1が要素数までの処理
            #バブル界面ベクトルを求める
            vector_x = xy_arr[i+1][0]-xy_arr[i][0]#x2-x1
            vector_y = xy_arr[i+1][1]-xy_arr[i][1]
            vector_arr[i, 0] = vector_x#vector_arrayに格納
            vector_arr[i, 1] = vector_y

            #隣接する輪郭の座標のどうしの中点を求める
            contours_center_x = (xy_arr[i+1][0]+xy_arr[i][0])/2#(x1+x2)/2
            contours_center_y = (xy_arr[i+1][1]+xy_arr[i][1])/2
            contours_center_arr[i, 0] = contours_center_x
            contours_center_arr[i, 1] = contours_center_y
            

        else:#iが要素数まで来たら最初と最後をつなげる
            #バブル界面ベクトルを求める
            vector_x = xy_arr[0][0]-xy_arr[i][0]#x2-x1
            vector_y = xy_arr[0][1]-xy_arr[i][1]
            vector_arr[i, 0] = vector_x#vector_arrayに格納
            vector_arr[i, 1] = vector_y

            #隣接する輪郭の座標のどうしの中点を求める
            contours_center_x = (xy_arr[0][0]+xy_arr[i][0])/2#(x1+x2)/2
            contours_center_y = (xy_arr[0][1]+xy_arr[i][1])/2
            contours_center_arr[i, 0] = contours_center_x
            contours_center_arr[i, 1] = contours_center_y
            
    return vector_arr, contours_center_arr


#バブル界面ベクトル，隣接輪郭座標中点を受け取るvector_arr=[x2-x1,y2-y1]，contours_center_arr=[(x1+x2)/2,(y1+y2)/2]
vector_arr1, contours_center_arr1 = InterfaceVelocity_and_ContousCenter(xy_arr1)
vector_arr2, contours_center_arr2 = InterfaceVelocity_and_ContousCenter(xy_arr2)



#バブル界面ベクトルの傾きと，法線ベクトルの傾き
def tyokusen(vector_arr, contours_center_arr):
#直線の方程式y=ax+b
    kaimen_gradient_intercept_arr = np.zeros((len(vector_arr), 4))#[a,b,x,y],(vector_arrの要素数　× 4の0埋めndarray作成)
    housen_gradient_intercept_arr = np.zeros((len(vector_arr), 4))#[a,b,x,y]=[傾きa，切片b, 法線がx=一定の場合のx, y]
    for i in range(0, len(vector_arr)):
        if vector_arr[i, 0] != 0:#(x2-x1)≠0の時，バブル界面ベクトルの傾き存在a≠∞，-∞
            kaimen_gradient = vector_arr[i, 1]/vector_arr[i, 0]#(y2-y1)/(x2-x1)
            
            if kaimen_gradient == 0:#界面ベクトルの傾き0の時,a=0
                #界面はy=一定，法線はx=一定
                #界面の直線y=(y1+y2)/2
                kaimen_gradient_intercept_arr[i, 0] = 0#a=0
                kaimen_gradient_intercept_arr[i, 1] = contours_center_arr[i, 1]#b=(y1+y2)/2

                #法線の直線x=(x1+x2)/2
                housen_gradient_intercept_arr[i, 2] = contours_center_arr[i, 0]#x=(x1+x2)/2

            else:#界面ベクトルの傾きが0以外なら，a≠∞，-∞, 0
                #界面はy=ax+b，法線はy=-1/a*x+b
                #界面の直線y=ax+b,a=(y2-y1)/(x2-x1), b=(y1+y2)/2
                kaimen_gradient_intercept_arr[i, 0] = kaimen_gradient
                kaimen_gradient_intercept_arr[i, 1] = contours_center_arr[i, 1]-kaimen_gradient*contours_center_arr[i, 0]
                #法線の直線y=-1/a*x+b
                housen_gradient_intercept_arr[i, 0] = -1/kaimen_gradient
                housen_gradient_intercept_arr[i, 1] = contours_center_arr[i, 1]-housen_gradient_intercept_arr[i, 0]*contours_center_arr[i, 0]

        else:#(x2-x1)=0の時，バブル界面ベクトルの傾き無限大a=∞，-∞
            #界面ベクトルはx=一定，法線はy=一定
            #界面の直線x=(x1+x2)/2
            kaimen_gradient_intercept_arr[i, 2] = contours_center_arr[i, 0]#x=(x1+x2)/2

            #法線の直線y=(y1+y2)/2
            housen_gradient_intercept_arr[i, 0] = 0#a=0
            housen_gradient_intercept_arr[i, 1] = contours_center_arr[i, 1]#b=(y1+y2)/2
            
            
    #[傾き，切片，x, y]  法線がx=一定となるときのxをarr中のxに保存      
    return kaimen_gradient_intercept_arr, housen_gradient_intercept_arr


#[傾き，切片，x, y] を取得
kaimen_gradient_intercept_arr1, housen_gradient_intercept_arr1 = tyokusen(vector_arr1, contours_center_arr1)
kaimen_gradient_intercept_arr2, housen_gradient_intercept_arr2 = tyokusen(vector_arr2, contours_center_arr2)

    


#frame1の法線とframe2の界面の直線の交点を求める
#kouten_arr=[x, y, d]([交点x, 交点y, 変位量d])
kouten_arr = np.zeros((len(kaimen_gradient_intercept_arr2), 3))#kaimen_gradient_intercept_arr1の要素数　× 3の0埋めndarray作成
min_kouten_arr = np.zeros((len(housen_gradient_intercept_arr1), 3))
min_kouten_arr[:, 2] = np.inf
for i in range(0, len(housen_gradient_intercept_arr1)):#法線の数(i本)だけ交点を求めるループ
    for j in range(0, len(kaimen_gradient_intercept_arr2)):#1本の法線に対して，j個のframe2の界面との交点を探すループ
        if housen_gradient_intercept_arr1[i, 2] != 0:#frame1の法線の式がx=一定でなければ
            #(a1-a2)≠0平行でないとき，
            if (housen_gradient_intercept_arr1[i,0]-kaimen_gradient_intercept_arr2[j, 0]) != 0:
                #i番目の法線に対するj個の交点の座標
                #交点のｘ座標x=(b2-b1)/(a1-a2)
                kouten_arr[j, 0] = (kaimen_gradient_intercept_arr2[j, 1]-housen_gradient_intercept_arr1[i, 1])\
                                        /(housen_gradient_intercept_arr1[i,0]-kaimen_gradient_intercept_arr2[j, 0])
                #交点のy座標y=ax+b 
                kouten_arr[j, 1] = housen_gradient_intercept_arr1[i, 0]*kouten_arr[j, 0]+housen_gradient_intercept_arr1[i, 1]
                
            else:#(a1-a2)=0平行のとき，変位を∞とする
                kouten_arr[j, 2] = np.inf
          
        else:#frame1法線の式がx=一定(中点のx座標の値)の場合
            if kaimen_gradient_intercept_arr1[i, 2] != 0:#frame2の界面の式がx=一定でなければ(交点あり)
                #x=(x1+x2)/2
                kouten_arr[j, 0] = housen_gradient_intercept_arr1[i, 2]
                #交点のy座標y=(a2)*x+(b2) 
                kouten_arr[j, 1] = kaimen_gradient_intercept_arr2[j, 0]*kouten_arr[j, 0]\
                                    +kaimen_gradient_intercept_arr2[j, 1]
                
            else:#frame2の界面の式がx=一定(交点なし)
                kouten_arr[j, 2] = np.inf
            
        #交点(x,y)とframe1の中点(Xc, Yc)との距離(変位量)d=√((x-Xc)^2+(y-Yc)^2)
        kouten_arr[j, 2] = math.sqrt((kouten_arr[j, 0]-contours_center_arr1[j, 0])**2\
                                     +(kouten_arr[j, 1]-contours_center_arr1[j, 1])**2)
    
    #dの最小値のインデックスを取得
    min_kouten = np.argmin(kouten_arr, axis=0)
    min_d_idx = min_kouten[2]
    #dの最小値とそれに対応した交点座標を格納[x, y, d]
    min_kouten_arr[i, 0] = kouten_arr[min_d_idx, 0]
    min_kouten_arr[i, 1] = kouten_arr[min_d_idx, 1]
    min_kouten_arr[i, 2] = kouten_arr[min_d_idx, 2]
  

plt.scatter(x_list1, y_list1, s=1)
plt.scatter(x_list2, y_list2, s=1)
plt.scatter(min_kouten_arr[:, 0], min_kouten_arr[:, 1])
plt.axes().set_aspect('equal', 'datalim')


### 試したこと

ここに問題に対して試したことを記載してください。

### 補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

fana

2020/11/17 01:54

（私には読めないので）コードを見ずに話の意味を問いますが… オレンジの輪郭上のある１点についてこの処理を考えたとき，この点から計算された法線Ｎとすると，緑の輪郭上の点の個数がm個であるときに，この法線Nに関して最大m個の「交点」が算出されると思うのですが，どのような基準でどれを選ぶんでしょう？

mk_taro5

2020/11/17 02:13

質問を読んでいただきありがとうございます。おっしゃる通り、frame1の1本の法線に対して、frame2の界面を示す直線m本との交点をｍ個求めます。その後、frame1から求めた隣接する２つの輪郭座標の中点と、交点との距離をｍ通り求め、その距離が最短となる交点を選びます。このようにして選択された交点が、frame1から求めた隣接する２つの輪郭座標の中点のframe2時における移動先であると考えています。

fana

2020/11/17 04:04

とりあえず何か単純なテストケースみたいなのを作って，ある１本の法線に対して求められた複数の「交点」群を表示等して確認してみてはいかがでしょうか．その中に正当な（狙い通りの）位置を示すものが存在するならば，選択部分のミスでしょうし，存在しないなら交点計算とかそこらへんに誤りがあるのだと判断できるように思います．

mk_taro5

2020/11/17 08:20

アドバイスありがとうございます．やってみたいと思います．