diff --git a/BuildGUI.py b/BuildGUI.py new file mode 100644 index 0000000..a7b20d5 --- /dev/null +++ b/BuildGUI.py @@ -0,0 +1,167 @@ +import numpy as np +import os +from PIL import Image +import cvui +from UtilProcess import convert_tiff2numpy, convert16to8bit, resize_square +import cv2 + +from UtilsGUI import DisplayClock +from SignalProcess import Roi2Signal +from ImageProcess import ImageProcess +from rrEstimation import rrEstimation +from timeit import default_timer as timer +from rqiCalculation import rqiCalculation +from yoloV3.detect import YOLO +from GLOBAL import * +from flirpy.camera.boson import Boson + + +WINDOW_NAME = "RespThermal" +time_for_analysis = 15 + +cycle = 60/RR + +dt = cycle/8 + + +class PositionInGUI: + def __init__(self, x, y, width=0, height=0, size=0.0): + self.x = x + self.y = y + self.w = width + self.h = height + self.size = size + + +def EstimateRespiratoryRate(yolo): + camera = Boson() + + # GUI上のパーツのPositionを定義 + p_fps_text = PositionInGUI(90, 15) + p_detection_frame = PositionInGUI(40, 60) + p_start_button = PositionInGUI(10, 10) + p_rr_text = PositionInGUI(80, 20, 0, 0, 0.7) + + # cvuiの初期化 + window_frame = np.zeros((400, 500, 3), np.uint8) + cvui.init(WINDOW_NAME) + + is_on = False + accum_time = 0 + curr_fps = 0 + fps = "FPS: ??" + rr_text = "Resp. rate = ??" + mouse_rr_text = "Resp. rate = ??" + prev_time = timer() + + respiration_count = 0 + + # init Signal class for each grid　 + signals_for_each_grid = [] + for i in range(NUM_WIDTH * NUM_HEIGHT): + signals_for_each_grid.append(Roi2Signal()) + + # init ImageProcess class + image_process_class = ImageProcess(dst_dir=None) + + estimated_rrs = [] + + while True: + window_frame[:] = (49, 52, 49) # フレームの色を指定(BGR) + + if cvui.button(window_frame, p_start_button.x, p_start_button.y, "Start"): + # ON/OFFを切り替える + is_on = not is_on + start_time = timer() + + # ボタンが押されたら、鼻の検出を開始 + if is_on: + curr_time = timer() + exec_time = curr_time - prev_time + prev_time = curr_time + accum_time = accum_time + exec_time + curr_fps = curr_fps + 1 + if accum_time > 1: + accum_time = accum_time - 1 + fps = curr_fps + curr_fps = 0 + + if curr_time-start_time > 60: + break + + if curr_time-start_time > cycle * respiration_count: + respiration_count = respiration_count + 1 + + # preprocess image + u16_image = camera.grab(0) + u16_image = np.stack([u16_image, u16_image, u16_image], axis=-1) + u8_image = convert16to8bit(u16_image) + u16_resized_image = resize_square(u16_image) + u8_resized_image = resize_square(u8_image) + + # detect face by YOLO v3 + detection_result, face_box = yolo.detect_face(Image.fromarray(u8_resized_image), + display_score=False) + detection_image = np.asarray(detection_result) + + # calculate ROI + face_height, face_width = (face_box[2] - face_box[0], face_box[3] - face_box[1]) + grid_height, grid_width = face_height // NUM_HEIGHT, face_width // NUM_WIDTH + + # define each grid + block_boxes = [] + for j in range(NUM_HEIGHT): + for k in range(NUM_WIDTH): + top, left = (grid_height * j + face_box[0], grid_width * k + face_box[1]) + bottom, right = (grid_height + top, grid_width + left) + block_boxes.append(np.array([top, left, bottom, right])) + + time = curr_time-start_time + + # calculate signal intensity for each grid + for j in range(NUM_WIDTH * NUM_HEIGHT): + image_process_class.load(u16_resized_image, u8_resized_image, + detection_image, block_boxes[j]) + signals_for_each_grid[j].append_signal(time, image_process_class.get_signal_intensity()) + + # 一定時間(15s)以上信号がたまったら，呼吸数を計算 + signal_size = signals_for_each_grid[0].get_size() + window_size = int(fps*time_for_analysis) + if fps > 7 and signal_size > window_size: + f_rrs = [] + t_rrs = [] + rqis = [] + d_rrs = [] + for i in range(NUM_WIDTH * NUM_HEIGHT): + respClass = rrEstimation(signals_for_each_grid[i], dst_dir=None, + sampling_rate=fps, window=False) + respClass.calculate_PSD() + + f_rrs.append(respClass.estimate_f_rr()) + t_rrs.append(respClass.estimate_t_rr(size=MV_SIZE, iteration=MV_ITER)) + d_rrs.append(OUTLIER if t_rrs[-1] == OUTLIER + else abs(t_rrs[-1] - f_rrs[-1])) + rqiClass = rqiCalculation(respClass, RQI_PARAM, d_rrs[-1], PENALTY_ID) + rqis.append(rqiClass.get_rqi(rqi_id=RQI_ID)) + + # select grid and rr + estimated_place = rqis.index(max(rqis)) + estimated_rr = f_rrs[estimated_place] + + if signal_size % 18 == 0: + rr_text = "Resp. rate = " + str(round(estimated_rr)) + " bpm" + estimated_rrs.append(estimated_rr) + + cvui.text(window_frame, p_fps_text.x, p_fps_text.y, "FPS:" + str(fps), p_fps_text.size) + cvui.image(window_frame, p_detection_frame.x, p_detection_frame.y, detection_image) + cvui.text(window_frame, p_rr_text.x, p_rr_text.y, rr_text, p_rr_text.size) + + cvui.imshow(WINDOW_NAME, window_frame) + + if cv2.waitKey(20) == 27: + break + + +if __name__ == '__main__': + EstimateRespiratoryRate(YOLO(**FLAGS)) + diff --git a/EstimateRRbyFaceRQI.py b/EstimateRRbyFaceRQI.py new file mode 100644 index 0000000..f3ac9d2 --- /dev/null +++ b/EstimateRRbyFaceRQI.py @@ -0,0 +1,163 @@ +import numpy as np +import os +from PIL import Image +from glob import glob +from LoadInputData import load_input_data +from UtilProcess import convert_tiff2numpy, convert16to8bit, resize_square +from SignalProcess import Roi2Signal +from ImageProcess import ImageProcess +from rrEstimation import rrEstimation +from rqiCalculation import rqiCalculation +from yoloV3.detect import YOLO +from GLOBAL import * + + +def estimateRR_usingRQI(_src_dir, _dst_dir, yolo): + # load dataset info + _gt_rr, thermal_fps, _ = load_input_data(os.path.join(_src_dir, "info.txt")) + + # init Signal class for each grid　 + signals_for_each_grid = [] + for i in range(NUM_WIDTH*NUM_HEIGHT): + signals_for_each_grid.append(Roi2Signal()) + + # init ImageProcess class + image_process_class = ImageProcess(dst_dir=_dst_dir) + + # load thermal image(16bit) + u16_images = convert_tiff2numpy(os.path.join(_src_dir, "video.tiff")) + + # ----------------------------------------------- + # calculate signal intensity + for i, u16_image in enumerate(u16_images): + # preprocess image + u16_image = np.stack([u16_image, u16_image, u16_image], axis=-1) + u8_image = convert16to8bit(u16_image) + u16_resized_image = resize_square(u16_image) + u8_resized_image = resize_square(u8_image) + + # detect face by YOLO v3 + detection_result, face_box = yolo.detect_face(Image.fromarray(u8_resized_image), + display_score=False) + detection_image = np.asarray(detection_result) + + # calculate ROI + face_height, face_width = (face_box[2] - face_box[0], face_box[3] - face_box[1]) + grid_height, grid_width = face_height // NUM_HEIGHT, face_width // NUM_WIDTH + + # define each grid + block_boxes = [] + for j in range(NUM_HEIGHT): + for k in range(NUM_WIDTH): + top, left = (grid_height * j + face_box[0], grid_width * k + face_box[1]) + bottom, right = (grid_height + top, grid_width + left) + block_boxes.append(np.array([top, left, bottom, right])) + + time = i/thermal_fps + + # calculate signal intensity for each grid + for j in range(NUM_WIDTH * NUM_HEIGHT): + image_process_class.load(u16_resized_image, u8_resized_image, + detection_image, block_boxes[j]) + signals_for_each_grid[j].append_signal(time, image_process_class.get_signal_intensity()) + + # ----------------------------------------------- + # estimate Respiratory Rate + f_rrs = [] + t_rrs = [] + rqis = [] + d_rrs = [] + # estimate RR in time domain, and frequency domain + # estimate rqis + for i in range(NUM_WIDTH*NUM_HEIGHT): + respClass = rrEstimation(signals_for_each_grid[i], dst_dir=_dst_dir, + sampling_rate=thermal_fps, window=False) + respClass.calculate_PSD() + + f_rrs.append(respClass.estimate_f_rr()) + t_rrs.append(respClass.estimate_t_rr(size=MV_SIZE, iteration=MV_ITER)) + d_rrs.append(OUTLIER if t_rrs[-1] == OUTLIER + else abs(t_rrs[-1]-f_rrs[-1])) + rqiClass = rqiCalculation(respClass, RQI_PARAM, d_rrs[-1], PENALTY_ID) + rqis.append(rqiClass.get_rqi(rqi_id=RQI_ID)) + respClass.write_signal2png(str(i)+'_signal.png', _dst_dir) + respClass.write_psd2png(str(i)+'_psd.png', _dst_dir) + + # select grid and rr + estimated_place = rqis.index(max(rqis)) + print(estimated_place) + estimated_rr = f_rrs[estimated_place] + + print(os.path.basename(_dst_dir), _gt_rr, estimated_rr) + # ---------------------------------------------------------------- + # save result + + return _gt_rr, estimated_rr, estimated_place + + +if __name__ == '__main__': + _yolo = YOLO(**FLAGS) + results = [] + src_dirs = glob(SRC_DIRS) + i = 0 + for src_dir in src_dirs: + base_name = os.path.basename(src_dir) + dst_dir = os.path.join(DST_DIRS, base_name) + + if 'ishikawa' not in dst_dir: + continue + + # print(base_name) + if not os.path.exists(dst_dir): + os.mkdir(dst_dir) + + gt_rr, rr, index = estimateRR_usingRQI(src_dir, dst_dir, _yolo) + abs_error = np.round(abs(gt_rr - rr), 2) + results.append([base_name, gt_rr, rr, index, abs_error]) + i += 1 + + print("base_name, GT, rr, index, AE") + AEs = [] + for result in results: + print(result[0], result[1], result[2], result[3], result[4]) + AEs.append(result[4]) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/EstimateRRbyNoseMouth.py b/EstimateRRbyNoseMouth.py new file mode 100644 index 0000000..af66440 --- /dev/null +++ b/EstimateRRbyNoseMouth.py @@ -0,0 +1,126 @@ +import numpy as np +import os +from PIL import Image +from glob import glob +from LoadInputData import load_input_data +from UtilProcess import convert_tiff2numpy, convert16to8bit, resize_square +from SignalProcess import Roi2Signal +from ImageProcess import ImageProcess +from rrEstimation import rrEstimation +from rqiCalculation import rqiCalculation +from yoloV3.detect import YOLO +from GLOBAL import * + + +def estimateRR_usingNM(_src_dir, _dst_dir, yolo): + # load dataset info + _gt_rr, thermal_fps, _ = load_input_data(os.path.join(_src_dir, "info.txt")) + + # init ImageProcess class + nose_roi = ImageProcess(dst_dir=_dst_dir) + mouth_roi = ImageProcess(dst_dir=_dst_dir) + + # init SignalClass + nose_signal = Roi2Signal() + mouth_signal = Roi2Signal() + + # load thermal image(16bit) + u16_images = convert_tiff2numpy(os.path.join(_src_dir, "video.tiff")) + + # ----------------------------------------------- + # calculate signal intensity + for i, u16_image in enumerate(u16_images): + # preprocess image + u16_image = np.stack([u16_image, u16_image, u16_image], axis=-1) + u8_image = convert16to8bit(u16_image) + u16_resized_image = resize_square(u16_image) + u8_resized_image = resize_square(u8_image) + + # detect nose and mouth by YOLO v3 + detection_result, nose_box, mouth_box, face_box = yolo.detect_fromConf(Image.fromarray(u8_resized_image)) + + detection_image = np.asarray(detection_result) + + time = i/thermal_fps + + # calculate signal intensity for each grid + nose_roi.load(u16_resized_image, u8_resized_image, detection_image, nose_box) + mouth_roi.load(u16_resized_image, u8_resized_image, detection_image, mouth_box) + + nose_signal.append_signal(time, nose_roi.get_signal_intensity()) + mouth_signal.append_signal(time, mouth_roi.get_signal_intensity()) + + # ----------------------------------------------- + # estimate Respiratory Rate + # estimate RR in time domain, and frequency domain + # estimate rqis + resp_nose = rrEstimation(nose_signal, dst_dir=_dst_dir, sampling_rate=thermal_fps, window=False) + resp_mouth = rrEstimation(mouth_signal, dst_dir=_dst_dir, sampling_rate=thermal_fps, window=False) + _nose_rr = resp_nose.estimate_f_rr() + _mouth_rr = resp_mouth.estimate_f_rr() + + return _gt_rr, _nose_rr, _mouth_rr + + +if __name__ == '__main__': + _yolo = YOLO(**FLAGS) + results = [] + src_dirs = glob(SRC_DIRS) + i = 0 + for src_dir in src_dirs: + base_name = os.path.basename(src_dir) + dst_dir = os.path.join(DST_DIRS, base_name) + # print(base_name) + if not os.path.exists(dst_dir): + os.mkdir(dst_dir) + + gt_rr, nose_rr, mouth_rr = estimateRR_usingNM(src_dir, dst_dir, _yolo) + abs_error = np.round(abs(gt_rr - rr), 2) + results.append([base_name, gt_rr, rr, index, abs_error]) + i += 1 + + print("base_name, GT, rr, index, AE") + AEs = [] + for result in results: + print(result[0], result[1], result[2], result[3], result[4]) + AEs.append(result[4]) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/GLOBAL.py b/GLOBAL.py new file mode 100644 index 0000000..fc5c1c4 --- /dev/null +++ b/GLOBAL.py @@ -0,0 +1,37 @@ +# in rrEstimation +LF_THRESH = 0.10 +HF_THRESH = 2 + +RATIO_BF = 0.25 +RATIO_HF = 0.1 + +# ------------------------------------------- +# in EstimateRRbyFaceRQI +RQI_ID = 6 +PENALTY_ID = 1 +NUM_WIDTH, NUM_HEIGHT = 4, 6 +MV_SIZE, MV_ITER = 5, 3 +OUTLIER = 100 +RQI_PARAM = 5 +FLAGS = {'image': False, 'input': './path2your_video', 'output': ''} +SRC_DIRS = r"I:\3.Data\RespData\yolo_src2\*" +DST_DIRS = r"I:\3.Data\RespData\yolo_dst_sqi4" +DST_DIRS = r"C:\Users\takah\OneDrive - 千葉大学\Sensors of Journal\Image_New" + +# -------------------------------------------- +# in YOLOv3 +CLASS_NAMES = ["face", "mouth", "nose"] +CLASS_COLORS = {"face": (255, 0, 0), "mouth": (0, 255, 0), "nose": (0, 0, 255)} +# MODEL_PATH = r"weights/weights.h5" +MODEL_PATH = r"I:\3.Data\RespData\weights\01_08_2\FLIR2.h5" +ANCHORS_PATH = r'model_data/yolo_anchors.txt' +CLASSES_PATH = r"model_data/classes.txt" +SCORE = 0.25 +IOU = 0.45 +MODEL_SIZE_IMAGE = (320, 320) +GPU_NUM = 1 + + + + + diff --git a/ImageProcess.py b/ImageProcess.py new file mode 100644 index 0000000..efe6eb6 --- /dev/null +++ b/ImageProcess.py @@ -0,0 +1,60 @@ +import numpy as np +import cv2 +import os + + +class ImageProcess: + def __init__(self, dst_dir): + self.__u16_images = [] + self.__u8_images = [] + self.__detection_images = [] + self.__mask_images = [] + + self.__boxes = [] + self.__dst_dir = dst_dir + # TODO リアルタイム推定の時などに使えるように可変にする + self.__frame_rate = 8.6 + + def load(self, image_u16, image_u8, detection_result, box): + self.__u16_images.append(image_u16) + self.__u8_images.append(image_u8) + self.__detection_images.append(detection_result) + self.__boxes.append(box) + + def extract_region(self, top, left, bottom, right): + interest_region = self.__u8_images[-1][top:bottom, left:right, :] + ret, mask = cv2.threshold(interest_region[:, :, 0], 0, 255, cv2.THRESH_OTSU) + mask = cv2.bitwise_not(mask) + return mask + + def get_signal_intensity(self): + if np.all(self.__boxes[-1] == 0): + raise Exception("検出が正常に行われませんでした") + mouse_box = self.__boxes[-1] + top, left, bottom, right = np.array([int(mouse_box[i]) + for i in range(mouse_box.shape[0])]) + + interest_region = self.__u16_images[-1][top:bottom, left:right, :] + intensity = np.mean(interest_region) + # mask = self.extract_region(top, left, bottom, right) + # interest_region = interest_region[:, :, 0] * mask + # intensity = np.sum(interest_region) / np.count_nonzero(interest_region > 0) + + return intensity + + def save_detection_result(self, file_path): + output_path = os.path.join(self.__dst_dir, file_path) + print(output_path) + size = (self.__detection_images[0].shape[1], self.__detection_images[0].shape[0]) + + fmt = cv2.VideoWriter_fourcc('m', 'p', '4', 'v') + writer = cv2.VideoWriter(output_path, fmt, self.__frame_rate, size) # ライター作成 + for frame in self.__detection_images: + writer.write(frame) # 画像を1フレーム分として書き込み + # ファイルを閉じる + writer.release() + + + + + diff --git a/LoadInputData.py b/LoadInputData.py new file mode 100644 index 0000000..15952ea --- /dev/null +++ b/LoadInputData.py @@ -0,0 +1,58 @@ +import numpy as np + + +def load_input_data(path): + """ + 実験データの読み込み + :param path: ぱす + :return: RR(Ground Truth), サーモグラフィのfps, Bitalinoのfps + """ + gt_rr = None + thermal_fps = None + bitalino_fps = None + + with open(path, 'r') as f: + s = f.read() + split_sentence = s.split('\n') + for i in range(len(split_sentence)): + temp = split_sentence[i].split('=') + if len(temp) < 2: + continue + if temp[0] in 'estimated fps': + thermal_fps = float(temp[1]) + elif temp[0] in 'revised_RR': + gt_rr = float(temp[1]) + elif temp[0] in 'Sampling Rate': + bitalino_fps = float(temp[1]) + + if (gt_rr is None) or (thermal_fps is None) or (bitalino_fps is None): + raise Exception("入力データセットが異なります") + + else: + return gt_rr, thermal_fps, bitalino_fps + + +def load_signals_fromCSV(csv_path): + csv_list = [] + with open(csv_path, 'r') as f: + line = f.readline() + while line: + if not line == '\n': + csv_list.append(line.replace('\n', '')) + line = f.readline() + + times = np.zeros(len(csv_list)) + signals = np.zeros(len(csv_list)) + + for i in range(len(csv_list)): + split_csv = csv_list[i].split(',') + times[i] = split_csv[0] + signals[i] = split_csv[1] + + return times, signals + + + + + + diff --git a/NonlinearFunction.py b/NonlinearFunction.py new file mode 100644 index 0000000..4f6ef00 --- /dev/null +++ b/NonlinearFunction.py @@ -0,0 +1,23 @@ +from math import e, tanh, atan, pi + + +def sigmoid(x): + y = 1/(1+e**(-x)) + return y + + +def arc_tan(x): + y = atan(x)/pi + 1/2 + return y + + +def step(x): + if x>=0: + return 1 + else : + return 0 + + +def tanh_01(x): + y = (tanh(x)+1)/2 + return y diff --git a/README.md b/README.md new file mode 100644 index 0000000..10bf5a2 --- /dev/null +++ b/README.md @@ -0,0 +1,19 @@ +# What is RespThermal : +- サーモグラフィ画像から呼吸数を推定するシステム + - 顔領域から呼吸数を推定 : ```EstimateRRbyFaceRQI.py``` + - 口・鼻領域から呼吸数を推定(Respiratory Quality Index) : ```EstimateRRbyNoseMouth.py``` + +- ```GLOBAL.py```でパラメタを管理. 適宜変更してください. +- 重み```weights.weights.h5```は共有ファイルからダウンロードしてください. +# Environment : +- 基本的には[keras-yolov3](https://github.com/qqwweee/keras-yolo3) が動く環境がベース. 自分はanacondaで環境構築. +- その他細かいライブラリはrequirement.yamlを参照. ただ別プロジェクト用のパッケージも入ってるのであしからず． + - サーモグラフィ用ライブラリ : ```pip install flirpy``` + + +# Demonstration : +- サーモグラフィカメラはBoson320を使用 +- リアルタイム推論モード : ```BuildGUI.py``` + +![デモ動画](Images/RespThermal.gif) + diff --git a/SignalProcess.py b/SignalProcess.py new file mode 100644 index 0000000..53b97bb --- /dev/null +++ b/SignalProcess.py @@ -0,0 +1,78 @@ +import numpy as np +from GLOBAL import * + + +class Roi2Signal: + def __init__(self): + self._times = [] + self._signals = [] + + def append_signal(self, time, signal): + """ + Detection結果から算出された信号強度を保存 + (Detection失敗により)信号強度が算出されなかった場合, 0-paddingするか前回の信号強度を参照 + :param time: 時間 + :param signal: 信号強度 + :return: None + """ + if len(self._times) == 0: + self._signals.append(signal if signal is not None else 0) + else: + self._signals.append(signal if signal is not None + else self._signals[-1]) + + self._times.append(time) + + def get_size(self): + return len(self._times) + + def load_signals(self, times, signals): + """ + 時系列の信号を読み込み + :param times: 時間データ + :param signals: 信号強度データ + :return: None + """ + self._times = times + self._signals = signals + + def get_times(self, window=False): + """ + 時間データをNumpyに変換してreturn + window=Falseの時すべてreturn, window=intの時窓の長さ分return + :param window: 窓の長さ + :return: 時間データ(Numpy) + """ + if window: + return np.array(self._times[-window:]) + else: + return np.array(self._times) + + def get_signal(self, window=False): + """ + 信号データをNumpyに変換してreturn + window=Falseの時すべてreturn, window=intの時窓の長さ分return + :param window: 窓の長さ + :return: 信号データ + """ + if window: + return np.array(self._signals[-window:]) + else: + return np.array(self._signals) + + def apply_moving_average_filter(self, num=5, window=False): + """ + 移動平均 + :param num:移動平均 + :param window: 窓の長さ + :return: 平滑化された信号 + """ + b = np.ones(num)/float(num) + + if window: + windowed_signals = np.array(self._signals[-window:]) + else: + windowed_signals = np.array(self._signals) + filtered_signals = np.convolve(windowed_signals, b, mode="same") + + return filtered_signals diff --git a/UtilProcess.py b/UtilProcess.py new file mode 100644 index 0000000..d167083 --- /dev/null +++ b/UtilProcess.py @@ -0,0 +1,37 @@ +import tifffile +import numpy as np +import cv2 + + +def convert_tiff2numpy(_src_path): + """[frame, height, width]の形であることに注意""" + np_images = tifffile.imread(_src_path) + # print(np_images.shape) + # print(np_images.dtype) + + return np_images + + +def convert16to8bit(images_u16): + maxV, minV = np.amax(images_u16), np.amin(images_u16) + # maxV, minV = (60000, 50000) + # print(maxV, minV) + alpha = 255.0 / (maxV - minV) + images_u8 = np.add(images_u16, -minV) + images_u8 = images_u8 * alpha + images_u8 = images_u8.astype(np.uint8) + return images_u8 + + +def resize_square(img, size=(320, 320)): + height, width, _ = img.shape + if height < width: + margin = (width-height) // 2 + reshaped_img = img[:, margin:-margin, :] + else: + margin = (height - width) // 2 + reshaped_img = img[margin:-margin, :, :] + + reshaped_img = cv2.resize(reshaped_img, dsize=size) + + return reshaped_img diff --git a/UtilsGUI.py b/UtilsGUI.py new file mode 100644 index 0000000..284f660 --- /dev/null +++ b/UtilsGUI.py @@ -0,0 +1,42 @@ +import cvui + + +inspiration3 = 'breath in 3' +inspiration2 = 'breath in 2' +inspiration1 = 'breath in 1' +inspiration0 = 'breath in 0' + +exhalation3 = 'breath out 3' +exhalation2 = 'breath out 2' +exhalation1 = 'breath out 1' +exhalation0 = 'breath out 0' + +inspiration_color = 0x00ffff +exhalation_color = 0xffd700 + + +def DisplayClock(displayed_time, window_frame, p_instruction_text, dt): + if displayed_time < dt: + cvui.text(window_frame, p_instruction_text.x, p_instruction_text.y, + inspiration0, p_instruction_text.size, inspiration_color) + elif displayed_time < dt * 2: + cvui.text(window_frame, p_instruction_text.x, p_instruction_text.y, + inspiration1, p_instruction_text.size, inspiration_color) + elif displayed_time < dt * 3: + cvui.text(window_frame, p_instruction_text.x, p_instruction_text.y, + inspiration2, p_instruction_text.size, inspiration_color) + elif displayed_time < dt * 4: + cvui.text(window_frame, p_instruction_text.x, p_instruction_text.y, + inspiration3, p_instruction_text.size, inspiration_color) + elif displayed_time < dt * 5: + cvui.text(window_frame, p_instruction_text.x, p_instruction_text.y, + exhalation0, p_instruction_text.size, exhalation_color) + elif displayed_time < dt * 6: + cvui.text(window_frame, p_instruction_text.x, p_instruction_text.y, + exhalation1, p_instruction_text.size, exhalation_color) + elif displayed_time < dt * 7: + cvui.text(window_frame, p_instruction_text.x, p_instruction_text.y, + exhalation2, p_instruction_text.size, exhalation_color) + else: + cvui.text(window_frame, p_instruction_text.x, p_instruction_text.y, + exhalation3, p_instruction_text.size, exhalation_color) \ No newline at end of file diff --git a/requirement.yaml b/requirement.yaml new file mode 100644 index 0000000..45fadef --- /dev/null +++ b/requirement.yaml @@ -0,0 +1,173 @@ +name: py36 +channels: + - conda-forge + - defaults +dependencies: + - _tflow_select=2.1.0=gpu + - absl-py=0.11.0=py36haa95532_0 + - argon2-cffi=20.1.0=py36he774522_1 + - asn1crypto=1.4.0=py_0 + - astor=0.8.1=py36_0 + - async_generator=1.10=py36h28b3542_0 + - attrs=20.3.0=pyhd3eb1b0_0 + - backcall=0.2.0=py_0 + - blas=1.0=mkl + - bleach=3.2.1=py_0 + - brotlipy=0.7.0=py36he774522_1000 + - ca-certificates=2020.6.20=hecda079_0 + - certifi=2020.6.20=py36hd36e781_2 + - cffi=1.14.3=py36h7a1dbc1_0 + - chardet=3.0.4=py36_1003 + - colorama=0.4.4=py_0 + - cryptography=2.3.1=py36h74b6da3_0 + - cudatoolkit=10.0.130=0 + - cudnn=7.6.5=cuda10.0_0 + - cycler=0.10.0=py36h009560c_0 + - cython=0.29.21=py36ha925a31_0 + - decorator=4.4.2=py_0 + - defusedxml=0.6.0=py_0 + - entrypoints=0.3=py36_0 + - enum34=1.1.10=py36h9f0ad1d_2 + - freetype=2.10.4=hd328e21_0 + - gast=0.4.0=py_0 + - grpcio=1.14.1=py36h5c4b210_0 + - h5py=2.10.0=py36h5e291fa_0 + - hdf5=1.10.4=h7ebc959_0 + - icc_rt=2019.0.0=h0cc432a_1 + - icu=58.2=ha925a31_3 + - idna=2.10=py_0 + - imagecodecs-lite=2019.12.3=py36h4f3e613_3 + - importlib-metadata=2.0.0=py_1 + - importlib_metadata=2.0.0=1 + - intel-openmp=2020.1=216 + - ipykernel=5.3.4=py36h5ca1d4c_0 + - ipython=7.16.1=py36h5ca1d4c_0 + - ipython_genutils=0.2.0=py36_0 + - jedi=0.17.2=py36_0 + - jinja2=2.11.2=py_0 + - jpeg=9d=he774522_0 + - json5=0.9.5=py_0 + - jsonschema=3.2.0=py_2 + - jupyter_client=6.1.7=py_0 + - jupyter_core=4.6.3=py36_0 + - jupyterlab=2.2.8=py_0 + - jupyterlab_pygments=0.1.2=py_0 + - jupyterlab_server=1.2.0=py_0 + - keras=2.2.4=0 + - keras-applications=1.0.8=py_1 + - keras-base=2.2.4=py36_0 + - keras-preprocessing=1.1.0=py_1 + - kiwisolver=1.3.0=py36hd77b12b_0 + - libblas=3.8.0=16_mkl + - libcblas=3.8.0=16_mkl + - liblapack=3.8.0=16_mkl + - liblapacke=3.8.0=16_mkl + - libpng=1.6.37=h2a8f88b_0 + - libprotobuf=3.13.0.1=h200bbdf_0 + - libsodium=1.0.18=h62dcd97_0 + - libtiff=4.1.0=h56a325e_1 + - libwebp=1.0.1=he774522_0 + - lz4-c=1.9.2=hf4a77e7_3 + - m2w64-gcc-libgfortran=5.3.0=6 + - m2w64-gcc-libs=5.3.0=7 + - m2w64-gcc-libs-core=5.3.0=7 + - m2w64-gmp=6.1.0=2 + - m2w64-libwinpthread-git=5.0.0.4634.697f757=2 + - markdown=3.3.2=py36_0 + - markupsafe=1.1.1=py36he774522_0 + - matplotlib=3.3.2=0 + - matplotlib-base=3.3.2=py36hba9282a_0 + - mistune=0.8.4=py36he774522_0 + - mkl=2020.1=216 + - mkl-service=2.3.0=py36hb782905_0 + - mkl_fft=1.2.0=py36h45dec08_0 + - mkl_random=1.1.1=py36h47e9c7a_0 + - msys2-conda-epoch=20160418=1 + - natsort=7.0.1=py_0 + - nbclient=0.5.1=py_0 + - nbconvert=6.0.7=py36_0 + - nbformat=5.0.8=py_0 + - nest-asyncio=1.4.1=py_0 + - notebook=6.1.4=py36_0 + - numpy=1.19.2=py36hadc3359_0 + - numpy-base=1.19.2=py36ha3acd2a_0 + - olefile=0.46=py36_0 + - opencv=4.0.1=py36hce2de41_201 + - openssl=1.0.2u=hfa6e2cd_0 + - packaging=20.4=py_0 + - pandas=1.1.1=py36ha925a31_0 + - pandoc=2.11=h9490d1a_0 + - pandocfilters=1.4.2=py36_1 + - parso=0.7.0=py_0 + - pathlib=1.0.1=py36h9f0ad1d_3 + - pickleshare=0.7.5=py36_0 + - pillow=8.0.1=py36h4fa10fc_0 + - pip=20.2.4=py36haa95532_0 + - prometheus_client=0.8.0=py_0 + - prompt-toolkit=3.0.8=py_0 + - protobuf=3.13.0.1=py36ha925a31_1 + - pycparser=2.20=py_2 + - pygments=2.7.2=pyhd3eb1b0_0 + - pyopenssl=19.0.0=py36_0 + - pyparsing=2.4.7=py_0 + - pyqt=5.6.0=py36ha878b3d_6 + - pyreadline=2.1=py36_1 + - pyrsistent=0.17.3=py36he774522_0 + - pysocks=1.7.1=py36_0 + - python=3.6.10=h9f7ef89_2 + - python-dateutil=2.8.1=py_0 + - python_abi=3.6=1_cp36m + - pytz=2020.1=py_0 + - pywin32=227=py36he774522_1 + - pywinpty=0.5.7=py36_0 + - pyyaml=5.3.1=py36he774522_1 + - pyzmq=19.0.2=py36ha925a31_1 + - qt=5.6.2=vc14h6f8c307_12 + - requests=2.24.0=py_0 + - scipy=1.5.2=py36h9439919_0 + - send2trash=1.5.0=py36_0 + - setuptools=49.6.0=py36_1 + - sip=4.18.1=py36h6538335_2 + - six=1.15.0=py_0 + - sqlite=3.33.0=h2a8f88b_0 + - tensorboard=1.14.0=py36he3c9ec2_0 + - tensorflow=1.14.0=gpu_py36h305fd99_0 + - tensorflow-base=1.14.0=gpu_py36h55fc52a_0 + - tensorflow-estimator=1.14.0=py_0 + - tensorflow-gpu=1.14.0=h0d30ee6_0 + - termcolor=1.1.0=py36_1 + - terminado=0.9.1=py36_0 + - testpath=0.4.4=py_0 + - tifffile=2019.7.26.2=py36_0 + - tk=8.6.10=he774522_0 + - tornado=6.0.4=py36he774522_1 + - traitlets=4.3.3=py36_0 + - urllib3=1.25.11=py_0 + - vc=14.1=h0510ff6_4 + - vs2015_runtime=14.16.27012=hf0eaf9b_3 + - wcwidth=0.2.5=py_0 + - webencodings=0.5.1=py36_1 + - werkzeug=1.0.1=py_0 + - wheel=0.35.1=py_0 + - win_inet_pton=1.1.0=py36_0 + - wincertstore=0.2=py36h7fe50ca_0 + - winpty=0.4.3=4 + - wrapt=1.12.1=py36he774522_1 + - xz=5.2.5=h62dcd97_0 + - yaml=0.2.5=he774522_0 + - zeromq=4.3.2=ha925a31_3 + - zipp=3.4.0=pyhd3eb1b0_0 + - zlib=1.2.11=h62dcd97_4 + - zstd=1.4.5=h04227a9_0 + - pip: + - cvui==2.7 + - flirpy==0.2.3 + - psutil==5.7.2 + - pycocotools==2.0 + - pydicom==2.0.0 + - pyftdi==0.51.2 + - pyserial==3.4 + - pyudev==0.22.0 + - pyusb==1.1.0 + - tqdm==4.50.2 +prefix: H:\Anaconda\envs\py36 diff --git a/rqiCalculation.py b/rqiCalculation.py new file mode 100644 index 0000000..4a99185 --- /dev/null +++ b/rqiCalculation.py @@ -0,0 +1,120 @@ +from rrEstimation import rrEstimation +import numpy as np +from NonlinearFunction import sigmoid, arc_tan, step, tanh_01 +from GLOBAL import * + + +class rqiCalculation: + def __init__(self, rr_class: rrEstimation, a, d_rr, penalty_id=0): + self._frequency = rr_class.get_frequency() + self._spectrum = rr_class.get_spectrum() + self.lp_index = np.abs(self._frequency - LF_THRESH).argmin() + self.hp_index = np.abs(self._frequency - HF_THRESH).argmin() + + self._bp_spectrum = self._spectrum[self.lp_index:self.hp_index] + self._hp_spectrum = self._spectrum[self.hp_index:] + self.d_rr = d_rr + + if penalty_id == 0: + self._penalty = sigmoid(a - d_rr) + elif penalty_id == 1: + self._penalty = tanh_01(a - d_rr) + elif penalty_id == 2: + self._penalty = step(a - d_rr) + elif penalty_id == 3: + self._penalty = arc_tan(a-d_rr) + else: + raise Exception("penalty関数の設定に誤りがあります") + + def get_rqi(self, rqi_id): + if rqi_id == 0: + return self.calculate_sqi() + elif rqi_id == 1: + return self.calculate_rqi1(ratio_bp=RATIO_BF) + elif rqi_id == 2: + return self.calculate_rqi2(ratio_bp=RATIO_BF) + elif rqi_id == 3: + return self.calculate_rqi3(ratio_bp=RATIO_BF, ration_hp=RATIO_HF) + elif rqi_id == 4: + return self.calculate_rqi4(ratio_bp=RATIO_BF, ratio_hp=RATIO_HF) + elif rqi_id == 5: + return self.calculate_rqi5() + elif rqi_id == 6: + return self.calculate_best_rqi(ratio_bp=RATIO_BF, ratio_hp=RATIO_HF) + else: + raise Exception("存在しないRQIが指定されました") + + def calculate_rqi1(self, ratio_bp=RATIO_BF): + # バンド帯に存在する大きいスペクトルに着目 + Fbp = np.sum(self._bp_spectrum > np.max(self._bp_spectrum) * ratio_bp) / len(self._bp_spectrum) + + # 高周波領域のノイズに着目 + Fhp = np.mean(self._hp_spectrum) / np.max(self._hp_spectrum) + + return np.round(1 - Fbp - Fhp, 2) + + def calculate_rqi2(self, ratio_bp=RATIO_BF): + # バンド帯に存在する大きいスペクトルに着目 + Fbp = np.sum(self._bp_spectrum > np.max(self._bp_spectrum) * ratio_bp) / len(self._bp_spectrum) + + # 高周波領域のノイズに着目 + Fhp = np.mean(self._hp_spectrum) / np.max(self._hp_spectrum) + + rqi = np.max(self._bp_spectrum) * np.round(1 - Fbp - Fhp, 2) + + return rqi + + def calculate_rqi3(self, ratio_bp=RATIO_BF, ration_hp=RATIO_HF): + + # バンド帯に存在する大きいスペクトルに着目 + Fbp = np.sum(self._bp_spectrum > np.max(self._bp_spectrum) * ratio_bp) / len(self._bp_spectrum) + + # 高周波領域のノイズに着目 + Fhp = np.sum(self._hp_spectrum > np.max(self._hp_spectrum) * ration_hp) / len(self._hp_spectrum) + + return np.round(1 - Fbp - Fhp, 2) + + # 専門ゼミのSQIを再現 + def calculate_sqi(self): + max_spectrum = np.max(self._spectrum) + + F1 = np.max(self._spectrum[self.hp_index:]) / max_spectrum + F2 = np.sum(self._spectrum[self.hp_index:] > max_spectrum * 0.1) / len(self._spectrum[index_2Hz:]) + F3 = (np.max(self._spectrum[:self.lp_index]) - + np.max(self._spectrum[self.lp_index:self.hp_index])) + F3 /= max_spectrum + F4 = (np.max(self._spectrum[:self.lp_index]) / np.max(self._spectrum[self.lp_index:self.hp_index])) + + if F4 > 2: + sqi = 1 - (F3 / 2 + (F1 + F2) / 4) + else: + sqi = 1 - (F1 + F2) / 2 + + return sqi + + def calculate_rqi4(self, ratio_bp=RATIO_BF, ratio_hp=RATIO_HF): + # バンド帯に存在する大きいスペクトルに着目 + Fbp = np.sum(self._bp_spectrum > np.max(self._bp_spectrum) * ratio_bp) / len(self._bp_spectrum) + + # 高周波領域のノイズに着目 + Fhp = np.sum(self._hp_spectrum > np.max(self._hp_spectrum) * ratio_hp) / len(self._hp_spectrum) + + return np.round(1 - (Fbp + Fhp)/2, 2) + + def calculate_best_rqi(self, ratio_bp=RATIO_BF, ratio_hp=RATIO_HF): + # バンド帯に存在する大きいスペクトルに着目 + Fbp = np.sum(self._bp_spectrum > np.max(self._bp_spectrum) * ratio_bp) / len(self._bp_spectrum) + # 高周波領域のノイズに着目 + Fhp = np.sum(self._hp_spectrum > np.max(self._hp_spectrum) * ratio_hp) / len(self._hp_spectrum) + + return self._penalty * np.round(1 - (Fbp + Fhp) / 2, 2) + + def calculate_rqi5(self, ration_bp=RATIO_BF, ration_hp=RATIO_HF): + return np.max(self._bp_spectrum) + + + + + + + diff --git a/rrEstimation.py b/rrEstimation.py new file mode 100644 index 0000000..085b0a4 --- /dev/null +++ b/rrEstimation.py @@ -0,0 +1,162 @@ +import cv2 +import numpy as np +from scipy.fftpack import fft +from scipy import signal +import matplotlib.pyplot as plt +from statistics import mean, median, variance, stdev +import matplotlib.pyplot as plt +import statistics +from SignalProcess import Roi2Signal +from GLOBAL import * +import seaborn as sns + +import os +from math import e, tanh, atan, pi +# sampling rate = 8.6, RR分解能 = 0.1にするためのpadding +FREQUENCY_RESOLUTION = 0.1 + + +class rrEstimation: + def __init__(self, __signal: Roi2Signal, dst_dir, sampling_rate, window=False): + self._signal = __signal.get_signal(window) + self._time = __signal.get_times(window) + self._dst_dir = dst_dir + self._sampling_rate = sampling_rate + + self._dt = self._time[-1] - self._time[-2] + self._N = len(self._time) + self._frequency = np.linspace(0, 1.0 / self._dt, self._N) + + self._spectrum = 0 + + # 分解能をあげるためにpaddingを設定 + self._n_padding = int(60*self._sampling_rate/FREQUENCY_RESOLUTION) + + def get_frequency(self): + return self._frequency + + def get_spectrum(self): + return self._spectrum + + def get_preprocessed_signal(self): + return signal.detrend(self._signal) + + # TODO 周波数解析の他クラス化 + def calculate_FFT(self): + preprocessed_signal = self.get_preprocessed_signal() + self._frequency = np.linspace(0, 1.0 / self._dt, self._N) + self._spectrum = np.abs(fft(preprocessed_signal) / (self._N / 2)) + + def calculate_PSD(self): + # TODO:PSDの前にナイキスト周波数以上の高周波を除去する必要あり? + + # preprocessed_signal = self.get_preprocessed_signal() + self. _frequency, self._spectrum = signal.welch(self._signal, 1 / self._dt, + detrend='linear', + nfft=self._n_padding) + + def calculate_PSD_Welch(self, window='boxcar', per_seg=2, det='linear'): + self._frequency, self._spectrum = signal.welch(self._signal, 1 / self._dt, + window=window, detrend=det, + nperseg=self._N // per_seg) + + def estimate_f_rr(self): + lp_index = np.abs(self._frequency - LF_THRESH).argmin() + hp_index = np.abs(self._frequency - HF_THRESH).argmin() + + # BP域最大スペクトルに対応する周波数を呼吸周波数とする + rr_index = np.argmax(self._spectrum[lp_index:hp_index]) + lp_index + + return round(self._frequency[rr_index] * 60, 2) + + def estimate_t_rr(self, size=5, iteration=3): + signalClass = Roi2Signal() + preprocessed_signal = self.get_preprocessed_signal() + signalClass.load_signals(self._time, preprocessed_signal) + + smoothed_signal = [] + # 平滑化 + for i in range(iteration): + smoothed_signal = signalClass.apply_moving_average_filter(num=size, window=False) + signalClass.load_signals(self._time, smoothed_signal) + + # zero-cross点の取得 + zero_cross_points = [] + for i in range(1, len(smoothed_signal)): + if smoothed_signal[i - 1] * smoothed_signal[i] < 0: + zero_cross_point = (self._time[i - 1] + self._time[i]) / 2 + zero_cross_points.append(zero_cross_point) + + # 呼吸周期候補の算出 + rr_candidate = [] + for i in range(1, len(zero_cross_points)): + periodic_time = (zero_cross_points[i] - zero_cross_points[i - 1]) * 2 + rr = 60 / periodic_time + # 算出されたrrがBP域の時のみ呼吸数の候補としてみなす + if (rr > LF_THRESH * 60) and (rr < HF_THRESH * 60): + rr_candidate.append(rr) + + # 呼吸数の算出 + # 0.1Hz~2Hzの場合のzero-cross点の数と合わない場合，エラー判定 + t_window = len(self._time)/self._sampling_rate + if (int(LF_THRESH * t_window * 2 + 0.5) < len(rr_candidate)) and (len(rr_candidate) < int(HF_THRESH * t_window * 2 + 0.5)): + rr = statistics.median(rr_candidate) + else: + rr = OUTLIER + + # print("t_rr :", rr) + return rr + + def write_signal2png(self, filename, dst_dir): + file_path = os.path.join(dst_dir, filename) + signalClass = Roi2Signal() + preprocessed_signal = self.get_preprocessed_signal() + signalClass.load_signals(self._time, preprocessed_signal) + smoothed_signal = signalClass.apply_moving_average_filter(num=5, window=False) + + plt.figure(figsize=(8, 6), dpi=200) + plt.xticks(fontsize=15) + plt.yticks(fontsize=15) + plt.plot(self._time[3:-3], smoothed_signal[3:-3]) + + plt.xlabel(r"$Time[s]$", style='italic', ha='center', fontsize=25, + fontname="Times New Roman") + plt.ylabel(r"$Signal$" + " " + r"$Intensity$", style='italic', ha='center', + fontsize=25, fontname="Times New Roman") + plt.tight_layout() + plt.savefig(os.path.join(file_path)) + + def write_psd2png(self, filename, dst_dir): + file_path = os.path.join(dst_dir, filename) + + plt.figure(figsize=(8, 6), dpi=200) + plt.xticks(fontsize=15) + plt.yticks(fontsize=15) + plt.plot(self._frequency, self._spectrum) + plt.xlabel(r"$Frequency[Hz]$", style='italic', ha='center', fontsize=25, + fontname="Times New Roman") + plt.ylabel(r"$Power$" + " " + r"$Spectrum$" + " " + r"$Density$", style='italic', + ha='center', fontsize=25, fontname="Times New Roman") + plt.tight_layout() + plt.savefig(os.path.join(file_path)) + + + + + + + + + + + + + + + + + + + + +