import {
getIdentityMat33,
REAL_WORLD_EPSILON,
REAL_WORLD_EXPONENT
} from '../math/matrix.js';
import {Point3D, Point} from '../math/point.js';
import {Vector3D} from '../math/vector.js';
import {Index} from '../math/index.js';
import {getBasicStats} from '../math/stats.js';
import {precisionRound} from '../utils/string.js';
import {logger} from '../utils/logger.js';
import {Size} from './size.js';
import {Spacing} from './spacing.js';
import {BooleanResult} from '../utils/result.js';
import {arrayMap} from '../utils/array.js';
import {Matrix33} from '../math/matrix.js';
/**
* 2D/3D Geometry class.
*/
export class Geometry {
/**
* Array of origins.
*
* @type {Point3D[]}
*/
#origins;
/**
* Data size.
*
* @type {Size}
*/
#size;
/**
* Data spacing.
*
* @type {Spacing}
*/
#spacing;
/**
* Local helper object for time points.
*
* @type {Object<string, Point3D[]>}
*/
#timeOrigins = {};
/**
* Initial time index.
*
* @type {number}
*/
#initialTime;
/**
* Data orientation.
*
* @type {Matrix33}
*/
#orientation = getIdentityMat33();
/**
* @param {Point3D[]} origins The object origins.
* @param {Size} size The object size.
* @param {Spacing} spacing The object spacing.
* @param {Matrix33} [orientation] The object orientation (3*3 matrix,
* default to 3*3 identity).
* @param {number} [time] Optional time index.
*/
constructor(origins, size, spacing, orientation, time) {
this.#origins = origins;
this.#size = size;
this.#spacing = spacing;
if (typeof time !== 'undefined') {
this.#initialTime = time;
this.#timeOrigins[time] = origins.slice();
}
// check input orientation
if (typeof orientation !== 'undefined') {
this.#orientation = orientation;
}
}
/**
* Clone the geometry.
*
* @returns {Geometry} A clone of this geometry.
*/
clone() {
return new Geometry(
this.#origins.slice(),
new Size(this.#size.getValues()),
new Spacing(this.#spacing.getValues()),
new Matrix33(this.#orientation.getValues()),
this.#initialTime
);
}
/**
* Get the time value that was passed at construction.
*
* @returns {number} The time value.
*/
getInitialTime() {
return this.#initialTime;
}
/**
* Set the initial time.
*
* @param {number} time The new time.
*/
setInitialTime(time) {
if (typeof this.#timeOrigins[this.#initialTime] !== 'undefined') {
delete this.#timeOrigins[this.#initialTime];
}
this.#initialTime = time;
this.#timeOrigins[time] = this.#origins.slice();
}
/**
* Get the total number of slices.
* Can be different from what is stored in the size object
* during a volume with time points creation process.
*
* @returns {number} The total count.
*/
getCurrentTotalNumberOfSlices() {
const keys = Object.keys(this.#timeOrigins);
if (keys.length === 0) {
return this.#origins.length;
}
let count = 0;
for (let i = 0; i < keys.length; ++i) {
count += this.#timeOrigins[keys[i]].length;
}
return count;
}
/**
* Check if a time point has associated slices.
*
* @param {number} time The time point to check.
* @returns {boolean} True if slices are present.
*/
hasSlicesAtTime(time) {
return typeof this.#timeOrigins[time] !== 'undefined';
}
/**
* Get the number of slices stored for time points preceding
* the input one.
*
* @param {number} time The time point to check.
* @returns {number|undefined} The count.
*/
getCurrentNumberOfSlicesBeforeTime(time) {
const keys = Object.keys(this.#timeOrigins);
if (keys.length === 0) {
return undefined;
}
let count = 0;
for (let i = 0; i < keys.length; ++i) {
const key = keys[i];
if (parseInt(key, 10) === time) {
break;
}
count += this.#timeOrigins[key].length;
}
return count;
}
/**
* Get the object origin.
* This should be the lowest origin to ease calculations (?).
*
* @returns {Point3D} The object origin.
*/
getOrigin() {
return this.#origins[0];
}
/**
* Get the object origins.
*
* @returns {Point3D[]} The object origins.
*/
getOrigins() {
return this.#origins;
}
/**
* Check if a point is in the origin list.
*
* @param {Point3D} point3D The point to check.
* @param {number} [tol] Optional number comparison tolerance
* defaults to Number.EPSILON.
* @returns {boolean} True if in list.
*/
includesOrigin(point3D, tol) {
for (let i = 0; i < this.#origins.length; ++i) {
if (this.#origins[i].isSimilar(point3D, tol)) {
return true;
}
}
return false;
}
/**
* Get the object size.
* Warning: the size comes as stored in DICOM, meaning that it could
* be oriented.
*
* @param {Matrix33} [viewOrientation] The view orientation (optional).
* @returns {Size} The object size.
*/
getSize(viewOrientation) {
let res = this.#size;
if (viewOrientation && typeof viewOrientation !== 'undefined') {
let values = getOrientedArray3D(
[
this.#size.get(0),
this.#size.get(1),
this.#size.get(2)
],
viewOrientation);
values = values.map(Math.abs);
res = new Size(values.concat(this.#size.getValues().slice(3)));
}
return res;
}
/**
* Calculate slice spacing from origins and replace current
* if needed.
*/
updateSliceSpacing() {
const geoSliceSpacing = getSliceGeometrySpacing(this.#origins);
// update local if needed
if (typeof geoSliceSpacing !== 'undefined') {
// 1 is default
if (this.#spacing.get(2) !== 1 &&
this.#spacing.get(2) !== geoSliceSpacing) {
logger.debug('Using geometric spacing ' + geoSliceSpacing +
' instead of tag spacing ' + this.#spacing.get(2));
}
const values = this.#spacing.getValues();
values[2] = geoSliceSpacing;
this.#spacing = new Spacing(values);
}
}
/**
* Get the object spacing.
* Warning: the spacing comes as stored in DICOM, meaning that it could
* be oriented.
*
* @param {Matrix33} [viewOrientation] The view orientation (optional).
* @returns {Spacing} The object spacing.
*/
getSpacing(viewOrientation) {
let res = this.#spacing;
if (viewOrientation && typeof viewOrientation !== 'undefined') {
let orientedValues = getOrientedArray3D(
[
this.#spacing.get(0),
this.#spacing.get(1),
this.#spacing.get(2)
],
viewOrientation);
orientedValues = orientedValues.map(Math.abs);
res = new Spacing(orientedValues);
}
return res;
}
/**
* Get the image spacing in real world.
*
* @returns {Spacing} The object spacing.
*/
getRealSpacing() {
// asOneAndZeros to not change spacing values...
return this.getSpacing(
this.#orientation.getInverse().asOneAndZeros()
);
}
/**
* Get the object orientation.
*
* @returns {Matrix33} The object orientation.
*/
getOrientation() {
return this.#orientation;
}
/**
* Get the slice position of a point in the current slice layout.
* Slice indices increase with decreasing origins (high index -> low origin),
* this simplified the handling of reconstruction since it means
* the displayed data is in the same 'direction' as the extracted data.
* As seen in the getOrigin method, the main origin is the lowest one.
* This implies that the index to world and reverse method do some flipping
* magic...
*
* @param {Point3D} point The point to evaluate.
* @param {number} [time] Optional time index.
* @returns {number} The slice index.
*/
getSliceIndex(point, time) {
// cannot use this.worldToIndex(point).getK() since
// we cannot guaranty consecutive slices...
let localOrigins = this.#origins;
if (typeof time !== 'undefined') {
localOrigins = this.#timeOrigins[time];
}
// find the closest origin
const closestOriginIndex = point.getClosest(localOrigins);
const closestOrigin = localOrigins[closestOriginIndex];
// direction between the input point and the closest origin
const pointDir = point.minus(closestOrigin);
// use third orientation matrix column as plane normal vector
const normal = new Vector3D(
this.#orientation.get(0, 2),
this.#orientation.get(1, 2),
this.#orientation.get(2, 2)
);
// codirectional vectors: above slice index
// oposite vectors: below slice index
const isCodirectional = normal.isCodirectional(pointDir);
const sliceIndex = isCodirectional
? closestOriginIndex + 1 : closestOriginIndex;
return sliceIndex;
}
/**
* Check if another geometry is compatible with this one, ie that
* slices with these geometries can be merged into one volume.
*
* @param {Geometry} rhs The geometry to check.
* @returns {BooleanResult} Result with success set to true if
* the geometry is compatible.
*/
canAppend(rhs) {
// check size
const rhsSize = rhs.getSize();
if (rhsSize.get(2) !== 1) {
return {
success: false,
message: 'Cannot append more than one slice'
};
}
const size = this.getSize();
if (size.get(0) !== rhsSize.get(0)) {
return {
success: false,
message: 'Cannot append a slice with different number of columns'
};
}
if (size.get(1) !== rhsSize.get(1)) {
return {
success: false,
message: 'Cannot append a slice with different number of rows'
};
}
// check orientation
if (!this.getOrientation().isSimilar(
rhs.getOrientation(), REAL_WORLD_EPSILON)) {
return {
success: false,
message: 'Cannot append a slice with different orientation'
};
}
// check origin
const canAppendOrigin =
this.canAppendOrigin(rhs.getOrigin(), rhs.getInitialTime());
if (!canAppendOrigin.success) {
return canAppendOrigin;
}
return new BooleanResult(true);
}
/**
* Check if an origin can be appended to this geometry.
*
* @param {Point3D} origin The origin to append.
* @param {number} [time] Optional time index.
* @returns {BooleanResult} Result with success set to true if
* the origin is compatible.
*/
canAppendOrigin(origin, time) {
// equal callback
const equalToOrigin = function (element) {
return element.equals(origin);
};
if (typeof time !== 'undefined' && time !== this.#initialTime) {
// check if not already in list
if (typeof this.#timeOrigins[time] !== 'undefined') {
const found = this.#timeOrigins[time].find(equalToOrigin);
if (typeof found !== 'undefined') {
return {
success: false,
message: 'Cannot append same time origin twice'
};
}
}
} else {
// check if not already in list
const found = this.#origins.find(equalToOrigin);
if (typeof found !== 'undefined') {
return {
success: false,
message: 'Cannot append same origin twice'
};
}
}
return new BooleanResult(true);
}
/**
* Append an origin to the geometry.
*
* @param {Point3D} origin The origin to append.
* @param {number} index The index at which to append.
* @param {number} [time] Optional time index.
*/
appendOrigin(origin, index, time) {
// check if possible
const canAppend = this.canAppendOrigin(origin, time);
if (!canAppend.success) {
throw new Error(canAppend.message);
}
// add in origin array
if (typeof time !== 'undefined') {
this.#timeOrigins[time].splice(index, 0, origin);
}
if (typeof time === 'undefined' || time === this.#initialTime) {
this.#origins.splice(index, 0, origin);
// increment second dimension
const values = this.#size.getValues();
values[2] += 1;
this.#size = new Size(values);
}
}
/**
* Append a frame to the geometry.
*
* @param {Point3D} origin The origin to append.
* @param {number} time Optional time index.
*/
appendFrame(origin, time) {
// add origin to list
this.#timeOrigins[time] = [origin];
// increment third dimension
const sizeValues = this.#size.getValues();
const spacingValues = this.#spacing.getValues();
if (sizeValues.length === 4) {
sizeValues[3] += 1;
} else {
sizeValues.push(2);
spacingValues.push(1);
}
this.#size = new Size(sizeValues);
this.#spacing = new Spacing(spacingValues);
}
/**
* Get a string representation of the geometry.
*
* @returns {string} The geometry as a string.
*/
toString() {
return 'Origin: ' + this.getOrigin() +
', Size: ' + this.getSize() +
', Spacing: ' + this.getSpacing() +
', Orientation: ' + this.getOrientation();
}
/**
* Check for equality.
*
* @param {Geometry} rhs The object to compare to.
* @returns {boolean} True if both objects are equal.
*/
equals(rhs) {
return rhs !== null &&
typeof rhs !== 'undefined' &&
this.getOrigin().equals(rhs.getOrigin()) &&
this.getSize().equals(rhs.getSize()) &&
this.getSpacing().equals(rhs.getSpacing()) &&
this.getOrientation().equals(rhs.getOrientation());
}
/**
* Check that a point is within bounds.
*
* @param {Point} point The point to check.
* @returns {boolean} True if the given coordinates are within bounds.
*/
isInBounds(point) {
return this.isIndexInBounds(this.worldToIndex(point));
}
/**
* Check that a index is within bounds.
*
* @param {Index} index The index to check.
* @param {number[]} [dirs] Optional list of directions to check.
* @returns {boolean} True if the given coordinates are within bounds.
*/
isIndexInBounds(index, dirs) {
return this.getSize().isInBounds(index, dirs);
}
/**
* Get the geometrical range, ie the minimum and maximum
* positions.
*
* @returns {Point[]} The min and max positions.
*/
getRange() {
const nDims = this.getSize().length();
const minValues = new Array(nDims);
minValues.fill(0);
const minIndex = new Index(minValues);
const maxIndex = new Index(this.getSize().getValues());
return [
this.indexToWorld(minIndex),
this.indexToWorld(maxIndex)
];
}
/**
* Convert an index into world coordinates.
*
* @param {Index} index The index to convert.
* @returns {Point} The corresponding point.
*/
indexToWorld(index) {
// apply spacing
// (spacing is oriented, apply before orientation)
const spacing = this.getSpacing();
const orientedPoint3D = new Point3D(
index.get(0) * spacing.get(0),
index.get(1) * spacing.get(1),
index.get(2) * spacing.get(2)
);
// de-orient
const point3D = this.getOrientation().multiplyPoint3D(orientedPoint3D);
// keep >3d values
const values = index.getValues();
const origin = this.getOrigin();
values[0] = origin.getX() + point3D.getX();
values[1] = origin.getY() + point3D.getY();
values[2] = origin.getZ() + point3D.getZ();
// return point
return new Point(values);
}
/**
* Convert a 3D point into world coordinates.
*
* @param {Point3D} point The 3D point to convert.
* @returns {Point3D} The corresponding world 3D point.
*/
pointToWorld(point) {
// apply spacing
// (spacing is oriented, apply before orientation)
const spacing = this.getSpacing();
const orientedPoint3D = new Point3D(
point.getX() * spacing.get(0),
point.getY() * spacing.get(1),
point.getZ() * spacing.get(2)
);
// de-orient
const point3D = this.getOrientation().multiplyPoint3D(orientedPoint3D);
// return point3D
const origin = this.getOrigin();
return new Point3D(
origin.getX() + point3D.getX(),
origin.getY() + point3D.getY(),
origin.getZ() + point3D.getZ()
);
}
/**
* Convert world coordinates into an index.
*
* @param {Point} point The point to convert.
* @returns {Index} The corresponding index.
*/
worldToIndex(point) {
// compensate for origin
// (origin is not oriented, compensate before orientation)
// TODO: use slice origin...
const origin = this.getOrigin();
const point3D = new Point3D(
point.get(0) - origin.getX(),
point.get(1) - origin.getY(),
point.get(2) - origin.getZ()
);
// orient
const orientedPoint3D =
this.getOrientation().getInverse().multiplyPoint3D(point3D);
// keep >3d values
const values = point.getValues();
// apply spacing and floor
// (precisionRound to avoid float number inaccuracies)
const spacing = this.getSpacing();
values[0] = Math.floor(precisionRound(
orientedPoint3D.getX() / spacing.get(0),
REAL_WORLD_EXPONENT
));
values[1] = Math.floor(precisionRound(
orientedPoint3D.getY() / spacing.get(1),
REAL_WORLD_EXPONENT
));
values[2] = Math.floor(precisionRound(
orientedPoint3D.getZ() / spacing.get(2),
REAL_WORLD_EXPONENT
));
// return index
return new Index(values);
}
/**
* Convert world coordinates into an point.
*
* @param {Point} point The world point to convert.
* @returns {Point3D} The corresponding point.
*/
worldToPoint(point) {
// compensate for origin
// (origin is not oriented, compensate before orientation)
const origin = this.getOrigin();
const point3D = new Point3D(
point.get(0) - origin.getX(),
point.get(1) - origin.getY(),
point.get(2) - origin.getZ()
);
// orient
const orientedPoint3D =
this.getOrientation().getInverse().multiplyPoint3D(point3D);
// keep >3d values
const values = point.getValues();
// apply spacing and round
const spacing = this.getSpacing();
values[0] = orientedPoint3D.getX() / spacing.get(0);
values[1] = orientedPoint3D.getY() / spacing.get(1);
values[2] = orientedPoint3D.getZ() / spacing.get(2);
// return index
return new Point3D(values[0], values[1], values[2]);
}
} // class Geometry
/**
* Get the oriented values of an input 3D array.
*
* @param {number[]} array3D The 3D array.
* @param {Matrix33} orientation The orientation 3D matrix.
* @returns {number[]} The values reordered according to the orientation.
*/
export function getOrientedArray3D(array3D, orientation) {
// values = orientation * orientedValues
// -> inv(orientation) * values = orientedValues
return orientation.getInverse().multiplyArray3D(array3D);
}
/**
* Get the raw values of an oriented input 3D array.
*
* @param {number[]} array3D The 3D array.
* @param {Matrix33} orientation The orientation 3D matrix.
* @returns {number[]} The values reordered to compensate the orientation.
*/
export function getDeOrientedArray3D(array3D, orientation) {
// values = orientation * orientedValues
return orientation.multiplyArray3D(array3D);
}
/**
* Get the slice spacing from the difference in the Z directions
* of input origins.
*
* @param {Point3D[]} origins An array of Point3D.
* @returns {number|undefined} The spacing.
*/
export function getSliceGeometrySpacing(origins) {
// check origins
if (origins.length <= 1) {
return;
}
const spacings = [];
for (let i = 0; i < origins.length - 1; ++i) {
const origin1 = origins[i];
const origin2 = origins[i + 1];
const sliceSpacing = origin1.getDistance(origin2);
if (sliceSpacing === 0) {
throw new Error('Zero slice spacing ' +
origin1.toString() + ' ' + origin2.toString());
}
spacings.push(sliceSpacing);
}
// use rounded mean value as spacing
const stats = getBasicStats(spacings);
const spacing = precisionRound(stats.mean, 4);
// warn if non constant
if (stats.stdDev > REAL_WORLD_EPSILON) {
logger.warn('Varying slice spacing, value: ' + spacing +
' (mean: ' + stats.mean +
', min: ' + stats.min +
', max: ' + stats.max +
', stdDev: ' + stats.stdDev + ')');
}
return spacing;
}
/**
* Merge two geometries into one using the merge size and
* smallest resolution.
* WARNING: needs to be called with the final geometries.
*
* @param {Geometry} geometry1 The first geometry.
* @param {Geometry} geometry2 The second geometry.
* @returns {Geometry} The merged geometry.
*/
export function mergeGeometries(geometry1, geometry2) {
const orientation = geometry1.getOrientation();
// check input
if (!orientation.isSimilar(geometry2.getOrientation(), REAL_WORLD_EPSILON)) {
throw new Error('Cannot merge geometries with different orientation');
}
const invOrientation = orientation.getInverse();
const minByIndex = function (array1, array2) {
return arrayMap(array1, array2, Math.min);
};
const maxByIndex = function (array1, array2) {
return arrayMap(array1, array2, Math.max);
};
const newSpacing = new Spacing(minByIndex(
geometry1.getSpacing().getValues(),
geometry2.getSpacing().getValues()
));
const deorient = function (item) {
return invOrientation.multiplyPoint3D(item.get3D());
};
const range1 = geometry1.getRange().map(deorient);
const range2 = geometry2.getRange().map(deorient);
const minRangeValues = minByIndex(
range1[0].getValues(),
range2[0].getValues()
);
const maxRangeValues = maxByIndex(
range1[1].getValues(),
range2[1].getValues()
);
let sizeValues = [];
for (let i = 0; i < minRangeValues.length; ++i) {
sizeValues.push(Math.round(
Math.abs(maxRangeValues[i] - minRangeValues[i]) / newSpacing.get(i)
));
}
// handle extra dims
const maxSize = maxByIndex(
geometry1.getSize().getValues(),
geometry2.getSize().getValues()
);
sizeValues = sizeValues.concat(maxSize.slice(3));
const newSize = new Size(sizeValues);
const newOrigins = [];
for (let i = 0; i < newSize.get(2); ++i) {
const values = minRangeValues.slice();
values[2] = minRangeValues[2] + i * newSpacing.get(2);
const origin = new Point3D(values[0], values[1], values[2]);
newOrigins.push(orientation.multiplyPoint3D(origin));
}
return new Geometry(
newOrigins,
newSize,
newSpacing,
orientation
);
}