ImageMarchingSquares

Introduction

vtkImageMarchingSquares - isocontour an image (or slice from a volume)

Given a specified contour value, generate isolines using the
Marching Squares algorithm (the 2D version of the 3D Marching Cubes
algorithm).

Public API

contourValues

Set/Get an array of isocontour values.

slice

Set/Get the IJK slice number of the input volume. By default the
slice = 0.

slicingMode

Set/Get the slicing mode (IJK) which determines the direction to slice. By
default the slicingMode = 2.

mergePoints

As lines forming the isolines are generated, indicate whether
conincident points are to be merged. Merging produces connected polylines
at the cost of additional memory and computation.

Source

caseTable.js
// ----------------------------------------------------------------------------
// Marching squares case functions (using lines to generate the 2D tessellation).
// For each case, a list of edge ids that form the triangles. A -1 marks the
// end of the list of edges. Edges are taken three at a time to generate
// triangle points.
// ----------------------------------------------------------------------------
const MARCHING_SQUARES_CASES = [
[-1, -1, -1, -1, -1] /* 0 */,
[0, 3, -1, -1, -1] /* 1 */,
[1, 0, -1, -1, -1] /* 2 */,
[1, 3, -1, -1, -1] /* 3 */,
[2, 1, -1, -1, -1] /* 4 */,
[0, 3, 2, 1, -1] /* 5 */,
[2, 0, -1, -1, -1] /* 6 */,
[2, 3, -1, -1, -1] /* 7 */,
[3, 2, -1, -1, -1] /* 8 */,
[0, 2, -1, -1, -1] /* 9 */,
[1, 0, 3, 2, -1] /* 10 */,
[1, 2, -1, -1, -1] /* 11 */,
[3, 1, -1, -1, -1] /* 12 */,
[0, 1, -1, -1, -1] /* 13 */,
[3, 0, -1, -1, -1] /* 14 */,
[-1, -1, -1, -1, -1] /* 15 */,
];

const EDGES = [
[0, 1],
[1, 3],
[2, 3],
[0, 2],
];

function getCase(index) {
return MARCHING_SQUARES_CASES[index];
}

// Define the four edges of the pixel by the following pairs of vertices
function getEdge(eid) {
return EDGES[eid];
}

// ----------------------------------------------------------------------------
// Static API
// ----------------------------------------------------------------------------
export default {
getCase,
getEdge,
};
index.js
import macro from 'vtk.js/Sources/macros';
import vtkPolyData from 'vtk.js/Sources/Common/DataModel/PolyData';
import vtkEdgeLocator from 'vtk.js/Sources/Common/DataModel/EdgeLocator';

import vtkCaseTable from './caseTable';

const { vtkErrorMacro, vtkDebugMacro } = macro;

// ----------------------------------------------------------------------------
// vtkImageMarchingSquares methods
// ----------------------------------------------------------------------------

function vtkImageMarchingSquares(publicAPI, model) {
/**
* Get the X,Y kernels based on the set slicing mode.
* @returns {[number, number]}
*/
function getKernels() {
let kernelX = 0; // default K slicing mode
let kernelY = 1;
if (model.slicingMode === 1) {
kernelX = 0;
kernelY = 2;
} else if (model.slicingMode === 0) {
kernelX = 1;
kernelY = 2;
}

return [kernelX, kernelY];
}

// Set our className
model.classHierarchy.push('vtkImageMarchingSquares');

/**
* Get the list of contour values.
* @returns {number[]}
*/
publicAPI.getContourValues = () => model.contourValues;

/**
* Set the list contour values.
* @param {number[]} cValues
*/
publicAPI.setContourValues = (cValues) => {
model.contourValues = cValues;
publicAPI.modified();
};

const ids = [];
const pixelScalars = [];
const pixelPts = [];
const edgeLocator = vtkEdgeLocator.newInstance();

/**
* Retrieve scalars and pixel coordinates.
* @param {Vector3} ijk origin of the pixel
* @param {Vector3} dims dimensions of the image
* @param {TypedArray} scalars list of scalar values
* @param {Vector3} increments IJK slice increments
* @param {number} kernelX index of the X element
* @param {number} kernelY index of the Y element
*/
publicAPI.getPixelScalars = (
ijk,
dims,
scalars,
increments,
kernelX,
kernelY
) => {
const [i, j, k] = ijk;

// First get the indices for the pixel
ids[0] = k * dims[1] * dims[0] + j * dims[0] + i; // i, j, k
ids[1] = ids[0] + increments[kernelX]; // i+1, j, k
ids[2] = ids[0] + increments[kernelY]; // i, j+1, k
ids[3] = ids[2] + increments[kernelX]; // i+1, j+1, k

// Now retrieve the scalars
for (let ii = 0; ii < 4; ++ii) {
pixelScalars[ii] = scalars[ids[ii]];
}
};

/**
* Retrieve pixel coordinates.
* @param {Vector3} ijk origin of the pixel
* @param {Vector3} origin origin of the image
* @param {Vector3} spacing spacing of the image
* @param {number} kernelX index of the X element
* @param {number} kernelY index of the Y element
*/
publicAPI.getPixelPoints = (ijk, origin, spacing, kernelX, kernelY) => {
const i = ijk[kernelX];
const j = ijk[kernelY];

// (i,i+1),(j,j+1),(k,k+1) - i varies fastest; then j; then k
pixelPts[0] = origin[kernelX] + i * spacing[kernelX]; // 0
pixelPts[1] = origin[kernelY] + j * spacing[kernelY];

pixelPts[2] = pixelPts[0] + spacing[kernelX]; // 1
pixelPts[3] = pixelPts[1];

pixelPts[4] = pixelPts[0]; // 2
pixelPts[5] = pixelPts[1] + spacing[kernelY];

pixelPts[6] = pixelPts[2]; // 3
pixelPts[7] = pixelPts[5];
};

/**
* Produce points and lines for the polydata.
* @param {number[]} cVal list of contour values
* @param {Vector3} ijk origin of the pixel
* @param {Vector3} dims dimensions of the image
* @param {Vector3} origin origin of the image
* @param {Vector3} spacing sapcing of the image
* @param {TypedArray} scalars list of scalar values
* @param {number[]} points list of points
* @param {number[]} lines list of lines
* @param {Vector3} increments IJK slice increments
* @param {number} kernelX index of the X element
* @param {number} kernelY index of the Y element
*/
publicAPI.produceLines = (
cVal,
ijk,
dims,
origin,
spacing,
scalars,
points,
lines,
increments,
kernelX,
kernelY
) => {
const k = ijk[model.slicingMode];
const CASE_MASK = [1, 2, 8, 4]; // case table is actually for quad
const xyz = [];
let pId;

publicAPI.getPixelScalars(ijk, dims, scalars, increments, kernelX, kernelY);

let index = 0;
for (let idx = 0; idx < 4; idx++) {
if (pixelScalars[idx] >= cVal) {
index |= CASE_MASK[idx]; // eslint-disable-line no-bitwise
}
}

const pixelLines = vtkCaseTable.getCase(index);
if (pixelLines[0] < 0) {
return; // don't get the pixel coordinates, nothing to do
}

publicAPI.getPixelPoints(ijk, origin, spacing, kernelX, kernelY);

const z = origin[model.slicingMode] + k * spacing[model.slicingMode];
for (let idx = 0; pixelLines[idx] >= 0; idx += 2) {
lines.push(2);
for (let eid = 0; eid < 2; eid++) {
const edgeVerts = vtkCaseTable.getEdge(pixelLines[idx + eid]);
pId = undefined;
if (model.mergePoints) {
pId = edgeLocator.isInsertedEdge(
ids[edgeVerts[0]],
ids[edgeVerts[1]]
)?.value;
}
if (pId === undefined) {
const t =
(cVal - pixelScalars[edgeVerts[0]]) /
(pixelScalars[edgeVerts[1]] - pixelScalars[edgeVerts[0]]);
const x0 = pixelPts.slice(edgeVerts[0] * 2, (edgeVerts[0] + 1) * 2);
const x1 = pixelPts.slice(edgeVerts[1] * 2, (edgeVerts[1] + 1) * 2);
xyz[kernelX] = x0[0] + t * (x1[0] - x0[0]);
xyz[kernelY] = x0[1] + t * (x1[1] - x0[1]);
xyz[model.slicingMode] = z;
pId = points.length / 3;
points.push(xyz[0], xyz[1], xyz[2]);

if (model.mergePoints) {
edgeLocator.insertEdge(ids[edgeVerts[0]], ids[edgeVerts[1]], pId);
}
}
lines.push(pId);
}
}
};

publicAPI.requestData = (inData, outData) => {
// implement requestData
const input = inData[0];

if (!input) {
vtkErrorMacro('Invalid or missing input');
return;
}

if (
model.slicingMode == null ||
model.slicingMode < 0 ||
model.slicingMode > 2
) {
vtkErrorMacro('Invalid or missing slicing mode');
return;
}

console.time('msquares');

// Retrieve output and volume data
const origin = input.getOrigin();
const spacing = input.getSpacing();
const dims = input.getDimensions();
const extent = input.getExtent();
const increments = input.computeIncrements(extent);
const scalars = input.getPointData().getScalars().getData();
const [kernelX, kernelY] = getKernels();

// Points - dynamic array
const points = [];

// Cells - dynamic array
const lines = [];

// Ensure slice is valid
let k = Math.round(model.slice);
if (k >= dims[model.slicingMode]) {
k = 0;
}

// Loop over all contour values, and then pixels, determine case and process
const ijk = [0, 0, 0];
ijk[model.slicingMode] = k;
for (let cv = 0; cv < model.contourValues.length; ++cv) {
for (let j = 0; j < dims[kernelY] - 1; ++j) {
ijk[kernelY] = j;
for (let i = 0; i < dims[kernelX] - 1; ++i) {
ijk[kernelX] = i;

publicAPI.produceLines(
model.contourValues[cv],
ijk,
dims,
origin,
spacing,
scalars,
points,
lines,
increments,
kernelX,
kernelY
);
}
}

edgeLocator.initialize();
}

// Update output
const polydata = vtkPolyData.newInstance();
polydata.getPoints().setData(new Float32Array(points), 3);
polydata.getLines().setData(new Uint32Array(lines));
outData[0] = polydata;

vtkDebugMacro('Produced output');
console.timeEnd('msquares');
};
}

// ----------------------------------------------------------------------------
// Object factory
// ----------------------------------------------------------------------------

const DEFAULT_VALUES = {
contourValues: [],
slicingMode: 2,
slice: 0,
mergePoints: false,
};

// ----------------------------------------------------------------------------

export function extend(publicAPI, model, initialValues = {}) {
Object.assign(model, DEFAULT_VALUES, initialValues);

// Make this a VTK object
macro.obj(publicAPI, model);

// Also make it an algorithm with one input and one output
macro.algo(publicAPI, model, 1, 1);

macro.setGet(publicAPI, model, ['slicingMode', 'slice', 'mergePoints']);

// Object specific methods
macro.algo(publicAPI, model, 1, 1);
vtkImageMarchingSquares(publicAPI, model);
}

// ----------------------------------------------------------------------------

export const newInstance = macro.newInstance(extend, 'vtkImageMarchingSquares');

// ----------------------------------------------------------------------------

export default { newInstance, extend };