ちょいめも

物理/Python/Cの雑記帳

エッジ像からMTFを計算

下のエッジ像の中心ラインだけ抽出して計算

f:id:choimemo:20211207003418j:plain
エッジ像
#第1引数:画像ファイルパス
import sys
import numpy as np
import matplotlib.pyplot as plt
from imageio import imread

#設定
pixel_pitch = 0.005 #mm 画素ピッチ lp/mmを算出するため


class edge2mtf():

	def __init__(self, load_img):
		self.img = imread(load_img).astype(np.float64)

	#ESF
	def ESF(self):
		center = int(self.img.shape[0] / 2) #画像縦方向中心
		self.esf = self.img[center] #中心の輝度値を抽出 本来は横ライン全て取得して疑似的にサンプリング増やした方が良いらしい
		self.esf_x = np.arange(self.esf.shape[0]) #x軸
		
	#LSF
	def LSF(self):
		self.lsf = np.diff(self.esf)
		self.lsf_x = np.arange(self.lsf.shape[0]) #x軸

	#MTF
	def MTF(self):
		self.otf = np.fft.fft(self.lsf)
		self.mtf = np.abs(self.otf) / np.abs(self.otf).max()
		self.mtf_x = np.fft.fftfreq(len(self.lsf_x)) #freq  x軸

		half_idx = int(len(self.mtf)/2)+1 #フーリエ変換は正負対象なので、正部分を抽出するためのインデックス
		self.mtf = self.mtf[0:half_idx]
		self.mtf_norm_x = self.mtf_x[0:half_idx] #pix/mm
		self.mtf_x = self.mtf_norm_x / pixel_pitch #lp/mm

	#グラフ描画
	def PlotGraph(self):

		self.ESF()
		self.LSF()
		self.MTF()

		fig, (esf_g, lsf_g, mtf_norm_g, mtf_g) = plt.subplots(nrows = 4, figsize=(5, 10), sharex = False)
		plt.subplots_adjust(wspace = 0.5, hspace = 0.5) # 余白を設定

		esf_g.plot(self.esf_x, self.esf)
		esf_g.set_title('ESF')
		esf_g.set_xlabel('x (px)')
		esf_g.set_ylabel('Brightness')
		esf_g.set_xlim(self.esf_x.min(), self.esf_x.max())
		
		lsf_g.plot(self.lsf_x, self.lsf)
		lsf_g.set_title('LSF')
		lsf_g.set_xlabel('x (px)')
		lsf_g.set_ylabel('LSF')
		lsf_g.set_xlim(self.lsf_x.min(), self.lsf_x.max())

		mtf_norm_g.plot(self.mtf_norm_x, self.mtf)
		mtf_norm_g.set_title('MTF Normalize')
		mtf_norm_g.set_xlabel('cycles / pixel')
		mtf_norm_g.set_ylabel('MTF')
		mtf_norm_g.set_xlim(self.mtf_norm_x.min(), self.mtf_norm_x.max())

		mtf_g.plot(self.mtf_x, self.mtf)
		mtf_g.set_title('MTF')
		mtf_g.set_xlabel('lp / mm')
		mtf_g.set_ylabel('MTF')
		mtf_g.set_xlim(self.mtf_x.min(), self.mtf_x.max())

		fig.show()
		

if __name__ == '__main__':

	#エラー処理なし
	argvs = sys.argv
	argc = len(argvs)

	img = edge2mtf(argvs[1])
	img.PlotGraph()


LSFのデータからも

#第1引数:LSFの入ったCSV
import sys
import numpy as np
import matplotlib.pyplot as plt

#設定
pixel_pitch = 0.005 #mm 画素ピッチ lp/mmを算出するため

class lsf2csv():

	#LSF読み込み
	def __init__(self, load_img):
		self.lsf = np.loadtxt(load_img).astype(np.float64)
		self.lsf_x = np.arange(self.lsf.shape[0]) #x軸

	#MTF
	def MTF(self):
		self.otf = np.fft.fft(self.lsf)
		self.mtf = np.abs(self.otf) / np.abs(self.otf).max()
		self.mtf_x = np.fft.fftfreq(len(self.lsf_x)) #freq  x軸

		half_idx = int(len(self.mtf)/2)+1 #フーリエ変換は正負対象なので、正部分を抽出するためのインデックス
		self.mtf = self.mtf[0:half_idx]
		self.mtf_norm_x = self.mtf_x[0:half_idx] #pix/mm
		self.mtf_x = self.mtf_norm_x / pixel_pitch #lp/mm

	#グラフ描画
	def PlotGraph(self):

		self.MTF()

		fig, (lsf_g, mtf_norm_g, mtf_g) = plt.subplots(nrows = 3, figsize=(5, 10), sharex = False)
		plt.subplots_adjust(wspace = 0.5, hspace = 0.5) # 余白を設定
		
		lsf_g.plot(self.lsf_x, self.lsf)
		lsf_g.set_title('LSF')
		lsf_g.set_xlabel('x (px)')
		lsf_g.set_ylabel('LSF')
		lsf_g.set_xlim(self.lsf_x.min(), self.lsf_x.max())

		mtf_norm_g.plot(self.mtf_norm_x, self.mtf)
		mtf_norm_g.set_title('MTF Normalize')
		mtf_norm_g.set_xlabel('cycles / pixel')
		mtf_norm_g.set_ylabel('MTF')
		mtf_norm_g.set_xlim(self.mtf_norm_x.min(), self.mtf_norm_x.max())

		mtf_g.plot(self.mtf_x, self.mtf)
		mtf_g.set_title('MTF')
		mtf_g.set_xlabel('lp / mm')
		mtf_g.set_ylabel('MTF')
		mtf_g.set_xlim(self.mtf_x.min(), self.mtf_x.max())

		fig.show()
		

if __name__ == '__main__':

	#エラー処理なし
	argvs = sys.argv
	argc = len(argvs)

	img = lsf2csv(argvs[1])
	img.PlotGraph()

あってんのかな?