（医学三维重建）MATLAB体绘制算法：最大密度投影（MIP）

by HPC_ZY

---

算法原理

1. 体绘制原理（大多数）
   用一矩形屏幕正对三维模型，从屏幕投射出M*N条平行光线，每条光线穿过模型。
   每条线会穿过很多像素块，就像这样：
   
   当每条线间隔大的时候，它们会穿过不同的像素块；而间隔很小（采样精度高）时，多条线也会穿过同一像素块，就像这样：
   
   最后我们单独处理每一条光线上的像素块，并将结果输出到屏幕上（光线发出的位置），就得到了体绘制结果。

2. MIP原理
   最大密度投影（Maximum Intensity Projection, MIP）：顾名思义，把最大的值投射出来。
   即，计算每一条光线穿过的像素块中灰度最大的，作为输出结果。
   举例：假如，
   1）我们有一个3x3x3的模型，数据如下图（右中）；
   2）屏幕大小（输出图像尺寸）为3x3，其坐标如下图（左上）；
   3）光线每前进1个单位，进行一次采样。
   
   我们从第一个位置发出一条射线，间隔取1，得到光线穿过的所有像素位置，如上图（右上）；
   然后可从模型对应位置，取出灰度值存入数组，如上图（右下）。由于模型没有第四层，所以（1,1,4）位置值为0。
   然后找出最大值，返回屏幕，如上图（左下）。

---

代码实现

注：为更好的展示算法原理，方便大家理解，以下代码未做任何优化，基本按照C语言的格式书写。

1. 数据准备
   关于model怎么生成，请查看三维体数据的生成

% 不得不说的事：
% 这个model是从医院拿的DICOM数据，然后事先读取进来的。
% 也不一定非得医学图像，你可以用任意你自己准备好的三维数据。
% 再简单一点讲，model就是一个M*M*D的数组，里面很很多数值，就能用来重建。
load model % 我的数据
model = im2double(model); % 模型归一化（这一步的目的是为了后面使用imshow，否则显示出来只有黑白色）
[mRows,mCols,mDims] = size(model); % 模型尺寸
mCenter = ([mRows,mCols,mDims]+1)/2; % 模型中心(因为下标从1开始)

值得注意的是，我这里使用的了model = im2double(model);, 目的是将数据归一化到0~1，后面使用imshow才能正常显示，否则显示出来只有黑白色

2. 参数设置

% 相机成像范围(默认等于尺寸中最大值)
cRows = max([mRows,mCols,mDims]);
cCols = max([mRows,mCols,mDims]);
cDims = max([mRows,mCols,mDims]);
% 相机成像分辨率（输出尺寸）
rRows =