How do I calculate the orientation in aorta_centerline.json used in the ImageCPRMapper example? I saw in another question the matrix obtained by combining Normal, Tangent and Bitangent, but I don’t know how to calculate these parameters. I first resampled the curve. Here is my code. Can you help me?
import { vec3, mat3, quat, mat4 } from 'gl-matrix';
/**
* Calculates the total length of a curve defined by points.
* @param {Array<vec3>} points - Array of 3D points (gl-matrix vec3).
* @returns {number} Total length of the curve.
*/
function calculateCurveLength(points) {
let length = 0;
for (let i = 1; i < points.length; i++) {
length += vec3.distance(points[i - 1], points[i]);
}
return length;
}
/**
* Resamples a curve to have points approximately equidistant.
* Simple linear interpolation implementation.
* @param {Array<vec3>} points - Array of 3D points (gl-matrix vec3).
* @param {number} desiredSpacing - The desired distance between points.
* @returns {Array<vec3>} Array of resampled 3D points.
*/
function resampleCurve(points, desiredSpacing) {
if (!points || points.length < 2 || desiredSpacing <= 0) {
return points ? [...points] : []; // Return copy or empty
}
const resampledPoints = [vec3.clone(points[0])];
let currentPoint = vec3.clone(points[0]);
let distanceAccumulated = 0;
let segmentIndex = 0;
const totalLength = calculateCurveLength(points);
if (totalLength === 0) {
return [vec3.clone(points[0])]; // Single point or all points coincident
}
const numSamples = Math.max(2, Math.round(totalLength / desiredSpacing) + 1);
const actualSpacing = totalLength / (numSamples - 1); // Adjust spacing to fit length
let distanceAlongCurve = actualSpacing; // Target distance for the next point
while (segmentIndex < points.length - 1 && resampledPoints.length < numSamples) {
const startPoint = points[segmentIndex];
const endPoint = points[segmentIndex + 1];
const segmentVector = vec3.subtract(vec3.create(), endPoint, startPoint);
const segmentLength = vec3.length(segmentVector);
if (segmentLength > 1e-6) {
// Avoid division by zero for coincident points
while (distanceAlongCurve <= distanceAccumulated + segmentLength) {
const distanceNeeded = distanceAlongCurve - distanceAccumulated;
const fraction = distanceNeeded / segmentLength;
const newPoint = vec3.scaleAndAdd(vec3.create(), startPoint, segmentVector, fraction);
resampledPoints.push(newPoint);
currentPoint = newPoint; // Update current position
if (resampledPoints.length >= numSamples) break;
distanceAlongCurve += actualSpacing; // Move to next target distance
}
}
distanceAccumulated += segmentLength;
segmentIndex++;
if (resampledPoints.length >= numSamples) break;
}
// Ensure the last point is included if slightly short due to precision
if (
resampledPoints.length < numSamples &&
vec3.distance(resampledPoints[resampledPoints.length - 1], points[points.length - 1]) > 1e-6
) {
resampledPoints.push(vec3.clone(points[points.length - 1]));
} else if (resampledPoints.length === numSamples && numSamples > 0) {
// Ensure the last point *is* the original last point for accuracy
vec3.copy(resampledPoints[resampledPoints.length - 1], points[points.length - 1]);
}
return resampledPoints;
}
/**
* Calculates the centroid of a set of points.
* @param {Array<vec3>} points - Array of 3D points (gl-matrix vec3).
* @returns {vec3} The centroid.
*/
function calculateCentroid(points) {
const centroid = vec3.create();
if (!points || points.length === 0) return centroid;
for (const p of points) {
vec3.add(centroid, centroid, p);
}
vec3.scale(centroid, centroid, 1.0 / points.length);
return centroid;
}
/**
* Computes the 4x4 transformation matrix for Curved Planar Reformat.
* This matrix defines the orientation and position of the straightened
* coordinate system in world space.
* Focuses on the 'straightening' mode (minimizing torsion).
*
* @param {Array<Array<number>>} originalCurvePointsWorld - Array of [x, y, z] points defining the curve.
* @param {Array<number>} sliceSizeMm - [width, height] of the desired slice in mm.
* @param {number} spacingAlongCurveMm - Desired spacing between slices along the straightened curve in mm.
* @param {number} [rotationDeg=0.0] - Optional rotation (in degrees) around the curve's axis (Z).
* @param {number} [resamplingSpacingFactor=5.0] - How much coarser the resampling for orientation calculation is compared to output spacing.
* @returns {mat4 | null} The 4x4 transformation matrix (gl-matrix mat4) or null if computation fails.
*/
export function computeStraighteningTransformMatrix(
originalCurvePointsWorld,
sliceSizeMm,
spacingAlongCurveMm,
rotationDeg = 0.0,
viewCamera,
resamplingSpacingFactor = 5.0 // Matches Python logic
) {
if (!originalCurvePointsWorld || originalCurvePointsWorld.length < 2) {
return null;
}
if (!sliceSizeMm || sliceSizeMm.length < 2) {
return null;
}
if (spacingAlongCurveMm <= 0) {
return null;
}
const originalCurveVec3 = originalCurvePointsWorld.map(p => vec3.fromValues(p[0], p[1], p[2]));
const resamplingCurveSpacing = spacingAlongCurveMm * resamplingSpacingFactor;
const resampledPoints = resampleCurve(originalCurveVec3, resamplingCurveSpacing);
if (resampledPoints.length < 2) {
return null;
}
const numResampledPoints = resampledPoints.length;
const curveStartPoint = resampledPoints[0];
const curveEndPoint = resampledPoints[numResampledPoints - 1];
const transformGridAxisZ = vec3.subtract(vec3.create(), curveEndPoint, curveStartPoint);
const curveLength = vec3.length(transformGridAxisZ);
if (curveLength < 1e-6) {
return null;
}
vec3.normalize(transformGridAxisZ, transformGridAxisZ);
const transformGridAxisY = vec3.copy(vec3.create(), viewCamera.viewPlaneNormal);
const transformGridAxisX = vec3.cross(vec3.create(), transformGridAxisZ, transformGridAxisY);
const worldUp = vec3.fromValues(0, 0, 1); // Try Z up first
let dotProd = Math.abs(vec3.dot(worldUp, transformGridAxisZ));
if (dotProd > 0.99) {
vec3.set(worldUp, 0, 1, 0);
dotProd = Math.abs(vec3.dot(worldUp, transformGridAxisZ));
if (dotProd > 0.99) {
vec3.set(worldUp, 1, 0, 0);
}
}
// Calculate X = normalize(cross(up, Z))
vec3.cross(transformGridAxisX, worldUp, transformGridAxisZ);
if (vec3.length(transformGridAxisX) < 1e-6) {
console.error(
"Failed to calculate X axis (Z and Up likely parallel). This shouldn't happen with the checks above."
);
// Fallback: try cross(Z, world_X_or_Y_that_is_not_up)
const fallbackUp = vec3.equals(worldUp, [0, 0, 1]) ? [0, 1, 0] : [0, 0, 1]; // Choose different axis
vec3.cross(transformGridAxisX, fallbackUp, transformGridAxisZ);
if (vec3.length(transformGridAxisX) < 1e-6) {
console.error('Critical failure to find orthogonal X axis.');
return null;
}
}
vec3.normalize(transformGridAxisX, transformGridAxisX);
// Calculate Y = normalize(cross(Z, X))
vec3.cross(transformGridAxisY, transformGridAxisZ, transformGridAxisX);
vec3.normalize(transformGridAxisY, transformGridAxisY); // Should already be normalized if Z,X are orthonormal
if (Math.abs(rotationDeg) > 1e-6) {
const rotationRad = rotationDeg * (Math.PI / 180.0);
const rotationMatrix = mat4.fromZRotation(mat4.create(), rotationRad);
vec3.transformMat4(transformGridAxisX, transformGridAxisX, rotationMatrix);
vec3.transformMat4(transformGridAxisY, transformGridAxisY, rotationMatrix);
}
const curveCenter = calculateCentroid(resampledPoints); // 重采样曲线的中心
const offsetX = vec3.scale(vec3.create(), transformGridAxisX, sliceSizeMm[0] / 2.0);
const offsetY = vec3.scale(vec3.create(), transformGridAxisY, sliceSizeMm[1] / 2.0);
const offsetZ = vec3.scale(vec3.create(), transformGridAxisZ, curveLength / 2.0); // 从开始/终点使用长度
// Calculate final origin: center - offsetX - offsetY - offsetZ
const transformGridOrigin = vec3.create();
vec3.subtract(transformGridOrigin, curveCenter, offsetX);
vec3.subtract(transformGridOrigin, transformGridOrigin, offsetY);
vec3.subtract(transformGridOrigin, transformGridOrigin, offsetZ);
const transformMatrix = mat4.create(); // Creates identity
console.log('transformMatrix', transformMatrix);
// Set columns directly
mat4.set(
transformMatrix,
transformGridAxisX[0],
transformGridAxisX[1],
transformGridAxisX[2],
0, // Col 0 = X axis
transformGridAxisY[0],
transformGridAxisY[1],
transformGridAxisY[2],
0, // Col 1 = Y axis
transformGridAxisZ[0],
transformGridAxisZ[1],
transformGridAxisZ[2],
0, // Col 2 = Z axis
0,
0,
0,
1 // Col 3 = Origin
);
return transformMatrix;
}