已有数据

现有大小为 137*299*348 的数字人体素数据，不同数字代表不同器官，0 代表空气，在文本文档中依次列出。

目的

将三维体素数据转换为STL网格模型，希望实现该数字人的可视化，便于判断数据质量以及为后续实验设计作参考，stl 文件还可以实现 3D打印。

实现方法

1. 参数配置与初始化

• 使用灵活的参数输入系统，支持缩放、分辨率、简化因子等多种处理选项

• 设置默认参数值，确保基本功能无需复杂配置即可使用

2. 数据读取与预处理

• 高效数据读取：采用缓冲区分块读取策略，适合处理大型文本格式体素数据

• 内存优化：使用uint8数据类型减少内存占用

3. 组织筛选与掩码创建

• 支持选择性提取特定组织类型（通过OrganSelect参数）

• 创建二值掩码标识目标体素区域

4. 数据降采样（可选）

• 三维体数据降采样处理，大幅提高处理速度

• 保持主要结构特征的同时减少数据量

5. 表面网格提取

• 使用MATLAB内置isosurface算法提取等值面

• 应用物理分辨率参数，确保输出模型尺寸准确

6. 网格后处理与优化

• 法线统一：确保所有面法线方向一致（朝外）

• 网格简化：按比例减少三角形数量，优化文件大小

• 平滑处理：迭代平滑算法改善表面质量

• Z轴方向修正：解决医学影像常见的朝向问题

7. 尺寸缩放与输出

• 支持各轴向独立缩放

• 高效二进制STL文件输出，包含进度显示

function voxelToSTL(inputFile, outputFile, options)
    % 参数设置
    arguments
        inputFile string
        outputFile string
        options.Scale (1,3) = [1, 1, 1]   % 缩放比例 (x,y,z)
        options.Resolution (1,1) = 4.84       % 体素实际尺寸(mm)
        options.SimplifyFactor (1,1) = 0.5 % 网格简化比例(0-1)
        options.SmoothIter (1,1) = 0       % 平滑迭代次数
        options.OrganSelect double = []    % 选择特定组织
        options.DownsampleFactor (1,1) = 1 % 降采样因子(1=不降采样)
        options.FlipZ (1,1) = true         % 修正Z轴方向(默认开启)
    end
    
    % 记录开始时间
    tic;
    fprintf('开始处理体素数据...\n');
    
    % 1. 快速读取文本数据
    [voxelData, totalVoxels] = readVoxelDataFast(inputFile);
    
    % 2. 创建组织掩码
    if isempty(options.OrganSelect)
        mask = voxelData > 0;
    else
        mask = ismember(voxelData, options.OrganSelect);
    end
    
    % 3. 降采样处理(如果需要)
    if options.DownsampleFactor > 1
        fprintf('降采样处理(因子: %d)...\n', options.DownsampleFactor);
        mask = downsampleVolume(mask, options.DownsampleFactor);
    end
    
    % 4. 提取表面网格(使用优化的方法)
    [f, v] = extractSurfaceOptimized(mask, options.Resolution);
    
    % 5. 修正Z轴方向(解决头朝下问题)
    if options.FlipZ
        fprintf('修正模型方向(Z轴翻转)...\n');
        v(:, 3) = max(v(:, 3)) - v(:, 3);  % 翻转Z轴
    end
    
    % 6. 网格优化
    if options.SimplifyFactor < 1
        fprintf('简化网格 (因子: %.1f)...\n', options.SimplifyFactor);
        [f, v] = reducepatch(f, v, options.SimplifyFactor);
    end
    
    if options.SmoothIter > 0
        fprintf('平滑网格 (迭代: %d)...\n', options.SmoothIter);
        v = smoothmesh(v, f, options.SmoothIter);
    end
    
    % 7. 应用缩放
    if any(options.Scale ~= 1)
        fprintf('应用缩放 (X:%.2f, Y:%.2f, Z:%.2f)...\n', ...
                options.Scale(1), options.Scale(2), options.Scale(3));
        v = v .* options.Scale;
    end
    
    % 8. 写入STL文件(使用优化的写入方法)
    fprintf('写入STL文件: %s...\n', outputFile);
    stlwrite_optimized(outputFile, f, v);
    
    % 计算处理时间
    elapsedTime = toc;
    fprintf('\n处理完成!\n');
    fprintf('总处理时间: %.2f 分钟\n', elapsedTime/60);
    fprintf('输入体素数: %s\n', formatNumber(totalVoxels));
    fprintf('输出三角形数: %s\n', formatNumber(size(f, 1)));
    fprintf('输出文件大小: %.2f MB\n', dir(outputFile).bytes/1e6);
end

%% 辅助函数
function [voxelData, totalVoxels] = readVoxelDataFast(inputFile)
    % 使用高效方法读取数据
    tic;
    fprintf('读取体素数据...\n');
    
    % 设置维度
    xSize = 137; % 每"逻辑行"的元素数量
    ySize = 299; % 行数
    zSize = 348; % 层数
    
    % 预分配三维数组
    voxelData = zeros(xSize, ySize, zSize, 'uint8');
    
    % 打开文件
    fid = fopen(inputFile, 'r');
    
    % 创建进度条
    h = waitbar(0, '读取体素数据...', 'Name', '处理进度');
    
    % 初始化计数器
    totalElementsRead = 0;
    currentBuffer = [];
    
    % 逐层读取数据
    for z = 1:zSize
        % 更新进度条
        if mod(z, 10) == 0
            waitbar(z/zSize, h, sprintf('处理中: %d/%d 层 (%.1f%%)', ...
                z, zSize, z/zSize*100));
        end
        
        % 逐行读取当前层
        for y = 1:ySize
            % 我们需要读取xSize(137)个元素
            neededElements = xSize;
            
            % 从缓冲区和文件中读取足够的数据
            while length(currentBuffer) < neededElements && ~feof(fid)
                % 读取一行数据
                line = fgetl(fid);
                
                % 将行数据转换为数值数组
                numbers = sscanf(line, '%d')';
                
                % 添加到缓冲区
                currentBuffer = [currentBuffer, numbers];
            end
            
            % 从缓冲区提取所需数量的元素
            if length(currentBuffer) >= neededElements
                % 提取137个元素
                elements = currentBuffer(1:neededElements);
                
                % 从缓冲区移除这些元素
                currentBuffer = currentBuffer(neededElements+1:end);
                
                % 填充到三维数组
                voxelData(:, y, z) = elements;
                
                % 更新计数器
                totalElementsRead = totalElementsRead + neededElements;
            else
                error('文件中的数据不足，无法填充整个三维数组');
            end
        end
    end
    
    % 关闭文件和进度条
    fclose(fid);
    close(h);
    
    % 验证读取的元素数量
    expectedElements = xSize * ySize * zSize;
    if totalElementsRead ~= expectedElements
        warning('读取的元素数量(%d)与预期(%d)不匹配', totalElementsRead, expectedElements);
    end
    
    % 计算总体素数
    totalVoxels = numel(voxelData);
    fprintf('读取完成! 耗时: %.2f 秒\n', toc);
    fprintf('读取了 %d 个元素，填充到 %dx%dx%d 数组中\n', ...
        totalElementsRead, xSize, ySize, zSize);
end

function [faces, vertices] = extractSurfaceOptimized(mask, resolution)
    % 使用优化的等值面提取方法
    tic;
    fprintf('提取表面网格...\n');
    
    % 创建网格坐标
    [X, Y, Z] = meshgrid(1:size(mask,2), 1:size(mask,1), 1:size(mask,3));
    
    % 应用实际分辨率
    X = X * resolution;
    Y = Y * resolution;
    Z = Z * resolution;
    
    % 使用优化的等值面提取方法
    [faces, vertices] = isosurface(X, Y, Z, mask, 0.5);
    
    % 确保法线方向一致
    faces = unifyMeshNormals(faces, vertices);
    
    fprintf('表面提取完成! 顶点数: %s, 面数: %s, 耗时: %.2f 秒\n', ...
        formatNumber(size(vertices, 1)), formatNumber(size(faces, 1)), toc);
end

function volDown = downsampleVolume(vol, factor)
    % 三维体数据降采样
    if factor <= 1
        volDown = vol;
        return;
    end
    
    % 计算新尺寸
    newSize = floor(size(vol) / factor);
    volDown = false(newSize);
    
    % 创建进度条
    h = waitbar(0, '降采样处理...', 'Name', '处理进度');
    
    % 降采样处理
    for z = 1:newSize(3)
        % 更新进度条
        if mod(z, 10) == 0
            waitbar(z/newSize(3), h, sprintf('降采样: %d/%d 层 (%.1f%%)', ...
                z, newSize(3), z/newSize(3)*100));
        end
        
        for y = 1:newSize(2)
            for x = 1:newSize(1)
                % 提取当前块
                block = vol((x-1)*factor+1:x*factor, ...
                           (y-1)*factor+1:y*factor, ...
                           (z-1)*factor+1:z*factor);
                
                % 如果块中有任何组织存在，则设置为组织
                volDown(x, y, z) = any(block(:));
            end
        end
    end
    
    % 关闭进度条
    close(h);
end

function faces = unifyMeshNormals(faces, vertices)
    % 确保所有面法线方向一致（向外）
    % 计算面法线
    v1 = vertices(faces(:,2),:) - vertices(faces(:,1),:);
    v2 = vertices(faces(:,3),:) - vertices(faces(:,1),:);
    faceNormals = cross(v1, v2, 2);
    
    % 计算中心点
    centers = (vertices(faces(:,1),:) + vertices(faces(:,2),:) + vertices(faces(:,3),:)) / 3;
    
    % 计算包围盒中心
    bboxCenter = mean([min(vertices); max(vertices)]);
    
    % 计算从面中心指向包围盒中心的向量
    toCenter = bboxCenter - centers;
    
    % 检查法线方向（如果法线与toCenter的夹角小于90度，则方向大致向外）
    dotProducts = sum(faceNormals .* toCenter, 2);
    
    % 反转方向错误的法线
    flipFaces = dotProducts > 0;
    faces(flipFaces, :) = faces(flipFaces, [1,3,2]);
end

function stlwrite_optimized(filename, faces, vertices)
    % 优化的STL写入函数（使用二进制格式）
    tic;
    fprintf('写入STL文件...\n');
    
    % 打开文件
    fid = fopen(filename, 'wb');
    
    % 写入80字节头文件
    header = sprintf('Generated from Voxel Data on %s', datestr(now));
    header = [header, repmat(' ', 1, 80-length(header))];
    fwrite(fid, header, 'char*1');
    
    % 写入三角形数量
    numFaces = size(faces, 1);
    fwrite(fid, numFaces, 'uint32');
    
    % 创建进度条
    h = waitbar(0, '写入STL文件...', 'Name', '处理进度');
    
    % 计算法线（批量计算提高性能）
    v1 = vertices(faces(:,2),:) - vertices(faces(:,1),:);
    v2 = vertices(faces(:,3),:) - vertices(faces(:,1),:);
    normals = cross(v1, v2, 2);
    normals = normals ./ vecnorm(normals, 2, 2); % 归一化
    
    % 写入每个三角形
    for i = 1:numFaces
        % 更新进度条
        if mod(i, 10000) == 0
            waitbar(i/numFaces, h, sprintf('写入: %d/%d 三角形 (%.1f%%)', ...
                i, numFaces, i/numFaces*100));
        end
        
        % 写入法线
        fwrite(fid, normals(i,:), 'float32');
        
        % 写入顶点
        fwrite(fid, vertices(faces(i,1),:), 'float32');
        fwrite(fid, vertices(faces(i,2),:), 'float32');
        fwrite(fid, vertices(faces(i,3),:), 'float32');
        
        % 属性字节计数 (设为0)
        fwrite(fid, 0, 'uint16');
    end
    
    % 关闭文件和进度条
    fclose(fid);
    close(h);
    
    fprintf('STL写入完成! 耗时: %.2f 秒\n', toc);
end

function str = formatNumber(num)
    % 格式化大数字显示
    if num >= 1e9
        str = sprintf('%.1fB', num/1e9);
    elseif num >= 1e6
        str = sprintf('%.1fM', num/1e6);
    elseif num >= 1e3
        str = sprintf('%.1fK', num/1e3);
    else
        str = num2str(num);
    end
end

数字人可视化效果图如下：