目的

影像工具的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   1const 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 spacinglet d = 0; // default linear dimensionlet e = 0; // default extent startlet 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 boundsoutOrigin[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影像一般为1const 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.NEARESTmodel.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 1if (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 1const 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;}// 将影像坐标系矩阵复制给transformif (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// 否则返回falselet isIdentity = publicAPI.isIdentityMatrix(transform)；// the outMatrix takes OutputData indices to OutputData coordinates,// the inMatrix takes InputData coordinates to InputData indicesfor (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       1inMatrix[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   1inMatrix[4 * 3 + i] = -inOrigin[i] / inSpacing[i];//   osx       0        0      0//    0      osy       0      0//    0        0       osz     0//    0        0        0      1outMatrix[4 * i + i] = outSpacing[i];//   osx       0        0      0//    0    osy       0      0//    0        0       osz     0//   oox      ooy      ooz     1outMatrix[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 inputif (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 / iszmat4.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，第四行的第三项不为1const 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值，所以元组维度为1const inComponents = inScalars.getNumberOfComponents(); // interpolator.GetNumberOfComponents();// 获取组件偏移量，不设置为0const 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;// 不设置为falseconst rescaleScalars =model.scalarShift !== 0.0 || model.scalarScale !== 1.0;// is nearest neighbor optimization possible?// 是否近邻优化，一般为trueconst 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数据一般为1const 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、spacingconst 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 doublelet floatPtr = null;if (!optimizeNearest) { |floatPtr = new Float64Array(inComponents * (outExt[1] - outExt[0]));}// 初始化背景色，255位，全为0const background = macro.newTypedArray(inputScalarType,model.backgroundColor);// set color for area outside of input volume extent// void *background;// vtkAllocBackgroundPixel(&amp;background,//    self->GetBackgroundColor(), scalarType, scalarSize, outComponents);// get various helper functions// 一般为falseconst 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// 起始ylet idY = outExt[2] - 1;// 起始zlet 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*512const outPtr0 = iter.getScalars(output, 0);let outPtrIndex = 0;// outTmp的数据长度为输出数据的对角线长度乘2，数据类型为scalarType，dicom数据一般为int16const 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*512for (; !iter.isAtEnd(); iter.nextSpan()) {// 每一行x轴需要处理的数据数量const span = iter.spanEndId() - iter.getId();// 处理的起始点位置outPtrIndex = iter.getId() * scalarSize * outComponents;// 如果是模版内的数据则条件为false，可以看如果为trueif (!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为falseif (!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;}// 是否透视一般为falseif (perspective) {// only do perspective if necessaryconst f = 1 / inPoint[3]; |inPoint[0] *= f;inPoint[1] *= f;inPoint[2] *= f;}// 优化转换，一般为falseif (optimizedTransform !== null) {// apply the AbstractTransform if there is onepublicAPI.applyTransform(optimizedTransform,inPoint,inOrigin,inInvSpacing);}// checkBoundsIJK函数为校验点inPoint是否在范围内if (model.interpolator.checkBoundsIJK(inPoint)) {// do the interpolationisInBounds = 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，后续都为waswasInBounds = idX > idXmin ? wasInBounds : isInBounds;} // write a segment to the output// 当isInBounds和wasInBounds保持一致说明遍历完才退出，则idX多了1const endIdX = idX - 1 - (isInBounds !== wasInBounds);// 获取处理过数据的数量const numpixels = endIdX - startIdX + 1;let n = 0;// 如果先前在范围内，一般都在if (wasInBounds) {// 一般为falseif (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);}// 将值放到outPtr0for (let i = 0; i < n; ++i) {outPtr0[outPtrIndex++] = outTmp[i];}// startIdX加上处理数据数量startIdX += numpixels;// wasInBound置为isInBounds// wasInBounds和isInBounds应该为了避免边界情况，但是目前不清楚这样做wasInBounds = isInBounds;}} else {// optimize for nearest-neighbor interpolationconst inPtrTmp0 = inPtr;const outPtrTmp = outPtr;// inInc为每个轴的增长数量，以(512, 512, 1)，则是1，512，512*512const inIncX = inInc[0] * inputScalarSize;const inIncY = inInc[1] * inputScalarSize;const inIncZ = inInc[2] * inputScalarSize;// 获取dimensionconst 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 pixelsstartIdX = iidX;isInBounds = true;const n = setpixels(outTmp,background,outComponents,startIdX - idXmin);// 将outTmp赋值给outPtr0for (let i = 0; i < n; ++i) {outPtr0[outPtrIndex++] = outTmp[i];}}// set the final index that was within input boundsendIdX = 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 byteslet oc = 0;do {outPtr0[outPtrIndex++] = inPtrTmp0[offset++];} while (++oc !== bytesPerPixel);break;}} } else if (isInBounds) {// leaving input bounds// 之前在范围内，现在不在呢，则说明到达边界，则退出循环break;}}// clear trailing out-of-bounds pixelsoutPtr = outPtrTmp;// 遍历完毕，没有遍历到的证明不在范围内，则用背景色填充const n = setpixels(outTmp,background,outComponents,idXmax - endIdX);// 填充到outPtr0for (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。

有不懂或者有问题的可私聊或评论。

vtk的requestData 、getIndexMatrix和vtkImageResliceExecute讲解相关推荐

vtk相机_C#开发PACS医学影像三维重建(一)使用VTK重建3D影像
VTK简介: VTK是一个开源的免费软件系统,主要用于三维计算机图形学.图像处理和可视化.Vtk是在面向对象原理的基础上设计和实现的,它的内核是用C++构建的. 因为使用C#语言开发,而VTK是C++ ...
基础矩阵，本质矩阵，单应性矩阵讲解
ORB-SLAM点云地图中相机的位姿初始化,无论算法工作在平面场景,还是非平面场景下,都能够完成初始化的工作.其中主要是使用了适用于平面场景的单应性矩阵H和适用于非平面场景的基础矩阵F,程序中通过一个 ...
win下使用QT添加VTK插件实现点云可视化GUI
摘要大家在做点云的时候经常会用到QT,但是我们需要使用QT做点云的可视化的时候又需要VTK,虽然我们在windows下安装PCL的时候就已经安装了VTK,由于跟着PCL安装的VTK是没有和QT联合 ...
01-从零开始学习VTK
1.从零开始学习VTK 可能在这之前你没有使用过VTK,甚至不知道VTK是什么东西.这里假定你没有一点VTK基础,但已经有了一些基本的C/C++编程基础,以及计算机图形学的理论知识储备,想使用VTK从 ...
VTK序列图像的读取
VTK序列图像的读取分类: VTK系列教程 2013-05-07 19:50 2673人阅读评论(2) 收藏举报 VTK序列图像读取医学图像处理的应用程序中,经常会碰到读取一个序列图像的操作. ...
5、VTK在图像处理中的应用
5.VTK在图像处理中的应用图像是VTK中一个非常重要的数据.数字图像广泛应用于工业生产.生物医学.媒体娱乐.地质.气象等重要领域,数字图像处理具有重要的应用价值.我们在掌握了VTK的基本知识后,这 ...
VTK修炼之道10：可视化管道的连接与执行
1.可视化管道综述 vtkProp; vtkAbstractMapper; vtkProperty; vtkCamera; vtkLight; vtkRenderer; vtkRenderW ...
VTK修炼之道6_仔细分析一个复杂程序
1.程序代码 #include <vtkAutoInit.h> VTK_MODULE_INIT(vtkRenderingOpenGL); / #include <vtkSmartPo ...
creator qt 设置换行方式_win下使用QT添加VTK插件实现点云可视化GUI
大家在做点云的时候经常会用到QT,但是我们需要使用QT做点云的可视化的时候又需要VTK,虽然我们在windows下安装PCL的时候就已经安装了VTK,由于跟着PCL安装的VTK是没有和QT联合编译的, ...

vtk的requestData 、getIndexMatrix和vtkImageResliceExecute讲解

目的