ImageStreamline

Introduction

vtkImageStreamline - integrate streamlines in a vtkImageData

vtkImageStreamline is a filter that generates streamlines from a vtkImageData
input over which a vector field is defined. This filter will look for vectors
(i.e. getVectors()) in the input. It will then integrate these vectors, using
Runge-Kutta 2, from a starting set of seeds defined by the points of the 2nd
input until a specified maximum number of steps is reached or until the
streamline leaves the domain.

The output will be a vtkPolyData which contains a polyline for each
streamline. Currently, this filter does not interpolate any input fields to
the points of the streamline.

Methods

computeNextStep

Argument Type Required Description
velArray Yes
image Yes
delT Number Yes
xyz Array. Yes

computeStructuredCoordinates

Argument Type Required Description
x Vector3 Yes
ijk Vector3 Yes
pcoords Vector3 Yes
extent Extent Yes
spacing Vector3 Yes
origin Vector3 Yes
bounds Bounds Yes

extend

Method used to decorate a given object (publicAPI+model) with vtkImageStreamline characteristics.

Argument Type Required Description
publicAPI Yes object on which methods will be bounds (public)
model Yes object on which data structure will be bounds (protected)
initialValues IImageStreamlineInitialValues No (default: {})

getIntegrationStep

Get the step length (delT) used during integration.

getMaximumNumberOfSteps

Get the number of steps to be used in the integration.

getVoxelIndices

Argument Type Required Description
ijk Vector3 Yes
dims Vector2 Yes
ids Array. Yes

interpolationFunctions

Argument Type Required Description
pcoords Vector3 Yes
sf Array. Yes

newInstance

Method used to create a new instance of vtkImageStreamline

Argument Type Required Description
initialValues IImageStreamlineInitialValues No for pre-setting some of its content

requestData

Argument Type Required Description
inData Yes
outData Yes

setIntegrationStep

Set the step length (delT) used during integration.

Argument Type Required Description
integrationStep Number Yes

setMaximumNumberOfSteps

Set the number of steps to be used in the integration.
Integration can terminal earlier if the streamline leaves the domain.

Argument Type Required Description
maximumNumberOfSteps Number Yes

streamIntegrate

Argument Type Required Description
velArray Yes
image vtkImageData Yes
seed Array. Yes
offset Number Yes

vectorAt

Argument Type Required Description
xyz Array. Yes
velArray Yes
image vtkImageData Yes
velAtArg Yes

Source

index.d.ts
import vtkImageData from "../../../Common/DataModel/ImageData";
import { vtkAlgorithm, vtkObject } from "../../../interfaces";
import { Bounds, Extent, Vector2, Vector3 } from "../../../types";

/**
*
*/
export interface IImageStreamlineInitialValues {
integrationStep?: number,
maximumNumberOfSteps?: number,
}

type vtkImageStreamlineBase = vtkObject & vtkAlgorithm;

export interface vtkImageStreamline extends vtkImageStreamlineBase {

/**
*
* @param velArray
* @param image
* @param {Number} delT
* @param {Number[]} xyz
*/
computeNextStep(velArray: any, image: any, delT: number, xyz: number[]): boolean;

/**
*
* @param {Vector3} x
* @param {Vector3} ijk
* @param {Vector3} pcoords
* @param {Extent} extent
* @param {Vector3} spacing
* @param {Vector3} origin
* @param {Bounds} bounds
*/
computeStructuredCoordinates(x: Vector3, ijk: Vector3, pcoords: Vector3, extent: Extent, spacing: Vector3, origin: Vector3, bounds: Bounds): boolean;

/**
* Get the step length (delT) used during integration.
*/
getIntegrationStep(): number;

/**
* Get the number of steps to be used in the integration.
*/
getMaximumNumberOfSteps(): number;

/**
*
* @param {Vector3} ijk
* @param {Vector2} dims
* @param {Number[]} ids
*/
getVoxelIndices(ijk: Vector3, dims: Vector2, ids: number[]): void;

/**
*
* @param {Vector3} pcoords
* @param {Number[]} sf
*/
interpolationFunctions(pcoords: Vector3, sf: number[]): void;
/**
*
* @param inData
* @param outData
*/
requestData(inData: any, outData: any): void;

/**
* Set the step length (delT) used during integration.
* @param {Number} integrationStep
*/
setIntegrationStep(integrationStep: number): boolean;

/**
* Set the number of steps to be used in the integration.
* Integration can terminal earlier if the streamline leaves the domain.
* @param {Number} maximumNumberOfSteps
*/
setMaximumNumberOfSteps(maximumNumberOfSteps: number): boolean;

/**
*
* @param velArray
* @param {vtkImageData} image
* @param {Number[]} seed
* @param {Number} offset
*/
streamIntegrate(velArray: any, image: vtkImageData, seed: number[], offset: number): any[];

/**
*
* @param {Number[]} xyz
* @param velArray
* @param {vtkImageData} image
* @param velAtArg
*/
vectorAt(xyz: number[], velArray: any, image: vtkImageData, velAtArg: any): boolean;
}

/**
* Method used to decorate a given object (publicAPI+model) with vtkImageStreamline characteristics.
*
* @param publicAPI object on which methods will be bounds (public)
* @param model object on which data structure will be bounds (protected)
* @param {IImageStreamlineInitialValues} [initialValues] (default: {})
*/
export function extend(publicAPI: object, model: object, initialValues?: IImageStreamlineInitialValues): void;

/**
* Method used to create a new instance of vtkImageStreamline
* @param {IImageStreamlineInitialValues} [initialValues] for pre-setting some of its content
*/
export function newInstance(initialValues?: IImageStreamlineInitialValues): vtkImageStreamline;


/**
* vtkImageStreamline - integrate streamlines in a vtkImageData
*
* vtkImageStreamline is a filter that generates streamlines from a vtkImageData
* input over which a vector field is defined. This filter will look for vectors
* (i.e. getVectors()) in the input. It will then integrate these vectors, using
* Runge-Kutta 2, from a starting set of seeds defined by the points of the 2nd
* input until a specified maximum number of steps is reached or until the
* streamline leaves the domain.
*
* The output will be a vtkPolyData which contains a polyline for each
* streamline. Currently, this filter does not interpolate any input fields to
* the points of the streamline.
*/
export declare const vtkImageStreamline: {
newInstance: typeof newInstance;
extend: typeof extend;
}
export default vtkImageStreamline;
index.js
import macro from 'vtk.js/Sources/macros';
import vtkPolyData from 'vtk.js/Sources/Common/DataModel/PolyData';

const { vtkErrorMacro } = macro;

// ----------------------------------------------------------------------------
// vtkImageStreamline methods
// ----------------------------------------------------------------------------

function vtkImageStreamline(publicAPI, model) {
// Set our className
model.classHierarchy.push('vtkImageStreamline');

const indices = new Int32Array(3);
const paramCoords = new Float32Array(3);
const weights = new Float32Array(8);
const voxelIndices = new Uint32Array(8);
const dimensions = new Uint32Array(3);
const velAt = new Float32Array(3);
const xtmp = new Float32Array(3);

publicAPI.interpolationFunctions = (pcoords, sf) => {
const r = pcoords[0];
const s = pcoords[1];
const t = pcoords[2];

const rm = 1.0 - r;
const sm = 1.0 - s;
const tm = 1.0 - t;

sf[0] = rm * sm * tm;
sf[1] = r * sm * tm;
sf[2] = rm * s * tm;
sf[3] = r * s * tm;
sf[4] = rm * sm * t;
sf[5] = r * sm * t;
sf[6] = rm * s * t;
sf[7] = r * s * t;
};

publicAPI.computeStructuredCoordinates = (
x,
ijk,
pcoords,
extent,
spacing,
origin,
bounds
) => {
// tolerance is needed for 2D data (this is squared tolerance)
const tol2 = 1e-12;
//
// Compute the ijk location
//
let isInBounds = true;
for (let i = 0; i < 3; i++) {
const d = x[i] - origin[i];
const doubleLoc = d / spacing[i];
// Floor for negative indexes.
ijk[i] = Math.floor(doubleLoc);
pcoords[i] = doubleLoc - ijk[i];

let tmpInBounds = false;
const minExt = extent[i * 2];
const maxExt = extent[i * 2 + 1];

// check if data is one pixel thick
if (minExt === maxExt) {
const dist = x[i] - bounds[2 * i];
if (dist * dist <= spacing[i] * spacing[i] * tol2) {
pcoords[i] = 0.0;
ijk[i] = minExt;
tmpInBounds = true;
}
} else if (ijk[i] < minExt) {
if (
(spacing[i] >= 0 && x[i] >= bounds[i * 2]) ||
(spacing[i] < 0 && x[i] <= bounds[i * 2 + 1])
) {
pcoords[i] = 0.0;
ijk[i] = minExt;
tmpInBounds = true;
}
} else if (ijk[i] >= maxExt) {
if (
(spacing[i] >= 0 && x[i] <= bounds[i * 2 + 1]) ||
(spacing[i] < 0 && x[i] >= bounds[i * 2])
) {
// make sure index is within the allowed cell index range
pcoords[i] = 1.0;
ijk[i] = maxExt - 1;
tmpInBounds = true;
}
} else {
tmpInBounds = true;
}

// clear isInBounds if out of bounds for this dimension
isInBounds = isInBounds && tmpInBounds;
}

return isInBounds;
};

publicAPI.getVoxelIndices = (ijk, dims, ids) => {
ids[0] = ijk[2] * dims[0] * dims[1] + ijk[1] * dims[0] + ijk[0];
ids[1] = ids[0] + 1; // i+1, j, k
ids[2] = ids[0] + dims[0]; // i, j+1, k
ids[3] = ids[2] + 1; // i+1, j+1, k
ids[4] = ids[0] + dims[0] * dims[1]; // i, j, k+1
ids[5] = ids[4] + 1; // i+1, j, k+1
ids[6] = ids[4] + dims[0]; // i, j+1, k+1
ids[7] = ids[6] + 1; // i+1, j+1, k+1
};

publicAPI.vectorAt = (xyz, velArray, image, velAtArg) => {
if (
!publicAPI.computeStructuredCoordinates(
xyz,
indices,
paramCoords,
image.getExtent(),
image.getSpacing(),
image.getOrigin(),
image.getBounds()
)
) {
return false;
}

publicAPI.interpolationFunctions(paramCoords, weights);
const extent = image.getExtent();
dimensions[0] = extent[1] - extent[0] + 1;
dimensions[1] = extent[3] - extent[2] + 1;
dimensions[2] = extent[5] - extent[4] + 1;
publicAPI.getVoxelIndices(indices, dimensions, voxelIndices);
velAtArg[0] = 0.0;
velAtArg[1] = 0.0;
velAtArg[2] = 0.0;
const vel = new Array(3);
for (let i = 0; i < 8; i++) {
velArray.getTuple(voxelIndices[i], vel);
for (let j = 0; j < 3; j++) {
velAtArg[j] += weights[i] * vel[j];
}
}

return true;
};

publicAPI.computeNextStep = (velArray, image, delT, xyz) => {
// This does Runge-Kutta 2

// Start with evaluating velocity @ initial point
if (!publicAPI.vectorAt(xyz, velArray, image, velAt)) {
return false;
}
// Now find the mid point
for (let i = 0; i < 3; i++) {
xtmp[i] = xyz[i] + (delT / 2.0) * velAt[i];
}
// Use the velocity @ that point to project
if (!publicAPI.vectorAt(xtmp, velArray, image, velAt)) {
return false;
}
for (let i = 0; i < 3; i++) {
xyz[i] += delT * velAt[i];
}

if (!publicAPI.vectorAt(xyz, velArray, image, velAt)) {
return false;
}

return true;
};

publicAPI.streamIntegrate = (velArray, image, seed, offset) => {
const retVal = [];

const maxSteps = model.maximumNumberOfSteps;
const delT = model.integrationStep;
const xyz = new Float32Array(3);
xyz[0] = seed[0];
xyz[1] = seed[1];
xyz[2] = seed[2];

const pointsBuffer = [];

let step = 0;
for (step = 0; step < maxSteps; step++) {
if (!publicAPI.computeNextStep(velArray, image, delT, xyz)) {
break;
}
for (let i = 0; i < 3; i++) {
pointsBuffer[3 * step + i] = xyz[i];
}
}

const pd = vtkPolyData.newInstance();

const points = new Float32Array(pointsBuffer);
retVal[0] = points;

pd.getPoints().setData(points, 3);

const npts = points.length / 3;
const line = new Uint32Array(npts + 1);
line[0] = npts;
for (let i = 0; i < npts; i++) {
line[i + 1] = i + offset;
}
retVal[1] = line;

pd.getLines().setData(line);

return retVal;
};

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

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

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

const seedPts = seeds.getPoints();
const nSeeds = seedPts.getNumberOfPoints();

let offset = 0;
const datas = [];
const vectors = input.getPointData().getVectors();
const point = [];
for (let i = 0; i < nSeeds; i++) {
seedPts.getTuple(i, point);
const retVal = publicAPI.streamIntegrate(vectors, input, point, offset);
offset += retVal[0].length / 3;
datas.push(retVal);
}

let cellArrayLength = 0;
let pointArrayLength = 0;
datas.forEach((data) => {
cellArrayLength += data[1].length;
pointArrayLength += data[0].length;
});
offset = 0;
let offset2 = 0;
const cellArray = new Uint32Array(cellArrayLength);
const pointArray = new Float32Array(pointArrayLength);
datas.forEach((data) => {
cellArray.set(data[1], offset);
offset += data[1].length;
pointArray.set(data[0], offset2);
offset2 += data[0].length;
});

const output = vtkPolyData.newInstance();
output.getPoints().setData(pointArray, 3);
output.getLines().setData(cellArray);

outData[0] = output;
};
}

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

const DEFAULT_VALUES = {
integrationStep: 1,
maximumNumberOfSteps: 1000,
};

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

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, 2, 1);

// Generate macros for properties
macro.setGet(publicAPI, model, ['integrationStep', 'maximumNumberOfSteps']);

// Object specific methods
vtkImageStreamline(publicAPI, model);
}

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

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

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

export default { newInstance, extend };