diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a047a94 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +.idea/ +__pycache__/ \ No newline at end of file diff --git a/DicomProcessor.py b/DicomProcessor.py new file mode 100644 index 0000000..44fc729 --- /dev/null +++ b/DicomProcessor.py @@ -0,0 +1,185 @@ +import pydicom +import glob +import sys, os +import numpy as np +import cv2 +import utils +import matplotlib.pyplot as plt +from tqdm import tqdm + + +class DicomProcessor: + + def __init__(self, dicom_dir, args): + self.args = args + self.dicom_dir = dicom_dir + self.dicom_file_list = glob.glob(os.path.join(dicom_dir, '*')) + self.dicom_file_list = sorted(self.dicom_file_list, key=lambda x: float(pydicom.read_file(x).SliceLocation), reverse=True) + self.dicom_size = len(self.dicom_file_list) + + if self.dicom_size == 0: + print("dicom file has no files") + raise + + sample_dfile = pydicom.read_file(self.dicom_file_list[0]) + sample_dfile2 = pydicom.read_file(self.dicom_file_list[1]) + self.dcm_wc = sample_dfile.WindowCenter + self.dcm_ww = sample_dfile.WindowWidth + print("original wc: {}".format(self.dcm_wc)) + print("original ww: {}".format(self.dcm_ww)) + self.dicom_row_num = sample_dfile.Rows + self.dicom_col_num = sample_dfile.Columns + + self.pixel_spacing_dx, self.pixel_spacing_dy = sample_dfile.PixelSpacing + self.slice_spacing = abs(sample_dfile.ImagePositionPatient[2] - sample_dfile2.ImagePositionPatient[2]) + self.image_position_patient = sample_dfile.ImagePositionPatient + #なぜか3DSlicer上ではSが-0.25されている + # self.image_position_patient[2] -= 2.5 + self.CT_for_imshow = self.load_initial_CT() + self.original_eso_radius = args.original_eso_radius + self.target_eso_radius = int(args.eso_radius / self.pixel_spacing_dx) + + def calc_ijk2LPS_mat(self, dst=None): + ijk_samples = self.calc_ijk_samples() + LPS_samples = self.calc_LPS_samples() + + Ls = LPS_samples[:, 0] + Ps = LPS_samples[:, 1] + Ss = LPS_samples[:, 2] + + row1 = np.linalg.solve(ijk_samples, Ls) + row2 = np.linalg.solve(ijk_samples, Ps) + row3 = np.linalg.solve(ijk_samples, Ss) + row4 = np.array([0, 0, 0, 1]) + + IJKtoLPS_mat = np.array([row1, row2, row3, row4]) + + if dst != None: + dst = IJKtoLPS_mat + + return IJKtoLPS_mat + + def calc_ijk_samples(self): + origin_k = len(self.dicom_file_list) - 1 + ijk = np.array([[0, 0, origin_k, 1], + [1, 0, origin_k, 1], + [0, 1, origin_k, 1], + [0, 0, origin_k - 1, 1]]) + + return ijk + + def calc_LPS_samples(self): + origin_L, origin_P, origin_S = self.image_position_patient + + LPS = np.array([[origin_L, origin_P, origin_S], + [origin_L + self.pixel_spacing_dx, origin_P, origin_S], + [origin_L, origin_P + self.pixel_spacing_dy, origin_S], + [origin_L, origin_P, origin_S - self.slice_spacing]]) + + return LPS + + def load_initial_CT(self): + CT_img = pydicom.read_file(self.dicom_file_list[0]) + CT_img = CT_img.pixel_array + + return self.row2uint8(CT_img) + + def row2uint8(self, CT_row, delete_0s=False): + if delete_0s: + CT_row = np.where(CT_row == 0, 0, CT_row - np.min(CT_row[CT_row != 0])) + + window_max = self.dcm_wc + self.dcm_ww / 2 + window_min = self.dcm_wc - self.dcm_ww / 2 + + CT_img = CT_row.astype(np.float) + CT_img[CT_img < window_min] = window_min + CT_img[window_max < CT_img] = window_max + + CT_img -= np.mean(CT_img) + CT_img = CT_img / (np.max(np.abs(CT_img)) + 1e-5) * 128.0 + CT_img -= np.mean(CT_img) + CT_img = np.clip(CT_img, 0, 255).astype(np.uint8) + CT_img = cv2.cvtColor(CT_img, cv2.COLOR_GRAY2BGR) + + return CT_img + + def adapt_window(self, array, window_low, window_high): + return np.clip(array, window_low, window_high) + + def get_CT_by_index(self, dicom_index): + # print(self.dicom_file_list[dicom_index]) + CT_row = pydicom.read_file(self.dicom_file_list[dicom_index]).pixel_array + + return self.row2uint8(CT_row) + + def calc_k_on_ijk_coordinates(self, dicom_index): + return self.dicom_size - dicom_index - 1 + + def make_distorted_dicom(self, ijk_eso_centers, img_eso_radius): + pad_eso_centers = self.liner_pad_list(ijk_eso_centers) + os.makedirs(self.args.out_dicom_dir, exist_ok=True) + if (self.dicom_col_num % 2 != 0) and (self.dicom_row_num % 2 != 0): + print('It may not be able to get correct dicom image. In detail, please ask to Yukiya') + raise + half_col_num = self.dicom_col_num / 2 + half_row_num = self.dicom_row_num / 2 + fixed_points = np.array([[half_col_num, self.dicom_row_num], [self.dicom_col_num, self.dicom_row_num], + [self.dicom_col_num, half_row_num], [self.dicom_col_num, 0], + [half_col_num, 0], [0, 0], + [0, half_row_num], [0, self.dicom_row_num]], dtype=np.float32) + for dicom_index, eso_center in enumerate(pad_eso_centers): + print("a") + file_name = os.path.basename(self.dicom_file_list[dicom_index]) + new_dicom_path = os.path.join(self.args.out_dicom_dir, file_name) + if eso_center is None: + ds = pydicom.dcmread(self.dicom_file_list[dicom_index]) + # d.save_as(new_dicom_path) + # d.save_as(self.dicom_file_list[dicom_index]) + ds.save_as(new_dicom_path) + else: + show_CT_img = self.get_CT_by_index(dicom_index) + cv2.imshow("src", show_CT_img) + original_radius_by_pixel = np.int(self.original_eso_radius / self.pixel_spacing_dx) + src_8neighbors = utils.get_8neighbors_with_radius(original_radius_by_pixel, eso_center) + target_8neighbors = utils.get_8neighbors_with_radius(self.target_eso_radius, eso_center) + + src_points = np.vstack((src_8neighbors, fixed_points)).reshape((1, -1, 2)) + target_points = np.vstack((target_8neighbors, fixed_points)).reshape((1, -1, 2)) + + matches = list() + for i in range(len(src_points[0, :, :])): + matches.append(cv2.DMatch(i, i, 0)) + + tps = cv2.createThinPlateSplineShapeTransformer() + tps.estimateTransformation(target_points, src_points, matches) + ds = pydicom.read_file(self.dicom_file_list[dicom_index]) + ds.decompress() + out = tps.warpImage(ds.pixel_array) + show_out = tps.warpImage(show_CT_img[:, :, 0]) + cv2.imshow("distorted", show_out) + + cv2.waitKey(1) + + ds.PixelData = out.tobytes() + ds.save_as(new_dicom_path) + + + def liner_pad_list(self, ijk_eso_centers): + pad_center_list = [] + Is_first = True + prev_index = 0 + for index, eso_center in enumerate(ijk_eso_centers): + if (eso_center is not None) and Is_first: + prev_index = index + Is_first = False + elif (eso_center is not None) and (Is_first is not True): + item_distance = index - prev_index + if item_distance != 1: + average_weight = 1 / item_distance + for new_index in range((prev_index + 1), index, 1): + pad_center = (ijk_eso_centers[prev_index] * (average_weight * (index - new_index)) + + ijk_eso_centers[index] * (average_weight * (new_index - prev_index))) + pad_center_list[new_index] = pad_center + prev_index = index + pad_center_list.append(eso_center) + return pad_center_list diff --git a/EsoOBJMaker.py b/EsoOBJMaker.py new file mode 100644 index 0000000..7b37635 --- /dev/null +++ b/EsoOBJMaker.py @@ -0,0 +1,212 @@ +import numpy as np +import sys, os +import math +from scipy import interpolate + + +class EsoOBJMaker: + + def __init__(self, ijk_eso_centers, img_eso_radius, ijk2LPS_mat, args): + self.args = args + self.ijk_eso_centers = self.preprocess_eso_centers(ijk_eso_centers) + self.ijk_eso_centers = sorted([x for x in ijk_eso_centers if x is not None], key=lambda item: item[2]) + self.img_eso_radius = img_eso_radius + self.ijk2LPS_mat = ijk2LPS_mat + + def preprocess_eso_centers(self, ijk_eso_centers): + ijk_eso_centers = self.pad_centers(ijk_eso_centers) + row_centers = sorted([x for x in ijk_eso_centers if x is not None], key=lambda item: item[2]) + spline_centers = self.calc_spline_centers(row_centers) + return spline_centers + + def pad_centers(self, ijk_eso_centers): + x_list = [] + i_list = [] + j_list = [] + k_list = [] + for cur_index, cur_center in enumerate(ijk_eso_centers): + if cur_center is not None: + x_list.append(cur_index) + i_list.append(cur_center[0]) + j_list.append(cur_center[1]) + k_list.append(cur_center[2]) + f_i = interpolate.interp1d(x_list, i_list, kind='cubic') + f_j = interpolate.interp1d(x_list, j_list, kind='cubic') + f_k = interpolate.interp1d(x_list, k_list, kind='linear') + # 中心点をスプライン補完 + prev_center = None + prev_index = 0 + for cur_index, cur_center in enumerate(ijk_eso_centers): + if cur_center is not None: + if prev_center is None: + prev_center = cur_center + prev_index = cur_index + else: + for i in range(prev_index + 1, cur_index): + ijk_eso_centers[i] = np.array([f_i(i), f_j(i), f_k(i)]) + prev_center = cur_center + return ijk_eso_centers + + + def calc_spline_centers(self, row_centers): + i_list = [x[0] for x in row_centers] + j_list = [x[1] for x in row_centers] + k_list = [x[2] for x in row_centers] + i_CS = interpolate.interp1d(k_list, i_list, kind='cubic') + j_CS = interpolate.interp1d(k_list, j_list, kind='cubic') + k_for_CS = np.linspace(k_list[0], k_list[len(k_list) - 1], self.args.spline_num * (len(k_list) - 1) + 1) + spline_ijk = [np.array([i_CS(k), j_CS(k), k]) for k in k_for_CS] + return spline_ijk + + def create_obj(self): + os.makedirs(self.args.output_dir, exist_ok=True) + eso_center_LPS = self.eso_centers_ijk2LPS() + ijk_eso_vertices = self.calc_ijk_eso_vertices() + LPS_eso_vertices = self.convert_ijk2LPS(ijk_eso_vertices) + eso_vn = self.calc_vn(LPS_eso_vertices) + with open(os.path.join(self.args.output_dir, 'Eso_center_only.obj'), 'w') as f: + self.write_head_obj_file(f) + self.write_center_obj_file(f, eso_center_LPS) + + # with open(os.path.join(self.args.output_dir, 'Eso_made_by_python.obj'), 'w') as f: + # self.write_head_obj_file(f) + # self.write_vertices_obj_file(f, LPS_eso_vertices) + # self.write_vertex_normals_obj(f, eso_vn) + # self.write_faces_obj_file(f, LPS_eso_vertices, eso_vn) + + def write_head_obj_file(self, f): + f.write('# n-lab program made by Yukiya OBJ file\n') + f.write('# chiba univ\n') + f.write('mtlib Eso_made_by_oython.mtl\n') + f.write('o Eso_made_by_python.obj\n') + + def write_vertices_obj_file(self, f, LPS_eso_vertices): + print('start writing vertices') + for circle_vertices in LPS_eso_vertices: + for model_vertx in circle_vertices: + f.write('v {:.6f} {:.6f} {:.6f}\n'.format(model_vertx[0], model_vertx[1], model_vertx[2])) + print('end writing vertices') + + def write_vertex_normals_obj(self, f, vn): + print('start writing vertex normals') + for circle_vn in vn: + for n in circle_vn: + f.write('vn {:.6f} {:.6f} {:.6f}\n'.format(n[0], n[1], n[2])) + print('end writing vertex normals') + + def write_faces_obj_file(self, f, LPS_eso_vertices, eso_vn): + print('start writing face') + circle_vertices_size = len(LPS_eso_vertices[0]) + vn_num = 1 + for level_index in range(len(LPS_eso_vertices) - 1): + level_init_num = level_index * circle_vertices_size + 1 + for vertex_index in range(circle_vertices_size): + if vertex_index == (circle_vertices_size - 1): + vertex_set = [level_init_num + vertex_index, level_init_num, + level_init_num + circle_vertices_size + vertex_index, + level_init_num + circle_vertices_size] + + else: + vertex_set = [level_init_num + vertex_index, level_init_num + vertex_index + 1, + level_init_num + circle_vertices_size + vertex_index, + level_init_num + circle_vertices_size + vertex_index + 1] + f.write('f {0}//{4} {1}//{4} {3}//{4} {2}//{4}\n'.format(vertex_set[0], vertex_set[1], + vertex_set[2], vertex_set[3], vn_num)) + vn_num += 1 + #f.write('f {0}//{3} {1}//{3} {2}//{3}\n'.format(vertex_set[0], vertex_set[2], + # vertex_set[3], vn_num)) + #vn_num += 1 + #f.write('f {0}//{3} {1}//{3} {2}//{3}\n'.format(vertex_set[0], vertex_set[1], + # vertex_set[3], vn_num)) + #vn_num += 1 + print('end writing face') + + + def calc_vn(self, LPS_eso_vertices): + circle_vertices_size = len(LPS_eso_vertices[0]) + eso_vn = [] + for level_index in range(len(LPS_eso_vertices) - 1): + vn_set = [] + for vertex_index in range(circle_vertices_size): + if vertex_index == (circle_vertices_size - 1): + vector1 = LPS_eso_vertices[level_index][vertex_index] - LPS_eso_vertices[level_index + 1][ + vertex_index] + vector2 = LPS_eso_vertices[level_index + 1][0] - LPS_eso_vertices[level_index + 1][ + vertex_index] + vn = np.cross(vector1, vector2) + vn /= np.sum(np.abs(vn)) + vn_set.append(vn) + #vector1 = LPS_eso_vertices[level_index][vertex_index] - LPS_eso_vertices[level_index + 1][ + # 0] + #vector2 = LPS_eso_vertices[level_index][0] - LPS_eso_vertices[level_index + 1][ + # 0] + #vn = np.cross(vector1, vector2) + #vn /= np.sum(np.abs(vn)) + #vn_set.append(vn) + else: + vector1 = LPS_eso_vertices[level_index][vertex_index] - LPS_eso_vertices[level_index + 1][ + vertex_index] + vector2 = LPS_eso_vertices[level_index + 1][vertex_index + 1] - LPS_eso_vertices[level_index + 1][ + vertex_index] + vn = np.cross(vector1, vector2) + vn /= np.sum(np.abs(vn)) + vn_set.append(vn) + #vector1 = LPS_eso_vertices[level_index][vertex_index] - LPS_eso_vertices[level_index + 1][ + # vertex_index + 1] + #vector2 = LPS_eso_vertices[level_index][vertex_index + 1] - LPS_eso_vertices[level_index + 1][ + # vertex_index + 1] + #vn = np.cross(vector1, vector2) + #vn /= np.sum(np.abs(vn)) + #vn_set.append(vn) + eso_vn.append(vn_set) + return eso_vn + + def calc_ijk_eso_vertices(self): + ijk_eso_vertices = [] + init_relative_radius = np.array([0.0, float(self.img_eso_radius)]) + step_theta = 360.0 / self.args.circle_divisions + for ijk_center in self.ijk_eso_centers: + tmp_ij = np.array([ijk_center[0], ijk_center[1]]) + theta = 0 + eso_circle_vertices = [] + for i in range(self.args.circle_divisions): + rotae_mat = np.array([[math.cos(math.radians(theta)), -math.sin(math.radians(theta))], + [math.sin(math.radians(theta)), math.cos(math.radians(theta))]]) + ij = tmp_ij + np.dot(rotae_mat, np.transpose(init_relative_radius)) + eso_circle_vertices.append(np.array([ij[0], ij[1], ijk_center[2]])) + theta += step_theta + ijk_eso_vertices.append(eso_circle_vertices) + return ijk_eso_vertices + + def convert_ijk2LPS(self, ijk_vertices): + LPS_eso_verices = [] + for ijk_circle_vertices in ijk_vertices: + LPS_circle_vertice = [] + for ijk_vertice in ijk_circle_vertices: + tmp = np.ones((4, 1), dtype=np.float) + tmp[:3, :] = np.transpose(ijk_vertice[np.newaxis, :]) + tmp = np.dot(self.ijk2LPS_mat, tmp) + LPS_circle_vertice.append(np.transpose(tmp)[0, :3]) + LPS_eso_verices.append(LPS_circle_vertice) + return LPS_eso_verices + + def eso_centers_ijk2LPS(self): + eso_center_LPS = list() + for center_ijk in self.ijk_eso_centers: + tmp = np.ones((4, 1), dtype=np.float) + tmp[:3, :] = np.transpose(center_ijk[np.newaxis, :]) + tmp = np.dot(self.ijk2LPS_mat, tmp) + eso_center_LPS.append(np.transpose(tmp)[0, :3]) + return eso_center_LPS + + def write_center_obj_file(self, f, LPS_eso_vertices): + print('start writing vertices') + num_vertices = 0 + for model_vertx in LPS_eso_vertices: + f.write('v {:.6f} {:.6f} {:.6f}\n'.format(model_vertx[0], model_vertx[1], model_vertx[2])) + num_vertices += 1 + print('end writing vertices') + for i in range(3, num_vertices): + f.write('f {0}//{1} {1}//{2} {2}//{1}\n'.format(i - 2, i - 1, i)) + # def write_eso_centers(self, f): + diff --git a/GUI.py b/GUI.py new file mode 100644 index 0000000..ef7e070 --- /dev/null +++ b/GUI.py @@ -0,0 +1,251 @@ +import numpy as np +import cv2 +from DicomProcessor import DicomProcessor +from EsoOBJMaker import EsoOBJMaker + +GUI_LOAD_NEXT = 1 +GUI_LOAD_PREV = 0 + + +class GUIController: + + def __init__(self, args): + self.args = args + self.Is_continue = True + self.Is_zoomed = False + self.dicom_processor = DicomProcessor(args.dicom_dir, args) + self.dicom_size = len(self.dicom_processor.dicom_file_list) + self.CT_img_for_show = self.dicom_processor.CT_for_imshow + self.cur_src_CT_img = self.CT_img_for_show + self.src_CT_size = self.cur_src_CT_img.shape + self.zoomed_CT_size = self.preprocess_zoomed_CT_size() + self.zoom_point = [0, 0] + + # menuボタン系 + self.target_color = (24, 185, 237) + self.non_target_color = (220, 220, 220) + + self.window_button_state = { + "pushed_wc": True, + "pushed_ww": False + } + self.wc_width = None + self.ww_width = None + self.state_board = None + self.button_height, self.button_width = 0, 0 + self.show_wc_str = "WC: {}".format(round(self.dicom_processor.dcm_wc)) + self.show_ww_str = "WW: {}".format(round(self.dicom_processor.dcm_ww)) + self.wc_str_color = self.target_color + self.ww_str_color = self.non_target_color + self.button_img = self.create_button_img() + + # まだ画像のピクセル幅がx, yで同じときのみに対応している． + self.img_eso_radius = args.eso_radius / self.dicom_processor.pixel_spacing_dx + self.set_init_window() + + #出力用変数 + self.dicom_index = 0 + self.ijk_eso_centers = [None] * self.dicom_size + + self.GUI_imshow() + + def set_init_window(self): + cv2.namedWindow('menu') + cv2.namedWindow('CT image') + + def preprocess_zoomed_CT_size(self): + zoomed_CT_size = [int(x / self.args.magnification_ratio) for x in self.CT_img_for_show.shape][:2] + if zoomed_CT_size[0] % 2 == 0: + zoomed_CT_size[0] -= 1 + if zoomed_CT_size[1] % 2 == 0: + zoomed_CT_size[1] -= 1 + return zoomed_CT_size + + def create_button_img(self): + buttons = ['back', 'next', 'export'] + button_size = (self.args.resized_button_width, self.args.resized_button_height) + for i, button_name in enumerate(buttons): + if i == 0: + button_img = cv2.imread('./button_img/' + button_name + '.png') + button_img = cv2.resize(button_img, button_size) + self.button_height, self.button_width, _ = button_img.shape + else: + tmp = cv2.imread('./button_img/' + button_name + '.png') + tmp = cv2.resize(tmp, button_size) + button_img = cv2.hconcat([button_img, tmp]) + + # 2021/1/12追加のウィンドウ処理処理用のボタン + wc_pushed = cv2.imread("./button_img/wc_pushed.png") + aspect_rate = self.args.resized_button_height / wc_pushed.shape[0] + wc_pushed = cv2.resize(wc_pushed, (round(aspect_rate * wc_pushed.shape[1]), round(aspect_rate * wc_pushed.shape[0]))) + + ww = cv2.imread("./button_img/ww.png") + aspect_rate = self.args.resized_button_height / ww.shape[0] + ww = cv2.resize(ww, (round(aspect_rate * ww.shape[1]), round(aspect_rate * ww.shape[0]))) + + width = button_img.shape[1] - wc_pushed.shape[1] - ww.shape[1] + self.state_board = np.zeros((ww.shape[0], width, 3)).astype(np.uint8) + cv2.putText(self.state_board, self.show_wc_str, (20, 30), cv2.FONT_HERSHEY_SIMPLEX, 1, self.wc_str_color, + 1, cv2.LINE_AA) + cv2.putText(self.state_board, self.show_ww_str, (20, 60), cv2.FONT_HERSHEY_SIMPLEX, 1, self.ww_str_color, + 1, cv2.LINE_AA) + + window_buttons = np.hstack([wc_pushed, ww, self.state_board]) + + button_img = np.vstack([button_img, window_buttons]) + self.wc_width = wc_pushed.shape[1] + self.ww_width = ww.shape[1] + + return button_img + + def change_cur_CT(self, mode): + if mode is GUI_LOAD_NEXT: + self.dicom_index += 1 + if mode is GUI_LOAD_PREV: + self.dicom_index -= 1 + print('index : {}'.format(self.dicom_index)) + self.CT_img_for_show = self.dicom_processor.get_CT_by_index(self.dicom_index) + self.cur_src_CT_img = self.CT_img_for_show + self.Is_zoomed = False + self.GUI_imshow() + + def GUI_imshow(self): + cv2.imshow("menu", self.button_img) + cv2.imshow("CT image", self.CT_img_for_show) + cv2.waitKey(1) + + def menu_callbacks(self, event, x, y, flags, param): + # click back + if event == cv2.EVENT_LBUTTONUP and (0 <= x and x <= self.button_width) and y <= self.button_height and self.dicom_index != 0: + self.change_cur_CT(GUI_LOAD_PREV) + + # click next + elif event == cv2.EVENT_LBUTTONUP and (self.button_width < x and x < (2 * self.button_width)) and y <= self.button_height and self.dicom_index != (self.dicom_size - 1): + self.change_cur_CT(GUI_LOAD_NEXT) + + # click export + elif event == cv2.EVENT_LBUTTONUP and ((2 * self.button_width) < x and x < (3 * self.button_width)) and y <= self.button_height: + print("click export") + # self.dicom_processor.make_distorted_dicom(self.ijk_eso_centers, self.img_eso_radius) + eso_obj_maker = EsoOBJMaker(self.ijk_eso_centers, self.img_eso_radius, + self.dicom_processor.calc_ijk2LPS_mat(), self.args) + eso_obj_maker.create_obj() + self.Is_continue = False + print('exit') + + # click wc + elif event == cv2.EVENT_LBUTTONUP and x < self.wc_width and self.button_height < y and not self.window_button_state["pushed_wc"]: + # この辺は関数か出来る． + replaced_button = cv2.resize(cv2.imread("./button_img/wc_pushed.png"), (self.wc_width, self.button_height)) + self.button_img[self.button_height:, :self.wc_width, :] = replaced_button + replaced_button = cv2.resize(cv2.imread("./button_img/ww.png"), (self.ww_width, self.button_height)) + self.button_img[self.button_height:, self.wc_width:(self.wc_width + self.ww_width), :] = replaced_button + self.window_button_state["pushed_wc"], self.wc_str_color = True, self.target_color + self.window_button_state["pushed_ww"], self.ww_str_color = False, self.non_target_color + self.update_state_board() + + # click ww + elif event == cv2.EVENT_LBUTTONUP and self.wc_width < self.wc_width + self.ww_width and self.button_height < y and not self.window_button_state["pushed_ww"]: + # この辺は関数か出来る． + replaced_button = cv2.resize(cv2.imread("./button_img/wc.png"), (self.wc_width, self.button_height)) + self.button_img[self.button_height:, :self.wc_width, :] = replaced_button + replaced_button = cv2.resize(cv2.imread("./button_img/ww_pushed.png"), (self.ww_width, self.button_height)) + self.button_img[self.button_height:, self.wc_width:(self.wc_width + self.ww_width), :] = replaced_button + self.window_button_state["pushed_wc"], self.wc_str_color = False, self.non_target_color + self.window_button_state["pushed_ww"], self.ww_str_color = True, self.target_color + self.update_state_board() + + def CT_image_callbacks(self, event, x, y, flags, param): + if event is cv2.EVENT_LBUTTONUP and self.Is_zoomed is True: + zoom_i_ratio, zoom_j_ratio = ((self.src_CT_size[0] / self.zoomed_CT_size[0]), + (self.src_CT_size[1] / self.zoomed_CT_size[1])) + cv2.circle(self.CT_img_for_show, (x, y), 5, (0, 0, 255), -1, cv2.LINE_AA) + cv2.circle(self.CT_img_for_show, (x, y), int(self.img_eso_radius * zoom_i_ratio), + (200, 50, 50),lineType=cv2.LINE_AA) + non_zoomed_i, non_zoomed_j = self.calc_non_zoomed_ij(x, y, zoom_i_ratio, zoom_j_ratio) + print('clicked i:{} j:{} k:{}'.format(non_zoomed_i, non_zoomed_j, + self.dicom_processor.calc_k_on_ijk_coordinates(self.dicom_index))) + self.ijk_eso_centers[self.dicom_index] = np.array([non_zoomed_i, non_zoomed_j, + self.dicom_processor.calc_k_on_ijk_coordinates(self.dicom_index)]) + elif event is cv2.EVENT_LBUTTONUP and self.Is_zoomed is False: + cv2.circle(self.CT_img_for_show, (x, y), 5, (0, 0, 255), -1, lineType=cv2.LINE_AA) + cv2.circle(self.CT_img_for_show, (x, y), int(self.img_eso_radius), (200, 50, 50), lineType=cv2.LINE_AA) + print('clicked i:{} j:{} k:{}'.format(x, y, self.dicom_processor.calc_k_on_ijk_coordinates(self.dicom_index))) + self.ijk_eso_centers[self.dicom_index] = np.array([x, y, self.dicom_processor.calc_k_on_ijk_coordinates(self.dicom_index)]) + elif event is cv2.EVENT_RBUTTONUP and self.Is_zoomed is False: + self.zoom_function(x, y) + elif event is cv2.EVENT_RBUTTONUP and self.Is_zoomed is True: + self.CT_img_for_show = self.cur_src_CT_img + self.Is_zoomed = False + + def calc_non_zoomed_ij(self, i, j, zoom_i_ratio, zoom_j_ratio): + non_zoomed_i, non_zoomed_j = (self.zoom_point[0] + (i / zoom_i_ratio), + self.zoom_point[1] + (j / zoom_j_ratio)) + return non_zoomed_i, non_zoomed_j + + def zoom_function(self, x, y): + self.zoom_point = [int(x - ((self.zoomed_CT_size[0] - 1) / 2)), int(y - ((self.zoomed_CT_size[1] - 1) / 2))] + if self.zoom_point[0] < 0: + self.zoom_point[0] = 0 + if self.src_CT_size[1] <= 0: + self.zoom_point[1] = 0 + tmp = self.cur_src_CT_img[self.zoom_point[1]:(self.zoom_point[1] + self.zoomed_CT_size[1]), + self.zoom_point[0]:(self.zoom_point[0] + self.zoomed_CT_size[0]), :] + self.CT_img_for_show = cv2.resize(tmp, self.src_CT_size[:2], interpolation=cv2.INTER_LANCZOS4) + self.Is_zoomed = True + + def key_function(self): + key = cv2.waitKey(50) + if key is ord('d') and self.dicom_index != (self.dicom_size - 1): + self.change_cur_CT(GUI_LOAD_NEXT) + elif key is ord('a') and self.dicom_index != 0: + self.change_cur_CT(GUI_LOAD_PREV) + + # 2021/1/12追加の + elif key is ord('w'): + if self.window_button_state["pushed_wc"]: + self.dicom_processor.dcm_wc += 1 + print("changed => wc: ", self.dicom_processor.dcm_wc) + self.CT_img_for_show = self.dicom_processor.get_CT_by_index(self.dicom_index) + self.update_state_board() + + elif self.window_button_state["pushed_ww"]: + self.dicom_processor.dcm_ww += 1 + print("changed => ww: ", self.dicom_processor.dcm_ww) + self.CT_img_for_show = self.dicom_processor.get_CT_by_index(self.dicom_index) + self.update_state_board() + + elif key is ord('x'): + if self.window_button_state["pushed_wc"]: + self.dicom_processor.dcm_wc -= 1 + print("changed => wc: ", self.dicom_processor.dcm_wc) + self.CT_img_for_show = self.dicom_processor.get_CT_by_index(self.dicom_index) + self.update_state_board() + + elif self.window_button_state["pushed_ww"]: + self.dicom_processor.dcm_ww -= 1 + print("changed => ww: ", self.dicom_processor.dcm_ww) + self.CT_img_for_show = self.dicom_processor.get_CT_by_index(self.dicom_index) + self.update_state_board() + + def update_state_board(self): + self.state_board = np.zeros_like(self.state_board).astype(np.uint8) + self.show_wc_str = "WC: {}".format(round(self.dicom_processor.dcm_wc)) + self.show_ww_str = "WW: {}".format(round(self.dicom_processor.dcm_ww)) + cv2.putText(self.state_board, self.show_wc_str, (20, 30), cv2.FONT_HERSHEY_SIMPLEX, 1, + self.wc_str_color, + 1, cv2.LINE_AA) + cv2.putText(self.state_board, self.show_ww_str, (20, 60), cv2.FONT_HERSHEY_SIMPLEX, 1, + self.ww_str_color, + 1, cv2.LINE_AA) + self.button_img[self.button_height:, (self.wc_width + self.ww_width):, :] = self.state_board + + def run(self): + cv2.setMouseCallback('menu', self.menu_callbacks) + cv2.setMouseCallback('CT image', self.CT_image_callbacks) + + while True: + self.key_function() + self.GUI_imshow() + if not self.Is_continue: + break diff --git a/button_img/back.png b/button_img/back.png new file mode 100644 index 0000000..432cf61 --- /dev/null +++ b/button_img/back.png Binary files differ diff --git a/button_img/export.png b/button_img/export.png new file mode 100644 index 0000000..75716e7 --- /dev/null +++ b/button_img/export.png Binary files differ diff --git a/button_img/next.png b/button_img/next.png new file mode 100644 index 0000000..604e426 --- /dev/null +++ b/button_img/next.png Binary files differ diff --git a/button_img/wc.png b/button_img/wc.png new file mode 100644 index 0000000..708a4a0 --- /dev/null +++ b/button_img/wc.png Binary files differ diff --git a/button_img/wc_pushed.png b/button_img/wc_pushed.png new file mode 100644 index 0000000..bc0f720 --- /dev/null +++ b/button_img/wc_pushed.png Binary files differ diff --git a/button_img/ww.png b/button_img/ww.png new file mode 100644 index 0000000..29987a5 --- /dev/null +++ b/button_img/ww.png Binary files differ diff --git a/button_img/ww_pushed.png b/button_img/ww_pushed.png new file mode 100644 index 0000000..cede782 --- /dev/null +++ b/button_img/ww_pushed.png Binary files differ diff --git a/main.py b/main.py new file mode 100644 index 0000000..6056cce --- /dev/null +++ b/main.py @@ -0,0 +1,22 @@ +import argparse +from GUI import GUIController +from DicomProcessor import DicomProcessor + +parser = argparse.ArgumentParser() +parser.add_argument("--dicom_dir", type=str, required=True, help="DICOMファイルへのパス") +parser.add_argument("--output_dir", type=str, default="./obj", help='OBJファイルの出力先を指定') +parser.add_argument("--out_dicom_dir", type=str, default="./extend_dicom_tmp", help='食道拡張を行ったdicomの出力先') +parser.add_argument("--resized_button_height", type=int, default=154, help="メニューボタンの高さ") +parser.add_argument("--resized_button_width", type=int, default=314, help="メニューボタンの幅") +parser.add_argument("--eso_radius", type=float, default=15.0, help='mm単位で食道モデルの半径を指定') +parser.add_argument("--circle_divisions", type=int, default=50, help='食道モデルの半径方向への分解能') +parser.add_argument("--spline_num", type=int, default=10, help='食道の中心点間を補完する点の数') +parser.add_argument("--magnification_ratio", type=int, default=4, help='ダブルクリック時の拡大倍率（偶数のほうがいいかも）') +parser.add_argument("--original_eso_radius", type=float, default=5, help="TPS用の画像中の食道半径(mm単位)") +parser.add_argument("--window_low", type=float, default=-500.0, help="ウィンドウサイズの下限") +parser.add_argument("--window_high", type=float, default=500.0, help="ウィンドウサイズの上限") +args = parser.parse_args() + +if __name__ == '__main__': + GUI = GUIController(args) + GUI.run() diff --git a/moduleTest2.py b/moduleTest2.py new file mode 100644 index 0000000..443859e --- /dev/null +++ b/moduleTest2.py @@ -0,0 +1,5 @@ +import cv2 + +a = cv2.imread("./button_img/ww_pushed.png") +cv2.imshow("a", a) +cv2.waitKey() diff --git a/module_test.py b/module_test.py new file mode 100644 index 0000000..4152b6d --- /dev/null +++ b/module_test.py @@ -0,0 +1,109 @@ +import numpy as np +import cv2 +import pydicom + +xy_s = [] +center_xy = [] +dicom_num = 16 + +def mouse_callbacks(event, x, y, flags, param): + if event == cv2.EVENT_LBUTTONUP: + xy = (x, y) + print("x: {} y: {}".format(x, y)) + xy_s.append(xy) + + if event == cv2.EVENT_RBUTTONUP: + print("RBUTTONUP") + center_xy.append((x, y)) + + +def write_circle(xy_s, mat): + for xy in xy_s: + cv2.circle(mat, xy, 4, (0, 0, 200), -1) + return mat + +a = "{:03}".format(dicom_num) +ds = pydicom.dcmread(r'C:\Users\Planck\Desktop\3Dsyokudo\Image{:03}'.format(dicom_num)) +CT_row = ds.pixel_array +CT_row = np.where(CT_row == 0, 0, CT_row - np.min(CT_row[CT_row != 0])) + + +tmp = np.array(255 * (CT_row / np.max(CT_row)), dtype=np.uint8) +CT_img = cv2.cvtColor(tmp, cv2.COLOR_GRAY2BGR) +show_CT = CT_img.copy() +cv2.imshow("a", show_CT) +cv2.setMouseCallback('a', mouse_callbacks) + +while True: + cv2.imshow("a", show_CT) + key = cv2.waitKey(1) + if key == ord('a'): + write_circle(xy_s, show_CT) + spacing_dx, spacing_dy = ds.PixelSpacing + radius15_int = np.int(15.0 / spacing_dx) + dicom_row_num = 512 + dicom_col_num = 512 + half_col_num = np.int(dicom_col_num / 2) + half_row_num = np.int(dicom_row_num / 2) + fixed_points = np.array([[half_col_num, dicom_row_num], [dicom_col_num, dicom_row_num], + [dicom_col_num, half_row_num], [dicom_col_num, 0], + [half_col_num, 0], [0, 0], + [0, half_row_num], [0, dicom_row_num]], dtype=np.float32) + center = center_xy[0] + after_points = [(center[0], center[1] + radius15_int), (center[0] + radius15_int, center[1]), + (center[0], center[1] - radius15_int), (center[0] - radius15_int, center[1])] + for p in after_points: + cv2.circle(show_CT, p, 4, (0, 200, 0), thickness=-1) + cv2.imwrite("before{}.png".format(dicom_num), show_CT) + xy_s = np.array(xy_s, dtype=np.float32) + after_points = np.array(after_points, dtype=np.float32) + src = np.vstack((xy_s, fixed_points)).reshape((1, -1, 2)) + target = np.vstack((after_points, fixed_points)).reshape((1, -1, 2)) + + matches = list() + for i in range(len(src[0, :, :])): + matches.append(cv2.DMatch(i, i, 0)) + + tps = cv2.createThinPlateSplineShapeTransformer() + tps.estimateTransformation(target, src, matches) + out = tps.warpImage(CT_img) + cv2.imshow("dist", out) + cv2.imwrite("after{}.png".format(dicom_num), out) + + + if key == ord('b'): + center = center_xy[0] + spacing_dx, spacing_dy = ds.PixelSpacing + src_radius = np.int(5 / spacing_dx) + xy_s = [(center[0], center[1] + src_radius), (center[0] + src_radius, center[1]), + (center[0], center[1] - src_radius), (center[0] - src_radius, center[1])] + write_circle(xy_s, show_CT) + radius15_int = np.int(15.0 / spacing_dx) + dicom_row_num = 512 + dicom_col_num = 512 + half_col_num = np.int(dicom_col_num / 2) + half_row_num = np.int(dicom_row_num / 2) + fixed_points = np.array([[half_col_num, dicom_row_num], [dicom_col_num, dicom_row_num], + [dicom_col_num, half_row_num], [dicom_col_num, 0], + [half_col_num, 0], [0, 0], + [0, half_row_num], [0, dicom_row_num]], dtype=np.float32) + + after_points = [(center[0], center[1] + radius15_int), (center[0] + radius15_int, center[1]), + (center[0], center[1] - radius15_int), (center[0] - radius15_int, center[1])] + for p in after_points: + cv2.circle(show_CT, p, 4, (0, 200, 0), thickness=-1) + cv2.imwrite("before{}.png".format(dicom_num), show_CT) + xy_s = np.array(xy_s, dtype=np.float32) + after_points = np.array(after_points, dtype=np.float32) + src = np.vstack((xy_s, fixed_points)).reshape((1, -1, 2)) + target = np.vstack((after_points, fixed_points)).reshape((1, -1, 2)) + + matches = list() + for i in range(len(src[0, :, :])): + matches.append(cv2.DMatch(i, i, 0)) + + tps = cv2.createThinPlateSplineShapeTransformer() + tps.estimateTransformation(target, src, matches) + out = tps.warpImage(CT_img) + cv2.imshow("dist", out) + cv2.imwrite("after{}.png".format(dicom_num), out) diff --git a/utils.py b/utils.py new file mode 100644 index 0000000..7034d1e --- /dev/null +++ b/utils.py @@ -0,0 +1,97 @@ +import numpy as np +import cv2 +import math + +def region_growing(bin_mat, seed): + new_points = [seed] + array_for_next = [[0, 1], [1, 0], [0, -1], [-1, 0]] + + mask = np.zeros_like(bin_mat, dtype=np.bool) + mask[seed[0], seed[1]] = True + + while len(new_points) != 0: + stack_point = [] + for candidate_point in new_points: + for array in array_for_next: + next = (candidate_point[0] + array[0], candidate_point[1] + array[1]) + print(bin_mat[next[0], next[1]]) + if (bin_mat[next[0], next[1]] == 0) and (mask[next[0], next[1]] != True): + mask[next[0], next[1]] = True + stack_point.append(next) + #cv2.imshow("region", 255 * mask.astype(np.uint8)) + #cv2.waitKey(100) + new_points = stack_point + + return mask + + +def get_8_neighbor_with_mask(mask, center_rc): + bin_mat = 255 * mask.astype(np.uint8) + test = 255 * bin_mat[:, :, np.newaxis].astype(np.uint8) + test = cv2.cvtColor(test, cv2.COLOR_GRAY2BGR) + + if bin_mat[center_rc[0], center_rc[1]] == 0: + raise + bin_shape = bin_mat.shape + edges_rc = [] + + # check vertical edge + for j in range(1, bin_shape[0], 1): + v_derivative = bin_mat[j, center_rc[1]] - bin_mat[j - 1, center_rc[1]] + if v_derivative != 0: + edges_rc.append(np.array([j, center_rc[1]])) + cv2.circle(test, (center_rc[1], j), 3, (0, 0, 200), -1) + + # check horizonal edge + for i in range(1, bin_shape[1], 1): + h_derivative = bin_mat[center_rc[0], i] - bin_mat[center_rc[0], i - 1] + if h_derivative != 0: + edges_rc.append(np.array([center_rc[0], i])) + cv2.circle(test, (i, center_rc[0]), 3, (200, 0, 0), -1) + + cv2.imshow("test", test) + cv2.waitKey(1) + return edges_rc + + +def get_8neighbors_with_radius(src_radius, center_rc): + # get_vertical points + TPS_src_xy = [] + init_relative_vector = np.array([0, src_radius]) + center_xy_array = np.array([center_rc[0], center_rc[1]]) + step_theta = 45 + theta = 0 + for i in range(8): + rotate_mat = np.array([[math.cos(math.radians(theta)), -math.sin(math.radians(theta))], + [math.sin(math.radians(theta)), math.cos(math.radians(theta))]]) + src_wall_point = center_xy_array + np.dot(rotate_mat, np.transpose(init_relative_vector)) + TPS_src_xy.append(np.array([src_wall_point[0], src_wall_point[1]], dtype=np.float32)) + theta += step_theta + TPS_src_xy = np.array(TPS_src_xy, dtype=np.float32) + return TPS_src_xy + + +def thin_plate_spline(sshape, tshape, src_mat, n_matches=8): + tps = cv2.createThinPlateSplineShapeTransformer() + + sshape = np.array(sshape, dtype=np.float32) + tshape = np.array(tshape, dtype=np.float32) + sshape = sshape.reshape((1, -1, 2)) + tshape = tshape.reshape((1, -1, 2)) + + matches = list() + for i in range(n_matches): + matches.append(cv2.DMatch(i, i, 0)) + + tps.estimateTransformation(sshape, tshape, matches) + + ret, tshape_ = tps.applyTransformation(sshape) + + tps.estimateTransformation(tshape, sshape, matches) + + distorted_mat = tps.warpImage(src_mat) + return distorted_mat + + +def convert_rc_to_cv2point(rc): + return [rc[1], rc[0]]