当前位置：首页 > java >正文

CASToR 生成的文件进行转换

java 2025/9/6 6:40:38

一代码功能

实现把下列CASToR生成的文件转换为NIFTI格式文件，然后用itk-snap软件进行打开

二源码

#!/usr/bin/env python
# -*- coding: utf-8 -*-
""" 
将CT/PET图像从INTERFILE格式(hdr/img)转换为NIfTI格式的脚本
可用于将my_images-castor-output目录下的所有hdr/img文件对转换为可用ITK-SNAP打开的NIfTI文件
"""import os
import struct
import numpy as np
import SimpleITK as sitk
import argparse
from glob import globdef parse_interfile_header(hdr_path):"""解析INTERFILE格式的HDR文件，提取图像元数据信息参数:hdr_path: HDR文件路径返回:包含图像元数据的字典"""header = {}try:with open(hdr_path, 'r', encoding='utf-8', errors='ignore') as f:for line in f:line = line.strip()if not line or line.startswith('!END OF INTERFILE'):break# 解析键值对if ':=' in line:parts = line.split(':=', 1)key = parts[0].strip().lower()value = parts[1].strip()# 提取有用的元数据if 'matrix size [1]' in key:header['matrix_size_x'] = int(value)elif 'matrix size [2]' in key:header['matrix_size_y'] = int(value)elif 'matrix size [3]' in key:header['matrix_size_z'] = int(value)elif 'number of bytes per pixel' in key:header['bytes_per_pixel'] = int(value)elif 'number format' in key:header['data_type'] = value.lower()elif 'name of data file' in key:header['data_file'] = valueelif 'imagedata byte order' in key:header['byte_order'] = value.lower()elif 'scaling factor (mm/pixel) [1]' in key:header['spacing_x'] = float(value)elif 'scaling factor (mm/pixel) [2]' in key:header['spacing_y'] = float(value)elif 'scaling factor (mm/pixel) [3]' in key:header['spacing_z'] = float(value)elif 'data rescale slope' in key:header['rescale_slope'] = float(value)elif 'data rescale offset' in key:header['rescale_offset'] = float(value)except Exception as e:print(f"警告: 解析HDR文件 '{hdr_path}' 时出错: {str(e)}")# 设置默认值（如果某些键不存在）if 'matrix_size_x' not in header:header['matrix_size_x'] = 220if 'matrix_size_y' not in header:header['matrix_size_y'] = 220if 'matrix_size_z' not in header:header['matrix_size_z'] = 161if 'bytes_per_pixel' not in header:header['bytes_per_pixel'] = 4if 'data_type' not in header:header['data_type'] = 'short float'if 'rescale_slope' not in header:header['rescale_slope'] = 1.0if 'rescale_offset' not in header:header['rescale_offset'] = 0.0if 'spacing_x' not in header:header['spacing_x'] = 0.25if 'spacing_y' not in header:header['spacing_y'] = 0.25if 'spacing_z' not in header:header['spacing_z'] = 0.15return headerdef get_numpy_dtype(data_type, bytes_per_pixel):"""根据INTERFILE中的数据类型和每像素字节数确定NumPy数据类型"""if 'float' in data_type:if bytes_per_pixel == 4:return np.float32elif bytes_per_pixel == 8:return np.float64elif 'short' in data_type and 'float' not in data_type:return np.int16elif 'long' in data_type:return np.int32elif 'unsigned' in data_type:if bytes_per_pixel == 2:return np.uint16elif bytes_per_pixel == 4:return np.uint32# 默认返回float32return np.float32def read_img_data(img_path, header):"""读取IMG文件中的二进制数据"""try:# 确定数据类型dtype = get_numpy_dtype(header['data_type'], header['bytes_per_pixel'])# 计算数据大小data_size = (header['matrix_size_z'], header['matrix_size_y'], header['matrix_size_x'])# 读取二进制数据with open(img_path, 'rb') as f:data = np.fromfile(f, dtype=dtype)# 重塑数据为3D数组try:data = data.reshape(data_size)except ValueError as e:print(f"警告: 数据重塑失败，尝试自动调整维度顺序。错误: {str(e)}")# 尝试不同的维度顺序排列possible_shapes = [(header['matrix_size_z'], header['matrix_size_y'], header['matrix_size_x']),(header['matrix_size_x'], header['matrix_size_y'], header['matrix_size_z']),(header['matrix_size_x'], header['matrix_size_z'], header['matrix_size_y'])]for shape in possible_shapes:try:if np.prod(shape) == data.size:data = data.reshape(shape)print(f"成功: 使用替代维度顺序 {shape}")breakexcept:continue# 如果所有尝试都失败，使用原始数据if data.ndim != 3:print("警告: 无法正确重塑数据维度，使用原始一维数据")# 创建一个最小的3D数组以避免SimpleITK错误data = np.expand_dims(data[:header['matrix_size_x']*header['matrix_size_y']*header['matrix_size_z']], axis=0)# 应用缩放因子if header['rescale_slope'] != 1.0 or header['rescale_offset'] != 0.0:data = data * header['rescale_slope'] + header['rescale_offset']return dataexcept Exception as e:print(f"错误: 读取IMG文件 '{img_path}' 时出错: {str(e)}")return Nonedef process_single_file(hdr_path, output_dir):"""处理单个hdr/img文件对的转换参数:hdr_path: HDR文件的完整路径output_dir: 输出NIfTI文件的目录"""try:# 获取文件名（不包含扩展名）base_name = os.path.splitext(os.path.basename(hdr_path))[0]input_dir = os.path.dirname(hdr_path)print(f"正在处理单个文件: {base_name}.hdr/.img")# 解析HDR文件header = parse_interfile_header(hdr_path)# 构建对应的img文件路径if 'data_file' in header and header['data_file']:# 处理HDR中指定的数据文件名可能包含路径的情况if os.path.isabs(header['data_file']):img_path = header['data_file']else:# 尝试在同一目录下查找img_path = os.path.join(input_dir, header['data_file'])else:img_path = os.path.join(input_dir, base_name + '.img')# 检查img文件是否存在if not os.path.exists(img_path):print(f"警告: 跳过 '{hdr_path}'，因为对应的img文件 '{img_path}' 不存在")return# 读取IMG数据img_data = read_img_data(img_path, header)if img_data is None:return# 创建SimpleITK图像image = sitk.GetImageFromArray(img_data)# 设置图像间距image.SetSpacing((header['spacing_x'], header['spacing_y'], header['spacing_z']))# 构建输出文件路径output_path = os.path.join(output_dir, base_name + '.nii.gz')# 保存为NIfTI格式sitk.WriteImage(image, output_path)print(f"已保存: {output_path}")print("单个文件转换完成！")except Exception as e:print(f"错误: 处理单个文件 '{hdr_path}' 时出错: {str(e)}")def convert_hdr_img_to_nifti(input_dir, output_dir=None):"""将指定目录下的所有hdr/img文件对转换为NIfTI格式参数:input_dir: 包含hdr/img文件的输入目录output_dir: 输出NIfTI文件的目录，默认与输入目录相同"""# 检查输入目录是否存在if not os.path.exists(input_dir):print(f"错误: 输入目录 '{input_dir}' 不存在")return# 设置输出目录if output_dir is None:output_dir = input_dir# 确保输出目录存在os.makedirs(output_dir, exist_ok=True)# 查找所有hdr文件hdr_files = glob(os.path.join(input_dir, '*.hdr'))if not hdr_files:print(f"警告: 在目录 '{input_dir}' 中未找到hdr文件")returnprint(f"找到 {len(hdr_files)} 个hdr文件，开始转换...")success_count = 0fail_count = 0for hdr_path in hdr_files:try:# 调用单个文件处理函数process_single_file(hdr_path, output_dir)success_count += 1except Exception as e:print(f"错误: 处理文件 '{hdr_path}' 时出错: {str(e)}")fail_count += 1print(f"转换完成！成功: {success_count} 个文件, 失败: {fail_count} 个文件")def main():"""主函数，处理命令行参数"""parser = argparse.ArgumentParser(description='将INTERFILE格式(hdr/img)医学图像转换为NIfTI格式',formatter_class=argparse.RawTextHelpFormatter,epilog='''使用示例:\n'''''' 1. 转换默认目录下的文件，自动创建输出目录:\n''''''    python convert_to_nifti.py\n'''''' 2. 指定输入目录:\n''''''    python convert_to_nifti.py -i /path/to/input\n'''''' 3. 指定输入和输出目录:\n''''''    python convert_to_nifti.py -i /path/to/input -o /path/to/output''')# 使用当前目录作为默认输入目录，提高脚本可移植性parser.add_argument('-i', '--input_dir', default=os.path.join(os.getcwd()),help='包含hdr/img文件的输入目录 (默认: 当前目录)')parser.add_argument('-o', '--output_dir', default=None,help='输出NIfTI文件的目录 (默认: 自动在输入目录下创建"nifti_output"目录)')# 添加单个文件转换模式parser.add_argument('-f', '--file', default=None,help='单个hdr文件的路径 (使用此参数时将忽略input_dir参数)')args = parser.parse_args()# 确保路径格式兼容不同操作系统if args.file:# 处理单个文件hdr_path = os.path.abspath(args.file)input_dir = os.path.dirname(hdr_path)output_dir = args.output_dir# 如果未指定输出目录，自动创建一个if output_dir is None:output_dir_name = "nifti_output"output_dir = os.path.join(input_dir, output_dir_name)print(f"未指定输出目录，将自动创建: {output_dir}")# 确保输出目录存在os.makedirs(output_dir, exist_ok=True)# 单独处理这个文件process_single_file(hdr_path, output_dir)else:# 处理整个目录input_dir = os.path.abspath(args.input_dir)output_dir = args.output_dir# 如果未指定输出目录，自动创建一个if output_dir is None:output_dir_name = "nifti_output"output_dir = os.path.join(input_dir, output_dir_name)print(f"未指定输出目录，将自动创建: {output_dir}")else:output_dir = os.path.abspath(output_dir)convert_hdr_img_to_nifti(input_dir, output_dir)if __name__ == "__main__":main()