このプログラムは、上から下まで撮影した写真の上下の写真のずれを取得し、グラフ化するプログラムである。

beampattern2.py

# -*- coding: utf-8 -*-
"""
Created on Wed Aug 17 15:27:51 2016

@author: ryousuke
"""


import numpy as np
import cv2
from array import array
#import re

import pandas as pd
from os.path import join, relpath
from glob import glob
import itertools
import ROOT
import math
import mytool


ROOT.std.__file__ = 'dummy_for_old_pyastro'



ROOT.gROOT.Reset()

path = "C:\\Users\\GTR\\Documents\\lab_log\\log\\H28_11\\20161117\\beam_info\\GTR_test\\img\\"

#img_lob = []

flag = 0
i = 0
base_u = 71
base_d = 83
fin = 157

myarea = 20

num = array('d')
b_x = array('d')
b_y = array('d')
b_xsum = array('d')
b_ysum = array('d')

zero = 0
num.append(float(zero))
b_x.append(float(zero))
b_y.append(float(zero))
b_xsum.append(float(zero))
b_ysum.append(float(zero))

while i < fin:
    
    if flag == 1:
        i = base_d
        print i
    
    if i == base_u:
        flag = 1
        print i
    
#上側の写真を読み取る場所
    if i < 10:
        img_u = cv2.imread(path + "mod64pl200{0}.png".format(i),0);
    elif i < 100:
        img_u = cv2.imread(path + "mod64pl20{0}.png".format(i),0);
    else:
        img_u = cv2.imread(path + "mod64pl2{0}.png".format(i),0);

#下側の写真を読み取る場所    
    s = i + 1
    if s == base_u + 1:
        s = base_d
    
    if s < 10:
        img_d = cv2.imread(path + "mod64pl200{0}.png".format(s),0);
    elif s < 100:
        img_d = cv2.imread(path + "mod64pl20{0}.png".format(s),0);
    else:
        img_d = cv2.imread(path + "mod64pl2{0}.png".format(s),0);
    

    if i == base_d:
        flag = 0 


    num.append(float(s))
#    img_lob.append(i)

#上側の画像から輝度値のまとまりを読み取る。
#    img_u = cv2.imread(path + "mod64pl2000.png",0);
    cont_u = mytool.funcContrust(img_u)
    bin_u = mytool.funcThreshold(cont_u,31,60,255)
#    con_u = mytool.mycontour(bin_u,myarea)
    contours_u = cv2.findContours(bin_u, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)    

#下側の画像から輝度値のまとまりを読み取る。    
#    img_d = cv2.imread(path + "mod64pl2000.png",0);
    cont_d = mytool.funcContrust(img_d)
    bin_d = mytool.funcThreshold(cont_d,31,60,255)
#    con_d = mytool.mycontour(bin_d,myarea)
    contours_d = cv2.findContours(bin_d, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)    


    mom_u = []
    mom_d = []
    
    for cu in contours_u[0]:
        mom = {}
        M = cv2.moments(cu)
        if M['m00'] != 0:
            cx = int(M['m10']/M['m00'])
            cy = int(M['m01']/M['m00'])
            mom['cx'] = int(cx)
            mom['cy'] = int(cy)
            mom_u.append(mom)
        else:
            cx = 0
            cy = 0

    for cd in contours_d[0]:
        mom = {}
        M = cv2.moments(cd)
        if M['m00'] != 0:
            cx = int(M['m10']/M['m00'])
            cy = int(M['m01']/M['m00'])
        
            mom['cx'] = int(cx)
            mom['cy'] = int(cy)
            mom_d.append(mom)
        else:
            cx = 0
            cy = 0

    c1 = ROOT.TCanvas('c1','Example with Formula',200,10,700,700)
    h2dydx = ROOT.TH2D("h2dydx","img:{0} to img:{1};x[pix];y[pix]".format(i,s),100, -50, 50, 100, -50, 50)

    for b_u,b_d in itertools.product(mom_u,mom_u):
#    print(b_u,b_d)
    
#もしかしたら、画像の中心を（0，0）になるようにしなければならないかも
#その場合は、絶対値（abs）を使用する必要あり    
        r = math.sqrt((b_u['cx'] - b_d['cx'])**2 + (b_u['cy'] - b_d['cy'])**2)
    
        if r > 100:
            continue
    
        dx = b_u['cx'] - b_d['cx']
        dy = b_u['cy'] - b_d['cy']
        h2dydx.Fill(dx,dy)

    print 'finish'
    
    x = h2dydx.GetMean(1)
    y = h2dydx.GetMean(2)
    sum_x = b_xsum[-1] + x
    sum_y = b_ysum[-1] + y
    
    b_x.append(float(x))
    b_y.append(float(y))
    b_xsum.append(float(sum_x))
    b_ysum.append(float(sum_y))            
    
    h2dydx.Draw('colz')
    c1.Print(path + 'Beampattern:img{0}_to_img{1}.png'.format(i,s))
    i += 1

#各パターンマッチのヒストグラムが作成できたら、その統計情報を使用して各画像のずれをグラフにする。

n = len(num)
#n_flo = float(n)

gr1 = ROOT.TGraph(n,num,b_x)
gr1.GetXaxis().SetLimits(0,fin + 10)
gr1.SetMaximum(10)
gr1.SetMinimum(-10)
gr1.SetTitle('bx')
gr1.SetMarkerStyle(7)
gr1.SetMarkerColor(1)
gr1.Draw('AP')
c1.Print(path + 'rusult:bx.png')


gr2 = ROOT.TGraph(n,num,b_y)
gr2.GetXaxis().SetLimits(0,fin + 10)
gr2.SetMaximum(10)
gr2.SetMinimum(-10)
gr2.SetTitle('by')
gr2.SetMarkerStyle(7)
gr2.SetMarkerColor(8)
gr2.Draw('AP')
c1.Print(path + 'rusult:by.png')

gr3 = ROOT.TGraph(n,num,b_xsum)
gr3.GetXaxis().SetLimits(0,fin + 10)
gr3.SetMaximum(50)
gr3.SetMinimum(-70)
gr3.SetTitle('bx_sum')
gr3.SetMarkerStyle(7)
gr3.SetMarkerColor(1)
gr3.Draw('AP')
c1.Print(path + 'rusult:bx_sum.png')

gr4 = ROOT.TGraph(n,num,b_ysum)
gr4.GetXaxis().SetLimits(0,fin + 10)
gr4.SetMaximum(50)
gr4.SetMinimum(-70)
gr4.SetTitle('by_sum')
gr4.SetMarkerStyle(7)
gr4.SetMarkerColor(8)
gr4.Draw('AP')
c1.Print(path + 'rusult:by_sum.png')