ちょいめも

物理/Python/Cの雑記帳

画像のグラフなどを任意関数でフィッティング

画像になってしまってるグラフから線の部分の座標を抽出して
任意の関数でフィッティングしたいと思ったのでPythonで適当に作成。

class Rfit():
	
	def __init__(self, load_path='', mmppx = 1):
		
		#同じフォルダでもフルパスでもファイル名と拡張子を返す
		#splitext:パス名を(filename,ext)に分割、basename:pathの末尾ファイル各部分(拡張子)を返す
		self.filename, self.ext = os.path.splitext(os.path.basename(load_path)) 
		#load_pathそのまま
		self.load_path = load_path 
		#同じフォルダのファイルをload_pathにしたときはdirpathは空、フルパスだと現在のフォルダをフルパスで表示
		self.load_up_path, self.f_name = os.path.split(load_path) 
		#現在のフォルダ 同じフォルダのファイルは空
		self.current_dirname = os.path.basename(self.load_up_path)

		self.mmppx = mmppx

		
	def Img2xy(self):
		#画像から曲線が描かれている(Rの)xy座標を抽出
		
		img = cv2.imread(self.load_path,0) #0:grayscale
		img = cv2.medianBlur(img,5) #メディアンフィルタを用いて画像の平滑化
		
		kernel = np.ones((10,10),np.uint8)
#		img = cv2.erode(img,kernel,iterations = 1) #圧縮フィルタ
#		img = cv2.dilate(img,kernel,iterations = 1) #膨張フィルタ
		img = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel) #オープニング処理 収縮の後に膨張
#		img = cv2.morphologyEx(img, cv2.MORPH_CLOSE, kernel) #膨張の後に収縮
		
		#ret,img = cv2.threshold(img,127,255,cv2.THRESH_BINARY)
		#img = cv2.adaptiveThreshold(img,255,cv2.ADAPTIVE_THRESH_MEAN_C,cv2.THRESH_BINARY,11,2)
		img = cv2.adaptiveThreshold(img,255,cv2.ADAPTIVE_THRESH_GAUSSIAN_C,cv2.THRESH_BINARY,11,2)
		img = cv2.GaussianBlur(img,(5,5),0)
		ret,img = cv2.threshold(img,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
		
#		2値化画像のチェック 何かキーを押すとウインドウを閉じる
#		cv2.imshow('image', img)
#		cv2.waitKey(0)
#		cv2.destroyAllWindows()

		xdata = [i for i in range(0,img.shape[0])]
		ydata = [img[i].argmin() for i in range(0,img.shape[0])] #arrayの一番左の要素番号
		
		img_inv = np.fliplr(img)
		ydata_inv = [img_inv[i].argmin() for i in range(0,img_inv.shape[0])] #arrayの一番右の要素番号
																
		# ===== fitting確認 =====
		if len(self.current_dirname) == 0:
			#np.savetxt(self.filename + '_xydata.csv', np.array([xdata, ydata]), delimiter = ',')
			np.savetxt(self.filename + '_value.csv', img, delimiter = ',')
			np.savetxt(self.filename + '_value_inv.csv', img_inv, delimiter = ',')
		else:
			#np.savetxt(self.load_up_path + '\\' + self.filename + '_xydata.csv', np.array([xdata, ydata]), delimiter = ',')
			np.savetxt(self.load_up_path + '\\' + self.filename + '_value.csv', img, delimiter = ',')
			np.savetxt(self.load_up_path + '\\' + self.filename + '_value_inv.csv', img_inv, delimiter = ',')
		
		self.xdata = xdata
		self.ydata = ydata
		self.ydata_inv = ydata_inv
		
		
	def Value2xy(self):
		#輝度値csvから曲線が描かれている(Rの)xy座標を取得
		
		value = np.loadtxt(self.load_path, delimiter=',')
		
		xdata = [i for i in range(0,value.shape[0])]
		ydata = [value[i].argmin() for i in range(0,value.shape[0])]
											
		value_inv = np.fliplr(value)
		ydata_inv = [value_inv[i].argmin() for i in range(0,value_inv.shape[0])] #arrayの一番右の要素番号

		self.xdata = xdata
		self.ydata = ydata
		self.ydata_inv = ydata_inv
	
		
	def Xy2r(self):
		#Rのxy座標を円でフィッティング
		
		#フィッティング(左側)
		cxe,cye,re = circlefitting(self.xdata,self.ydata)
		quad_optimal, covariance = scipy.optimize.curve_fit(quad, self.xdata, self.ydata) #得られたデータから2次関数でFITしてRの正負判定に使用
		#print(quad_optimal[0])
		if quad_optimal[0] > 0: # R>0
			fit_r2 = re*self.mmppx
			circle_optimal, circle_covariance = scipy.optimize.curve_fit(circle_m, self.xdata, self.ydata,np.array([re, cxe, cye]))
			fit_r = circle_optimal[0]*self.mmppx
			r_sigma = np.sqrt(circle_covariance[0][0])*self.mmppx
			print(fit_r)
			#print(fit_r2)
			print(r_sigma) #Rのフィッティング標準偏差
			
		else: # R<0
			fit_r2 = -re*self.mmppx
			circle_optimal, circle_covariance = scipy.optimize.curve_fit(circle_p, self.xdata, self.ydata,np.array([re, cxe, cye]))
			fit_r = -circle_optimal[0]*self.mmppx
			r_sigma = np.sqrt(circle_covariance[0][0])*self.mmppx
			print(fit_r)
			#print(fit_r2)
			print(r_sigma) #Rのフィッティング標準偏差
			
		#反転してフィッティング(右側)
		cxe_inv,cye_inv,re_inv = circlefitting(self.xdata,self.ydata_inv)
		quad_optimal_inv, covariance_inv = scipy.optimize.curve_fit(quad, self.xdata, self.ydata_inv)
		if quad_optimal_inv[0] > 0: # R<0
			fit_r2_inv = -re_inv*self.mmppx
			circle_optimal_inv, circle_covariance_inv = scipy.optimize.curve_fit(circle_m, self.xdata, self.ydata_inv, np.array([re_inv, cxe_inv, cye_inv]))
			fit_r_inv = -circle_optimal_inv[0]*self.mmppx
			r_sigma_inv = np.sqrt(circle_covariance_inv[0][0])*self.mmppx
			print(fit_r_inv)
			#print(fit_r2_inv)
			print(r_sigma_inv) #Rのフィッティング標準偏差
			
		else: # R>0
			fit_r2_inv = re*self.mmppx
			circle_optimal_inv, circle_covariance_inv = scipy.optimize.curve_fit(circle_p, self.xdata, self.ydata_inv, np.array([re_inv, cxe_inv, cye_inv]))
			fit_r_inv = circle_optimal_inv[0]*self.mmppx
			r_sigma_inv = np.sqrt(circle_covariance_inv[0][0])*self.mmppx
			print(fit_r_inv)
			#print(fit_r2_inv)
			print(r_sigma_inv) #Rのフィッティング標準偏差

		# ===== file,r,mmppx書き込み =====
		if len(self.current_dirname) == 0:
			if not os.path.isfile('_Rfit.csv'):
				f = open('_Rfit.csv', 'w', newline='')
				f.write('FileName,R_left,R_left_sigma,R_right,R_right_sigma,mm/px\n')
			else:
				f = open('_Rfit.csv', 'a', newline='')
		else:
			if not os.path.isfile(self.load_up_path + '\\_' + self.current_dirname + '_Rfit.csv'):
				f = open(self.load_up_path + '\\_' + self.current_dirname + '_Rfit.csv', 'w', newline = '')
				f.write('FileName,R_left,R_left_sigma,R_right,R_right_sigma,mm/px\n')
			else:
				f = open(self.load_up_path + '\\_' + self.current_dirname + '_Rfit.csv', 'a', newline = '')
				
		csvWriter = csv.writer(f)
		csvWriter.writerow([self.filename, fit_r, r_sigma, fit_r_inv, r_sigma_inv, self.mmppx])
		f.close()

if __name__ == '__main__':
	
	argvs = sys.argv
	argc = len(argvs)
	
	if argc == 1:
		sys.exit()
	elif argc > 3:
		sys.exit()
	
	if re.search('.jpg|.jpeg',argvs[1].lower()):
		if argc == 2:
			img = Rfit(argvs[1])
		else:
			img = Rfit(argvs[1],float(argvs[2]))

		img.Img2xy()
		img.Xy2r()
	
	elif re.search('.csv',argvs[1].lower()):
		if argc == 2:
			value = Rfit(argvs[1])
		else:
			value = Rfit(argvs[1],float(argvs[2]))

		value.Value2xy()
		value.Xy2r()
		
	else:
		jpg_list = glob.glob('%s/*.jp*' % argvs[1])
		
		if len(jpg_list) != 0:
			for i in jpg_list:
				if argc == 2:
					img = Rfit(i)
				else:
					img = Rfit(i,float(argvs[2]))

				img.Img2xy()
				img.Xy2r()
		else:
			print('Error : Input file/folder is not supported.')
			sys.exit()