ImageData

Introduction

vtkImageData is a data object that is a concrete implementation of vtkDataSet.
vtkImageData represents a geometric structure that is a topological and geometrical regular array of points. Examples include volumes (voxel data) and pixmaps.
All vtkDataSet functions are inherited.

getExtent(), setExtent(array)

The extent of a dataset is a set of 6 integers. It says what the first and last pixel indices are in each of the three directions. E.g.

extent = [ i_min, i_max, j_min, j_max, k_min, k_max ]

getDimensions(), setDimensions(array[i,j,k]), setDimensions(i,j,k)

Get the dimensions of the image data as an array of 3 points defining the i,j,k dimensions of the data.
Set the values of the extent, from 0 to (i-1), etc.

getBounds()

The Bounds of a vtkImage are returned as pairs of world coordinates

[x_min, x_max, y_min, y_max, z_min, z_max]

these are calculated from the Extent, Origin, and Spacing, defined through

bounds[6] =
[
  i_min*Spacing[0] + Origin[0], i_max*Spacing[0] + Origin[0],
  j_min*Spacing[1] + Origin[1], j_max*Spacing[1] + Origin[1],
  k_min*Spacing[2] + Origin[2], k_max*Spacing[2] + Origin[2]
];

You can’t directly set the bounds. First you need to decide how many pixels across your image will be (i.e. what the extent should be), and then you must find the origin and spacing that will produce the bounds that you need from the extent that you have. This is simple algebra.

In general, always set the extent to start at zero, e.g. [0, 9, 0, 9, 0, 9] for a 10x10x10 image. Calling setDimensions(10,10,10) does exactly the same thing as setExtent(0,9,0,9,0,9) but you should always do the latter to be explicit about where your extent starts.

getDirection(), setDirection(array[9]), setDirection(…array)

Direction is a mat3 matrix corresponding to the axes directions in world coordinates for the I, J, K axes of the image. Direction must form an orthonormal basis. setDirection can be called with an array of length 9, or each part provided as individual arguments.

getNumberOfCells()

Standard vtkDataSet API method.

getNumberOfPoints()

Determine the number of points composing the dataset.

getPoint(index)

Returns the world position of a data point. Index is the point’s index in the 1D data array.

indexToWorld(in, out), indexToWorldVec3(vin, vout)

Converts the input index vector [i,j,k] to world values [x,y,z]. Modifies the out vector array in place, but also returns it.
Can be sped up by providing gl-matrix vec3 objects directly to indexToWorldVec3

worldToIndex(in, out), worldToIndexVec3(vin, vout)

Converts the input world vector [x,y,z] to approximate index values [i,j,k]. Should be rounded to integers before attempting to access the index. Modifies the out vector array in place, but also returns it.
Can be sped up by providing gl-matrix vec3 objects directly to worldToIndexVec3

getIndexToWorld(), getWorldToIndex()

Returns the mat4 matrices used to convert between world and index. worldToIndex is the inverse matrix of indexToWorld. Both are made with Float64Array.

indexToWorldBounds(in[6], out[6]?)

Calculate the corresponding world bounds for the given index bounds [i_min, i_max, j_min, j_max, k_min, k_max]. Modifies out in place if provided, or returns a new array.

worldToIndexBounds(in[6], out[6]?)

Calculate the corresponding index bounds for the given world bounds [x_min, x_max, y_min, y_max, z_min, z_max]. Modifies out in place if provided, or returns a new array.

getCenter()

Returns an [x,y,z] location of the center of the imageData.

getOffsetIndexFromWorld(vec3)

Returns the data array index for the point at the provided world position.

getScalarValueFromWorld(vec3, comp)

Returns the scalar value for the point at the provided world position, or NaN if the world bounds are outside the volumeData bounds. comp is the scalar component index, for multi-component scalar data.

computeHistogram(worldBounds[6], voxelFunc?)

Returns an object with { minimum, maximum, average, variance, sigma } of the imageData points found within the provided worldBounds.

voxelFunc(index, bounds) is an optional function that is called with the [i,j,k] index and index bounds, expected to return truthy if the data point should be counted in the histogram, and falsey if not.

(internal) extentToBounds(extent)

Returns a bounds array from a given Extent, useful if you need to calculate the world bounds of a subset of the imageData’s data.

(internal) computeTransforms()

Calculates the indexToWorld and worldToIndex conversion matrices from the origin, direction, and spacing. Shouldn’t need to call this as it is handled internally, and updated whenever the vtkImageData is modified.

(internal) computeIncrements(extent, numberOfComponents = 1)

Returns an array[3] of values to multiply an [i,j,k] index to convert into the actual data array index, from the provided extent. numberOfComponents should match the Scalar components.

(internal) computeOffsetIndex([i, j, k])

Converts an [i,j,k] index to the flat data array index. Returns NaN if any of the i,j,k bounds are outside the data Extent.

Source

index.js

import macro from 'vtk.js/Sources/macro';
import * as vtkMath from 'vtk.js/Sources/Common/Core/Math';
import vtkBoundingBox from 'vtk.js/Sources/Common/DataModel/BoundingBox';
import vtkDataSet from 'vtk.js/Sources/Common/DataModel/DataSet';
import vtkStructuredData from 'vtk.js/Sources/Common/DataModel/StructuredData';
import { StructuredType } from 'vtk.js/Sources/Common/DataModel/StructuredData/Constants';
import { vec3, mat3, mat4 } from 'gl-matrix';

const { vtkErrorMacro } = macro;

// ----------------------------------------------------------------------------
// vtkImageData methods
// ----------------------------------------------------------------------------

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

  publicAPI.setExtent = (...inExtent) => {
    if (model.deleted) {
      vtkErrorMacro('instance deleted - cannot call any method');
      return false;
    }

    const extentArray = inExtent.length === 1 ? inExtent[0] : inExtent;

    if (extentArray.length !== 6) {
      return false;
    }

    let changeDetected = false;
    model.extent.forEach((item, index) => {
      if (item !== extentArray[index]) {
        if (changeDetected) {
          return;
        }
        changeDetected = true;
      }
    });

    if (changeDetected) {
      model.extent = extentArray.slice();
      model.dataDescription = vtkStructuredData.getDataDescriptionFromExtent(
        model.extent
      );
      publicAPI.modified();
    }
    return changeDetected;
  };

  publicAPI.setDimensions = (...dims) => {
    let i;
    let j;
    let k;

    if (model.deleted) {
      vtkErrorMacro('instance deleted - cannot call any method');
      return;
    }

    if (dims.length === 1) {
      const array = dims[0];
      i = array[0];
      j = array[1];
      k = array[2];
    } else if (dims.length === 3) {
      i = dims[0];
      j = dims[1];
      k = dims[2];
    } else {
      vtkErrorMacro('Bad dimension specification');
      return;
    }

    publicAPI.setExtent(0, i - 1, 0, j - 1, 0, k - 1);
  };

  publicAPI.getDimensions = () => [
    model.extent[1] - model.extent[0] + 1,
    model.extent[3] - model.extent[2] + 1,
    model.extent[5] - model.extent[4] + 1,
  ];

  publicAPI.getNumberOfCells = () => {
    const dims = publicAPI.getDimensions();
    let nCells = 1;

    for (let i = 0; i < 3; i++) {
      if (dims[i] === 0) {
        return 0;
      }
      if (dims[i] > 1) {
        nCells *= dims[i] - 1;
      }
    }

    return nCells;
  };

  publicAPI.getNumberOfPoints = () => {
    const dims = publicAPI.getDimensions();
    return dims[0] * dims[1] * dims[2];
  };

  publicAPI.getPoint = (index) => {
    const dims = publicAPI.getDimensions();
    const ijk = vec3.fromValues(0, 0, 0);
    const coords = [0, 0, 0];

    if (dims[0] === 0 || dims[1] === 0 || dims[2] === 0) {
      vtkErrorMacro('Requesting a point from an empty image.');
      return null;
    }

    switch (model.dataDescription) {
      case StructuredType.EMPTY:
        return null;

      case StructuredType.SINGLE_POINT:
        break;

      case StructuredType.X_LINE:
        ijk[0] = index;
        break;

      case StructuredType.Y_LINE:
        ijk[1] = index;
        break;

      case StructuredType.Z_LINE:
        ijk[2] = index;
        break;

      case StructuredType.XY_PLANE:
        ijk[0] = index % dims[0];
        ijk[1] = index / dims[0];
        break;

      case StructuredType.YZ_PLANE:
        ijk[1] = index % dims[1];
        ijk[2] = index / dims[1];
        break;

      case StructuredType.XZ_PLANE:
        ijk[0] = index % dims[0];
        ijk[2] = index / dims[0];
        break;

      case StructuredType.XYZ_GRID:
        ijk[0] = index % dims[0];
        ijk[1] = (index / dims[0]) % dims[1];
        ijk[2] = index / (dims[0] * dims[1]);
        break;

      default:
        vtkErrorMacro('Invalid dataDescription');
        break;
    }

    const vout = vec3.create();
    publicAPI.indexToWorldVec3(ijk, vout);
    vec3.copy(coords, vout);
    return coords;
  };

  // vtkCell *GetCell(vtkIdType cellId) VTK_OVERRIDE;
  // void GetCell(vtkIdType cellId, vtkGenericCell *cell) VTK_OVERRIDE;
  // void GetCellBounds(vtkIdType cellId, double bounds[6]) VTK_OVERRIDE;
  // virtual vtkIdType FindPoint(double x, double y, double z)
  // {
  //   return this->vtkDataSet::FindPoint(x, y, z);
  // }
  // vtkIdType FindPoint(double x[3]) VTK_OVERRIDE;
  // vtkIdType FindCell(
  //   double x[3], vtkCell *cell, vtkIdType cellId, double tol2,
  //   int& subId, double pcoords[3], double *weights) VTK_OVERRIDE;
  // vtkIdType FindCell(
  //   double x[3], vtkCell *cell, vtkGenericCell *gencell,
  //   vtkIdType cellId, double tol2, int& subId,
  //   double pcoords[3], double *weights) VTK_OVERRIDE;
  // vtkCell *FindAndGetCell(double x[3], vtkCell *cell, vtkIdType cellId,
  //                                 double tol2, int& subId, double pcoords[3],
  //                                 double *weights) VTK_OVERRIDE;
  // int GetCellType(vtkIdType cellId) VTK_OVERRIDE;
  // void GetCellPoints(vtkIdType cellId, vtkIdList *ptIds) VTK_OVERRIDE
  //   {vtkStructuredData::GetCellPoints(cellId,ptIds,this->DataDescription,
  //                                     this->GetDimensions());}
  // void GetPointCells(vtkIdType ptId, vtkIdList *cellIds) VTK_OVERRIDE
  //   {vtkStructuredData::GetPointCells(ptId,cellIds,this->GetDimensions());}
  // void ComputeBounds() VTK_OVERRIDE;
  // int GetMaxCellSize() VTK_OVERRIDE {return 8;}; //voxel is the largest

  publicAPI.getBounds = () => publicAPI.extentToBounds(model.extent);

  publicAPI.extentToBounds = (ex) => {
    // prettier-ignore
    const corners = [
      ex[0], ex[2], ex[4],
      ex[1], ex[2], ex[4],
      ex[0], ex[3], ex[4],
      ex[1], ex[3], ex[4],
      ex[0], ex[2], ex[5],
      ex[1], ex[2], ex[5],
      ex[0], ex[3], ex[5],
      ex[1], ex[3], ex[5]];

    const idx = vec3.fromValues(corners[0], corners[1], corners[2]);
    const vout = vec3.create();
    publicAPI.indexToWorldVec3(idx, vout);
    const bounds = [vout[0], vout[0], vout[1], vout[1], vout[2], vout[2]];
    for (let i = 3; i < 24; i += 3) {
      vec3.set(idx, corners[i], corners[i + 1], corners[i + 2]);
      publicAPI.indexToWorldVec3(idx, vout);
      if (vout[0] < bounds[0]) {
        bounds[0] = vout[0];
      }
      if (vout[1] < bounds[2]) {
        bounds[2] = vout[1];
      }
      if (vout[2] < bounds[4]) {
        bounds[4] = vout[2];
      }
      if (vout[0] > bounds[1]) {
        bounds[1] = vout[0];
      }
      if (vout[1] > bounds[3]) {
        bounds[3] = vout[1];
      }
      if (vout[2] > bounds[5]) {
        bounds[5] = vout[2];
      }
    }

    return bounds;
  };

  // Internal, shouldn't need to call this manually.
  publicAPI.computeTransforms = () => {
    const trans = vec3.fromValues(
      model.origin[0],
      model.origin[1],
      model.origin[2]
    );
    mat4.fromTranslation(model.indexToWorld, trans);

    model.indexToWorld[0] = model.direction[0];
    model.indexToWorld[1] = model.direction[1];
    model.indexToWorld[2] = model.direction[2];

    model.indexToWorld[4] = model.direction[3];
    model.indexToWorld[5] = model.direction[4];
    model.indexToWorld[6] = model.direction[5];

    model.indexToWorld[8] = model.direction[6];
    model.indexToWorld[9] = model.direction[7];
    model.indexToWorld[10] = model.direction[8];

    const scale = vec3.fromValues(
      model.spacing[0],
      model.spacing[1],
      model.spacing[2]
    );
    mat4.scale(model.indexToWorld, model.indexToWorld, scale);

    mat4.invert(model.worldToIndex, model.indexToWorld);
  };

  //
  // The direction matrix is a 3x3 basis for the I, J, K axes
  // of the image. The rows of the matrix correspond to the
  // axes directions in world coordinates. Direction must
  // form an orthonormal basis, results are undefined if
  // it is not.
  //
  publicAPI.setDirection = (...args) => {
    if (model.deleted) {
      vtkErrorMacro('instance deleted - cannot call any method');
      return false;
    }

    let array = args;
    // allow an array passed as a single arg.
    if (
      array.length === 1 &&
      (Array.isArray(array[0]) ||
        array[0].constructor === Float32Array ||
        array[0].constructor === Float64Array)
    ) {
      array = array[0];
    }

    if (array.length !== 9) {
      throw new RangeError('Invalid number of values for array setter');
    }
    let changeDetected = false;
    model.direction.forEach((item, index) => {
      if (item !== array[index]) {
        if (changeDetected) {
          return;
        }
        changeDetected = true;
      }
    });

    if (changeDetected) {
      for (let i = 0; i < 9; ++i) {
        model.direction[i] = array[i];
      }
      publicAPI.modified();
    }
    return true;
  };

  // this is the fast version, requires vec3 arguments
  publicAPI.indexToWorldVec3 = (vin, vout) => {
    vec3.transformMat4(vout, vin, model.indexToWorld);
    return vout;
  };

  // slow version for generic arrays
  publicAPI.indexToWorld = (ain, aout = []) => {
    const vin = vec3.fromValues(ain[0], ain[1], ain[2]);
    const vout = vec3.create();
    vec3.transformMat4(vout, vin, model.indexToWorld);
    vec3.copy(aout, vout);
    return aout;
  };

  // this is the fast version, requires vec3 arguments
  publicAPI.worldToIndexVec3 = (vin, vout) => {
    vec3.transformMat4(vout, vin, model.worldToIndex);
    return vout;
  };

  // slow version for generic arrays
  publicAPI.worldToIndex = (ain, aout = []) => {
    const vin = vec3.fromValues(ain[0], ain[1], ain[2]);
    const vout = vec3.create();
    vec3.transformMat4(vout, vin, model.worldToIndex);
    vec3.copy(aout, vout);
    return aout;
  };

  publicAPI.indexToWorldBounds = (bin, bout = []) => {
    const in1 = [0, 0, 0];
    const in2 = [0, 0, 0];
    vtkBoundingBox.computeCornerPoints(in1, in2, bin);

    const out1 = [0, 0, 0];
    const out2 = [0, 0, 0];
    vec3.transformMat4(out1, in1, model.indexToWorld);
    vec3.transformMat4(out2, in2, model.indexToWorld);

    vtkMath.computeBoundsFromPoints(out1, out2, bout);

    return bout;
  };

  publicAPI.worldToIndexBounds = (bin, bout = []) => {
    const in1 = [0, 0, 0];
    const in2 = [0, 0, 0];
    vtkBoundingBox.computeCornerPoints(in1, in2, bin);

    const out1 = [0, 0, 0];
    const out2 = [0, 0, 0];
    vec3.transformMat4(out1, in1, model.worldToIndex);
    vec3.transformMat4(out2, in2, model.worldToIndex);

    vtkMath.computeBoundsFromPoints(out1, out2, bout);

    return bout;
  };

  // Make sure the transform is correct
  publicAPI.onModified(publicAPI.computeTransforms);
  publicAPI.computeTransforms();

  publicAPI.getCenter = () => {
    const bounds = publicAPI.getBounds();
    const center = [];

    for (let i = 0; i < 3; i++) {
      center[i] = (bounds[2 * i + 1] + bounds[2 * i]) / 2;
    }

    return center;
  };

  publicAPI.computeHistogram = (worldBounds, voxelFunc = null) => {
    const bounds = [0, 0, 0, 0, 0, 0];
    publicAPI.worldToIndexBounds(worldBounds, bounds);

    const point1 = [0, 0, 0];
    const point2 = [0, 0, 0];
    vtkBoundingBox.computeCornerPoints(point1, point2, bounds);

    vtkMath.roundVector(point1, point1);
    vtkMath.roundVector(point2, point2);

    const dimensions = publicAPI.getDimensions();

    vtkMath.clampVector(
      point1,
      [0, 0, 0],
      [dimensions[0] - 1, dimensions[1] - 1, dimensions[2] - 1],
      point1
    );
    vtkMath.clampVector(
      point2,
      [0, 0, 0],
      [dimensions[0] - 1, dimensions[1] - 1, dimensions[2] - 1],
      point2
    );

    const yStride = dimensions[0];
    const zStride = dimensions[0] * dimensions[1];

    const pixels = publicAPI
      .getPointData()
      .getScalars()
      .getData();

    let maximum = -Infinity;
    let minimum = Infinity;
    let sumOfSquares = 0;
    let isum = 0;
    let inum = 0;

    for (let z = point1[2]; z <= point2[2]; z++) {
      for (let y = point1[1]; y <= point2[1]; y++) {
        let index = point1[0] + y * yStride + z * zStride;
        for (let x = point1[0]; x <= point2[0]; x++) {
          if (!voxelFunc || voxelFunc([x, y, z], bounds)) {
            const pixel = pixels[index];

            if (pixel > maximum) maximum = pixel;
            if (pixel < minimum) minimum = pixel;
            sumOfSquares += pixel * pixel;
            isum += pixel;
            inum += 1;
          }

          ++index;
        }
      }
    }

    const average = inum > 0 ? isum / inum : 0;
    const variance = sumOfSquares - average * average;
    const sigma = Math.sqrt(variance);

    return {
      minimum,
      maximum,
      average,
      variance,
      sigma,
    };
  };

  // TODO: use the unimplemented `vtkDataSetAttributes` for scalar length, that is currently also a TODO (GetNumberOfComponents).
  // Scalar data could be tuples for color information?
  publicAPI.computeIncrements = (extent, numberOfComponents = 1) => {
    const increments = [];
    let incr = numberOfComponents;

    // Calculate array increment offsets
    // similar to c++ vtkImageData::ComputeIncrements
    for (let idx = 0; idx < 3; ++idx) {
      increments[idx] = incr;
      incr *= extent[idx * 2 + 1] - extent[idx * 2] + 1;
    }
    return increments;
  };

  /**
   * @param {Number[]} index the localized `[i,j,k]` pixel array position. Float values will be rounded.
   * @return {Number} the corresponding flattened index in the scalar array
   */
  publicAPI.computeOffsetIndex = ([i, j, k]) => {
    const extent = publicAPI.getExtent();
    const numberOfComponents = publicAPI
      .getPointData()
      .getScalars()
      .getNumberOfComponents();
    const increments = publicAPI.computeIncrements(extent, numberOfComponents);
    // Use the array increments to find the pixel index
    // similar to c++ vtkImageData::GetArrayPointer
    // Math.floor to catch "practically 0" e^-15 scenarios.
    return Math.floor(
      (Math.round(i) - extent[0]) * increments[0] +
        (Math.round(j) - extent[2]) * increments[1] +
        (Math.round(k) - extent[4]) * increments[2]
    );
  };

  /**
   * @param {Number[]} xyz the [x,y,z] Array in world coordinates
   * @return {Number|NaN} the corresponding pixel's index in the scalar array
   */
  publicAPI.getOffsetIndexFromWorld = (xyz) => {
    const extent = publicAPI.getExtent();
    const index = publicAPI.worldToIndex(xyz);

    // Confirm indexed i,j,k coords are within the bounds of the volume
    for (let idx = 0; idx < 3; ++idx) {
      if (index[idx] < extent[idx * 2] || index[idx] > extent[idx * 2 + 1]) {
        vtkErrorMacro(
          `GetScalarPointer: Pixel ${index} is not in memory. Current extent = ${extent}`
        );
        return NaN;
      }
    }

    // Assumed the index here is within 0 <-> scalarData.length, but doesn't hurt to check upstream
    return publicAPI.computeOffsetIndex(index);
  };
  /**
   * @param {Number[]} xyz the [x,y,z] Array in world coordinates
   * @param {Number?} comp the scalar component index for multi-component scalars
   * @return {Number|NaN} the corresponding pixel's scalar value
   */
  publicAPI.getScalarValueFromWorld = (xyz, comp = 0) => {
    const numberOfComponents = publicAPI
      .getPointData()
      .getScalars()
      .getNumberOfComponents();
    if (comp < 0 || comp >= numberOfComponents) {
      vtkErrorMacro(
        `GetScalarPointer: Scalar Component ${comp} is not within bounds. Current Scalar numberOfComponents: ${numberOfComponents}`
      );
      return NaN;
    }
    const offsetIndex = publicAPI.getOffsetIndexFromWorld(xyz);
    if (Number.isNaN(offsetIndex)) {
      // VTK Error Macro will have been tripped already, no need to do it again,
      return offsetIndex;
    }

    return publicAPI
      .getPointData()
      .getScalars()
      .getComponent(offsetIndex, comp);
  };
}

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

const DEFAULT_VALUES = {
  direction: null, // a mat3
  indexToWorld: null, // a mat4
  worldToIndex: null, // a mat4
  spacing: [1.0, 1.0, 1.0],
  origin: [0.0, 0.0, 0.0],
  extent: [0, -1, 0, -1, 0, -1],
  dataDescription: StructuredType.EMPTY,
};

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

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

  // Inheritance
  vtkDataSet.extend(publicAPI, model, initialValues);

  if (!model.direction) {
    model.direction = mat3.create();
  } else if (Array.isArray(model.direction)) {
    const dvals = model.direction.slice(0);
    model.direction = mat3.create();
    for (let i = 0; i < 9; ++i) {
      model.direction[i] = dvals[i];
    }
  }

  model.indexToWorld = new Float64Array(16);
  model.worldToIndex = new Float64Array(16);

  // Set/Get methods
  macro.get(publicAPI, model, ['direction', 'indexToWorld', 'worldToIndex']);
  macro.setGetArray(publicAPI, model, ['origin', 'spacing'], 3);
  macro.getArray(publicAPI, model, ['extent'], 6);

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

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

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

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

export default { newInstance, extend };