目的
影像工具的mpr和mip函数的实现逻辑与上述两个函数相关，两个函数同时包含影像中最重要的坐标系和矩阵知识，记录以防以后忘记。同时为了彻底理清上述两个函数，以确保最后实现的功能没有任何问题。
先前知识
世界坐标系
世界坐标系是真实物理条件下的坐标系，单位坐标距离为1
影像坐标系
相对于世界坐标系而产生的坐标系，x、y、z轴的表示为与世界坐标系x、y、z轴之间的cos值，原点值为(0,0,0)，单位坐标距离由dicom影像的属性计算得到。
下面所说的input和output是两个影像坐标系
input的轴方向、原点、单位距离由dicom影像的属性得到。
output一般是我们转动和移动input的坐标系后得到的坐标系，所以轴方向是我们转动input轴后的方向，原点方向是我们移动input原点后的坐标，spacing是由我们根据转动方向和input的单位长度计算得到的。
requestData
requestData(inData, outData) {
// implement requestData
// 获取输入的影像数据，这是一个对象，包含对影像数据各个参数获取和设置的方法
const input = inData[0];
if (!input) { |
vtkErrorMacro('Invalid or missing input'); |
return;
}
// Retrieve output and volume data
// 获取影像的原点位置，一般为(0, 0, 0)
const origin = input.getOrigin();
// 获取影像数据每个像素间距，一个像素代表一个坐标，像素间距代表单位坐标间距
const inSpacing = input.getSpacing();
// 获取输入影像数据的维度，如x轴有多少个像素，y轴有多少个像素，z轴有多少个像素
const dims = input.getDimensions();
// 影像数据特有属性
const inScalars = input.getPointData().getScalars();
// 获取像素数范围，代表影像数据整个坐标系大小
const inWholeExt = [0, dims[0] - 1, 0, dims[1] - 1, 0, dims[2] - 1]; |
// 初始化输出数据的属性
const outOrigin = [0, 0, 0];
const outSpacing = [1, 1, 1];
const outWholeExt = [0, 0, 0, 0, 0, 0];
const outDims = [0, 0, 0];
let matrix = null;
// resliceAxes就是我们输入的四位矩阵，也可以说切割坐标系，按照这个坐标系对影像数据进行切割
// 如果不存在resliceAxes就会初始化一个标准坐标系(和世界坐标系一致)
// cosx1 cosx2 cosx3 0 切割影像坐标系的x轴与世界坐标系x、y、z轴的cos值
// cosy1 cosy2 cosy3 0 切割影像坐标系的y轴与世界坐标系x、y、z轴的cos值
// cosz1 cosz2 cosz3 0 切割影像坐标系的z轴与世界坐标系x、y、z轴的cos值
// x y z 1 切割影像坐标系的原点值(为在世界坐标系下的坐标)
if (model.resliceAxes) {
matrix = model.resliceAxes;
} else {
matrix = mat4.identity(new Float64Array(16));
}
// 将矩阵反转，变成下述这样
// cosx1 cosy1 cosz1 x
// cosx2 cosy2 cosz2 y
// cosx3 cosy3 cosz3 z
// 0 0 0 1
const imatrix = new Float64Array(16);
mat4.invert(imatrix, matrix);
// 获取输入影像数据的中心点在世界坐标系下的坐标
const inCenter = [
origin[0] + 0.5 * (inWholeExt[0] + inWholeExt[1]) * inSpacing[0],
origin[1] + 0.5 * (inWholeExt[2] + inWholeExt[3]) * inSpacing[1],
origin[2] + 0.5 * (inWholeExt[4] + inWholeExt[5]) * inSpacing[2],
];
let maxBounds = null;
// 是否采用vtk内部方法对输入数据进行裁剪获取数据范围
// 算法逻辑为
// 1.遍历输入数据的8个点坐标(三维数据的8个顶点)
// 2.获取8个顶点数据在世界坐标系下的坐标
// 3.世界坐标系下的坐标映射到resliceAxes坐标系下的坐标
// 4.比较输入数据的bounds和resliceAxes坐标系下的维度
if (model.autoCropOutput) {
maxBounds = publicAPI.getAutoCroppedOutputBounds(input);
}
// 计算输出数据的spacing、dimension、center点、范围
// 循环每一层依次处理影像坐标系的x、y、z轴
for (let i = 0; i < 3; i++) {
let s = 0; // default output spacing
let d = 0; // default linear dimension
let e = 0; // default extent start
let c = 0; // transformed center-of-volume
// 如果将transformInputSampling设置为true，采用vtk的算法计算输出数据的各个值
// 否则就采用输入数据的属性
if (model.transformInputSampling) {
let r = 0.0;
// 遍历影像坐标系每个轴与世界坐标系x、y、z轴之间的关系
for (let j = 0; j < 3; j++) {
// 计算输入影像数据的中心点在输出影像坐标系下的坐标
c += imatrix[4 * j + i] * (inCenter[j] - matrix[4 * 3 + j]);
// 利用输入数据的信息计算输出数据的信息，通过与世界坐标系轴夹脚cos值计算
const tmp = matrix[4 * i + j] * matrix[4 * i + j];
s += tmp * Math.abs(inSpacing[j]);
d += tmp * (inWholeExt[2 * j + 1] - inWholeExt[2 * j]) * Math.abs(inSpacing[j]);
e += tmp * inWholeExt[2 * j];
r += tmp;
}
s /= r;
d /= r * Math.sqrt(r);
e /= r;
} else {
c = inCenter[i];
s = inSpacing[i];
d = (inWholeExt[2 * i + 1] - inWholeExt[2 * i]) * s;
e = inWholeExt[2 * i];
}
// 通过方法对输出数据的属性赋值优先级最高，优先采用该值
// 否则就采用计算得到的值
if (model.outputSpacing == null) {
outSpacing[i] = s;
} else {
outSpacing[i] = model.outputSpacing[i];
}
// 如果我们只需要2维数据，超过2维的都会置为0，例如我们只需要二维数据，其outWholeExt为[xmin,xmax,ymin,ymax,0,0]
// 如果我们输入了维数，outWholeExt会依据我们输入的维数而进行设置
// 如果我们没有输入维数，会依据我们是否采用了自动裁剪功能，如果自动裁剪了会使用算法裁剪输出数据维度
if (i >= model.outputDimensionality) {
outWholeExt[2 * i] = 0;
outWholeExt[2 * i + 1] = 0;
} else if (model.outputExtent == null) {
// 自动裁剪主要是用到了内部边界检测算法算出边界坐标
if (model.autoCropOutput) {
d = maxBounds[2 * i + 1] - maxBounds[2 * i];
}
outWholeExt[2 * i] = Math.round(e);
outWholeExt[2 * i + 1] = Math.round(outWholeExt[2 * i] + Math.abs(d / outSpacing[i]));
} else {
outWholeExt[2 * i] = model.outputExtent[2 * i];
outWholeExt[2 * i + 1] = model.outputExtent[2 * i + 1];
}
// 获取输出数据的origin位置，为输出数据的原点坐标，一般位于视图层的左下角，依据情况决定
// 和维数处理相同
// 需要注意这里的origin和矩阵中的第四行原点值不同，矩阵原点就是影像坐标系的原点位置，影像坐标系的原点位置
// 不一定在影像数据的左下角，可能在影像数据中的任意位置，origin就是影像坐标系的左下角位置，origin坐标是参
// 照影像坐标系原点的坐标，不一定是(0,0,0)
if (i >= model.outputDimensionality) {
outOrigin[i] = 0;
} else if (model.outputOrigin == null) {
if (model.autoCropOutput) {
// set origin so edge of extent is edge of bounds
outOrigin[i] = maxBounds[2 * i] - outWholeExt[2 * i] * outSpacing[i];
} else {
// center new bounds over center of input bounds
// 这里需要注意的一点是乘了一个outSpacing，这里不是影像坐标(有真实距离)，也不是世界坐标(没有进行坐标转换)
// 这里可以称作是带有真实距离的影像坐标，如果不是用自动裁剪，影像中心点为体影像在影像坐标系中的映射点
// 原点位置就是影像中心点减去影像长度的一半
outOrigin[i] = c - 0.5 * (outWholeExt[2 * i] + outWholeExt[2 * i + 1]) * outSpacing[i];
}
} else {
outOrigin[i] = model.outputOrigin[i];
}
// 记录影像维数
outDims[i] = outWholeExt[2 * i + 1] - outWholeExt[2 * i] + 1;
}
// 获取输出数据的数据类型
let dataType = inScalars.getDataType();
if (model.outputScalarType) {
dataType = model.outputScalarType;
}
// 每个像素右几个值确定，例如rgb为3个，rgba为4个，dicom影像一般为1
const numComponents = input
.getPointData()
.getScalars()
.getNumberOfComponents(); // or s.numberOfComponents;
// 创建输出数据
const outScalarsData = macro.newTypedArray(
dataType,
outDims[0] * outDims[1] * outDims[2] * numComponents
);
// 构建vtk的数据结构
const outScalars = vtkDataArray.newInstance({
name: 'Scalars',
values: outScalarsData,
numberOfComponents: numComponents,
});
// Update output
// 创建vtk数据对象
const output = vtkImageData.newInstance();
// 设置数据维数
output.setDimensions(outDims);
// 设置数据原点值
output.setOrigin(outOrigin);
// 设置输出数据单位坐标的距离
output.setSpacing(outSpacing);
// 设置数据结构类型
output.getPointData().setScalars(outScalars);
// 获取input和output的转换函数，这里最好自己手动推一下
publicAPI.getIndexMatrix(input, output);
// 获取插值模式
let interpolationMode = model.interpolationMode;
model.usePermuteExecute = false;
if (model.optimization) {
if (
optimizedTransform == null &&
model.slabSliceSpacingFraction === 1.0 &&
model.interpolator.isSeparable() &&
publicAPI.isPermutationMatrix(indexMatrix)
) {
model.usePermuteExecute = true;
if (publicAPI.canUseNearestNeighbor(indexMatrix, outWholeExt)) {
interpolationMode = InterpolationMode.NEAREST;
}
}
}
// 设置插值模式，插值在图像处理很多，因为例如数据只有26*512，用户设置的输出维度必须是512*512
// 这个时候就可以通过插值将26*512的数据变成512*512，默认为InterpolationMode.NEAREST
model.interpolator.setInterpolationMode(interpolationMode);
// 设置坐标超出border范围时如何处理
let borderMode = ImageBorderMode.CLAMP;
borderMode = model.wrap ? ImageBorderMode.REPEAT : borderMode;
borderMode = model.mirror ? ImageBorderMode.MIRROR : borderMode;
model.interpolator.setBorderMode(borderMode);
const mintol = 7.62939453125e-6;
const maxtol = 2.0 * 2147483647;
let tol = 0.5 * model.border;
tol = borderMode === ImageBorderMode.CLAMP ? tol : maxtol;
tol = tol > mintol ? tol : mintol;
// 设置容差
model.interpolator.setTolerance(tol);
// 使用input对model的属性scalar、spacing、origin、extent初始化
model.interpolator.initialize(input);
// 输出数据的属性得到后，该函数就是计算其真正的值
publicAPI.vtkImageResliceExecute(input, output);
model.interpolator.releaseData();
outData[0] = output;
}
getIndexMatrix(input, output)
获取input和output之间的转换矩阵，矩阵前三行是output影像坐标系的单位坐标距离转化为input影像坐标系的单位坐标距离，第四行是output的原点坐标在input影像坐标系下的坐标。
getIndexMatrix(input, output) {
// first verify that we have to update the matrix
// 初始化indexMatrix
// 1 0 0 0
// 0 1 0 0
// 0 0 1 0
// 0 0 0 1
if (indexMatrix === null) {
indexMatrix = mat4.identity(new Float64Array(16));
}
// 获取输入和输出的各个属性
// 我们假设值方便大家理解
// (iox, ioy, ioz)
const inOrigin = input.getOrigin();
// (isx, isy, isz)
const inSpacing = input.getSpacing();
// (oox, ooy, ooz)
const outOrigin = output.getOrigin();
// (osx, osy, osz)
const outSpacing = output.getSpacing();
// 初始化
// 1 0 0 0
// 0 1 0 0
// 0 0 1 0
// 0 0 0 1
const transform = mat4.identity(new Float64Array(16));
const inMatrix = mat4.identity(new Float64Array(16));
const outMatrix = mat4.identity(new Float64Array(16));
// 如果optimizedTransform存在则重置
if (optimizedTransform) {
optimizedTransform = null;
}
// 将影像坐标系矩阵复制给transform
if (model.resliceAxes) {
mat4.copy(transform, model.resliceAxes);
}
if (model.resliceTransform) {
// TODO
}
// check to see if we have an identity matrix
// 如果是
// 1 0 0 0
// 0 1 0 0
// 0 0 1 0
// 0 0 0 1
// 返回true
// 否则返回false
let isIdentity = publicAPI.isIdentityMatrix(transform)；
// the outMatrix takes OutputData indices to OutputData coordinates,
// the inMatrix takes InputData coordinates to InputData indices
for (let i = 0; i < 3; i++) {
if (
(optimizedTransform === null &&
(inSpacing[i] !== outSpacing[i] || inOrigin[i] !== outOrigin[i])) ||
(optimizedTransform !== null &&
(outSpacing[i] !== 1.0 || outOrigin[i] !== 0.0))
) {
isIdentity = false;
}
// 1/isx 0 0 0
// 0 1/isy 0 0
// 0 0 1/isz 0
// 0 0 0 1
inMatrix[4 * i + i] = 1.0 / inSpacing[i];
// 1/isx 0 0 0
// 0 1/isy 0 0
// 0 0 1/isz 0
// -iox/isx -ioy/isy -ioz/isz 1
inMatrix[4 * 3 + i] = -inOrigin[i] / inSpacing[i];
// osx 0 0 0
// 0 osy 0 0
// 0 0 osz 0
// 0 0 0 1
outMatrix[4 * i + i] = outSpacing[i];
// osx 0 0 0
// 0 osy 0 0
// 0 0 osz 0
// oox ooy ooz 1
outMatrix[4 * 3 + i] = outOrigin[i];
}
// 进入执行
if (!isIdentity) {
// transform.PreMultiply();
// transform.Concatenate(outMatrix);
// outMatrix ✖️ transform
// cosx1 * osx cosx2 * osx cosx3 * osx 0
// cosy1 * osy cosy2 * osy cosy3 * osy 0
// cosz1 * osz cosz2 * osz cosz3 * osz 0
// x y z 1
// output的origin的坐标是自带spacing的
// 第四行为origin在世界坐标系下的坐标
// x = x + oox * cosx1 + ooy * cosy1 + ooz * cosz1
// y = y + oox * cosx2 + ooy * cosy2 + ooz * cosz2
// z = z + oox * cosx3 + ooy * cosy3 + ooz * cosz3
// 上述目的是获取坐标系以output的origin为原点的转换矩阵
// vtk内部处理是以origin为初始点进行开始遍历的
mat4.multiply(transform, transform , outMatrix);
// the optimizedTransform requires data coords, not
// index coords, as its input
if (optimizedTransform == null) {
// transform->PostMultiply();
// transform->Concatenate(inMatrix);
// tramsform * inMatrix
// cosx1 * osx / isx cosx2 * osx / isy cosx3 * osx / isz 0
// cosy1 * osy / isx cosy2 * osy / isy cosy3 * osy / isz 0
// cosz1 * osz / isx cosz2 * osz / isy cosz3 * osz / isz 0
// x y z 1
// 第四行为origin在世界坐标系下的坐标
// 假定output坐标系下以origin为原点的点(x1,y1,z1,1)
// 世界坐标系下的x为(x1 * cosx1 * osx + y1 * cosy1 * osy + z1 * cosz1 * osz + x)
// 转化为input坐标系下以origin为原点的x为
// ((x1 * cosx1 * osx + y1 * cosy1 * osy + z1 * cosz1 * osz + x) - iox) / isx
// 看到这个x表达式其实可以看到transform的影子
// x = x / isx - iox / isx
// y = y / isy - ioy / isy
// z = z / isz - ioz / isz
mat4.multiply(transform, inMatrix, transform);
}
}
mat4.copy(indexMatrix, transform);
// indexMartrix就算处理完毕了，大家看到这里可能一脸懵逼，这么做为了啥？
// 如果把整个切割多看几遍应该会有眉目，切割是需要获取output数据，output数据哪里来？当然是从input中来了
// 因为要获取output数据，肯定遍历output数据中每一个位置，然后计算该位置对应到input数据中哪个位置，
// 然后将该位置的input数据赋值给output，这个时候懂了没？还没懂？
// 这个indexMatrix的作用就是将output坐标系中的坐标映射到input坐标系中，利用indexMatrix矩阵方便映射
return indexMatrix;
}
vtkImageResliceExecute(input, output)
依据input的spacing、原点坐标、x、y、z轴与世界坐标系的x、y、z轴的cos值和output的spacing、原点坐标、x、y、z轴与世界坐标系的x、y、z轴的cos值获取output的切面数据。
vtkImageResliceExecute(input, output) {
// const outDims = output.getDimensions();
// 获取数据结构
const inScalars = input.getPointData().getScalars();
const outScalars = output.getPointData().getScalars();
// 获取输出数据
let outPtr = outScalars.getData();
// 获取输出维度范围
const outExt = output.getExtent();
// 获取转化
const newmat = indexMatrix;
const outputStencil = null;
// multiple samples for thick slabs
// 获取mip的层数
const nsamples = Math.max(model.slabNumberOfSlices, 1);
// spacing between slab samples (as a fraction of slice spacing)
// 获取mip的层厚
const slabSampleSpacing = model.slabSliceSpacingFraction;
// check for perspective transformation
// 主要判断前三行的第三项不为0，第四行的第三项不为1
const perspective = publicAPI.isPerspectiveMatrix(newmat);
// extra scalar info for nearest-neighbor optimization
// 获取输入数据
let inPtr = inScalars.getData();
const inputScalarSize = 1;
// inScalars.getElementComponentSize();
// inScalars.getDataTypeSize();
// 获取输入数据类型
const inputScalarType = inScalars.getDataType();
// 获取每像素组件维度，元组是一个像素的表示维度
// dicom数据一般为ct值，所以元组维度为1
const inComponents = inScalars.getNumberOfComponents();
// interpolator.GetNumberOfComponents();
// 获取组件偏移量，不设置为0
const componentOffset = model.interpolator.getComponentOffset();
// 获取对于出界位置的处理方式
const borderMode = model.interpolator.getBorderMode();
// 获取输入数据维度
const inDims = input.getDimensions();
// 获取输入数据维度范围
const inExt = [0, inDims[0] - 1, 0, inDims[1] - 1, 0, inDims[2] - 1];
// interpolator->GetExtent();
// 输入数据每个维度数据数量
const inInc = [0, 0, 0];
// x轴每项数据数量
inInc[0] = inScalars.getNumberOfComponents();
// y轴每行数据数量---x轴
inInc[1] = inInc[0] * inDims[0];
// z轴每层数据数量---x,y平面
inInc[2] = inInc[1] * inDims[1];
// 输入数据总维度
const fullSize = inDims[0] * inDims[1] * inDims[2];
// 没使用到
if (componentOffset > 0 && componentOffset + inComponents < inInc[0]) {
inPtr = inPtr.subarray(inputScalarSize * componentOffset);
}
// 同样获取插值方式
let interpolationMode = InterpolationMode.NEAREST;
if (model.interpolator.isA('vtkImageInterpolator')) {
interpolationMode = model.interpolator.getInterpolationMode();
}
const convertScalars = null;
// 不设置为false
const rescaleScalars =
model.scalarShift !== 0.0 || model.scalarScale !== 1.0;
// is nearest neighbor optimization possible?
// 是否近邻优化，一般为true
const optimizeNearest =
interpolationMode === InterpolationMode.NEAREST &&
borderMode === ImageBorderMode.CLAMP &&
!(
optimizedTransform != null ||
perspective ||
convertScalars != null ||
rescaleScalars
) &&
inputScalarType === outScalars.getDataType() &&
fullSize === inScalars.getNumberOfTuples() &&
model.border === true &&
nsamples <= 1;
// get pixel information
// 获取输出数据的类型
const scalarType = outScalars.getDataType();
const scalarSize = 1;
// outScalars.getElementComponentSize()
// outScalars.scalarSize;
// 获取单项数据数量，dicom数据一般为1
const outComponents = outScalars.getNumberOfComponents();
// break matrix into a set of axes plus an origin
// (this allows us to calculate the transform Incrementally)
// 获取output坐标系下x、y、z轴坐标对input坐标系下的x、y、z轴的转化
// origin为input坐标系下的遍历起始点
const xAxis = [0, 0, 0, 0];
const yAxis = [0, 0, 0, 0];
const zAxis = [0, 0, 0, 0];
const origin = [0, 0, 0, 0];
for (let i = 0; i < 4; ++i) {
xAxis[i] = newmat[4 * 0 + i];
yAxis[i] = newmat[4 * 1 + i];
zAxis[i] = newmat[4 * 2 + i];
origin[i] = newmat[4 * 3 + i];
}
// get the input origin and spacing for conversion purposes
// 获取输入数据的origin、spacing
const inOrigin = model.interpolator.getOrigin();
const inSpacing = model.interpolator.getSpacing();
const inInvSpacing = [
1.0 / inSpacing[0],
1.0 / inSpacing[1],
1.0 / inSpacing[2],
];
// allocate an output row of type double
let floatPtr = null;
if (!optimizeNearest) { |
floatPtr = new Float64Array(inComponents * (outExt[1] - outExt[0]));
}
// 初始化背景色，255位，全为0
const background = macro.newTypedArray(
inputScalarType,
model.backgroundColor
);
// set color for area outside of input volume extent
// void *background;
// vtkAllocBackgroundPixel(&background,
// self->GetBackgroundColor(), scalarType, scalarSize, outComponents);
// get various helper functions
// 一般为false
const forceClamping =
interpolationMode > InterpolationMode.LINEAR ||
(nsamples > 1 && model.slabMode === SlabMode.SUM);
// 获取转化像素方法
const convertpixels = publicAPI.getConversionFunc(
inputScalarType,
scalarType,
model.scalarShift,
model.scalarScale,
forceClamping
);
// 获取设置像素的函数
// outComponents为1采用fill，如果不为1则使用for循环赋值
const setpixels = publicAPI.getSetPixelsFunc(
scalarType,
scalarSize,
outComponents,
outPtr
);
// 获取mip的形式，max最大，min最小，mean平均的函数
const composite = publicAPI.getCompositeFunc(
model.slabMode,
model.slabTrapezoidIntegration
);
// create some variables for when we march through the data
// 起始y
let idY = outExt[2] - 1;
// 起始z
let idZ = outExt[4] - 1;
// 内部使用的中间变量
const inPoint0 = [0.0, 0.0, 0.0, 0.0];
const inPoint1 = [0.0, 0.0, 0.0, 0.0];
// create an iterator to march through the data
// 获取数据迭代器
const iter = vtkImagePointDataIterator.newInstance();
// 利用output和outExt来初始化迭代器
iter.initialize(output, outExt, model.stencil, null);
// 获取输出数据，如果输出为(512, 512, 0)，则数据长度为512*512
const outPtr0 = iter.getScalars(output, 0);
let outPtrIndex = 0;
// outTmp的数据长度为输出数据的对角线长度乘2，数据类型为scalarType，dicom数据一般为int16
const outTmp = macro.newTypedArray(
scalarType,
vtkBoundingBox.getDiagonalLength(outExt) * outComponents * 2
);
// 初始化数据
const interpolatedPtr = new Float64Array(inComponents * nsamples);
const interpolatedPoint = new Float64Array(inComponents);
// iter.isAtEnd()是否遍历完所有输出数据
// 如果为mpr，假设输出数据的维度为(512, 512, 0)，第三项为0的原因为mpr为二维数据
// iter.isAtEnd()里面会有一个id为正在处理数据的位置，初始位置一般为0，判断id是否为end，end一般为512*512
for (; !iter.isAtEnd(); iter.nextSpan()) {
// 每一行x轴需要处理的数据数量
const span = iter.spanEndId() - iter.getId();
// 处理的起始点位置
outPtrIndex = iter.getId() * scalarSize * outComponents;
// 如果是模版内的数据则条件为false，可以看如果为true
if (!iter.isInStencil()) {
// clear any regions that are outside the stencil
// 模版外的数据将背景色进行填充
const n = setpixels(outTmp, background, outComponents, span);
for (let i = 0; i < n; ++i) {
outPtr0[outPtrIndex++] = outTmp[i];
}
} else {
// get output index, and compute position in input image
// 获取处理的数据的位置
const outIndex = iter.getIndex();
// if Z index increased, then advance position along Z axis
// outIndex如果是增加则代表在往z轴方向进行遍历
if (outIndex[2] > idZ) {
idZ = outIndex[2];
inPoint0[0] = origin[0] + idZ * zAxis[0];
inPoint0[1] = origin[1] + idZ * zAxis[1];
inPoint0[2] = origin[2] + idZ * zAxis[2];
inPoint0[3] = origin[3] + idZ * zAxis[3];
idY = outExt[2] - 1;
}
// if Y index increased, then advance position along Y axis
// outIndex如果是增长则代表往y轴方向进行遍历
if (outIndex[1] > idY) {
idY = outIndex[1];
inPoint1[0] = inPoint0[0] + idY * yAxis[0];
inPoint1[1] = inPoint0[1] + idY * yAxis[1];
inPoint1[2] = inPoint0[2] + idY * yAxis[2];
inPoint1[3] = inPoint0[3] + idY * yAxis[3];
}
// 上述是在output层进行从y轴和z后方向进行遍历
// march through one row of the output image
// idXmin为x的最小点，idMax为x的最大点
const idXmin = outIndex[0];
const idXmax = idXmin + span - 1;
// 如果是mip设置了层数，则optimizeNearest为false
if (!optimizeNearest) {
// 上一次是否在数据范围内
let wasInBounds = 1;
// 当前是否在数据范围内
let isInBounds = 1;
// 遍历起始点
let startIdX = idXmin;
// 遍历当前点
let idX = idXmin;
const tmpPtr = floatPtr;
// 影像数据下表
let pixelIndex = 0;
// 遍历x轴
while (startIdX <= idXmax) {
for (; idX <= idXmax && isInBounds === wasInBounds; idX++) {
// 对x轴进行遍历
const inPoint2 = [
inPoint1[0] + idX * xAxis[0],
inPoint1[1] + idX * xAxis[1],
inPoint1[2] + idX * xAxis[2],
inPoint1[3] + idX * xAxis[3],
];
// 中间变量
const inPoint3 = [0, 0, 0, 0];
let inPoint = inPoint2;
isInBounds = false;
let interpolatedPtrIndex = 0;
// 对mip的z轴进行遍历
for (let sample = 0; sample < nsamples; ++sample) {
if (nsamples > 1) {
let s = sample - 0.5 * (nsamples - 1);
s *= slabSampleSpacing;
// 对z轴移动
inPoint3[0] = inPoint2[0] + s * zAxis[0];
inPoint3[1] = inPoint2[1] + s * zAxis[1];
inPoint3[2] = inPoint2[2] + s * zAxis[2];
inPoint3[3] = inPoint2[3] + s * zAxis[3];
inPoint = inPoint3;
}
// 是否透视一般为false
if (perspective) {
// only do perspective if necessary
const f = 1 / inPoint[3]; |
inPoint[0] *= f;
inPoint[1] *= f;
inPoint[2] *= f;
}
// 优化转换，一般为false
if (optimizedTransform !== null) {
// apply the AbstractTransform if there is one
publicAPI.applyTransform(
optimizedTransform,
inPoint,
inOrigin,
inInvSpacing
);
}
// checkBoundsIJK函数为校验点inPoint是否在范围内
if (model.interpolator.checkBoundsIJK(inPoint)) {
// do the interpolation
isInBounds = 1;
// 获取input数据中inPoint位置的像素值存到interpolatedPoint中
model.interpolator.interpolateIJK(inPoint, interpolatedPoint);
// 将interpolatedPoint的像素存放到interpolatedPtr中
for (let i = 0; i < inComponents; ++i) {
interpolatedPtr[interpolatedPtrIndex++] = interpolatedPoint[i];
}
}
}
if (interpolatedPtrIndex > inComponents) {
// 采用对应的mip策略处理interpolatedPtr，值存放到interpolatedPtr[0]
composite(
interpolatedPtr,
inComponents,
interpolatedPtrIndex / inComponents
);
}
// 最后将值存放到tmpPtr中，注意tmpPtr和floatPtr引用同一份数据
for (let i = 0; i < inComponents; ++i) {
tmpPtr[pixelIndex++] = interpolatedPtr[i];
}
// set "was in" to "is in" if first pixel
// 从逻辑上看，第一次处理为is，后续都为was
wasInBounds = idX > idXmin ? wasInBounds : isInBounds;
}
// write a segment to the output
// 当isInBounds和wasInBounds保持一致说明遍历完才退出，则idX多了1
const endIdX = idX - 1 - (isInBounds !== wasInBounds);
// 获取处理过数据的数量
const numpixels = endIdX - startIdX + 1;
let n = 0;
// 如果先前在范围内，一般都在
if (wasInBounds) {
// 一般为false
if (outputStencil) {
outputStencil.insertNextExtent(startIdX, endIdX, idY, idZ);
}
// 是否存在缩放的像素
if (rescaleScalars) {
publicAPI.rescaleScalars(
floatPtr,
inComponents,
idXmax - idXmin + 1,
model.scalarShift,
model.scalarScale
);
}
// 如果有转换进行转换，将floatPtr值放入到outTmp中
if (convertScalars) {
convertScalars(
floatPtr.subarray(startIdX * inComponents),
outTmp,
inputScalarType,
inComponents,
numpixels,
startIdX,
idY,
idZ
);
n = numpixels * outComponents * scalarSize;
} else {
n = convertpixels(
outTmp,
floatPtr.subarray(startIdX * inComponents),
outComponents,
numpixels
);
}
} else {
n = setpixels(outTmp, background, outComponents, numpixels);
}
// 将值放到outPtr0
for (let i = 0; i < n; ++i) {
outPtr0[outPtrIndex++] = outTmp[i];
}
// startIdX加上处理数据数量
startIdX += numpixels;
// wasInBound置为isInBounds
// wasInBounds和isInBounds应该为了避免边界情况，但是目前不清楚这样做
wasInBounds = isInBounds;
}
} else {
// optimize for nearest-neighbor interpolation
const inPtrTmp0 = inPtr;
const outPtrTmp = outPtr;
// inInc为每个轴的增长数量，以(512, 512, 1)，则是1，512，512*512
const inIncX = inInc[0] * inputScalarSize;
const inIncY = inInc[1] * inputScalarSize;
const inIncZ = inInc[2] * inputScalarSize;
// 获取dimension
const inExtX = inExt[1] - inExt[0] + 1;
const inExtY = inExt[3] - inExt[2] + 1;
const inExtZ = inExt[5] - inExt[4] + 1;
// 遍历起始点
let startIdX = idXmin;
let endIdX = idXmin - 1;
let isInBounds = false;
// 单位像素的偏移
const bytesPerPixel = inputScalarSize * inComponents;
// 开始遍历，从x轴开始
for (let iidX = idXmin; iidX <= idXmax; iidX++) {
// 获取遍历点对应于input坐标系下的位置
const inPoint = [
inPoint1[0] + iidX * xAxis[0],
inPoint1[1] + iidX * xAxis[1],
inPoint1[2] + iidX * xAxis[2],
];
// 四舍五入获取对应的位置
const inIdX = vtkInterpolationMathRound(inPoint[0]) - inExt[0];
const inIdY = vtkInterpolationMathRound(inPoint[1]) - inExt[2];
const inIdZ = vtkInterpolationMathRound(inPoint[2]) - inExt[4];
// 如果在范围内则进行赋值
if (
inIdX >= 0 &&
inIdX < inExtX &&
inIdY >= 0 &&
inIdY < inExtY &&
inIdZ >= 0 &&
inIdZ < inExtZ
) {
// 如果之前不在范围内，现在在，则将之前的所有像素都置为背景色
if (!isInBounds) {
// clear leading out-of-bounds pixels
startIdX = iidX;
isInBounds = true;
const n = setpixels(
outTmp,
background,
outComponents,
startIdX - idXmin
);
// 将outTmp赋值给outPtr0
for (let i = 0; i < n; ++i) {
outPtr0[outPtrIndex++] = outTmp[i];
}
}
// set the final index that was within input bounds
endIdX = iidX;
// perform nearest-neighbor interpolation via pixel copy
// inPtrTmp0是一个一维数组，这个计算其实相当于三位数组偏移量变成一维数组偏移量
let offset = inIdX * inIncX + inIdY * inIncY + inIdZ * inIncZ;
// when memcpy is used with a constant size, the compiler will
// optimize away the function call and use the minimum number
// of instructions necessary to perform the copy
// 根据每个像素的长度使用不同方法赋值，可能为1个像素值或rgb值
switch (bytesPerPixel) {
case 1:
outPtr0[outPtrIndex++] = inPtrTmp0[offset];
break;
case 2:
case 3:
case 4:
case 8:
case 12:
case 16:
for (let i = 0; i < bytesPerPixel; ++i) {
outPtr0[outPtrIndex++] = inPtrTmp0[offset + i];
}
break;
default: {
// TODO: check bytes
let oc = 0;
do {
outPtr0[outPtrIndex++] = inPtrTmp0[offset++];
} while (++oc !== bytesPerPixel);
break;
}
}
} else if (isInBounds) {
// leaving input bounds
// 之前在范围内，现在不在呢，则说明到达边界，则退出循环
break;
}
}
// clear trailing out-of-bounds pixels
outPtr = outPtrTmp;
// 遍历完毕，没有遍历到的证明不在范围内，则用背景色填充
const n = setpixels(
outTmp,
background,
outComponents,
idXmax - endIdX
);
// 填充到outPtr0
for (let i = 0; i < n; ++i) {
outPtr0[outPtrIndex++] = outTmp[i];
}
if (outputStencil && endIdX >= startIdX) {
outputStencil.insertNextExtent(startIdX, endIdX, idY, idZ);
}
}
}
}
}
总结
上述没有使用到插值，看了下最近提交记录好像有了线性插值，代码还没细看，看了之后补充下。
我基于上述代码进行了二次开发，里面有些功能没有用到我自己删除了，只列举两个优化点：
- 上述代码每次循环会用一个内部变量获取这次循环获取到的像素值再作用到output，可以删除直接作用到out。
- mip的nsamples是确定的，s * zAxis[0]这个值是可以缓存而避免每次循环重新计算
我自己内部实现了一个三线性插值的方法，理论挺简单的，大家可以根据情况自己实现。
后面有时间会补充webgl和webgpu版本的mip。
有不懂或者有问题的可私聊或评论。
3493
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
