一 代码功能
实现把下列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 glob
def 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'] = value
elif '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'] = 220
if 'matrix_size_y' not in header:
header['matrix_size_y'] = 220
if 'matrix_size_z' not in header:
header['matrix_size_z'] = 161
if 'bytes_per_pixel' not in header:
header['bytes_per_pixel'] = 4
if 'data_type' not in header:
header['data_type'] = 'short float'
if 'rescale_slope' not in header:
header['rescale_slope'] = 1.0
if 'rescale_offset' not in header:
header['rescale_offset'] = 0.0
if 'spacing_x' not in header:
header['spacing_x'] = 0.25
if 'spacing_y' not in header:
header['spacing_y'] = 0.25
if 'spacing_z' not in header:
header['spacing_z'] = 0.15
return header
def get_numpy_dtype(data_type, bytes_per_pixel):
"""
根据INTERFILE中的数据类型和每像素字节数确定NumPy数据类型
"""
if 'float' in data_type:
if bytes_per_pixel == 4:
return np.float32
elif bytes_per_pixel == 8:
return np.float64
elif 'short' in data_type and 'float' not in data_type:
return np.int16
elif 'long' in data_type:
return np.int32
elif 'unsigned' in data_type:
if bytes_per_pixel == 2:
return np.uint16
elif bytes_per_pixel == 4:
return np.uint32
# 默认返回float32
return np.float32
def 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}")
break
except:
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 data
except Exception as e:
print(f"错误: 读取IMG文件 '{img_path}' 时出错: {str(e)}")
return None
def 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文件")
return
print(f"找到 {len(hdr_files)} 个hdr文件,开始转换...")
success_count = 0
fail_count = 0
for hdr_path in hdr_files:
try:
# 调用单个文件处理函数
process_single_file(hdr_path, output_dir)
success_count += 1
except Exception as e:
print(f"错误: 处理文件 '{hdr_path}' 时出错: {str(e)}")
fail_count += 1
print(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()