读取CT切片、进行三维重建并将结果导出为STL文件格式。
function ct_to_stl_reconstruction(dicomDir, outputSTL, threshold, smoothingFactor)% CT切片三维重建与STL导出% 输入参数:% dicomDir - DICOM文件目录路径% outputSTL - 输出的STL文件名% threshold - 阈值用于分割(默认-500)% smoothingFactor - 平滑因子(默认1.5)% 设置默认参数if nargin < 3threshold = -500; % 适用于骨骼分割的典型阈值endif nargin < 4smoothingFactor = 1.5; % 平滑因子end% 检查输出文件扩展名if ~contains(outputSTL, '.stl')outputSTL = [outputSTL, '.stl'];endfprintf('开始处理CT切片数据...\n');% 步骤1: 读取DICOM文件fprintf('读取DICOM文件...\n');[volume, spatialInfo] = readDicomSeries(dicomDir);% 步骤2: 预处理体积数据fprintf('预处理体积数据...\n');processedVolume = preprocessVolume(volume);% 步骤3: 分割感兴趣区域fprintf('分割感兴趣区域...\n');segmentedVolume = segmentVolume(processedVolume, threshold);% 步骤4: 三维重建fprintf('进行三维重建...\n');[faces, vertices] = createSurfaceMesh(segmentedVolume, spatialInfo, smoothingFactor);% 步骤5: 导出为STL文件fprintf('导出为STL文件: %s...\n', outputSTL);exportToSTL(faces, vertices, outputSTL);% 步骤6: 可视化结果fprintf('可视化结果...\n');visualizeResults(volume, segmentedVolume, faces, vertices);fprintf('处理完成!\n');
endfunction [volume, spatialInfo] = readDicomSeries(dicomDir)% 读取DICOM系列文件% 获取DICOM文件列表dicomFiles = dir(fullfile(dicomDir, '*.dcm'));if isempty(dicomFiles)error('在指定目录中未找到DICOM文件');end% 读取第一个文件获取信息firstFile = fullfile(dicomDir, dicomFiles(1).name);info = dicominfo(firstFile);% 确定图像尺寸rows = info.Rows;cols = info.Columns;numSlices = length(dicomFiles);% 初始化体积数据volume = zeros(rows, cols, numSlices, 'int16');% 读取所有切片for i = 1:numSlicesfilename = fullfile(dicomDir, dicomFiles(i).name);volume(:, :, i) = dicomread(filename);% 显示进度if mod(i, round(numSlices/10)) == 0fprintf('已读取 %d/%d 个切片\n', i, numSlices);endend% 存储空间信息spatialInfo.PixelSpacing = info.PixelSpacing;if isfield(info, 'SliceThickness')spatialInfo.SliceThickness = info.SliceThickness;elsespatialInfo.SliceThickness = 1.0; % 默认值endspatialInfo.ImagePositionPatient = info.ImagePositionPatient;
endfunction processedVolume = preprocessVolume(volume)% 预处理体积数据% 转换为双精度浮点数processedVolume = double(volume);% 应用高斯滤波减少噪声fprintf('应用高斯滤波...\n');processedVolume = imgaussfilt3(processedVolume, 1);% 对比度增强fprintf('增强对比度...\n');minVal = min(processedVolume(:));maxVal = max(processedVolume(:));processedVolume = (processedVolume - minVal) / (maxVal - minVal);
endfunction segmentedVolume = segmentVolume(volume, threshold)% 分割感兴趣区域% 二值化分割segmentedVolume = volume > threshold;% 形态学操作去除小噪声点fprintf('应用形态学操作...\n');segmentedVolume = bwareaopen(segmentedVolume, 100); % 移除小区域% 填充空洞segmentedVolume = imfill(segmentedVolume, 'holes');
endfunction [faces, vertices] = createSurfaceMesh(volume, spatialInfo, smoothingFactor)% 创建表面网格fprintf('生成等值面...\n');% 计算体素间距pixelSpacing = spatialInfo.PixelSpacing;sliceThickness = spatialInfo.SliceThickness;% 使用移动立方体算法生成等值面[faces, vertices] = isosurface(volume, 0.5);% 调整顶点坐标以反映实际物理尺寸vertices(:, 1) = vertices(:, 1) * pixelSpacing(1);vertices(:, 2) = vertices(:, 2) * pixelSpacing(2);vertices(:, 3) = vertices(:, 3) * sliceThickness;% 应用网格平滑if smoothingFactor > 0fprintf('平滑网格...\n');vertices = smoothMesh(vertices, faces, smoothingFactor);end% 减少面片数量（可选，用于大型数据集）fprintf('简化网格...\n');[faces, vertices] = reducepatch(faces, vertices, 0.5); % 减少50%的面片
endfunction smoothedVertices = smoothMesh(vertices, faces, factor)% 平滑网格% 创建邻接矩阵numVertices = size(vertices, 1);adjacency = sparse([faces(:,1); faces(:,2); faces(:,3)], ...[faces(:,2); faces(:,3); faces(:,1)], ...1, numVertices, numVertices);% 确保对称adjacency = max(adjacency, adjacency');% 计算拉普拉斯平滑smoothedVertices = vertices;for i = 1:size(vertices, 2)column = vertices(:, i);for iter = 1:3 % 迭代次数neighborAvg = (adjacency * column) ./ max(1, sum(adjacency, 2));column = column + factor * (neighborAvg - column);endsmoothedVertices(:, i) = column;end
endfunction exportToSTL(faces, vertices, filename)% 导出为STL文件% 确保顶点是单精度（STL标准要求）vertices = single(vertices);% 创建STL文件stlwrite(filename, faces, vertices);fprintf('STL文件已保存: %s\n', filename);
endfunction visualizeResults(originalVolume, segmentedVolume, faces, vertices)% 可视化结果figure('Position', [100, 100, 1200, 500]);% 显示原始CT切片subplot(1, 3, 1);midSlice = round(size(originalVolume, 3) / 2);imshow(originalVolume(:, :, midSlice), []);title('原始CT切片');colormap(gray);% 显示分割结果subplot(1, 3, 2);imshow(segmentedVolume(:, :, midSlice));title('分割结果');% 显示3D重建结果subplot(1, 3, 3);patch('Faces', faces, 'Vertices', vertices, ...'FaceColor', [0.8, 0.8, 1.0], 'EdgeColor', 'none');axis equal;grid on;title('3D重建');xlabel('X (mm)');ylabel('Y (mm)');zlabel('Z (mm)');view(3);camlight;lighting gouraud;% 添加旋转按钮uicontrol('Style', 'pushbutton', 'String', '旋转', ...'Position', [20, 20, 60, 30], ...'Callback', @(src,evt) rotate3d on);
end% 辅助函数：STL写入（如果MATLAB版本较旧，可能没有stlwrite函数）
function stlwrite(filename, faces, vertices)% 写入STL文件（二进制格式）% 打开文件fid = fopen(filename, 'wb');if fid == -1error('无法创建文件: %s', filename);end% 写入80字节的头部header = sprintf('MATLAB-generated STL from CT data. Created: %s', datestr(now));header = [header, zeros(1, 80 - length(header), 'uint8')];fwrite(fid, header, 'uint8');% 写入面片数量numFaces = size(faces, 1);fwrite(fid, numFaces, 'uint32');% 写入每个面片for i = 1:numFaces% 获取当前面片的顶点v1 = vertices(faces(i, 1), :);v2 = vertices(faces(i, 2), :);v3 = vertices(faces(i, 3), :);% 计算法向量normal = cross(v2 - v1, v3 - v1);normal = normal / norm(normal);% 写入法向量fwrite(fid, normal, 'single');% 写入顶点fwrite(fid, v1, 'single');fwrite(fid, v2, 'single');fwrite(fid, v3, 'single');% 写入属性字节计数（通常为0）fwrite(fid, 0, 'uint16');end% 关闭文件fclose(fid);
end% 使用示例函数
function example_usage()% 示例用法% 设置DICOM目录和输出文件名dicomDir = 'path/to/your/dicom/files'; % 替换为您的DICOM文件目录outputSTL = 'ct_reconstruction.stl';% 调用主函数ct_to_stl_reconstruction(dicomDir, outputSTL);% 或者使用自定义参数% ct_to_stl_reconstruction(dicomDir, outputSTL, -400, 2.0);
end
使用说明
-
准备工作：
- 确保您的MATLAB安装了Image Processing Toolbox
- 准备一组DICOM格式的CT切片，放在一个目录中
-
基本用法：
% 设置DICOM文件目录和输出文件名 dicomDir = 'C:\CT_Data\Patient1\'; outputFile = 'bone_reconstruction.stl';% 运行三维重建 ct_to_stl_reconstruction(dicomDir, outputFile); -
高级用法：
% 使用自定义参数 dicomDir = 'C:\CT_Data\Patient1\'; outputFile = 'bone_reconstruction.stl'; threshold = -400; % 分割阈值，调整以捕获不同组织 smoothing = 2.0; % 平滑因子，值越大表面越平滑ct_to_stl_reconstruction(dicomDir, outputFile, threshold, smoothing); -
参数说明：
threshold: 用于分割的HU值阈值。骨骼通常使用400-1000，软组织使用较低的值smoothingFactor: 控制表面平滑程度，值越大表面越平滑
功能特点
-
完整的处理流程：
- DICOM文件读取与解析
- 体积数据预处理（滤波、对比度增强）
- 基于阈值的图像分割
- 三维表面重建（移动立方体算法）
- STL文件导出
-
可视化功能：
- 原始CT切片显示
- 分割结果展示
- 3D重建模型可视化
-
可调参数：
- 分割阈值可调整以捕获不同组织
- 平滑因子控制表面光滑度
- 网格简化选项减少文件大小
参考代码 对CT切片进行三维重建，并把三维数据导出为STL www.youwenfan.com/contentcnk/64396.html
扩展功能
如果需要进一步扩展此代码，可以考虑：
- 添加多种组织的同时分割和重建
- 实现更先进的分割算法（如区域生长、水平集等）
- 添加测量功能（如体积计算、距离测量）
- 支持其他3D文件格式导出（如OBJ、PLY等）
这个实现提供了一个完整的从CT切片到3D打印模型的流程，适用于医学影像处理、生物医学工程和3D打印应用。
)
