ImageData

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;
}

if (!inExtent || inExtent.length !== 6) {
return;
}

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

if (changeDetected) {
model.extent = [].concat(inExtent);
model.dataDescription = vtkStructuredData.getDataDescriptionFromExtent(
model.extent
);
publicAPI.modified();
}
};

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;
};

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 = 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);
};

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;
};

// 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);
};

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

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;
};

// 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);
};

// 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,
};
};
}

// ----------------------------------------------------------------------------
// 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 };