需求
仅实现肋骨、脊柱上的效果，血管太复杂了一般要自动计算勾勒出沿着血管的路径才行，不然在大几百张图层上选择路径不是用户期望的。
项目背景
react + cornerstone的一个web影像渲染引擎
实现大体是以下几个步骤：
1.通过vtk实现mpr重建图像，操作面为横断位
2.自定义cornerstone-tools工具CPRTool，实现路径点/线的绘制
3.根据路径点将每一段路径截取对应的切面的pixeldata
4.拼接所有切面
5.将拼接后的pixel数据构造出一个CornerStoneImage，供cornerstone渲染canvas
step1
通过vtkjs获得多平面重建图像这个不赘述了，官网有相关的example[https://kitware.github.io/vtk-js/examples/ResliceCursorWidget.html]
首先，操作面是重建后的横断面，我们可以想象一下，相当于一把刀从横断的第一层切到最后一层（如图），所以曲面重建后图像的高度是固定了的
视角（原谅我不会作图）
step2
不用cornerstone及其tools这套框架的跳过这段就成，这边的目的是获得canvas上绘制的路径点
自定以路径绘制的cornerstone-tools工具CRPTools，可以参考官网提供的Probe工具[https://tools.cornerstonejs.org/examples/tools/probe.html] Probe工具落下的点其实就类比我们的路径点，只要使用drawLines多绘制一下中间的线段就行。
这边有一个问题，用户期望的是所有层上都可以看到完整路径，切换层只是为了方便找到准确位置选点，但是cornerstone-tools这儿toolState的机制默认是以imageId为key的，表现就是这些点只能在绘制的那层看到。解决办法就是自定义toolStateManager，cornerstone-tools提供setElementToolStateManager方法，去源码里把imageId的manager复制过来重写一个seriesInstanceUidStateManager.js，add和get方法都需要修改。
function addElementToolState(element, toolName, data) {
const enabledElement = external.cornerstone.getEnabledElement(element);
// 原先这边是imageId，现在改成和series关联的id
//这样一整个序列每张图的这个id都是一致的，路径就可以在所有image上完整显示了
if (!enabledElement.uuid) {
return;
}
addImageIdToolState(enabledElement.uuid, toolName, data);
}
得到每个点的image坐标之后，需要转成三维坐标才能给vtk使用，cornerstone-tools提供方法：imagePointToPatientPoint。
const v3 = imagePointToPatientPoint(eventData.currentPoints.image, imagePlane);
step3
接下来就是通过路径点来抽取每一spacing的像素值拼接起来。首先将得到的空间坐标转换成vtk world坐标系：
// Image Volume图像坐标系转换到VTK Volume 图像坐标系
projectToVtkVolumePosition(imagePosition) {
const pos = [imagePosition[0] - this.firstImageMatrix[0][3], imagePosition[1] - this.firstImageMatrix[1][3], imagePosition[2] - this.firstImageMatrix[2][3]];
let result = [0, 0, 0];
for (let i = 0; i < 3; i++) {
let vectorA = i === 2 ? this.zVector : [this.firstImageMatrix[0][i], this.firstImageMatrix[1][i], this.firstImageMatrix[2][i]];
let a = pos[0] * vectorA[0] + pos[1] * vectorA[1] + pos[2] * vectorA[2];
let b = Math.sqrt(vectorA[0] * vectorA[0] + vectorA[1] * vectorA[1] + vectorA[2] * vectorA[2]);
result[i] = a / (Math.abs(b) < 0.0001 ? 1 : b);
}
let t = !result ? imagePosition : result;
return [t[0] / this.pixelSpacingX, t[1] / this.pixelSpacingY, t[2] / this.spaceBetweenSlice];
}
//转换到world坐标系 这边的xyz就是step2中的v3
this.imageData.indexToWorldVec3(this.projectToVtkVolumePosition([x, y, z]), vout);
假设我们现在设置了两个点AB，连成第一条路径，得到空间坐标(xa,ya,za)和(xb,yb,zb)，结合上面横断面的图我们可以想象，其实要做的就是这段路径经过了多少“层”，|xb-xa|和|yb-ya|这两个方向的差值比较一下，x方向的差值比较大的话就以x来取切割点，反之用y取切割点，按该方向的pixelSpacing来分割：
取点
其实最好是做成拟合曲线来取点（贝塞尔曲线），会更加准确。这边我们就先用直线来取一下（代码如下），返回的数组就是AB上所有我们需要截取pixeldata数据的点的集合_points，这些点其实就是需要重建的center，在这边z是不敏感的，因为取的是冠状面或矢状面。
// 伪代码
const _points = linePointCalc(...)
// 每段路径直线上取点 begin和end就是A B
linePointCalc(beginVolume, endVolume, pixelSpacing, direction) {
let res = []
if (direction === 'y') {
let stepLength = Math.abs(endVolume.y - beginVolume.y);
//end点取不准 但是如果接近一个spacing的话 就当成一个点计算
const sign = (endVolume.y > beginVolume.y ? 1 : -1);
if ((stepLength % pixelSpacing !== 0)) {
if ((stepLength % pixelSpacing) > (pixelSpacing / 2)) {
stepLength = (lodash.ceil(stepLength / pixelSpacing)) * pixelSpacing;
} else {
stepLength = (lodash.floor(stepLength / pixelSpacing)) * pixelSpacing;
}
}
for (let i = 0; i <= stepLength; i += pixelSpacing) {
const yPos = beginVolume.y + i * sign;
const xPos = beginVolume.x + (yPos - beginVolume.y) * (endVolume.x - beginVolume.x) / (endVolume.y - beginVolume.y);
res.push({ x: xPos, y: yPos})
}
} else {
// x方向同理
}
return res;
}
遍历_points，我们要截取的就是每个point作为center的重建面，这个点从z为0到imageHeight(前面说过一刀下去就是整个横断的高度)的所有pixeldata，每个像素值都通过双三次插值计算，基于4*4的矩阵，所以要取这个点的前后两层重建面，具体看下面代码，getCoronalData做的是设置mpr center，获取mpr重建后冠状面的imagedata，不细说了。同时要记录一下minPixelValue和maxPixelValue，为了后续cornerstoneImage的构造。
ps: 写这个文章的时候突然意识到直接用point作为center重建的面，将点做双线性插值应该也是一样的，没必要写的这么复杂，但是写都写了，这边的代码还是先按三次卷积来写吧
// 这边按 y方向diff较大来做例子
// 曲面重建图像的宽度
let width = 0;
let resultValues = [];
for (let i = 0; i < _points.length; i++) {
const currentX = _points[i].x;
const currentY = _points[i].y;
const positionX = currentX / pixelSpacingX + 1;
// y前后推移两层 每行取值x前后推移2个
const prevY = lodash.floor(currentY / pixelSpacingY) * pixelSpacingY;
const layer1 = getCoronalData(currentX, (prevY - pixelSpacingY) < 0 ? 0 : prevY - 1);
const layer2 = getCoronalData(currentX, prevY);
const layer3 = getCoronalData(currentX, prevY + pixelSpacingY);
const layer4 = getCoronalData(currentX, prevY + 2 * pixelSpacingY);
// 按高度遍历 截取那一行的pixel
for (let j = 0; j < imageHeight; j++) {
// 双三次插值
const matrix = this.getCubicMatrix(j, positionX, layer1, layer2, layer3, layer4, imageWidth)
const pointPixelData = lodash.round(this.cubicInterpolation(coordinateX, coordinateY, matrix));
resultValues[j].push(pointPixelData);
// 同时要记录minPixelValue maxPixelValue
..........
}
}
同时把双三次插值这块实现的代码也贴上来：
// 组装双三次插值的矩阵 currentX是实际空间坐标的值 coordinateY是按pixeldata arr中的坐标
getCubicMatrix(z, coordinateIndex, layer1, layer2, layer3, layer4, imageWidth) {
const start = z * imageWidth;
const end = (z + 1) * imageWidth;
const line1 = layer1.slice(start, end);
const line2 = layer2.slice(start, end);
const line3 = layer3.slice(start, end);
const line4 = layer4.slice(start, end);
const index = lodash.floor(coordinateIndex);
const pos1 = index - 1;
const pos2 = index;
const pos3 = index + 1;
const pos4 = index + 2;
let matrix = [[], [], [], []];
matrix[0][0] = line1[pos1];
matrix[0][1] = line2[pos1];
matrix[0][2] = line3[pos1];
matrix[0][3] = line4[pos1];
matrix[1][0] = line1[pos2];
matrix[1][1] = line2[pos2];
matrix[1][2] = line3[pos2];
matrix[1][3] = line4[pos2];
matrix[2][0] = line1[pos3];
matrix[2][1] = line2[pos3];
matrix[2][2] = line3[pos3];
matrix[2][3] = line4[pos3];
matrix[3][0] = line1[pos4];
matrix[3][1] = line2[pos4];
matrix[3][2] = line3[pos4];
matrix[3][3] = line4[pos4];
return matrix;
}
/**
* 双三次插值
* @param x 当前点在矩阵中的x坐标
* @param y 当前点在矩阵中的y坐标
* @param data 4*4矩阵数据
* @returns {number}
*/
cubicInterpolation(x, y, data) {
let w = 4;
let h = 4;
// 计算x和y方向的最近的4*4的坐标和权重
let wcx = this.getCubicWeight(x);
let wcy = this.getCubicWeight(y);
let wx = wcx.weight;
let wy = wcy.weight;
// 坐标
let xs = wcx.coordinate;
let ys = wcy.coordinate;
let val = 0;
// 遍历4*4的点，根据权重相加
for (let j = 0; j < 4; j++) {
let py = ys[j];
py = py < 0 ? 0 : py;
py = py > h - 1 ? h - 1 : py;
for (let i = 0; i < 4; i++) {
let px = xs[i];
px = px < 0 ? 0 : px;
px = px > w - 1 ? w - 1 : px;
// 该点的值
let dv = data[py][px];
// 该点的权重
let w_x = wx[i];
let w_y = wy[j];
// 根据加权加起来
val += (dv * w_x * w_y);
}
}
return val;
}
step4
拼接其实在上一步的代码里面已经体现了，就是resultValues这个字段的处理，比如一段路径按spacing可以切分出10个点，那这段路径最终呈现的重建图像就是10 * Height，我们图像最终显示的pixelData的数组是按行排列下去的，比如一个4宽度、3高度的图像的pixelData其实是这样的：
[
1,2,3,4,
1,2,3,4,
1,2,3,4
]
那对应上面的resultValues[j]其实是获得了这样一个数组：
[
[1,2,3,4],
[1,2,3,4],
[1,2,3,4]
]
实际的多段路径就是height * point数量的数组：
[
[1,2,3,4.....], 长度为point数量
[1,2,3,4.....],
[1,2,3,4.....]
....
长度为height
]
然后lodash.flattenDeep(resultValues)把数据展开，就是我们期望的pixeldata了
step5
最后拼装一个cornerstoneImage给cornerstone渲染，需要把上面得到的pixeldata再Int16Array一下，然后min/maxPixelValue不能搞错，不然会影响到generateLut那边。因为实现了一下路径点回退的功能，所以发现退到width为0的时候，cornerstone那边的render不知道是哪步有问题，一旦为0之后，再drawImage或者updateImage都不会显示图像了，所以这边hack了一下最小也给1
const cornerStoneImage = {
height: coronalData.imageHeight,
rows: coronalData.imageHeight,
/* hack: width为0的时候cornerstone那边的render机制会出问题 后续都显示不出来 所以默认给1*/
width: width || 1,
columns: width,
pixelData: __resultValues,
getPixelData() {
return __resultValues;
},
imageId: 'cprgenerate',
columnPixelSpacing: sagittalData.pixelSpacingX,
rowPixelSpacing: sagittalData.pixelSpacingY,
slope: 1,
windowCenter: coronalData.windowCenter,
windowWidth: coronalData.windowLevel,
sizeInBytes: coronalData.imageHeight * width,
color: false,
invert: false,
intercept: coronalData.intercept,
minPixelValue: minPixelValue || 0,
maxPixelValue: maxPixelValue || 0
}
最后
以上是我的实现思路，下图是之前截的大体效果，左侧是横断的薄层，滚动才能看到完整胸骨走向，没仔细选点所以前面有几个比较大的断层可以无视，但后面那些马赛克纹路比较明显的，这是比较困惑我的地方，想要得到更平滑的重建图或许还是得做拟合曲线，去拟合左侧选点贴合在骨骼实际路径中心点的位置，如果大家有好的办法或者例子，麻烦告诉我，感谢~
效果图
参考：
CT后处理技巧-CPR https://mp.weixin.qq.com/s/nwRsx3nSowDvBAfoCYAvhw
js双三次插值 http://events.jianshu.io/p/3eec04ef39a1
网友评论