ImageReslice

Introduction

vtkImageReslice - Reslices a volume along a new set of axes

vtkImageReslice is the swiss-army-knife of image geometry filters:
It can permute, rotate, flip, scale, resample, deform, and pad image
data in any combination with reasonably high efficiency. Simple
operations such as permutation, resampling and padding are done
with similar efficiently to the specialized vtkImagePermute,
vtkImageResample, and vtkImagePad filters. There are a number of
tasks that vtkImageReslice is well suited for:

1. Application of simple rotations, scales, and translations to
   an image. It is often a good idea to use vtkImageChangeInformation
   to center the image first, so that scales and rotations occur around
   the center rather than around the lower-left corner of the image.
2. Extraction of slices from an image volume. The method
   SetOutputDimensionality(2) is used to specify that want to output a
   slice rather than a volume. You can use both the resliceAxes and the
   resliceTransform at the same time, in order to extract slices from a
   volume that you have applied a transformation to.

Usage

Provide the input to the filter via the standard
SetInput(Data/Connection) methods.

const imageReslice = vtkImageReslice.newInstance();
imageReslice.setInputData(imageData);
imageReslice.setOutputDimensionality(2);
const axes = mat4.identity(new Float64Array(16));
mat4.rotateX(axes, axes, 45 * Math.PI / 180);
imageReslice.setResliceAxes(axes);
imageReslice.setOutputScalarType('Uint16Array');
imageReslice.setScalarScale(65535 / 255);

const obliqueSlice = imageReslice.getOutputData();

Public API

ResliceAxes (set/get)

This method is used to set up the axes for the output voxels.
The output Spacing, Origin, and Extent specify the locations
of the voxels within the coordinate system defined by the axes.
The ResliceAxes are used most often to permute the data, e.g.
to extract ZY or XZ slices of a volume as 2D XY images.
The first column of the matrix specifies the x-axis
vector (the fourth element must be set to zero), the second
column specifies the y-axis, and the third column the
z-axis. The fourth column is the origin of the
axes (the fourth element must be set to one).

OutputDimensionality (set/get)

Force the dimensionality of the output to either 1, 2,
3 or 0 (default: 3). If the dimensionality is 2D, then
the Z extent of the output is forced to (0,0) and the Z
origin of the output is forced to 0.0 (i.e. the output
extent is confined to the xy plane). If the dimensionality
is 1D, the output extent is confined to the x axis.
For 0D, the output extent consists of a single voxel at
(0,0,0).

OutputOrigin (set/get)

Set the origin for the output data. The default output origin
is the input origin permuted through the ResliceAxes.

OutputSpacing (set/get)

Set the voxel spacing for the output data. The default output
spacing is the input spacing permuted through the ResliceAxes.

OutputExtent (set/get)

Set the extent for the output data. The default output extent
is the input extent permuted through the ResliceAxes.

OutputScalarType (set/get)

Set the scalar type of the output to be different from the input.
The default value is null, which means that the input scalar type will be
used to set the output scalar type. Otherwise, this must be set to one
of the following types: VtkDataTypes.CHAR, VtkDataTypes.SIGNED_CHAR,
VtkDataTypes.UNSIGNED_CHAR, VtkDataTypes.SHORT, VtkDataTypes.UNSIGNED_SHORT,
VtkDataTypes.INT, VtkDataTypes.UNSIGNED_INT, VtkDataTypes.FLOAT or
VtkDataTypes.DOUBLE. Other types are not permitted. If the output type
is an integer type, the output will be rounded and clamped to the limits of
the type.

See the documentation for vtkDataArray::getDataType() for additional data type settings.

ScalarShift (set/get)

Set a value to add to all the output voxels.
After a sample value has been interpolated from the input image, the
equation u = (v + ScalarShift)*ScalarScale will be applied to it before
it is written to the output image. The result will always be clamped to
the limits of the output data type.

ScalarScale (set/get)

Set multiplication factor to apply to all the output voxels.
After a sample value has been interpolated from the input image, the
equation u = (v + ScalarShift)*ScalarScale will be applied to it before
it is written to the output image. The result will always be clamped to
the limits of the output data type.

Wrap (set/get)

Turn on wrap-pad feature (default: false).

Mirror (set/get)

Turn on mirror-pad feature (default: false). This will override the wrap-pad.

Border (set/get)

Extend the apparent input border by a half voxel (default: On).
This changes how interpolation is handled at the borders of the
input image: if the center of an output voxel is beyond the edge
of the input image, but is within a half voxel width of the edge
(using the input voxel width), then the value of the output voxel
is calculated as if the input’s edge voxels were duplicated past
the edges of the input.
This has no effect if Mirror or Wrap are on.

BackgroundColor (set/get)

Set the background color (for multi-component images).

TransformInputSampling (set/get)

Specify whether to transform the spacing, origin and extent
of the Input (or the InformationInput) according to the
direction cosines and origin of the ResliceAxes before applying
them as the default output spacing, origin and extent
(default: true).

AutoCropOutput (set/get)

Turn this on if you want to guarantee that the extent of the
output will be large enough to ensure that none of the
data will be cropped (default: false).

Source

Constants.js

export const SlabMode = {
  MIN: 0,
  MAX: 1,
  MEAN: 2,
  SUM: 3,
};

export default {
  SlabMode,
};

index.js

import { vec4, mat4 } from 'gl-matrix';

import macro, { vtkWarningMacro } from 'vtk.js/Sources/macros';
import vtkDataArray from 'vtk.js/Sources/Common/Core/DataArray';
import vtkMath from 'vtk.js/Sources/Common/Core/Math';
import vtkMatrixBuilder from 'vtk.js/Sources/Common/Core/MatrixBuilder';
import { VtkDataTypes } from 'vtk.js/Sources/Common/Core/DataArray/Constants';
import vtkBoundingBox from 'vtk.js/Sources/Common/DataModel/BoundingBox';
import vtkImageData from 'vtk.js/Sources/Common/DataModel/ImageData';
import vtkImageInterpolator from 'vtk.js/Sources/Imaging/Core/ImageInterpolator';
import vtkImagePointDataIterator from 'vtk.js/Sources/Imaging/Core/ImagePointDataIterator';
import {
  ImageBorderMode,
  InterpolationMode,
} from 'vtk.js/Sources/Imaging/Core/AbstractImageInterpolator/Constants';
import {
  vtkInterpolationMathFloor,
  vtkInterpolationMathRound,
  vtkInterpolationMathClamp,
} from 'vtk.js/Sources/Imaging/Core/AbstractImageInterpolator/InterpolationInfo';
import Constants from 'vtk.js/Sources/Imaging/Core/ImageReslice/Constants';

const { SlabMode } = Constants;

const { vtkErrorMacro } = macro;

// ----------------------------------------------------------------------------
// vtkImageReslice methods
// ----------------------------------------------------------------------------

function vtkImageReslice(publicAPI, model) {
  // Set our className
  model.classHierarchy.push('vtkImageReslice');
  const superClass = { ...publicAPI };

  let indexMatrix = null;
  let optimizedTransform = null;

  function getImageResliceSlabTrap(tmpPtr, inComponents, sampleCount, f) {
    const n = sampleCount - 1;
    for (let i = 0; i < inComponents; i += 1) {
      let result = tmpPtr[i] * 0.5;
      for (let j = 1; j < n; j += 1) {
        result += tmpPtr[i + j * inComponents];
      }
      result += tmpPtr[i + n * inComponents] * 0.5;
      tmpPtr[i] = result * f;
    }
  }

  function getImageResliceSlabSum(tmpPtr, inComponents, sampleCount, f) {
    for (let i = 0; i < inComponents; i += 1) {
      let result = tmpPtr[i];
      for (let j = 1; j < sampleCount; j += 1) {
        result += tmpPtr[i + j * inComponents];
      }
      tmpPtr[i] = result * f;
    }
  }

  function getImageResliceCompositeMinValue(tmpPtr, inComponents, sampleCount) {
    for (let i = 0; i < inComponents; i += 1) {
      let result = tmpPtr[i];
      for (let j = 1; j < sampleCount; j += 1) {
        result = Math.min(result, tmpPtr[i + j * inComponents]);
      }
      tmpPtr[i] = result;
    }
  }

  function getImageResliceCompositeMaxValue(tmpPtr, inComponents, sampleCount) {
    for (let i = 0; i < inComponents; i += 1) {
      let result = tmpPtr[i];
      for (let j = 1; j < sampleCount; j += 1) {
        result = Math.max(result, tmpPtr[i + j * inComponents]);
      }
      tmpPtr[i] = result;
    }
  }

  function getImageResliceCompositeMeanValue(
    tmpPtr,
    inComponents,
    sampleCount
  ) {
    const f = 1.0 / sampleCount;
    getImageResliceSlabSum(tmpPtr, inComponents, sampleCount, f);
  }

  function getImageResliceCompositeMeanTrap(tmpPtr, inComponents, sampleCount) {
    const f = 1.0 / (sampleCount - 1);
    getImageResliceSlabTrap(tmpPtr, inComponents, sampleCount, f);
  }

  function getImageResliceCompositeSumValue(tmpPtr, inComponents, sampleCount) {
    const f = 1.0;
    getImageResliceSlabSum(tmpPtr, inComponents, sampleCount, f);
  }

  function getImageResliceCompositeSumTrap(tmpPtr, inComponents, sampleCount) {
    const f = 1.0;
    getImageResliceSlabTrap(tmpPtr, inComponents, sampleCount, f);
  }

  publicAPI.getMTime = () => {
    let mTime = superClass.getMTime();
    if (model.resliceTransform) {
      mTime = Math.max(mTime, model.resliceTransform.getMTime());
    }
    return mTime;
  };

  publicAPI.setResliceAxes = (resliceAxes) => {
    if (!model.resliceAxes) {
      model.resliceAxes = mat4.identity(new Float64Array(16));
    }

    if (!mat4.exactEquals(model.resliceAxes, resliceAxes)) {
      mat4.copy(model.resliceAxes, resliceAxes);

      publicAPI.modified();
      return true;
    }
    return null;
  };

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

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

    // console.time('reslice');

    // Retrieve output and volume data
    const origin = input.getOrigin();
    const inSpacing = input.getSpacing();
    const dims = input.getDimensions();
    const inScalars = input.getPointData().getScalars();
    const inWholeExt = [0, dims[0] - 1, 0, dims[1] - 1, 0, dims[2] - 1];

    const outOrigin = [0, 0, 0];
    const outSpacing = [1, 1, 1];
    const outWholeExt = [0, 0, 0, 0, 0, 0];
    const outDims = [0, 0, 0];

    const matrix = mat4.identity(new Float64Array(16));
    if (model.resliceAxes) {
      mat4.multiply(matrix, matrix, model.resliceAxes);
    }
    const imatrix = new Float64Array(16);
    mat4.invert(imatrix, matrix);

    const inCenter = [
      origin[0] + 0.5 * (inWholeExt[0] + inWholeExt[1]) * inSpacing[0],
      origin[1] + 0.5 * (inWholeExt[2] + inWholeExt[3]) * inSpacing[1],
      origin[2] + 0.5 * (inWholeExt[4] + inWholeExt[5]) * inSpacing[2],
    ];

    let maxBounds = null;
    if (model.autoCropOutput) {
      maxBounds = publicAPI.getAutoCroppedOutputBounds(input);
    }

    for (let i = 0; i < 3; i++) {
      let s = 0; // default output spacing
      let d = 0; // default linear dimension
      let e = 0; // default extent start
      let c = 0; // transformed center-of-volume

      if (model.transformInputSampling) {
        let r = 0.0;
        for (let j = 0; j < 3; j++) {
          c += imatrix[4 * j + i] * (inCenter[j] - matrix[4 * 3 + j]);
          const tmp = matrix[4 * i + j] * matrix[4 * i + j];
          s += tmp * Math.abs(inSpacing[j]);
          d +=
            tmp *
            (inWholeExt[2 * j + 1] - inWholeExt[2 * j]) *
            Math.abs(inSpacing[j]);
          e += tmp * inWholeExt[2 * j];
          r += tmp;
        }
        s /= r;
        d /= r * Math.sqrt(r);
        e /= r;
      } else {
        c = inCenter[i];
        s = inSpacing[i];
        d = (inWholeExt[2 * i + 1] - inWholeExt[2 * i]) * s;
        e = inWholeExt[2 * i];
      }

      if (model.outputSpacing == null) {
        outSpacing[i] = s;
      } else {
        outSpacing[i] = model.outputSpacing[i];
      }

      if (i >= model.outputDimensionality) {
        outWholeExt[2 * i] = 0;
        outWholeExt[2 * i + 1] = 0;
      } else if (model.outputExtent == null) {
        if (model.autoCropOutput) {
          d = maxBounds[2 * i + 1] - maxBounds[2 * i];
        }
        outWholeExt[2 * i] = Math.round(e);
        outWholeExt[2 * i + 1] = Math.round(
          outWholeExt[2 * i] + Math.abs(d / outSpacing[i])
        );
      } else {
        outWholeExt[2 * i] = model.outputExtent[2 * i];
        outWholeExt[2 * i + 1] = model.outputExtent[2 * i + 1];
      }

      if (i >= model.outputDimensionality) {
        outOrigin[i] = 0;
      } else if (model.outputOrigin == null) {
        if (model.autoCropOutput) {
          // set origin so edge of extent is edge of bounds
          outOrigin[i] = maxBounds[2 * i] - outWholeExt[2 * i] * outSpacing[i];
        } else {
          // center new bounds over center of input bounds
          outOrigin[i] =
            c -
            0.5 * (outWholeExt[2 * i] + outWholeExt[2 * i + 1]) * outSpacing[i];
        }
      } else {
        outOrigin[i] = model.outputOrigin[i];
      }
      outDims[i] = outWholeExt[2 * i + 1] - outWholeExt[2 * i] + 1;
    }

    let dataType = inScalars.getDataType();
    if (model.outputScalarType) {
      dataType = model.outputScalarType;
    }

    const numComponents = input
      .getPointData()
      .getScalars()
      .getNumberOfComponents(); // or s.numberOfComponents;

    const outScalarsData = macro.newTypedArray(
      dataType,
      outDims[0] * outDims[1] * outDims[2] * numComponents
    );
    const outScalars = vtkDataArray.newInstance({
      name: 'Scalars',
      values: outScalarsData,
      numberOfComponents: numComponents,
    });

    // Update output
    const output = vtkImageData.newInstance();
    output.setDimensions(outDims);
    output.setOrigin(outOrigin);
    output.setSpacing(outSpacing);
    if (model.outputDirection) {
      output.setDirection(model.outputDirection);
    }
    output.getPointData().setScalars(outScalars);

    publicAPI.getIndexMatrix(input, output);

    let interpolationMode = model.interpolationMode;
    model.usePermuteExecute = false;
    if (model.optimization) {
      if (
        optimizedTransform == null &&
        model.slabSliceSpacingFraction === 1.0 &&
        model.interpolator.isSeparable() &&
        publicAPI.isPermutationMatrix(indexMatrix)
      ) {
        model.usePermuteExecute = true;
        if (publicAPI.canUseNearestNeighbor(indexMatrix, outWholeExt)) {
          interpolationMode = InterpolationMode.NEAREST;
        }
      }
    }
    model.interpolator.setInterpolationMode(interpolationMode);

    let borderMode = ImageBorderMode.CLAMP;
    borderMode = model.wrap ? ImageBorderMode.REPEAT : borderMode;
    borderMode = model.mirror ? ImageBorderMode.MIRROR : borderMode;
    model.interpolator.setBorderMode(borderMode);

    const mintol = 7.62939453125e-6;
    const maxtol = 2.0 * 2147483647;
    let tol = 0.5 * model.border;
    tol = borderMode === ImageBorderMode.CLAMP ? tol : maxtol;
    tol = tol > mintol ? tol : mintol;
    model.interpolator.setTolerance(tol);

    model.interpolator.initialize(input);

    publicAPI.vtkImageResliceExecute(input, output);

    model.interpolator.releaseData();

    outData[0] = output;

    // console.timeEnd('reslice');
  };

  publicAPI.vtkImageResliceExecute = (input, output) => {
    // const outDims = output.getDimensions();
    const inScalars = input.getPointData().getScalars();
    const outScalars = output.getPointData().getScalars();
    let outPtr = outScalars.getData();
    const outExt = output.getExtent();
    const newmat = indexMatrix;
    const outputStencil = null;

    // multiple samples for thick slabs
    const nsamples = Math.max(model.slabNumberOfSlices, 1);

    // spacing between slab samples (as a fraction of slice spacing).
    const slabSampleSpacing = model.slabSliceSpacingFraction;

    // check for perspective transformation
    const perspective = publicAPI.isPerspectiveMatrix(newmat);

    // extra scalar info for nearest-neighbor optimization
    let inPtr = inScalars.getData();
    const inputScalarSize = 1; // inScalars.getElementComponentSize(); // inScalars.getDataTypeSize();
    const inputScalarType = inScalars.getDataType();
    const inComponents = inScalars.getNumberOfComponents(); // interpolator.GetNumberOfComponents();
    const componentOffset = model.interpolator.getComponentOffset();
    const borderMode = model.interpolator.getBorderMode();
    const inDims = input.getDimensions();
    const inExt = [0, inDims[0] - 1, 0, inDims[1] - 1, 0, inDims[2] - 1]; // interpolator->GetExtent();
    const inInc = [0, 0, 0];
    inInc[0] = inScalars.getNumberOfComponents();
    inInc[1] = inInc[0] * inDims[0];
    inInc[2] = inInc[1] * inDims[1];

    const fullSize = inDims[0] * inDims[1] * inDims[2];
    if (componentOffset > 0 && componentOffset + inComponents < inInc[0]) {
      inPtr = inPtr.subarray(inputScalarSize * componentOffset);
    }

    let interpolationMode = InterpolationMode.NEAREST;
    if (model.interpolator.isA('vtkImageInterpolator')) {
      interpolationMode = model.interpolator.getInterpolationMode();
    }

    const convertScalars = null;
    const rescaleScalars =
      model.scalarShift !== 0.0 || model.scalarScale !== 1.0;

    // is nearest neighbor optimization possible?
    const optimizeNearest =
      interpolationMode === InterpolationMode.NEAREST &&
      borderMode === ImageBorderMode.CLAMP &&
      !(
        optimizedTransform != null ||
        perspective ||
        convertScalars != null ||
        rescaleScalars
      ) &&
      inputScalarType === outScalars.getDataType() &&
      fullSize === inScalars.getNumberOfTuples() &&
      model.border === true &&
      nsamples <= 1;

    // get pixel information
    const scalarType = outScalars.getDataType();
    const scalarSize = 1; // outScalars.getElementComponentSize() // outScalars.scalarSize;
    const outComponents = outScalars.getNumberOfComponents();

    // break matrix into a set of axes plus an origin
    // (this allows us to calculate the transform Incrementally)
    const xAxis = [0, 0, 0, 0];
    const yAxis = [0, 0, 0, 0];
    const zAxis = [0, 0, 0, 0];
    const origin = [0, 0, 0, 0];
    for (let i = 0; i < 4; ++i) {
      xAxis[i] = newmat[4 * 0 + i];
      yAxis[i] = newmat[4 * 1 + i];
      zAxis[i] = newmat[4 * 2 + i];
      origin[i] = newmat[4 * 3 + i];
    }

    // get the input origin and spacing for conversion purposes
    const inOrigin = model.interpolator.getOrigin();
    const inSpacing = model.interpolator.getSpacing();
    const inInvSpacing = [
      1.0 / inSpacing[0],
      1.0 / inSpacing[1],
      1.0 / inSpacing[2],
    ];

    // allocate an output row of type double
    let floatPtr = null;
    if (!optimizeNearest) {
      floatPtr = new Float64Array(
        inComponents * (outExt[1] - outExt[0] + nsamples)
      );
    }

    const background = macro.newTypedArray(
      inputScalarType,
      model.backgroundColor
    );

    // set color for area outside of input volume extent
    // void *background;
    // vtkAllocBackgroundPixel(&background,
    //    self->GetBackgroundColor(), scalarType, scalarSize, outComponents);

    // get various helper functions
    const forceClamping =
      interpolationMode > InterpolationMode.LINEAR ||
      (nsamples > 1 && model.slabMode === SlabMode.SUM);
    const convertpixels = publicAPI.getConversionFunc(
      inputScalarType,
      scalarType,
      model.scalarShift,
      model.scalarScale,
      forceClamping
    );
    const setpixels = publicAPI.getSetPixelsFunc(
      scalarType,
      scalarSize,
      outComponents,
      outPtr
    );
    const composite = publicAPI.getCompositeFunc(
      model.slabMode,
      model.slabTrapezoidIntegration
    );

    // create some variables for when we march through the data
    let idY = outExt[2] - 1;
    let idZ = outExt[4] - 1;
    const inPoint0 = [0.0, 0.0, 0.0, 0.0];
    const inPoint1 = [0.0, 0.0, 0.0, 0.0];

    // create an iterator to march through the data
    const iter = vtkImagePointDataIterator.newInstance();
    iter.initialize(output, outExt, model.stencil, null);
    const outPtr0 = iter.getScalars(output, 0);
    let outPtrIndex = 0;
    const outTmp = macro.newTypedArray(
      scalarType,
      vtkBoundingBox.getDiagonalLength(outExt) * outComponents * 2
    );

    const interpolatedPtr = new Float64Array(inComponents * nsamples);
    const interpolatedPoint = new Float64Array(inComponents);

    for (; !iter.isAtEnd(); iter.nextSpan()) {
      const span = iter.spanEndId() - iter.getId();
      outPtrIndex = iter.getId() * scalarSize * outComponents;

      if (!iter.isInStencil()) {
        // clear any regions that are outside the stencil
        const n = setpixels(outTmp, background, outComponents, span);
        for (let i = 0; i < n; ++i) {
          outPtr0[outPtrIndex++] = outTmp[i];
        }
      } else {
        // get output index, and compute position in input image
        const outIndex = iter.getIndex();

        // if Z index increased, then advance position along Z axis
        if (outIndex[2] > idZ) {
          idZ = outIndex[2];
          inPoint0[0] = origin[0] + idZ * zAxis[0];
          inPoint0[1] = origin[1] + idZ * zAxis[1];
          inPoint0[2] = origin[2] + idZ * zAxis[2];
          inPoint0[3] = origin[3] + idZ * zAxis[3];
          idY = outExt[2] - 1;
        }

        // if Y index increased, then advance position along Y axis
        if (outIndex[1] > idY) {
          idY = outIndex[1];
          inPoint1[0] = inPoint0[0] + idY * yAxis[0];
          inPoint1[1] = inPoint0[1] + idY * yAxis[1];
          inPoint1[2] = inPoint0[2] + idY * yAxis[2];
          inPoint1[3] = inPoint0[3] + idY * yAxis[3];
        }

        // march through one row of the output image
        const idXmin = outIndex[0];
        const idXmax = idXmin + span - 1;

        if (!optimizeNearest) {
          let wasInBounds = 1;
          let isInBounds = 1;
          let startIdX = idXmin;
          let idX = idXmin;
          const tmpPtr = floatPtr;
          let pixelIndex = 0;

          while (startIdX <= idXmax) {
            for (; idX <= idXmax && isInBounds === wasInBounds; idX++) {
              const inPoint2 = [
                inPoint1[0] + idX * xAxis[0],
                inPoint1[1] + idX * xAxis[1],
                inPoint1[2] + idX * xAxis[2],
                inPoint1[3] + idX * xAxis[3],
              ];

              const inPoint3 = [0, 0, 0, 0];
              let inPoint = inPoint2;
              isInBounds = false;

              let interpolatedPtrIndex = 0;
              for (let sample = 0; sample < nsamples; ++sample) {
                if (nsamples > 1) {
                  let s = sample - 0.5 * (nsamples - 1);
                  s *= slabSampleSpacing;
                  inPoint3[0] = inPoint2[0] + s * zAxis[0];
                  inPoint3[1] = inPoint2[1] + s * zAxis[1];
                  inPoint3[2] = inPoint2[2] + s * zAxis[2];
                  inPoint3[3] = inPoint2[3] + s * zAxis[3];
                  inPoint = inPoint3;
                }

                if (perspective) {
                  // only do perspective if necessary
                  const f = 1 / inPoint[3];
                  inPoint[0] *= f;
                  inPoint[1] *= f;
                  inPoint[2] *= f;
                }

                if (optimizedTransform !== null) {
                  // apply the AbstractTransform if there is one
                  publicAPI.applyTransform(
                    optimizedTransform,
                    inPoint,
                    inOrigin,
                    inInvSpacing
                  );
                }

                if (model.interpolator.checkBoundsIJK(inPoint)) {
                  // do the interpolation
                  isInBounds = 1;
                  model.interpolator.interpolateIJK(inPoint, interpolatedPoint);
                  for (let i = 0; i < inComponents; ++i) {
                    interpolatedPtr[interpolatedPtrIndex++] =
                      interpolatedPoint[i];
                  }
                }
              }

              if (interpolatedPtrIndex > inComponents) {
                composite(
                  interpolatedPtr,
                  inComponents,
                  interpolatedPtrIndex / inComponents
                );
              }
              for (let i = 0; i < inComponents; ++i) {
                tmpPtr[pixelIndex++] = interpolatedPtr[i];
              }

              // set "was in" to "is in" if first pixel
              wasInBounds = idX > idXmin ? wasInBounds : isInBounds;
            }

            // write a segment to the output
            const endIdX = idX - 1 - (isInBounds !== wasInBounds);
            const numpixels = endIdX - startIdX + 1;

            let n = 0;
            if (wasInBounds) {
              if (outputStencil) {
                outputStencil.insertNextExtent(startIdX, endIdX, idY, idZ);
              }

              if (rescaleScalars) {
                publicAPI.rescaleScalars(
                  floatPtr,
                  inComponents,
                  idXmax - idXmin + 1,
                  model.scalarShift,
                  model.scalarScale
                );
              }

              if (convertScalars) {
                convertScalars(
                  floatPtr.subarray(startIdX * inComponents),
                  outTmp,
                  inputScalarType,
                  inComponents,
                  numpixels,
                  startIdX,
                  idY,
                  idZ
                );
                n = numpixels * outComponents * scalarSize;
              } else {
                n = convertpixels(
                  outTmp,
                  floatPtr.subarray(startIdX * inComponents),
                  outComponents,
                  numpixels
                );
              }
            } else {
              n = setpixels(outTmp, background, outComponents, numpixels);
            }
            for (let i = 0; i < n; ++i) {
              outPtr0[outPtrIndex++] = outTmp[i];
            }

            startIdX += numpixels;
            wasInBounds = isInBounds;
          }
        } else {
          // optimize for nearest-neighbor interpolation
          const inPtrTmp0 = inPtr;
          const outPtrTmp = outPtr;

          const inIncX = inInc[0] * inputScalarSize;
          const inIncY = inInc[1] * inputScalarSize;
          const inIncZ = inInc[2] * inputScalarSize;

          const inExtX = inExt[1] - inExt[0] + 1;
          const inExtY = inExt[3] - inExt[2] + 1;
          const inExtZ = inExt[5] - inExt[4] + 1;

          let startIdX = idXmin;
          let endIdX = idXmin - 1;
          let isInBounds = false;
          const bytesPerPixel = inputScalarSize * inComponents;

          for (let iidX = idXmin; iidX <= idXmax; iidX++) {
            const inPoint = [
              inPoint1[0] + iidX * xAxis[0],
              inPoint1[1] + iidX * xAxis[1],
              inPoint1[2] + iidX * xAxis[2],
            ];

            const inIdX = vtkInterpolationMathRound(inPoint[0]) - inExt[0];
            const inIdY = vtkInterpolationMathRound(inPoint[1]) - inExt[2];
            const inIdZ = vtkInterpolationMathRound(inPoint[2]) - inExt[4];

            if (
              inIdX >= 0 &&
              inIdX < inExtX &&
              inIdY >= 0 &&
              inIdY < inExtY &&
              inIdZ >= 0 &&
              inIdZ < inExtZ
            ) {
              if (!isInBounds) {
                // clear leading out-of-bounds pixels
                startIdX = iidX;
                isInBounds = true;
                const n = setpixels(
                  outTmp,
                  background,
                  outComponents,
                  startIdX - idXmin
                );
                for (let i = 0; i < n; ++i) {
                  outPtr0[outPtrIndex++] = outTmp[i];
                }
              }
              // set the final index that was within input bounds
              endIdX = iidX;

              // perform nearest-neighbor interpolation via pixel copy
              let offset = inIdX * inIncX + inIdY * inIncY + inIdZ * inIncZ;

              // when memcpy is used with a constant size, the compiler will
              // optimize away the function call and use the minimum number
              // of instructions necessary to perform the copy
              switch (bytesPerPixel) {
                case 1:
                  outPtr0[outPtrIndex++] = inPtrTmp0[offset];
                  break;
                case 2:
                case 3:
                case 4:
                case 8:
                case 12:
                case 16:
                  for (let i = 0; i < bytesPerPixel; ++i) {
                    outPtr0[outPtrIndex++] = inPtrTmp0[offset + i];
                  }
                  break;
                default: {
                  // TODO: check bytes
                  let oc = 0;
                  do {
                    outPtr0[outPtrIndex++] = inPtrTmp0[offset++];
                  } while (++oc !== bytesPerPixel);
                  break;
                }
              }
            } else if (isInBounds) {
              // leaving input bounds
              break;
            }
          }

          // clear trailing out-of-bounds pixels
          outPtr = outPtrTmp;
          const n = setpixels(
            outTmp,
            background,
            outComponents,
            idXmax - endIdX
          );
          for (let i = 0; i < n; ++i) {
            outPtr0[outPtrIndex++] = outTmp[i];
          }

          if (outputStencil && endIdX >= startIdX) {
            outputStencil.insertNextExtent(startIdX, endIdX, idY, idZ);
          }
        }
      }
    }
  };
  /**
   * The transform matrix supplied by the user converts output coordinates
   * to input coordinates.
   * To speed up the pixel lookup, the following function provides a
   * matrix which converts output pixel indices to input pixel indices.
   * This will also concatenate the ResliceAxes and the ResliceTransform
   * if possible (if the ResliceTransform is a 4x4 matrix transform).
   * If it does, this->OptimizedTransform will be set to nullptr, otherwise
   * this->OptimizedTransform will be equal to this->ResliceTransform.
   * @param {vtkPolyData} input
   * @param {vtkPolyData} output
   * @returns
   */
  publicAPI.getIndexMatrix = (input, output) => {
    // first verify that we have to update the matrix
    if (indexMatrix === null) {
      indexMatrix = mat4.identity(new Float64Array(16));
    }

    const inOrigin = input.getOrigin();
    const inSpacing = input.getSpacing();
    const inDirection = input.getDirection();
    const outOrigin = output.getOrigin();
    const outSpacing = output.getSpacing();

    const transform = mat4.identity(new Float64Array(16));
    const inMatrix = mat4.identity(new Float64Array(16));
    const outMatrix = mat4.identity(new Float64Array(16));

    if (optimizedTransform) {
      optimizedTransform = null;
    }

    if (model.resliceAxes) {
      mat4.copy(transform, model.resliceAxes);
    }
    if (model.resliceTransform) {
      if (model.resliceTransform.isA('vtkHomogeneousTransform')) {
        // transform->PostMultiply();
        // transform->Concatenate(
        // mat4.multiply(transform, transform, model.resliceTransform.getMatrix());
        mat4.multiply(transform, model.resliceTransform.getMatrix(), transform);
      } else {
        // TODO
        vtkWarningMacro('Non homogeneous transform have not yet been ported');
      }
    }

    if (!vtkMath.isIdentity3x3(inDirection)) {
      const imageTransform = vtkMatrixBuilder
        .buildFromRadian()
        .translate(inOrigin[0], inOrigin[1], inOrigin[2])
        .multiply3x3(inDirection)
        .translate(-inOrigin[0], -inOrigin[1], -inOrigin[2]);
      mat4.multiply(transform, imageTransform.getMatrix(), transform);
    }

    // check to see if we have an identity matrix
    let isIdentity = vtkMath.isIdentity(transform);

    // the outMatrix takes OutputData indices to OutputData coordinates,
    // the inMatrix takes InputData coordinates to InputData indices
    for (let i = 0; i < 3; i++) {
      if (
        (optimizedTransform == null &&
          (inSpacing[i] !== outSpacing[i] || inOrigin[i] !== outOrigin[i])) ||
        (optimizedTransform != null &&
          (outSpacing[i] !== 1.0 || outOrigin[i] !== 0.0))
      ) {
        isIdentity = false;
      }
      inMatrix[4 * i + i] = 1.0 / inSpacing[i];
      inMatrix[4 * 3 + i] = -inOrigin[i] / inSpacing[i];
      outMatrix[4 * i + i] = outSpacing[i];
      outMatrix[4 * 3 + i] = outOrigin[i];
    }

    if (!isIdentity) {
      // transform.PreMultiply();
      // transform.Concatenate(outMatrix);
      mat4.multiply(transform, transform, outMatrix);

      // the optimizedTransform requires data coords, not
      // index coords, as its input
      if (optimizedTransform == null) {
        // transform->PostMultiply();
        // transform->Concatenate(inMatrix);
        mat4.multiply(transform, inMatrix, transform);
      }
    }

    mat4.copy(indexMatrix, transform);

    return indexMatrix;
  };

  publicAPI.getAutoCroppedOutputBounds = (input) => {
    const inOrigin = input.getOrigin();
    const inSpacing = input.getSpacing();
    const inDirection = input.getDirection();
    const dims = input.getDimensions();
    const inWholeExt = [0, dims[0] - 1, 0, dims[1] - 1, 0, dims[2] - 1];

    const matrix = new Float64Array(16);
    if (model.resliceAxes) {
      mat4.invert(matrix, model.resliceAxes);
    } else {
      mat4.identity(matrix);
    }
    let transform = null;
    if (model.resliceTransform) {
      transform = model.resliceTransform.getInverse();
    }
    let imageTransform = null;
    if (!vtkMath.isIdentity3x3(inDirection)) {
      imageTransform = vtkMatrixBuilder
        .buildFromRadian()
        .translate(inOrigin[0], inOrigin[1], inOrigin[2])
        .multiply3x3(inDirection)
        .translate(-inOrigin[0], -inOrigin[1], -inOrigin[2])
        .invert()
        .getMatrix();
    }

    const bounds = [
      Number.MAX_VALUE,
      -Number.MAX_VALUE,
      Number.MAX_VALUE,
      -Number.MAX_VALUE,
      Number.MAX_VALUE,
      -Number.MAX_VALUE,
    ];

    const point = [0, 0, 0, 0];
    for (let i = 0; i < 8; ++i) {
      point[0] = inOrigin[0] + inWholeExt[i % 2] * inSpacing[0];
      point[1] =
        inOrigin[1] + inWholeExt[2 + (Math.floor(i / 2) % 2)] * inSpacing[1];
      point[2] =
        inOrigin[2] + inWholeExt[4 + (Math.floor(i / 4) % 2)] * inSpacing[2];
      point[3] = 1.0;

      if (imageTransform) {
        vec4.transformMat4(point, point, imageTransform);
      }

      if (model.resliceTransform) {
        transform.transformPoint(point, point);
      }

      vec4.transformMat4(point, point, matrix);

      const f = 1.0 / point[3];
      point[0] *= f;
      point[1] *= f;
      point[2] *= f;

      for (let j = 0; j < 3; ++j) {
        if (point[j] > bounds[2 * j + 1]) {
          bounds[2 * j + 1] = point[j];
        }
        if (point[j] < bounds[2 * j]) {
          bounds[2 * j] = point[j];
        }
      }
    }
    return bounds;
  };

  publicAPI.getDataTypeMinMax = (dataType) => {
    switch (dataType) {
      case 'Int8Array':
        return { min: -128, max: 127 };
      case 'Int16Array':
        return { min: -32768, max: 32767 };
      case 'Uint16Array':
        return { min: 0, max: 65535 };
      case 'Int32Array':
        return { min: -2147483648, max: 2147483647 };
      case 'Uint32Array':
        return { min: 0, max: 4294967295 };
      case 'Float32Array':
        return { min: -1.2e38, max: 1.2e38 };
      case 'Float64Array':
        return { min: -1.2e38, max: 1.2e38 };
      case 'Uint8Array':
      case 'Uint8ClampedArray':
      default:
        return { min: 0, max: 255 };
    }
  };

  publicAPI.clamp = (outPtr, inPtr, numscalars, n, min, max) => {
    const count = n * numscalars;
    for (let i = 0; i < count; ++i) {
      outPtr[i] = vtkInterpolationMathClamp(inPtr[i], min, max);
    }
    return count;
  };

  publicAPI.convert = (outPtr, inPtr, numscalars, n) => {
    const count = n * numscalars;
    for (let i = 0; i < count; ++i) {
      outPtr[i] = Math.round(inPtr[i]);
    }
    return count;
  };

  publicAPI.getConversionFunc = (
    inputType,
    dataType,
    scalarShift,
    scalarScale,
    forceClamping
  ) => {
    let useClamping = forceClamping;
    if (
      dataType !== VtkDataTypes.FLOAT &&
      dataType !== VtkDataTypes.DOUBLE &&
      !forceClamping
    ) {
      const inMinMax = publicAPI.getDataTypeMinMax(inputType);
      let checkMin = (inMinMax.min + scalarShift) * scalarScale;
      let checkMax = (inMinMax.max + scalarShift) * scalarScale;
      const outMinMax = publicAPI.getDataTypeMinMax(dataType);
      const outputMin = outMinMax.min;
      const outputMax = outMinMax.max;
      if (checkMin > checkMax) {
        const tmp = checkMax;
        checkMax = checkMin;
        checkMin = tmp;
      }
      useClamping = checkMin < outputMin || checkMax > outputMax;
    }

    if (
      useClamping &&
      dataType !== VtkDataTypes.FLOAT &&
      dataType !== VtkDataTypes.DOUBLE
    ) {
      const minMax = publicAPI.getDataTypeMinMax(dataType);
      const clamp = (outPtr, inPtr, numscalars, n) =>
        publicAPI.clamp(outPtr, inPtr, numscalars, n, minMax.min, minMax.max);
      return clamp;
    }
    return publicAPI.convert;
  };

  publicAPI.set = (outPtr, inPtr, numscalars, n) => {
    const count = numscalars * n;
    for (let i = 0; i < n; ++i) {
      outPtr[i] = inPtr[i];
    }
    return count;
  };

  publicAPI.set1 = (outPtr, inPtr, numscalars, n) => {
    outPtr.fill(inPtr[0], 0, n);
    return n;
  };

  publicAPI.getSetPixelsFunc = (dataType, dataSize, numscalars, dataPtr) =>
    numscalars === 1 ? publicAPI.set1 : publicAPI.set;

  publicAPI.getCompositeFunc = (slabMode, slabTrapezoidIntegration) => {
    let composite = null;
    // eslint-disable-next-line default-case
    switch (slabMode) {
      case SlabMode.MIN:
        composite = getImageResliceCompositeMinValue;
        break;
      case SlabMode.MAX:
        composite = getImageResliceCompositeMaxValue;
        break;
      case SlabMode.MEAN:
        if (slabTrapezoidIntegration) {
          composite = getImageResliceCompositeMeanTrap;
        } else {
          composite = getImageResliceCompositeMeanValue;
        }
        break;
      case SlabMode.SUM:
        if (slabTrapezoidIntegration) {
          composite = getImageResliceCompositeSumTrap;
        } else {
          composite = getImageResliceCompositeSumValue;
        }
        break;
    }

    return composite;
  };

  publicAPI.applyTransform = (newTrans, inPoint, inOrigin, inInvSpacing) => {
    inPoint[3] = 1;
    vec4.transformMat4(inPoint, inPoint, newTrans);
    inPoint[0] -= inOrigin[0];
    inPoint[1] -= inOrigin[1];
    inPoint[2] -= inOrigin[2];
    inPoint[0] *= inInvSpacing[0];
    inPoint[1] *= inInvSpacing[1];
    inPoint[2] *= inInvSpacing[2];
  };

  publicAPI.rescaleScalars = (
    floatData,
    components,
    n,
    scalarShift,
    scalarScale
  ) => {
    const m = n * components;
    for (let i = 0; i < m; ++i) {
      floatData[i] = (floatData[i] + scalarShift) * scalarScale;
    }
  };

  publicAPI.isPermutationMatrix = (matrix) => {
    for (let i = 0; i < 3; i++) {
      if (matrix[4 * i + 3] !== 0) {
        return false;
      }
    }
    if (matrix[4 * 3 + 3] !== 1) {
      return false;
    }
    for (let j = 0; j < 3; j++) {
      let k = 0;
      for (let i = 0; i < 3; i++) {
        if (matrix[4 * j + i] !== 0) {
          k++;
        }
      }
      if (k !== 1) {
        return 0;
      }
    }
    return 1;
  };

  // TODO: to move in vtkMath and add tolerance
  publicAPI.isIdentityMatrix = (matrix) => {
    for (let i = 0; i < 4; ++i) {
      for (let j = 0; j < 4; ++j) {
        if ((i === j ? 1.0 : 0.0) !== matrix[4 * j + i]) {
          return false;
        }
      }
    }
    return true;
  };

  publicAPI.isPerspectiveMatrix = (matrix) =>
    matrix[4 * 0 + 3] !== 0 ||
    matrix[4 * 1 + 3] !== 0 ||
    matrix[4 * 2 + 3] !== 0 ||
    matrix[4 * 3 + 3] !== 1;

  publicAPI.canUseNearestNeighbor = (matrix, outExt) => {
    // loop through dimensions
    for (let i = 0; i < 3; i++) {
      let j;
      for (j = 0; j < 3; j++) {
        if (matrix[4 * j + i] !== 0) {
          break;
        }
      }
      if (j >= 3) {
        return 0;
      }
      let x = matrix[4 * j + i];
      let y = matrix[4 * 3 + i];
      if (outExt[2 * j] === outExt[2 * j + 1]) {
        y += x * outExt[2 * i];
        x = 0;
      }
      const fx = vtkInterpolationMathFloor(x, 0).error;
      const fy = vtkInterpolationMathFloor(y, 0).error;
      if (fx !== 0 || fy !== 0) {
        return 0;
      }
    }
    return 1;
  };
}

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

const DEFAULT_VALUES = {
  transformInputSampling: true,
  autoCropOutput: false,
  outputDimensionality: 3,
  outputSpacing: null, // automatically computed if null
  outputOrigin: null, // automatically computed if null
  outputDirection: null, // identity if null
  outputExtent: null, // automatically computed if null
  outputScalarType: null,
  wrap: false, // don't wrap
  mirror: false, // don't mirror
  border: true, // apply a border
  interpolationMode: InterpolationMode.NEAREST, // only NEAREST supported so far
  slabMode: SlabMode.MIN,
  slabTrapezoidIntegration: false,
  slabNumberOfSlices: 1,
  slabSliceSpacingFraction: 1,
  optimization: false, // not supported yet
  scalarShift: 0, // for rescaling the data
  scalarScale: 1,
  backgroundColor: [0, 0, 0, 0],
  resliceAxes: null,
  // resliceTransform: null,
  interpolator: vtkImageInterpolator.newInstance(),
  usePermuteExecute: false, // no supported yet
};

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

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, [
    'outputDimensionality',
    'outputScalarType',
    'scalarShift',
    'scalarScale',
    'transformInputSampling',
    'autoCropOutput',
    'wrap',
    'mirror',
    'border',
    'interpolationMode',
    'resliceTransform',
    'slabMode',
    'slabTrapezoidIntegration',
    'slabNumberOfSlices',
    'slabSliceSpacingFraction',
  ]);

  macro.setGetArray(publicAPI, model, ['outputOrigin', 'outputSpacing'], 3);
  macro.setGetArray(publicAPI, model, ['outputExtent'], 6);
  macro.setGetArray(publicAPI, model, ['outputDirection'], 9);
  macro.setGetArray(publicAPI, model, ['backgroundColor'], 4);

  macro.get(publicAPI, model, ['resliceAxes']);

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

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

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

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

export default { newInstance, extend, ...Constants };