CT切片三维重建与STL导出MATLAB实现

读取CT切片、进行三维重建并将结果导出为STL文件格式。

function ct_to_stl_reconstruction(dicomDir, outputSTL, threshold, smoothingFactor)
    % CT切片三维重建与STL导出
    % 输入参数:
    %   dicomDir - DICOM文件目录路径
    %   outputSTL - 输出的STL文件名
    %   threshold - 阈值用于分割(默认-500)
    %   smoothingFactor - 平滑因子(默认1.5)
    
    % 设置默认参数
    if nargin < 3
        threshold = -500; % 适用于骨骼分割的典型阈值
    end
    if nargin < 4
        smoothingFactor = 1.5; % 平滑因子
    end
    
    % 检查输出文件扩展名
    if ~contains(outputSTL, '.stl')
        outputSTL = [outputSTL, '.stl'];
    end
    
    fprintf('开始处理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');
end

function [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:numSlices
        filename = fullfile(dicomDir, dicomFiles(i).name);
        volume(:, :, i) = dicomread(filename);
        
        % 显示进度
        if mod(i, round(numSlices/10)) == 0
            fprintf('已读取 %d/%d 个切片\n', i, numSlices);
        end
    end
    
    % 存储空间信息
    spatialInfo.PixelSpacing = info.PixelSpacing;
    if isfield(info, 'SliceThickness')
        spatialInfo.SliceThickness = info.SliceThickness;
    else
        spatialInfo.SliceThickness = 1.0; % 默认值
    end
    spatialInfo.ImagePositionPatient = info.ImagePositionPatient;
end

function 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);
end

function segmentedVolume = segmentVolume(volume, threshold)
    % 分割感兴趣区域
    
    % 二值化分割
    segmentedVolume = volume > threshold;
    
    % 形态学操作去除小噪声点
    fprintf('应用形态学操作...\n');
    segmentedVolume = bwareaopen(segmentedVolume, 100); % 移除小区域
    
    % 填充空洞
    segmentedVolume = imfill(segmentedVolume, 'holes');
end

function [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 > 0
        fprintf('平滑网格...\n');
        vertices = smoothMesh(vertices, faces, smoothingFactor);
    end
    
    % 减少面片数量（可选，用于大型数据集）
    fprintf('简化网格...\n');
    [faces, vertices] = reducepatch(faces, vertices, 0.5); % 减少50%的面片
end

function 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);
        end
        smoothedVertices(:, i) = column;
    end
end

function exportToSTL(faces, vertices, filename)
    % 导出为STL文件
    
    % 确保顶点是单精度（STL标准要求）
    vertices = single(vertices);
    
    % 创建STL文件
    stlwrite(filename, faces, vertices);
    
    fprintf('STL文件已保存: %s\n', filename);
end

function 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 == -1
        error('无法创建文件: %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

使用说明

1. 准备工作：

   • 确保您的MATLAB安装了Image Processing Toolbox
   • 准备一组DICOM格式的CT切片，放在一个目录中

2. 基本用法：

   % 设置DICOM文件目录和输出文件名
   dicomDir = 'C:\CT_Data\Patient1\';
   outputFile = 'bone_reconstruction.stl';
   
   % 运行三维重建
   ct_to_stl_reconstruction(dicomDir, outputFile);
   

3. 高级用法：

   % 使用自定义参数
   dicomDir = 'C:\CT_Data\Patient1\';
   outputFile = 'bone_reconstruction.stl';
   threshold = -400;      % 分割阈值，调整以捕获不同组织
   smoothing = 2.0;       % 平滑因子，值越大表面越平滑
   
   ct_to_stl_reconstruction(dicomDir, outputFile, threshold, smoothing);
   

4. 参数说明：

   • threshold: 用于分割的HU值阈值。骨骼通常使用400-1000，软组织使用较低的值
   • smoothingFactor: 控制表面平滑程度，值越大表面越平滑

功能特点

1. 完整的处理流程：

   • DICOM文件读取与解析
   • 体积数据预处理（滤波、对比度增强）
   • 基于阈值的图像分割
   • 三维表面重建（移动立方体算法）
   • STL文件导出

2. 可视化功能：

   • 原始CT切片显示
   • 分割结果展示
   • 3D重建模型可视化

3. 可调参数：

   • 分割阈值可调整以捕获不同组织
   • 平滑因子控制表面光滑度
   • 网格简化选项减少文件大小

参考代码 对CT切片进行三维重建，并把三维数据导出为STL www.youwenfan.com/contentcnk/64396.html

扩展功能

如果需要进一步扩展此代码，可以考虑：

1. 添加多种组织的同时分割和重建
2. 实现更先进的分割算法（如区域生长、水平集等）
3. 添加测量功能（如体积计算、距离测量）
4. 支持其他3D文件格式导出（如OBJ、PLY等）

这个实现提供了一个完整的从CT切片到3D打印模型的流程，适用于医学影像处理、生物医学工程和3D打印应用。