diff --git a/onlylian/oar4.py b/onlylian/oar4.py index 425c087..73f6984 100755 --- a/onlylian/oar4.py +++ b/onlylian/oar4.py @@ -125,6 +125,9 @@ LABEL_MAPPING: Dict[int, str] = { 5: "Right_Optic_Nerve", 6: "Left_Optic_Nerve" } + +# MAX_OAR = max(LABEL_MAPPING.keys()) +max_existing_label = max(LABEL_MAPPING.keys()) # 設定日誌 logging.basicConfig( @@ -346,7 +349,7 @@ def process_multiple_tumors(image_array): return result_img # 為每個腫瘤生成唯一的標籤值，從最大現有標籤值開始 - max_existing_label = max(LABEL_MAPPING.keys()) +# max_existing_label = max(LABEL_MAPPING.keys()) tumor_colors = generate_distinct_colors(num_tumors) # 更新標籤映射和顏色映射 @@ -374,37 +377,49 @@ def create_structure_masks(image: sitk.Image) -> Dict[str, np.ndarray]: return masks -def register(ct_fixed, moving, out_moving): +MODELS = ( + 'rigid', + 'affine', + 'joint', +) - with tempfile.TemporaryDirectory() as tmp: - fixed = os.path.join(tmp, 'clipped.nii.gz') - - ct = sf.load_volume(ct_fixed) - clipped = ct.clip(0, 80) - clipped.save(fixed) +def register(ct_fixed, moving): + + out_dir = tempfile.mkdtemp(dir=Path(__file__).resolve().parent/'0') - FREESURFER_HOME = str(Path(__file__).resolve().parent/'mri_synthmorph') + fixed = os.path.join(out_dir, 'clipped.nii.gz') + ct = sf.load_volume(ct_fixed) + clipped = ct.clip(0, 80) + clipped.save(fixed) + + FREESURFER_HOME = Path(__file__).resolve().parent/'mri_synthmorph' + my_env = os.environ.copy() + my_env["FREESURFER_HOME"] = str(FREESURFER_HOME) + + mri_synthmorph = str(FREESURFER_HOME/'mri_synthmorph') + + for m in MODELS: + + transform_extension = 'lta' if m in ('affine','rigid') else 'nii.gz' + args = [ - str(Path(__file__).resolve().parent/'mri_synthmorph/mri_synthmorph'), -# 'register', - '-m', 'affine', -# '-m', 'rigid', - '-o', out_moving, + mri_synthmorph, + '-m', m, + '-o', os.path.join(out_dir, '%s.nii.gz'%m), + + '-O', os.path.join(out_dir, 'out_fixed-%s.nii.gz'%m), + '-t', os.path.join(out_dir, 'moving_to_fixed-%s.%s'% (m, transform_extension)), + '-T', os.path.join(out_dir, 'fixed_to_moving-%s.%s'% (m, transform_extension)), + '-g', moving, fixed, ] + logger.info(f" 執行registration 指令：{' '.join(args)}") + result = subprocess.run(args, env=my_env) - my_env = os.environ.copy() - my_env["FREESURFER_HOME"] = str(Path(__file__).resolve().parent/'mri_synthmorph') - - -# logger.info(f"執行推論指令：{' '.join(predict_cmd)}") - result = subprocess.run( - args, - env=my_env, - ) + return out_dir ''' Namespace(command='register', extent=256, fixed='0/clipped.nii.gz', gpu=True, header_only=False, hyper=0.5, init=None, inverse=None, mid_space=False, model='affine', moving='0/t1c.nii.gz', out_dir=None, out_fixed=None, out_moving='0/out-aff.nii.gz', steps=7, threads=None, trans='0/trans.lta', verbose=False, weights=None) @@ -455,8 +470,8 @@ def register_synthmorph(ct_fixed, moving, out_moving): arg=argparse.Namespace(**default) FREESURFER_HOME = str(Path(__file__).resolve().parent/'mri_synthmorph') - print(arg) - print(FREESURFER_HOME) +# print(arg) +# print(FREESURFER_HOME) os.environ["FREESURFER_HOME"] = str(Path(__file__).resolve().parent/'mri_synthmorph') registration.register(arg) @@ -476,11 +491,7 @@ def inference(DCM_CT, DCM_MR): ROOT_DIR = os.path.dirname(os.path.dirname(DCM_CT)) NII_DIR = os.path.join(ROOT_DIR, 'nii') INPUT_DIR = os.path.join(ROOT_DIR, 'input') - OUTPUT_DIR = os.path.join(ROOT_DIR, 'output') - - # 設定 RTSTRUCT 輸出路徑 - head, tail = os.path.split(DCM_CT) - rtss_file = os.path.join(head, tail+'-rtss.dcm') + OUTPUT_DIR = os.path.join(ROOT_DIR, 'output') # 清理並建立工作目錄 for dir_path in [NII_DIR, INPUT_DIR, OUTPUT_DIR]: @@ -519,26 +530,26 @@ def inference(DCM_CT, DCM_MR): raise FileNotFoundError(f"找不到 CT 的 NIfTI 檔案，目錄內容：{os.listdir(NII_DIR)}") if not NII_MR: raise FileNotFoundError(f"找不到 MR 的 NIfTI 檔案，目錄內容：{os.listdir(NII_DIR)}") - + + + logger.info("Registration of %s and %s..."%(NII_CT, NII_MR)) + out_dir = register(NII_CT, NII_MR) + # 準備輸入檔案 basename = os.path.basename(NII_MR) old = '_'+basename.split('_')[-1] - input_ct = os.path.join(INPUT_DIR, basename.replace(old, '_0000.nii.gz')) - shutil.copy(NII_CT, input_ct) + for m in MODELS: + input_ct = os.path.join(INPUT_DIR, basename.replace(old, '_%s_0000.nii.gz'%m)) + shutil.copy(NII_CT, input_ct) - input_mr = os.path.join(INPUT_DIR, basename.replace(old, '_0001.nii.gz')) - - logger.info("Registration of %s and %s..."%(NII_CT, NII_MR)) - register(NII_CT, NII_MR, input_mr) + input_mr = os.path.join(INPUT_DIR, basename.replace(old, '_%s_0001.nii.gz'%m)) + shutil.copy(os.path.join(out_dir, '%s.nii.gz'%m), input_mr) + +# shutil.rmtree(out_dir) - output_file = os.path.join(OUTPUT_DIR, basename.replace(old, '.nii.gz')) - - basename_ct = os.path.basename(NII_CT) - old_ct = '_'+basename_ct.split('_')[-1] - label_file = os.path.join(NII_DIR, basename_ct.replace(old_ct, '.label.nii.gz')) # basename = os.path.basename(NII_MR) # old = '_'+basename.split('_')[-1] @@ -597,129 +608,133 @@ def inference(DCM_CT, DCM_MR): except subprocess.CalledProcessError as e: logger.error(f"模型推論失敗：\n輸出：{e.output}\n錯誤：{e.stderr}") raise RuntimeError(f"模型推論失敗：{e.stderr}") - - if not os.path.exists(output_file): - raise FileNotFoundError(f"找不到模型輸出檔案：{output_file}") - - logger.info("開始執行影像後處理...") - - try: - # 讀取預測結果 - predicted_image = sitk.ReadImage(output_file) - predicted_image = sitk.DICOMOrient(predicted_image) - predicted_array = sitk.GetArrayFromImage(predicted_image) - # 執行後處理 - processed_array = process_symmetric_structures(predicted_array) - processed_array = process_multiple_tumors(processed_array) - - processed_image = sitk.GetImageFromArray(processed_array) - processed_image.CopyInformation(predicted_image) - - # 暫時保存後處理結果 - processed_output = os.path.join(OUTPUT_DIR, 'processed_' + os.path.basename(output_file)) - sitk.WriteImage(processed_image, processed_output) - logger.info(f"後處理結果已保存至：{processed_output}") - - logger.info("開始執行影像配準...") -# reg_transform(NII_CT, NII_MR, processed_output, label_file) - shutil.copy(processed_output, label_file) - - except Exception as e: - logger.error(f"影像處理失敗：{str(e)}") - raise RuntimeError("後處理或配準失敗") - - if not os.path.exists(label_file): - raise FileNotFoundError(f"找不到配準後的標籤檔案：{label_file}") - - logger.info("開始建立 RTSTRUCT...") - - # 讀取 CT 影像序列 - try: - reader = sitk.ImageSeriesReader() - dicom_names = reader.GetGDCMSeriesFileNames(DCM_CT) - reader.SetFileNames(dicom_names) - reader.MetaDataDictionaryArrayUpdateOn() - reader.LoadPrivateTagsOn() - image = reader.Execute() - except Exception as e: - logger.error(f"讀取 CT 系列失敗：{str(e)}") - raise RuntimeError("無法讀取 CT 影像系列") - - # 讀取並重新取樣配準後的預測結果 - try: - final_image = sitk.ReadImage(label_file) - final_image = sitk.Resample( - final_image, - image, - sitk.Transform(), - sitk.sitkNearestNeighbor - ) - except Exception as e: - logger.error(f"重新取樣失敗：{str(e)}") - raise RuntimeError("預測結果重新取樣失敗") - - # 建立 RTSTRUCT - try: + rtss_combination = [] - images = [] - images.append(IntensityImage(image, Modality('CT'))) -# images.append(IntensityImage(sitk.ReadImage(input_mr), Modality('MR'))) + for m in MODELS: + + basename_ct = os.path.basename(NII_CT) + old_ct = '_'+basename_ct.split('_')[-1] + label_file = os.path.join(NII_DIR, basename_ct.replace(old_ct, '_%s.label.nii.gz'%m)) + + output_file = os.path.join(OUTPUT_DIR, basename.replace(old, '_%s.nii.gz'%m)) + + if not os.path.exists(output_file): + raise FileNotFoundError(f"找不到模型輸出檔案：{output_file}") + + logger.info("開始執行影像後處理...") + + try: + # 讀取預測結果 + predicted_image = sitk.ReadImage(output_file) + predicted_image = sitk.DICOMOrient(predicted_image) + predicted_array = sitk.GetArrayFromImage(predicted_image) + + # 執行後處理 + processed_array = process_symmetric_structures(predicted_array) + processed_array = process_multiple_tumors(processed_array) + + processed_image = sitk.GetImageFromArray(processed_array) + processed_image.CopyInformation(predicted_image) + + # 暫時保存後處理結果 + processed_output = os.path.join(OUTPUT_DIR, 'processed_' + os.path.basename(output_file)) + sitk.WriteImage(processed_image, processed_output) + logger.info(f"後處理結果已保存至：{processed_output}") + + logger.info("開始執行影像配準...") + # reg_transform(NII_CT, NII_MR, processed_output, label_file) + shutil.copy(processed_output, label_file) + + except Exception as e: + logger.error(f"影像處理失敗：{str(e)}") + raise RuntimeError("後處理或配準失敗") + + if not os.path.exists(label_file): + raise FileNotFoundError(f"找不到配準後的標籤檔案：{label_file}") + + logger.info("開始建立 RTSTRUCT...") + + # 讀取 CT 影像序列 + try: + reader = sitk.ImageSeriesReader() + dicom_names = reader.GetGDCMSeriesFileNames(DCM_CT) + reader.SetFileNames(dicom_names) + reader.MetaDataDictionaryArrayUpdateOn() + reader.LoadPrivateTagsOn() + image = reader.Execute() + except Exception as e: + logger.error(f"讀取 CT 系列失敗：{str(e)}") + raise RuntimeError("無法讀取 CT 影像系列") + + # 讀取並重新取樣配準後的預測結果 + try: + final_image = sitk.ReadImage(label_file) + final_image = sitk.Resample( + final_image, + image, + sitk.Transform(), + sitk.sitkNearestNeighbor + ) + except Exception as e: + logger.error(f"重新取樣失敗：{str(e)}") + raise RuntimeError("預測結果重新取樣失敗") + + # 建立 RTSTRUCT + try: + + images = [] + images.append(IntensityImage(image, Modality('CT'))) + # images.append(IntensityImage(sitk.ReadImage(input_mr), Modality('MR'))) + + for k, v in LABEL_MAPPING.items(): + images.append(SegmentationImage(final_image==k, v)) + + dcm_crawler = ps_io.SubjectDicomCrawler(DCM_CT, + modality_extractor=ExampleModalityExtractor(), + ) + dicom_series_info = dcm_crawler.execute() + + subject = Subject(dicom_series_info[0].get_patient_id(), images) + + reference_modality = 'CT' + + use_3d_conversion = True + if use_3d_conversion: + conv_conf = ps_io.RTSSConverter3DConfiguration() + else: + conv_conf = ps_io.RTSSConverter2DConfiguration() + meta_data = ps_io.RTSSMetaData( + series_description = '%s by %s'%(m, type(conv_conf).__name__) + ) + + converter = ps_io.SubjectToRTSSConverter( + subject, + dicom_series_info, + reference_modality, + conv_conf, + meta_data, + # colors = tuple(LABEL_COLORS.values()), + ) + + rtss = converter.convert() + + # 設定 RTSTRUCT 輸出路徑 + head, tail = os.path.split(DCM_CT) + rtss_file = os.path.join(head, tail+'_%s-rtss.dcm'%m) + + # Write the DICOM-RTSS to a separate subject directory + # and include the DICOM files crawled before + # Note: If you want to output just a subset of the + # original DICOM files you may use additional selectors + + output_dir_path = os.path.join(ROOT_DIR, 'pyradise') + rtss_filename = os.path.basename(rtss_file) - for k, v in LABEL_MAPPING.items(): - images.append(SegmentationImage(final_image==k, v)) + rtss_combination.append((rtss_filename, rtss)) - dcm_crawler = ps_io.SubjectDicomCrawler(DCM_CT, - modality_extractor=ExampleModalityExtractor(), - ) - dicom_series_info = dcm_crawler.execute() - - subject = Subject(dicom_series_info[0].get_patient_id(), images) - - reference_modality = 'CT' - - use_3d_conversion = True - if use_3d_conversion: - conv_conf = ps_io.RTSSConverter3DConfiguration() - else: - conv_conf = ps_io.RTSSConverter2DConfiguration() - meta_data = ps_io.RTSSMetaData( - series_description = '%s by onlylian'%type(conv_conf).__name__) - - converter = ps_io.SubjectToRTSSConverter( - subject, - dicom_series_info, - reference_modality, - conv_conf, - meta_data, -# colors = tuple(LABEL_COLORS.values()), - ) - - rtss = converter.convert() - - - # Write the DICOM-RTSS to a separate subject directory - # and include the DICOM files crawled before - # Note: If you want to output just a subset of the - # original DICOM files you may use additional selectors - - output_dir_path = os.path.join(ROOT_DIR, 'pyradise') - rtss_filename = os.path.basename(rtss_file) - rtss_combination = ((rtss_filename, rtss),) - print(output_dir_path) - print(rtss_combination) - -# shutil.rmtree(output_dir_path, ignore_errors=True) - os.makedirs(output_dir_path, exist_ok=True) - - writer = ps_io.DicomSeriesSubjectWriter() - writer.write(rtss_combination, - output_dir_path, -# subject.get_name(), - None, - dicom_series_info,) - - shutil.copy(os.path.join(output_dir_path, rtss_filename), rtss_file) + +# rtss_combination = ((rtss_filename, rtss),) # rtstruct = RTStructBuilder.create_new(dicom_series_path=DCM_CT) @@ -748,16 +763,30 @@ def inference(DCM_CT, DCM_MR): # logger.info(f"儲存 RTSTRUCT：{rtss_file}") # rtstruct.save(rtss_file) - except Exception as e: - logger.error(f"RTSTRUCT 生成失敗：{str(e)}") - raise RuntimeError("無法生成或儲存 RTSTRUCT") - - return rtss_file + except Exception as e: + logger.error(f"RTSTRUCT 生成失敗：{str(e)}") + raise RuntimeError("無法生成或儲存 RTSTRUCT") + +# print(output_dir_path) +# print(rtss_combination) + +# shutil.rmtree(output_dir_path, ignore_errors=True) + os.makedirs(output_dir_path, exist_ok=True) + + writer = ps_io.DicomSeriesSubjectWriter() + writer.write(rtss_combination, + output_dir_path, +# subject.get_name(), + None, + dicom_series_info,) except Exception as e: logger.error(f"推論過程發生錯誤：{str(e)}", exc_info=True) raise + return rtss_combination, output_dir_path + + def SendDCM(fp): """傳送 DICOM 檔案到 PACS""" logger.info(f"準備傳送 DICOM 檔案：{fp}") @@ -803,6 +832,62 @@ def SendDCM(fp): logger.error(f"DICOM 傳送過程發生錯誤：{str(e)}", exc_info=True) raise + +def SendDCMs(filelist): + """傳送 DICOM 檔案到 PACS""" + logger.info(f"準備傳送 DICOM 檔案：{' '.join(filelist)}") + debug_logger() + + dslist = [] + + for fp in filelist: + + if not os.path.exists(fp): + raise FileNotFoundError(f"找不到 DICOM 檔案：{fp}") + + try: + ds = dcmread(fp) + dslist.append(ds) + except Exception as e: + logger.error(f"讀取 DICOM 檔案失敗：{str(e)}") + raise RuntimeError("無法讀取 DICOM 檔案") + + + + try: + ae = AE() + ae.ae_title = 'OUR_STORE_SCP' + ae.add_requested_context(RTStructureSetStorage) + + + logger.info("正在連接 PACS 伺服器...") + try: + assoc = ae.associate("172.16.40.36", 104, ae_title='N1000_STORAGE') + except Exception as e: + logger.error(f"PACS 連接失敗：{str(e)}") + raise RuntimeError("無法連接 PACS 伺服器") + + if assoc.is_established: + try: + for ds in dslist: + status = assoc.send_c_store(ds) + if status: + logger.info(f'DICOM 傳送成功，狀態碼: 0x{status.Status:04x}') + else: + raise RuntimeError('傳送失敗：連接逾時或收到無效回應') + except Exception as e: + logger.error(f"DICOM 傳送失敗：{str(e)}") + raise + finally: + assoc.release() + else: + raise RuntimeError('傳送失敗：無法建立連接') + + except Exception as e: + logger.error(f"DICOM 傳送過程發生錯誤：{str(e)}", exc_info=True) + raise + + def main(): """主程式進入點""" if len(sys.argv) < 3: @@ -817,8 +902,18 @@ def main(): start_time = time.time() try: - rtss_file = inference(sys.argv[1], sys.argv[2]) - SendDCM(rtss_file) + + rtss_combination, output_dir_path = inference(sys.argv[1], sys.argv[2]) + + RTSS = [] + for rtss_filename, rtss in rtss_combination: + rtss_file = os.path.join(output_dir_path, rtss_filename) + shutil.copy(rtss_file, Path(sys.argv[1]).resolve().parent) + RTSS.append(rtss_file) +# SendDCM(rtss_file) + + SendDCMs(RTSS) + logger.info(f"處理完成，總耗時：{time.time() - start_time:.2f} 秒") except Exception as e: logger.error(f"處理過程發生錯誤：{str(e)}", exc_info=True)