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
`velArray`		Yes
`image`		Yes
`delT`	Number	Yes
`xyz`	Array.	Yes

computeStructuredCoordinates

Argument	Type	Required
`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
`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
`velArray`		Yes
`image`	vtkImageData	Yes
`seed`	Array.	Yes
`offset`	Number	Yes

vectorAt

Argument	Type	Required
`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 };