/* Copyright (c) 2001-2015, David A. Clunie DBA Pixelmed Publishing. All rights reserved. */ package org.weasis.dicom.codec.geometry; import java.awt.geom.Point2D; import java.util.ArrayList; import java.util.List; import javax.vecmath.Matrix3d; import javax.vecmath.Point3d; import javax.vecmath.Tuple3d; import javax.vecmath.Vector3d; import org.weasis.core.api.gui.util.MathUtil; /** * <p> * An abstract class that provides the basis for posting the position of specified slices and volumes on (usually * orthogonal) localizer images. * </p> * * <p> * This base class provides the interface, common storage and various utility methods, and specific methods of * performing the localization operation are performed by concrete sub-classes, instantiated through a factory class. * </p> * * <p> * Typically this would b used as follows, to get outlines in the form of a vector of shapes whose coordinates are those * of the localizer image: * </p> * * <pre> * GeometryOfSlice localizerGeometry = new GeometryOfSliceFromAttributeList(localizerAttributeList); * GeometryOfSlice postImageGeometry = new GeometryOfSliceFromAttributeList(postImageAttributeList); * LocalizerPoster localizerPoster = LocalizerPosterFactory.getLocalizerPoster(false, false); * localizerPoster.setLocalizerGeometry(localizerGeometry); * Vector shapes = localizerPoster.getOutlineOnLocalizerForThisGeometry(postImageGeometry); * </pre> * * @see com.pixelmed.geometry.LocalizerPosterFactory * @see com.pixelmed.geometry.GeometryOfSlice * @see com.pixelmed.dicom.GeometryOfSliceFromAttributeList * * @author dclunie */ public abstract class LocalizerPoster { protected Vector3d localizerRow; protected Vector3d localizerColumn; protected Vector3d localizerNormal; protected Point3d localizerTLHC; protected Tuple3d localizerVoxelSpacing; // row spacing (between centers of adjacent rows), then column spacing, // then slice spacing protected double[] localizerVoxelSpacingArray; protected Tuple3d localizerDimensions; // number of rows, then number of columns, then number of slices protected double[] localizerDimensionsArray; protected Matrix3d rotateIntoLocalizerSpace; public LocalizerPoster(Vector3d row, Vector3d column, Point3d tlhc, Tuple3d voxelSpacing, Tuple3d dimensions) { localizerRow = row; localizerColumn = column; localizerTLHC = tlhc; localizerVoxelSpacing = voxelSpacing; localizerDimensions = dimensions; doCommonConstructorStuff(); } public LocalizerPoster(GeometryOfSlice geometry) { localizerRow = geometry.getRow(); localizerColumn = geometry.getColumn(); localizerTLHC = geometry.getTLHC(); localizerVoxelSpacing = geometry.getVoxelSpacing(); localizerDimensions = geometry.getDimensions(); doCommonConstructorStuff(); } /** * <p> * Check that the row and column vectors are unit vectors and are orthogonal. * </p> * * @param row * the row direction cosines * @param column * the column direction cosines * @throws IllegalArgumentException * thrown if not */ public static void validateDirectionCosines(Vector3d row, Vector3d column) { if (Math.abs(row.lengthSquared() - 1) > 0.001) { throw new IllegalArgumentException("Row not a unit vector"); //$NON-NLS-1$ } if (Math.abs(column.lengthSquared() - 1) > 0.001) { throw new IllegalArgumentException("Column not a unit vector"); //$NON-NLS-1$ } if (row.dot(column) > 0.005) { // dot product should be cos(90)=0 if orthogonal throw new IllegalArgumentException("Row and column vectors are not orthogonal"); //$NON-NLS-1$ } } /** * <p> * Get the corners of a slice in the 3D coordinate space of that slice. * </p> * * @param row * the direction of the row as X, Y and Z components (direction cosines, unit vector) LPH+ * @param column * the direction of the column as X, Y and Z components (direction cosines, unit vector) LPH+ * @param originalTLHC * the position of the top left hand corner of the slice as a point (X, Y and Z) LPH+ * @param voxelSpacing * the row and column spacing and the slice interval * @param dimensions * the row and column dimensions and 1 for the third dimension * @return an array of four points that are the tlhc,trhc, brhc, blhc of the slice */ public static Point3d[] getCornersOfSourceRectangleInSourceSpace(Vector3d row, Vector3d column, Point3d originalTLHC, Tuple3d voxelSpacing, Tuple3d dimensions) { validateDirectionCosines(row, column); double[] spacingArray = new double[3]; voxelSpacing.get(spacingArray); double[] dimensionsArray = new double[3]; dimensions.get(dimensionsArray); Vector3d distanceAlongRow = new Vector3d(row); distanceAlongRow.scale((dimensionsArray[1]/* cols */) * spacingArray[1]/* between cols */); Vector3d distanceAlongColumn = new Vector3d(column); distanceAlongColumn.scale((dimensionsArray[0]/* rows */) * spacingArray[0]/* between rows */); // Build a square to project with 4 corners TLHC, TRHC, BRHC, BLHC ... Point3d tlhc = new Point3d(originalTLHC); // otherwise original TLHC gets changed later on Point3d trhc = new Point3d(tlhc); trhc.add(distanceAlongRow); Point3d blhc = new Point3d(tlhc); blhc.add(distanceAlongColumn); Point3d brhc = new Point3d(tlhc); brhc.add(distanceAlongRow); brhc.add(distanceAlongColumn); return new Point3d[] { tlhc, trhc, brhc, blhc }; } /** * <p> * Get the corners of a volume in the 3D coordinate space of that volume. * </p> * * @param row * the direction of the row as X, Y and Z components (direction cosines, unit vector) LPH+ * @param column * the direction of the column as X, Y and Z components (direction cosines, unit vector) LPH+ * @param originalTLHC * the position of the top left hand corner of the slice as a point (X, Y and Z) LPH+ * @param voxelSpacing * the row and column spacing and the slice interval * @param sliceThickness * the slice thickness * @param dimensions * the row and column dimensions and number of frames for the third dimension * @return an array of eight points that are the tlhcT, trhcT, brhcT, blhcT, tlhcB, trhcB, brhcB, blhcB of the * volume */ public static Point3d[] getCornersOfSourceCubeInSourceSpace(Vector3d row, Vector3d column, Point3d originalTLHC, Tuple3d voxelSpacing, double sliceThickness, Tuple3d dimensions) { validateDirectionCosines(row, column); Vector3d normal = new Vector3d(); normal.cross(row, column); // the normal to the plane is the cross product of the row and column double[] spacingArray = new double[3]; voxelSpacing.get(spacingArray); double[] dimensionsArray = new double[3]; dimensions.get(dimensionsArray); Vector3d distanceAlongRow = new Vector3d(row); distanceAlongRow.scale((dimensionsArray[1]/* cols */) * spacingArray[1]/* between cols */); Vector3d distanceAlongColumn = new Vector3d(column); distanceAlongColumn.scale((dimensionsArray[0]/* rows */) * spacingArray[0]/* between rows */); Vector3d distanceAlongNormal = new Vector3d(normal); distanceAlongNormal.scale((dimensionsArray[2] / 2) * sliceThickness); // divide by two ... half on either side // of center // Build the "top" square to project with 4 corners TLHC, TRHC, BRHC, BLHC ... Point3d tlhcT = new Point3d(originalTLHC); // otherwise original TLHC gets changed later on tlhcT.add(distanceAlongNormal); Point3d trhcT = new Point3d(tlhcT); trhcT.add(distanceAlongRow); Point3d blhcT = new Point3d(tlhcT); blhcT.add(distanceAlongColumn); Point3d brhcT = new Point3d(tlhcT); brhcT.add(distanceAlongRow); brhcT.add(distanceAlongColumn); // Build the "bottom" square to project with 4 corners TLHC, TRHC, BRHC, BLHC ... Point3d tlhcB = new Point3d(originalTLHC); // otherwise original TLHC gets changed later on tlhcB.sub(distanceAlongNormal); Point3d trhcB = new Point3d(tlhcB); trhcB.add(distanceAlongRow); Point3d blhcB = new Point3d(tlhcB); blhcB.add(distanceAlongColumn); Point3d brhcB = new Point3d(tlhcB); brhcB.add(distanceAlongRow); brhcB.add(distanceAlongColumn); return new Point3d[] { tlhcT, trhcT, brhcT, blhcT, tlhcB, trhcB, brhcB, blhcB }; } /** * <p> * Transform a point into the "viewport" defined by the localizer that we are an instance of. * </p> * * @param point * a 3D point to be transformed * @return a new, transformed point */ protected Point3d transformPointFromSourceSpaceIntoLocalizerSpace(Point3d point) { Point3d newPoint = new Point3d(point); // do not overwrite the supplied point newPoint.sub(localizerTLHC); // move everything to origin of the target localizer rotateIntoLocalizerSpace.transform(newPoint); return newPoint; } /** * <p> * Project a point in localizer 3D space onto the plane of the localizer by ignoring the Z component, and return the * X and Y coordinates as image-TLHC relative column and row offsets. * </p> * * <p> * Will return sub-pixel values ranging from 0.5 to 0.5 less than the maximum dimensions of the image, which allows * points at the very edges of the rendered image to be drawn (e.g. a column of 0.5 will draw at the extreme left * and a column of 255.5 will draw at the extreme right of a 256 pixel wide image (whereas 256.0 will not, though * 0.0 will). * </p> * * @param point * the point in 3D localizer space, the Z coordinate of which is ignored * @return an array of 2 values in which the column, then the row location on the image is returned */ protected Point2D.Double transformPointInLocalizerPlaneIntoImageSpace(Point3d point) { // number of rows double scaleSubPixelHeightOfColumn = (localizerDimensionsArray[0] - 1) / localizerDimensionsArray[0]; // number of cols double scaleSubPixelWidthOfRow = (localizerDimensionsArray[1] - 1) / localizerDimensionsArray[1]; /* * NB. x is the column; use as height number of rows spacing between rows * * NB. y is the row; use as width number of cols * spacing between cols */ return new Point2D.Double(point.x / localizerVoxelSpacingArray[1] * scaleSubPixelHeightOfColumn + 0.5, point.y / localizerVoxelSpacingArray[0] * scaleSubPixelWidthOfRow + 0.5); } protected List<Point2D.Double> drawOutlineOnLocalizer(List<Point3d> corners) { if (corners != null && !corners.isEmpty()) { Point3d[] cornersArray = new Point3d[corners.size()]; return drawOutlineOnLocalizer(corners.toArray(cornersArray)); } return null; } protected List<Point2D.Double> drawOutlineOnLocalizer(Point3d[] corners) { ArrayList<Point2D.Double> shapes = new ArrayList<>(); for (int i = 0; i < corners.length; ++i) { shapes.add(transformPointInLocalizerPlaneIntoImageSpace(corners[i])); } return shapes; } protected Point3d intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(double[] a, double[] b) { double[] u = new double[3]; // be careful not to divide by zero when slope infinite (and unnecessary, since multiplicand is then zero) // Z of unknown point is zero u[1] = (MathUtil.isEqual(b[2], a[2])) ? a[1] : (b[1] - a[1]) / (b[2] - a[2]) * (0 - a[2]) + a[1]; u[0] = (MathUtil.isEqual(b[1], a[1])) ? a[0] : (b[0] - a[0]) / (b[1] - a[1]) * (u[1] - a[1]) + a[0]; u[2] = 0; return new Point3d(u); } protected Point3d intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(Point3d aP, Point3d bP) { double[] a = new double[3]; aP.get(a); double[] b = new double[3]; bP.get(b); return intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(a, b); } protected List<Point2D.Double> drawLinesBetweenAnyPointsWhichIntersectPlaneWhereZIsZero(Point3d[] corners) { int size = corners.length; double[] thisArray = new double[3]; double[] nextArray = new double[3]; ArrayList<Point3d> intersections = new ArrayList<>(); for (int i = 0; i < size; ++i) { int next = (i == size - 1) ? 0 : i + 1; corners[i].get(thisArray); double thisZ = thisArray[2]; corners[next].get(nextArray); double nextZ = nextArray[2]; if ((thisZ <= 0 && nextZ >= 0) || (thisZ >= 0 && nextZ <= 0)) { Point3d intersection = intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(thisArray, nextArray); intersections.add(intersection); } } return !intersections.isEmpty() ? drawOutlineOnLocalizer(intersections) : null; } protected static boolean classifyCornersIntoEdgeCrossingZPlane(Point3d startCorner, Point3d endCorner) { double[] startArray = new double[3]; double[] endArray = new double[3]; startCorner.get(startArray); double startZ = startArray[2]; endCorner.get(endArray); double endZ = endArray[2]; return (startZ <= 0 && endZ >= 0) || (startZ >= 0 && endZ <= 0); } protected List<Point3d> getIntersectionsOfCubeWithZPlane(Point3d[] corners) { ArrayList<Point3d> intersections = new ArrayList<>(4); // the check and traversal order are very dependent on the order of the // corners which are: tlhcT, trhcT, brhcT, blhcT, tlhcB, trhcB, brhcB, blhcB // as established in LocalizerPoster.getCornersOfSourceCubeInSourceSpace() // traverse each of the (three) possibilities for which opposite edges intersect the Z plane ... // 0,1 2,3 4,5 6,7 // 0,3 1,2 4,7 5,6 // 0,4 1,5 2,6 3,7 if (classifyCornersIntoEdgeCrossingZPlane(corners[0], corners[1]) && classifyCornersIntoEdgeCrossingZPlane(corners[2], corners[3]) && classifyCornersIntoEdgeCrossingZPlane(corners[4], corners[5]) && classifyCornersIntoEdgeCrossingZPlane(corners[6], corners[7])) { intersections.add(intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(corners[0], corners[1])); intersections.add(intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(corners[2], corners[3])); intersections.add(intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(corners[6], corners[7])); intersections.add(intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(corners[4], corners[5])); } else if (classifyCornersIntoEdgeCrossingZPlane(corners[0], corners[3]) && classifyCornersIntoEdgeCrossingZPlane(corners[1], corners[2]) && classifyCornersIntoEdgeCrossingZPlane(corners[4], corners[7]) && classifyCornersIntoEdgeCrossingZPlane(corners[5], corners[6])) { intersections.add(intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(corners[0], corners[3])); intersections.add(intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(corners[1], corners[2])); intersections.add(intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(corners[5], corners[6])); intersections.add(intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(corners[4], corners[7])); } else if (classifyCornersIntoEdgeCrossingZPlane(corners[0], corners[4]) && classifyCornersIntoEdgeCrossingZPlane(corners[1], corners[5]) && classifyCornersIntoEdgeCrossingZPlane(corners[2], corners[6]) && classifyCornersIntoEdgeCrossingZPlane(corners[3], corners[7])) { intersections.add(intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(corners[0], corners[4])); intersections.add(intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(corners[1], corners[5])); intersections.add(intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(corners[2], corners[6])); intersections.add(intersectLineBetweenTwoPointsWithPlaneWhereZIsZero(corners[3], corners[7])); } return intersections; } protected void doCommonConstructorStuff() { validateDirectionCosines(localizerRow, localizerColumn); localizerNormal = new Vector3d(); localizerNormal.cross(localizerRow, localizerColumn); // the normal to the plane is the cross product of the row // and column localizerVoxelSpacingArray = new double[3]; localizerVoxelSpacing.get(localizerVoxelSpacingArray); localizerDimensionsArray = new double[3]; localizerDimensions.get(localizerDimensionsArray); rotateIntoLocalizerSpace = new Matrix3d(); rotateIntoLocalizerSpace.setRow(0, localizerRow); rotateIntoLocalizerSpace.setRow(1, localizerColumn); rotateIntoLocalizerSpace.setRow(2, localizerNormal); } /** * <p> * Get the shapes on the localizer of the specified slice. * </p> * * @param row * the unit vector (direction cosine) of the row direction in the DICOM coordinate system * @param column * the unit vector (direction cosine) of the row direction in the DICOM coordinate system * @param tlhc * the position in the DICOM coordinate system of the center of the top left hand corner pixel of the * slice * @param voxelSpacing * the row and column and slice interval in mm * @param sliceThickness * the slice thickness in mm * @param dimensions * the number of rows and columns and slices * @return vector of shapes {@link java.awt.Shape java.awt.Shape} to be drawn in the localizer row and column * coordinates */ public abstract List<Point2D.Double> getOutlineOnLocalizerForThisGeometry(Vector3d row, Vector3d column, Point3d tlhc, Tuple3d voxelSpacing, double sliceThickness, Tuple3d dimensions); /** * <p> * Get the shape on the localizer of a zero-thickness slice specified by the geometry of a 2D rectangle. * </p> * * @param geometry * @return vector of shapes {@link java.awt.Shape java.awt.Shape} to be drawn in the localizer row and column * coordinates */ public final List<Point2D.Double> getOutlineOnLocalizerForThisGeometry(GeometryOfSlice geometry) { return getOutlineOnLocalizerForThisGeometry(geometry.getRow(), geometry.getColumn(), geometry.getTLHC(), geometry.getVoxelSpacing(), geometry.getSliceThickness(), geometry.getDimensions()); } }