表面から8枚の画像を積層し、beamパターンを取得した後に、そのbeamパターンと各画像のずれを集計するプログラムである。

beampattern3.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 = []
addcount = 8
flag = 0
i = 0
base_u = 71
base_d = 83
fin = 157

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))

img = cv2.imread(path + "mod64pl2000.png",0);
beam = np.zeros_like(img)

while i < addcount:
    
    img_b = cv2.imread(path + "mod64pl200{0}.png".format(i),0);
#    cv2.imshow('img_b',img_b)
    cont_b = mytool.funcContrust(img_b)
    bin_b = mytool.funcThreshold(cont_b,31,60,1)
#    cv2.imshow('bin_b',bin_b * 255)
    beam = cv2.add(beam,bin_b)
#    cv2.imshow('beam',beam * 255)
    i += 1

threshold = 4
max_bri = 255
ret,thre = cv2.threshold(beam, threshold, max_bri, cv2.THRESH_BINARY)
contours_b = cv2.findContours(thre, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
#cv2.imshow('thre',thre)
#cv2.waitKey(0)
#cv2.destroyAllWindows()


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

        if i == base_d:
            flag = 0 

        num.append(float(i))

#下側の画像から輝度値のまとまりを読み取る。    
        cont_p = mytool.funcContrust(img_pattern)
        bin_p = mytool.funcThreshold(cont_p,31,60,255)
        contours_p = cv2.findContours(bin_p, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)    

        mom_b = []
        mom_p = []
    
        for cb in contours_b[0]:
            mom = {}
            M = cv2.moments(cb)
            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_b.append(mom)
            else:
                cx = 0
                cy = 0

        for cp in contours_p[0]:
            mom = {}
            M = cv2.moments(cp)
            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_p.append(mom)
            else:
                cx = 0
                cy = 0

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

        for b_b,b_p in itertools.product(mom_b,mom_p):
#    print(b_u,b_d)
    
#もしかしたら、画像の中心を（0，0）になるようにしなければならないかも
#その場合は、絶対値（abs）を使用する必要あり    
            r = math.sqrt((b_b['cx'] - b_p['cx'])**2 + (b_b['cy'] - b_p['cy'])**2)
    
            if r > 100:
                continue
    
            dx = b_b['cx'] - b_p['cx']
            dy = b_b['cy'] - b_p['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_3:beam_to_img{0}.png'.format(i))
        i += 1

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

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

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


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

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

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