对CT切片进行三维重建，并把三维数据导出为STL

jllllyuz

319人浏览 · 2025-12-19 16:08:16

jllllyuz · 2025-12-19 16:08:16 发布

一、文件列表

main_ct2stl.m % 一键：读片→重建→导出
readCTstack.m % 读取 DICOM/PNG 堆栈
volumeReconstruct.m % 体素重建（HU 值）
marchingCubesSTL.m % 手写 Marching Cubes + STL 写入
plotVolume.m % 3D 可视化（isosurface）

二、流程

① 读片 → 体素矩阵 V(x,y,z)
② 阈值分割（HU 值可调，默认 200）
③ Marching Cubes → 三角网格
④ STL 二进制写入 → 文件输出

三、核心代码

主程序

%% 0 参数
clear; clc; close all
dataFolder = 'CT_slices';     % DICOM 或 PNG 文件夹
isoHU      = 200;             % 骨阈值（HU）
outFile    = 'bone_model.stl';% 导出 STL

%% 1 读片堆栈
[V, spacing] = readCTstack(dataFolder);  % spacing=[dx,dy,dz] mm
fprintf('体素尺寸: %d×%d×%d, 体素间距: %.2f mm\n', size(V), spacing(1));

%% 2 三维重建（保持 HU）
V = volumeReconstruct(V, spacing);

%% 3 Marching Cubes → 三角网格
[F, Vtx] = marchingCubesSTL(V, isoHU, spacing);
fprintf('导出三角面片: %d\n', size(F,1));

%% 4 STL 写入
writeSTL(outFile, Vtx, F, 'binary');
fprintf('已导出: %s\n', outFile);

%% 5 可视化
plotVolume(V, isoHU, spacing);

读片堆栈（DICOM/PNG 自适应）

function [V, spacing] = readCTstack(folder)
files = dir(fullfile(folder, '*.dcm'));
if isempty(files), files = dir(fullfile(folder, '*.png')); end
N = length(files);
% 读第一张拿尺寸
I0 = imread(fullfile(folder, files(1).name));
V = zeros(size(I0,1), size(I0,2), N, 'single');
for k = 1:N
    I = imread(fullfile(folder, files(k).name));
    if size(I,3)==3, I = rgb2gray(I); end
    V(:,:,k) = single(I);
end
% 体素间距（DICOM 优先读 PixelSpacing）
info = imfinfo(fullfile(folder, files(1).name));
if isfield(info, 'PixelSpacing')
    spacing = str2num(info.PixelSpacing(1));  % mm
else
    spacing = [1 1 1];  % 默认 1 mm
end
spacing(3) = 1;  % 层厚暂设 1 mm（可改）
end

体素重建（HU 校准，可选）

function V = volumeReconstruct(V, spacing)
% 简单线性 HU 映射（PNG 0-255 → HU -1000~+1000）
if max(V(:)) > 255, return; end  % 已是 HU
V = (V/255)*2000 - 1000;        % -1000 ~ +1000 HU
end

手写 Marching Cubes + STL 写入

function [F, Vtx] = marchingCubesSTL(V, iso, spacing)
[nx,ny,nz] = size(V);
Vtx = []; F = [];
% 简易 MC 表（15 种三角配置）
edgeTable = uint16([...]);  % 256 种→15 三角
for i = 1:nx-1, for j = 1:ny-1, for k = 1:nz-1
    cube = V(i:i+1, j:j+1, k:k+1);
    config = sum(cube > iso, 'all');
    if config==0 || config==8, continue; end
    % 顶点坐标（mm）
    v0 = [(i-1)*spacing(1), (j-1)*spacing(2), (k-1)*spacing(3)];
    v1 = v0 + [spacing(1),0,0]; v2 = v0 + [0,spacing(2),0];
    v3 = v0 + [spacing(1),spacing(2),0]; v4 = v0 + [0,0,spacing(3)];
    v5 = v0 + [spacing(1),0,spacing(3)]; v6 = v0 + [0,spacing(2),spacing(3)];
    v7 = v0 + spacing;
    % 插值顶点（等值面）
    verts = zeros(12,3);
    if bitand(edgeTable(config+1), 1),  verts(1,:)  = v0 + (iso-cube(1))/(cube(2)-cube(1))*(v1-v0); end
    if bitand(edgeTable(config+1), 2),  verts(2,:)  = v0 + (iso-cube(2))/(cube(3)-cube(2))*(v2-v0); end
    ... % 其余 10 条边同理（省略）
    % 三角面
    tri = triTable(config+1,:);  tri(tri==0) = [];
    tri = reshape(tri,3,[]);
    base = size(Vtx,1);
    Vtx = [Vtx; verts(tri,:)];
    F   = [F;  base + (1:size(tri,1))'];
end, end, end
end

STL 二进制写入（零依赖）

function writeSTL(filename, Vtx, F, mode)
if strcmp(mode,'binary')
    fid = fopen(filename, 'wb');
    fwrite(fid, uint8(0), 'uint8', 80);          % 80 字节头
    fwrite(fid, uint32(size(F,1)), 'uint32');    % 面数
    for i = 1:size(F,1)
        fwrite(fid, [0 0 0], 'float32');         % 法向量（可填零）
        fwrite(fid, Vtx(F(i,1),:), 'float32');
        fwrite(fid, Vtx(F(i,2),:), 'float32');
        fwrite(fid, Vtx(F(i,3),:), 'float32');
        fwrite(fid, uint16(0), 'uint16');        % 属性字（0）
    end
    fclose(fid);
else
    % ASCII 模式（略）
end
end