diff --git a/CHANGELOG.md b/CHANGELOG.md index f69c77b..17424e4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,20 @@ # Terrasaur Changelog +## July 30, 2025 - v1.1.0 + +- Updates to existing tools + - AdjustShapeModelToOtherShapeModel + - fix intersection bug + - CreateSBMTStructure + - new options: -flipX, -flipY, -spice, -date, -observer, -target, -cameraFrame + - ValidateNormals + - new option: -fast to only check for overhangs if center and normal point in opposite directions + +- New tools: + - FacetInfo: Print info about a facet + - PointCloudOverlap: Find points in a point cloud which overlap a reference point cloud + - TriAx: Generate a triaxial ellipsoid in ICQ format + ## April 28, 2025 - v1.0.1 - Add MIT license to repository and source code diff --git a/doc/dist/Terrasaur-2025.03.03-e1a0e14-src.tar.gz b/doc/dist/Terrasaur-2025.03.03-e1a0e14-src.tar.gz deleted file mode 100644 index fd5ddba..0000000 Binary files a/doc/dist/Terrasaur-2025.03.03-e1a0e14-src.tar.gz and /dev/null differ diff --git a/doc/dist/Terrasaur-2025.03.03-e1a0e14_*.tar.gz b/doc/dist/Terrasaur-2025.03.03-e1a0e14_*.tar.gz deleted file mode 100644 index 0743b5d..0000000 Binary files a/doc/dist/Terrasaur-2025.03.03-e1a0e14_*.tar.gz and /dev/null differ diff --git a/doc/index.rst b/doc/index.rst index 4a18e52..1acee31 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -16,7 +16,7 @@ The Terrasaur package requires Java 21 or later. Some freely available versions Download ~~~~~~~~ -Binary packages for use on Mac OS X and Linux are available at ... +Binary packages for use on Mac OS X and Linux are available at `GitHub `__. We have not tried using the softare on Microsoft Windows, but users may try the Linux package with the `Windows Subsystem for Linux `__. diff --git a/doc/make.bat b/doc/make.bat index 922152e..2119f51 100644 --- a/doc/make.bat +++ b/doc/make.bat @@ -1,35 +1,35 @@ -@ECHO OFF - -pushd %~dp0 - -REM Command file for Sphinx documentation - -if "%SPHINXBUILD%" == "" ( - set SPHINXBUILD=sphinx-build -) -set SOURCEDIR=. -set BUILDDIR=_build - -if "%1" == "" goto help - -%SPHINXBUILD% >NUL 2>NUL -if errorlevel 9009 ( - echo. - echo.The 'sphinx-build' command was not found. Make sure you have Sphinx - echo.installed, then set the SPHINXBUILD environment variable to point - echo.to the full path of the 'sphinx-build' executable. Alternatively you - echo.may add the Sphinx directory to PATH. - echo. - echo.If you don't have Sphinx installed, grab it from - echo.http://sphinx-doc.org/ - exit /b 1 -) - -%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% -goto end - -:help -%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% - -:end -popd +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=. +set BUILDDIR=_build + +if "%1" == "" goto help + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/doc/tools/ColorSpots.rst b/doc/tools/ColorSpots.rst new file mode 100644 index 0000000..15485b2 --- /dev/null +++ b/doc/tools/ColorSpots.rst @@ -0,0 +1,70 @@ +.. _ColorSpots: + +########## +ColorSpots +########## + +ColorSpots takes as input a shape model and a file containing (x, y, z, value), +or (lat, lon, value). It writes out the mean and standard deviation of values +within a specified range for each facet. + +.. include:: ../toolDescriptions/ColorSpots.txt + :literal: + +******** +Examples +******** + +Download the :download:`Apophis<./support_files/apophis_g_15618mm_rdr_obj_0000n00000_v001.obj>` +shape model and the :download:`info<./support_files/xyzrandom.txt>` file containing +cartesian coordinates and a random value. + +Run ColorSpots: + +:: + + ColorSpots -obj apophis_g_15618mm_rdr_obj_0000n00000_v001.obj -xyz \ + -info xyzrandom.txt -outFile apophis_value_at_vertex.csv -noWeight \ + -allFacets -additionalFields n -searchRadius 0.015 -writeVertices + +The first few lines of apophis_value_at_vertex.csv look like: + +:: + + % head apophis_value_at_vertex.csv + 0.000000e+00, 0.000000e+00, 1.664960e-01, -3.805764e-02, 5.342315e-01, 4.000000e+01 + 1.589500e-02, 0.000000e+00, 1.591030e-01, 6.122849e-02, 6.017192e-01, 5.000000e+01 + 7.837000e-03, 1.486800e-02, 1.591670e-01, -6.072964e-03, 5.220682e-01, 5.700000e+01 + -7.747000e-03, 1.506300e-02, 1.621040e-01, 9.146163e-02, 5.488631e-01, 4.900000e+01 + -1.554900e-02, 0.000000e+00, 1.657970e-01, -8.172811e-03, 5.270302e-01, 3.400000e+01 + -7.982000e-03, -1.571100e-02, 1.694510e-01, -2.840524e-02, 5.045911e-01, 3.900000e+01 + 8.060000e-03, -1.543300e-02, 1.655150e-01, 3.531959e-02, 5.464390e-01, 4.900000e+01 + 3.179500e-02, 0.000000e+00, 1.515820e-01, -1.472434e-02, 5.967265e-01, 5.400000e+01 + 2.719700e-02, 1.658200e-02, 1.508930e-01, -9.050683e-03, 5.186966e-01, 4.700000e+01 + 1.554100e-02, 2.901300e-02, 1.530770e-01, -7.053547e-02, 4.980369e-01, 7.000000e+01 + +The columns are: + +.. list-table:: ColorSpots Vertex output + :header-rows: 1 + + * - Column + - Value + * - 1 + - X + * - 2 + - Y + * - 3 + - Z + * - 4 + - mean value in region + * - 5 + - standard deviation in region + * - 6 + - number of points in region + +.. figure:: images/ColorSpots-n.png + :alt: Number of points in region at each vertex + + This image shows the number of points in the region at each vertex. + diff --git a/doc/tools/CompareOBJ.rst b/doc/tools/CompareOBJ.rst index a0ad06c..6e95d46 100644 --- a/doc/tools/CompareOBJ.rst +++ b/doc/tools/CompareOBJ.rst @@ -23,7 +23,7 @@ Local Model Comparison Download the :download:`reference<./support_files/EVAL20_wtr.obj>` and :download:`comparison<./support_files/EVAL20.obj>` shape models. You can view them in a tool such as -`ParaView`. +`ParaView `__. .. figure:: images/CompareOBJ_local_1.png @@ -32,8 +32,8 @@ shape models. You can view them in a tool such as Run CompareOBJ to find the optimal transform to align the comparison with the reference: :: - CompareOBJ -computeOptimalRotationAndTranslation -model F3H-1/EVAL20.obj \ - -reference F3H-1/EVAL20_wtr.obj -computeVerticalError verticalError.txt \ + CompareOBJ -computeOptimalRotationAndTranslation -model EVAL20.obj \ + -reference EVAL20_wtr.obj -computeVerticalError verticalError.txt \ -saveOptimalShape optimal.obj -savePlateDiff plateDiff.txt -savePlateIndex plateIndex.txt The screen output is @@ -77,7 +77,7 @@ model for comparison: :: - ShapeFormatConverter -input Bennu/Bennu49k.obj -output BennuComparison.obj \ + ShapeFormatConverter -input Bennu49k.obj -output BennuComparison.obj \ -rotate 5,0,0,1 -translate 0.01,-0.01,0.01 This rotates the shape model by 5 degrees about the z axis and then translates @@ -94,11 +94,11 @@ Run CompareOBJ to find the optimal transform to align the comparison with the re CompareOBJ -computeOptimalRotationAndTranslation \ -model BennuComparison.obj \ - -reference Bennu/Bennu49k.obj \ - -computeVerticalError CompareOBJ/terrasaur-verticalError.txt \ - -saveOptimalShape CompareOBJ/terrasaur-optimal.obj \ - -savePlateDiff CompareOBJ/terrasaur-plateDiff.txt \ - -savePlateIndex CompareOBJ/terrasaur-plateIndex.txt + -reference Bennu49k.obj \ + -computeVerticalError terrasaur-verticalError.txt \ + -saveOptimalShape terrasaur-optimal.obj \ + -savePlateDiff terrasaur-plateDiff.txt \ + -savePlateIndex terrasaur-plateIndex.txt The screen output is diff --git a/doc/tools/PointCloudOverlap.rst b/doc/tools/PointCloudOverlap.rst new file mode 100644 index 0000000..e4f05a8 --- /dev/null +++ b/doc/tools/PointCloudOverlap.rst @@ -0,0 +1,38 @@ +.. _PointCloudOverlap: + +################# +PointCloudOverlap +################# + +***** +Usage +***** + +.. include:: ../toolDescriptions/PointCloudOverlap.txt + :literal: + +******** +Examples +******** + +Download the :download:`reference<./support_files/EVAL20_wtr.obj>` and :download:`comparison<./support_files/EVAL20.obj>` +shape models. You can view them in a tool such as +`ParaView `__. + +.. figure:: images/PointCloudOverlap_1.png + + This image shows the reference (pink) and input (blue) shape models. + +Run PointCloudOverlap: + +:: + + PointCloudOverlap -inputFile EVAL20.obj -referenceFile EVAL20_wtr.obj -outputFile overlap.vtk + +Note that OBJ is supported as an input format but not as an output format. + + +.. figure:: images/PointCloudOverlap_2.png + + The points in white are those in the input model that overlap the reference. + diff --git a/doc/tools/TriAx.rst b/doc/tools/TriAx.rst new file mode 100644 index 0000000..9fba689 --- /dev/null +++ b/doc/tools/TriAx.rst @@ -0,0 +1,43 @@ +.. _TriAx: + +===== +TriAx +===== + +TriAx is an implementation of the SPC tool TRIAX, which generates a triaxial ellipsoid in ICQ format. + +***** +Usage +***** + +.. include:: ../toolDescriptions/TriAx.txt + :literal: + +******* +Example +******* + +Generate an ellipsoid with dimensions 10, 8, 6, with q = 8. + +:: + + TriAx -A 10 -B 8 -C 6 -Q 8 -output triax.icq -saveOBJ + +The following ellipsoid is generated: + +.. container:: figures-row + + .. figure:: images/TriAx_X.png + :alt: looking down from the +X direction + + looking down from the +X direction + + .. figure:: images/TriAx_Y.png + :alt: looking down from the +Y direction + + looking down from the +Y direction + + .. figure:: images/TriAx_Z.png + :alt: looking down from the +Z direction + + looking down from the +Z direction \ No newline at end of file diff --git a/doc/tools/images/PointCloudOverlap_1.png b/doc/tools/images/PointCloudOverlap_1.png new file mode 100644 index 0000000..376b0ae Binary files /dev/null and b/doc/tools/images/PointCloudOverlap_1.png differ diff --git a/doc/tools/images/PointCloudOverlap_2.png b/doc/tools/images/PointCloudOverlap_2.png new file mode 100644 index 0000000..7f7a0b7 Binary files /dev/null and b/doc/tools/images/PointCloudOverlap_2.png differ diff --git a/doc/tools/images/TriAx_X.png b/doc/tools/images/TriAx_X.png new file mode 100644 index 0000000..f70b214 Binary files /dev/null and b/doc/tools/images/TriAx_X.png differ diff --git a/doc/tools/images/TriAx_Y.png b/doc/tools/images/TriAx_Y.png new file mode 100644 index 0000000..cc12a94 Binary files /dev/null and b/doc/tools/images/TriAx_Y.png differ diff --git a/doc/tools/images/TriAx_Z.png b/doc/tools/images/TriAx_Z.png new file mode 100644 index 0000000..7599c8e Binary files /dev/null and b/doc/tools/images/TriAx_Z.png differ diff --git a/mkPackage.bash b/mkPackage.bash index f61581a..f1655c8 100755 --- a/mkPackage.bash +++ b/mkPackage.bash @@ -89,7 +89,7 @@ function build_jar() { function make_scripts() { - classes=$(jar tf "${scriptPath}"/target/${packageName}.jar | grep $appSrcDir | grep -v '\$' | grep -v "package-info" | grep class) + classes=$(jar tf ${scriptPath}/target/${packageName}.jar | grep $appSrcDir | grep -v '\$' | grep -v "package-info" | grep -v "Immutable" | grep class) for class in $classes; do base=$(basename "$class" ".class") diff --git a/pom.xml b/pom.xml index 0da05f2..a2caa9e 100644 --- a/pom.xml +++ b/pom.xml @@ -1,5 +1,6 @@ + xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" + xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 https://maven.apache.org/xsd/maven-4.0.0.xsd"> 4.0.0 edu.jhuapl.ses.srn Terrasaur @@ -55,25 +56,29 @@ com.mycila license-maven-plugin 5.0.0 - - - -
com/mycila/maven/plugin/license/templates/MIT.txt
- - **/*.java - - - - - - - check-license - verify - - check - - - + + + +
com/mycila/maven/plugin/license/templates/MIT.txt
+ + **/*.java + + + 3rd-party/**/*.java + support-libraries/3rd-party/**/*.java + + + + + + + check-license + verify + + check + + + @@ -126,7 +131,8 @@ maven-surefire-plugin 3.5.3 - -Djava.library.path=${project.basedir}/3rd-party/${env.ARCH}/spice/JNISpice/lib + + -Djava.library.path=${project.basedir}/3rd-party/${env.ARCH}/spice/JNISpice/lib diff --git a/src/main/java/terrasaur/apps/AdjustShapeModelToOtherShapeModel.java b/src/main/java/terrasaur/apps/AdjustShapeModelToOtherShapeModel.java index 1312ff7..502a689 100644 --- a/src/main/java/terrasaur/apps/AdjustShapeModelToOtherShapeModel.java +++ b/src/main/java/terrasaur/apps/AdjustShapeModelToOtherShapeModel.java @@ -25,7 +25,6 @@ package terrasaur.apps; import java.io.File; import java.nio.charset.Charset; import java.util.*; - import org.apache.commons.cli.CommandLine; import org.apache.commons.cli.Option; import org.apache.commons.cli.Options; @@ -33,14 +32,14 @@ import org.apache.commons.io.FileUtils; import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; +import spice.basic.Plane; +import spice.basic.Vector3; import terrasaur.smallBodyModel.BoundingBox; import terrasaur.smallBodyModel.SmallBodyModel; import terrasaur.templates.TerrasaurTool; import terrasaur.utils.Log4j2Configurator; import terrasaur.utils.NativeLibraryLoader; import terrasaur.utils.PolyDataUtil; -import spice.basic.Plane; -import spice.basic.Vector3; import vtk.vtkGenericCell; import vtk.vtkPoints; import vtk.vtkPolyData; @@ -76,7 +75,7 @@ public class AdjustShapeModelToOtherShapeModel implements TerrasaurTool { return TerrasaurTool.super.fullDescription(options, header, ""); } - private static Options defineOptions() { + private static Options defineOptions() { Options options = TerrasaurTool.defineOptions(); options.addOption( Option.builder("from") @@ -182,6 +181,7 @@ public class AdjustShapeModelToOtherShapeModel implements TerrasaurTool { Vector3D origin = new Vector3D(0., 0., 0.); for (int i = 0; i < numberPoints; ++i) { points.GetPoint(i, p); + Vector3D thisPoint = new Vector3D(p); Vector3D lookDir; @@ -198,18 +198,20 @@ public class AdjustShapeModelToOtherShapeModel implements TerrasaurTool { } Vector3D lookPt = lookDir.scalarMultiply(diagonalLength); - lookPt = lookPt.add(origin); + lookPt = lookPt.add(thisPoint); List intersections = new ArrayList<>(); for (vtksbCellLocator cellLocator : cellLocators) { double[] intersectPoint = new double[3]; - // trace ray from the lookPt to the origin - first intersection is the farthest intersection - // from the origin + // trace ray from thisPoint to the lookPt - Assume cell intersection is the closest one if + // there are multiple? + // NOTE: result should return 1 in case of intersection but doesn't sometimes. + // Use the norm of intersection point to test for intersection instead. int result = cellLocator.IntersectWithLine( + thisPoint.toArray(), lookPt.toArray(), - origin.toArray(), tol, t, intersectPoint, @@ -219,38 +221,32 @@ public class AdjustShapeModelToOtherShapeModel implements TerrasaurTool { cell); Vector3D intersectVector = new Vector3D(intersectPoint); - if (fitPlane || localModel) { - // NOTE: result should return 1 in case of intersection but doesn't sometimes. - // Use the norm of intersection point to test for intersection instead. - - NavigableMap pointsMap = new TreeMap<>(); - if (intersectVector.getNorm() > 0) { - pointsMap.put(origin.subtract(intersectVector).getNorm(), intersectVector); - } - - lookPt = lookDir.scalarMultiply(-diagonalLength); - lookPt = lookPt.add(origin); - result = - cellLocator.IntersectWithLine( - lookPt.toArray(), - origin.toArray(), - tol, - t, - intersectPoint, - pcoords, - subId, - cell_id, - cell); - - intersectVector = new Vector3D(intersectPoint); - if (intersectVector.getNorm() > 0) { - pointsMap.put(origin.subtract(intersectVector).getNorm(), intersectVector); - } - - if (!pointsMap.isEmpty()) intersections.add(pointsMap.get(pointsMap.firstKey())); - } else { - if (result > 0) intersections.add(intersectVector); + NavigableMap pointsMap = new TreeMap<>(); + if (intersectVector.getNorm() > 0) { + pointsMap.put(thisPoint.subtract(intersectVector).getNorm(), intersectVector); } + + // look in the other direction + lookPt = lookDir.scalarMultiply(-diagonalLength); + lookPt = lookPt.add(thisPoint); + result = + cellLocator.IntersectWithLine( + thisPoint.toArray(), + lookPt.toArray(), + tol, + t, + intersectPoint, + pcoords, + subId, + cell_id, + cell); + + intersectVector = new Vector3D(intersectPoint); + if (intersectVector.getNorm() > 0) { + pointsMap.put(thisPoint.subtract(intersectVector).getNorm(), intersectVector); + } + + if (!pointsMap.isEmpty()) intersections.add(pointsMap.get(pointsMap.firstKey())); } if (intersections.isEmpty()) throw new Exception("Error: no intersections at all"); diff --git a/src/main/java/terrasaur/apps/CreateSBMTStructure.java b/src/main/java/terrasaur/apps/CreateSBMTStructure.java index be9ea69..f0f44e7 100644 --- a/src/main/java/terrasaur/apps/CreateSBMTStructure.java +++ b/src/main/java/terrasaur/apps/CreateSBMTStructure.java @@ -60,8 +60,7 @@ public class CreateSBMTStructure implements TerrasaurTool { @Override public String fullDescription(Options options) { - String header = - "This tool creates an SBMT ellipse file from a set of point on an image."; + String header = "This tool creates an SBMT ellipse file from a set of points on an image."; String footer = ""; return TerrasaurTool.super.fullDescription(options, header, footer); } @@ -129,20 +128,31 @@ public class CreateSBMTStructure implements TerrasaurTool { private static Options defineOptions() { Options options = TerrasaurTool.defineOptions(); + options.addOption( + Option.builder("flipX") + .desc("If present, negate the X coordinate of the input points.") + .build()); + options.addOption( + Option.builder("flipY") + .desc("If present, negate the Y coordinate of the input points.") + .build()); options.addOption( Option.builder("input") .required() .hasArg() .desc( - """ -Required. Name or input file. This is a text file with a pair of pixel coordinates per line. The pixel +""" +Required. Name or input file. This is a text file with a pair of pixel coordinates (X and Y) per line. The pixel coordinates are offsets from the image center. For example: # My test file -627.51274 876.11775 -630.53612 883.55992 -626.3499 881.46681 +89.6628 285.01 +97.8027 280.126 +95.0119 285.01 +-13.8299 323.616 +-1.9689 331.756 +-11.7367 330.826 Empty lines or lines beginning with # are ignored. @@ -161,11 +171,38 @@ axis and the third is a location for the semi-minor axis.""") .hasArg() .desc("Required. Name of output file.") .build()); + options.addOption( + Option.builder("spice") + .hasArg() + .desc( + "If present, name of metakernel to read. Other required options with -spice are -date, -observer, -target, and -cameraFrame.") + .build()); + options.addOption( + Option.builder("date") + .hasArgs() + .desc("Only used with -spice. Date of image (e.g. 2022 SEP 26 23:11:12.649).") + .build()); + options.addOption( + Option.builder("observer") + .hasArg() + .desc("Only used with -spice. Observing body (e.g. DART)") + .build()); + options.addOption( + Option.builder("target") + .hasArg() + .desc("Only used with -spice. Target body (e.g. DIMORPHOS).") + .build()); + options.addOption( + Option.builder("cameraFrame") + .hasArg() + .desc("Only used with -spice. Camera frame (e.g. DART_DRACO).") + .build()); options.addOption( Option.builder("sumFile") .required() .hasArg() - .desc("Required. Name of sum file to read.") + .desc( + "Required. Name of sum file to read. This is still required with -spice, but only used as a template to create a new sum file.") .build()); return options; } @@ -179,7 +216,7 @@ axis and the third is a location for the semi-minor axis.""") Map startupMessages = defaultOBJ.startupMessages(cl); for (MessageLabel ml : startupMessages.keySet()) - logger.info(String.format("%s %s", ml.label, startupMessages.get(ml))); + logger.info("{} {}", ml.label, startupMessages.get(ml)); NativeLibraryLoader.loadSpiceLibraries(); NativeLibraryLoader.loadVtkLibraries(); @@ -195,6 +232,8 @@ axis and the third is a location for the semi-minor axis.""") } RangeFromSumFile rfsf = new RangeFromSumFile(sumFile, polyData); + boolean flipX = cl.hasOption("flipX"); + boolean flipY = cl.hasOption("flipY"); List intercepts = new ArrayList<>(); List lines = FileUtils.readLines(new File(cl.getOptionValue("input")), Charset.defaultCharset()); @@ -204,6 +243,9 @@ axis and the third is a location for the semi-minor axis.""") int ix = (int) Math.round(Double.parseDouble(parts[0])); int iy = (int) Math.round(Double.parseDouble(parts[1])); + if (flipX) ix *= -1; + if (flipY) iy *= -1; + Map.Entry entry = rfsf.findIntercept(ix, iy); long cellID = entry.getKey(); if (cellID > -1) intercepts.add(entry.getValue()); @@ -216,12 +258,13 @@ axis and the third is a location for the semi-minor axis.""") // p1 and p2 define the long axis of the ellipse Vector3D p1 = intercepts.get(i); - Vector3D p2 = intercepts.get(i+1); + Vector3D p2 = intercepts.get(i + 1); // p3 lies on the short axis - Vector3D p3 = intercepts.get(i+2); + Vector3D p3 = intercepts.get(i + 2); - SBMTEllipseRecord record = createRecord(i/3, String.format("Ellipse %d", i/3), p1, p2, p3); + SBMTEllipseRecord record = + createRecord(i / 3, String.format("Ellipse %d", i / 3), p1, p2, p3); records.add(record); } diff --git a/src/main/java/terrasaur/apps/FacetInfo.java b/src/main/java/terrasaur/apps/FacetInfo.java new file mode 100644 index 0000000..b4d141f --- /dev/null +++ b/src/main/java/terrasaur/apps/FacetInfo.java @@ -0,0 +1,258 @@ +/* + * The MIT License + * Copyright © 2025 Johns Hopkins University Applied Physics Laboratory + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + * THE SOFTWARE. + */ +package terrasaur.apps; + +import java.io.PrintWriter; +import java.util.Arrays; +import java.util.Map; +import java.util.NavigableSet; +import java.util.TreeSet; +import org.apache.commons.cli.CommandLine; +import org.apache.commons.cli.Option; +import org.apache.commons.cli.Options; +import org.apache.commons.math3.geometry.euclidean.threed.SphericalCoordinates; +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.logging.log4j.LogManager; +import org.apache.logging.log4j.Logger; +import org.immutables.value.Value; +import terrasaur.templates.TerrasaurTool; +import terrasaur.utils.CellInfo; +import terrasaur.utils.NativeLibraryLoader; +import terrasaur.utils.PolyDataStatistics; +import terrasaur.utils.PolyDataUtil; +import vtk.vtkIdList; +import vtk.vtkOBBTree; +import vtk.vtkPolyData; + +public class FacetInfo implements TerrasaurTool { + + private static final Logger logger = LogManager.getLogger(); + + /** + * This doesn't need to be private, or even declared, but you might want to if you have other + * constructors. + */ + private FacetInfo() {} + + @Override + public String shortDescription() { + return "Print info about a facet."; + } + + @Override + public String fullDescription(Options options) { + String header = "Prints information about facet(s)."; + String footer = + """ + +This tool prints out facet center, normal, angle between center and +normal, and other information about the specified facet(s)."""; + + return TerrasaurTool.super.fullDescription(options, header, footer); + } + + private vtkPolyData polyData; + private vtkOBBTree searchTree; + private Vector3D origin; + + private FacetInfo(vtkPolyData polyData) { + this.polyData = polyData; + PolyDataStatistics stats = new PolyDataStatistics(polyData); + origin = new Vector3D(stats.getCentroid()); + + logger.info("Origin is at {}", origin); + + logger.info("Creating search tree"); + searchTree = new vtkOBBTree(); + searchTree.SetDataSet(polyData); + searchTree.SetTolerance(1e-12); + searchTree.BuildLocator(); + } + + /** + * @param cellId id of this cell + * @return Set of neighboring cells (ones which share a vertex with this one) + */ + private NavigableSet neighbors(long cellId) { + NavigableSet neighborCellIds = new TreeSet<>(); + + vtkIdList vertexIdlist = new vtkIdList(); + CellInfo.getCellInfo(polyData, cellId, vertexIdlist); + + vtkIdList facetIdlist = new vtkIdList(); + for (long i = 0; i < vertexIdlist.GetNumberOfIds(); i++) { + long vertexId = vertexIdlist.GetId(i); + polyData.GetPointCells(vertexId, facetIdlist); + } + for (long i = 0; i < facetIdlist.GetNumberOfIds(); i++) { + long id = facetIdlist.GetId(i); + if (id == cellId) continue; + neighborCellIds.add(id); + } + + return neighborCellIds; + } + + @Value.Immutable + public abstract static class FacetInfoLine { + + public abstract long index(); + + public abstract Vector3D radius(); + + public abstract Vector3D normal(); + + /** + * @return facets between this and origin + */ + public abstract NavigableSet interiorIntersections(); + + /** + * @return facets between this and infinity + */ + public abstract NavigableSet exteriorIntersections(); + + public static String getHeader() { + return "# Index, " + + "Center Lat (deg), " + + "Center Lon (deg), " + + "Radius, " + + "Radial Vector, " + + "Normal Vector, " + + "Angle between radial and normal (deg), " + + "facets between this and origin, " + + "facets between this and infinity"; + } + + public String toCSV() { + + SphericalCoordinates spc = new SphericalCoordinates(radius()); + + StringBuilder sb = new StringBuilder(); + + sb.append(String.format("%d, ", index())); + sb.append(String.format("%.4f, ", 90 - Math.toDegrees(spc.getPhi()))); + sb.append(String.format("%.4f, ", Math.toDegrees(spc.getTheta()))); + sb.append(String.format("%.6f, ", spc.getR())); + sb.append( + String.format("%.6f %.6f %.6f, ", radius().getX(), radius().getY(), radius().getZ())); + sb.append( + String.format("%.6f %.6f %.6f, ", normal().getX(), normal().getY(), normal().getZ())); + sb.append(String.format("%.3f, ", Math.toDegrees(Vector3D.angle(radius(), normal())))); + sb.append(String.format("%d", interiorIntersections().size())); + if (!interiorIntersections().isEmpty()) { + for (long id : interiorIntersections()) sb.append(String.format(" %d", id)); + } + sb.append(", "); + sb.append(String.format("%d", exteriorIntersections().size())); + if (!exteriorIntersections().isEmpty()) { + for (long id : exteriorIntersections()) sb.append(String.format(" %d", id)); + } + return sb.toString(); + } + } + + private FacetInfoLine getFacetInfoLine(long cellId) { + CellInfo ci = CellInfo.getCellInfo(polyData, cellId, new vtkIdList()); + + vtkIdList cellIds = new vtkIdList(); + searchTree.IntersectWithLine(origin.toArray(), ci.center().toArray(), null, cellIds); + + // count up all crossings of the surface between the origin and the facet. + NavigableSet insideIds = new TreeSet<>(); + for (long j = 0; j < cellIds.GetNumberOfIds(); j++) { + if (cellIds.GetId(j) == cellId) continue; + insideIds.add(cellIds.GetId(j)); + } + + Vector3D infinity = ci.center().scalarMultiply(1e9); + + cellIds = new vtkIdList(); + searchTree.IntersectWithLine(infinity.toArray(), ci.center().toArray(), null, cellIds); + + // count up all crossings of the surface between the infinity and the facet. + NavigableSet outsideIds = new TreeSet<>(); + for (long j = 0; j < cellIds.GetNumberOfIds(); j++) { + if (cellIds.GetId(j) == cellId) continue; + outsideIds.add(cellIds.GetId(j)); + } + + return ImmutableFacetInfoLine.builder() + .index(cellId) + .radius(ci.center()) + .normal(ci.normal()) + .interiorIntersections(insideIds) + .exteriorIntersections(outsideIds) + .build(); + } + + private static Options defineOptions() { + Options options = TerrasaurTool.defineOptions(); + options.addOption( + Option.builder("facet") + .required() + .hasArgs() + .desc("Facet(s) to query. Separate multiple indices with whitespace.") + .build()); + options.addOption( + Option.builder("obj").required().hasArg().desc("Shape model to validate.").build()); + options.addOption( + Option.builder("output").required().hasArg().desc("CSV file to write.").build()); + return options; + } + + public static void main(String[] args) { + TerrasaurTool defaultOBJ = new FacetInfo(); + + Options options = defineOptions(); + + CommandLine cl = defaultOBJ.parseArgs(args, options); + + Map startupMessages = defaultOBJ.startupMessages(cl); + for (MessageLabel ml : startupMessages.keySet()) + logger.info("{} {}", ml.label, startupMessages.get(ml)); + + NativeLibraryLoader.loadVtkLibraries(); + + try { + vtkPolyData polydata = PolyDataUtil.loadShapeModel(cl.getOptionValue("obj")); + FacetInfo app = new FacetInfo(polydata); + try (PrintWriter pw = new PrintWriter(cl.getOptionValue("output"))) { + pw.println(FacetInfoLine.getHeader()); + for (long cellId : + Arrays.stream(cl.getOptionValues("facet")).mapToLong(Long::parseLong).toArray()) { + pw.println(app.getFacetInfoLine(cellId).toCSV()); + + NavigableSet neighbors = app.neighbors(cellId); + for (long neighborCellId : neighbors) + pw.println(app.getFacetInfoLine(neighborCellId).toCSV()); + } + } + logger.info("Wrote {}", cl.getOptionValue("output")); + } catch (Exception e) { + logger.error(e.getLocalizedMessage(), e); + } + + logger.info("Finished."); + } +} diff --git a/src/main/java/terrasaur/apps/PointCloudFormatConverter.java b/src/main/java/terrasaur/apps/PointCloudFormatConverter.java index cfb7954..b53246b 100644 --- a/src/main/java/terrasaur/apps/PointCloudFormatConverter.java +++ b/src/main/java/terrasaur/apps/PointCloudFormatConverter.java @@ -305,8 +305,8 @@ public class PointCloudFormatConverter implements TerrasaurTool { } } else { if (halfSize < 0 || groundSampleDistance < 0) { - System.out.printf( - "Must supply -halfSize and -groundSampleDistance for %s output\n", + logger.error( + "Must supply -halfSize and -groundSampleDistance for {} output", outFormat); return; } diff --git a/src/main/java/terrasaur/apps/PointCloudOverlap.java b/src/main/java/terrasaur/apps/PointCloudOverlap.java new file mode 100644 index 0000000..911885c --- /dev/null +++ b/src/main/java/terrasaur/apps/PointCloudOverlap.java @@ -0,0 +1,417 @@ +/* + * The MIT License + * Copyright © 2025 Johns Hopkins University Applied Physics Laboratory + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + * THE SOFTWARE. + */ +package terrasaur.apps; + +import java.awt.geom.Path2D; +import java.awt.geom.Point2D; +import java.io.File; +import java.io.IOException; +import java.util.ArrayList; +import java.util.Collection; +import java.util.List; +import java.util.Map; +import nom.tam.fits.FitsException; +import org.apache.commons.cli.CommandLine; +import org.apache.commons.cli.Option; +import org.apache.commons.cli.Options; +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; +import org.apache.commons.math3.geometry.euclidean.twod.hull.MonotoneChain; +import org.apache.logging.log4j.LogManager; +import org.apache.logging.log4j.Logger; +import picante.math.coords.LatitudinalVector; +import spice.basic.SpiceException; +import spice.basic.Vector3; +import terrasaur.enums.FORMATS; +import terrasaur.templates.TerrasaurTool; +import terrasaur.utils.NativeLibraryLoader; +import terrasaur.utils.VectorStatistics; +import terrasaur.utils.tessellation.StereographicProjection; +import vtk.vtkPoints; + +public class PointCloudOverlap implements TerrasaurTool { + + private static final Logger logger = LogManager.getLogger(); + + @Override + public String shortDescription() { + return "Find points in a point cloud which overlap a reference point cloud."; + } + + @Override + public String fullDescription(Options options) { + String header = ""; + String footer = + "\nThis program finds all points in the input point cloud which overlap the points in a reference point cloud.\n\n" + + "Supported input formats are ASCII, BINARY, L2, OBJ, and VTK. Supported output formats are ASCII, BINARY, L2, and VTK. " + + "ASCII format is white spaced delimited x y z coordinates. BINARY files must contain double precision x y z coordinates.\n\n" + + "A plane is fit to the reference point cloud and all points in each cloud are projected onto this plane. Any point in the " + + "projected input cloud which falls within the outline of the projected reference cloud is considered to be overlapping."; + return TerrasaurTool.super.fullDescription(options, header, footer); + } + + private Path2D.Double polygon; + + private StereographicProjection proj; + + /*- Useful for debugging + private vtkPoints ref2DPoints; + private vtkPoints input2DPoints; + private vtkPolyData polygonPolyData; + private vtkCellArray polygonCells; + private vtkPoints polygonPoints; + private vtkDoubleArray polygonSuccessArray; + */ + + public PointCloudOverlap(Collection refPoints) { + if (refPoints != null) { + VectorStatistics vStats = new VectorStatistics(); + for (double[] pt : refPoints) vStats.add(new Vector3(pt)); + + Vector3D centerXYZ = vStats.getMean(); + + proj = + new StereographicProjection( + new LatitudinalVector(1, centerXYZ.getDelta(), centerXYZ.getAlpha())); + + createRefPolygon(refPoints); + } + } + + public StereographicProjection getProjection() { + return proj; + } + + private void createRefPolygon(Collection refPoints) { + List stereographicPoints = new ArrayList<>(); + for (double[] refPt : refPoints) { + Vector3D point3D = new Vector3D(refPt); + Point2D point = proj.forward(point3D.getDelta(), point3D.getAlpha()); + stereographicPoints.add(new Vector2D(point.getX(), point.getY())); + } + + /*- + ref2DPoints = new vtkPoints(); + input2DPoints = new vtkPoints(); + polygonPolyData = new vtkPolyData(); + polygonCells = new vtkCellArray(); + polygonPoints = new vtkPoints(); + polygonSuccessArray = new vtkDoubleArray(); + polygonSuccessArray.SetName("success") + polygonPolyData.SetPoints(polygonPoints); + polygonPolyData.SetLines(polygonCells); + polygonPolyData.GetCellData().AddArray(polygonSuccessArray); + + for (Vector2D refPoint : refPoints) + ref2DPoints.InsertNextPoint(refPoint.getX(), refPoint.getY(), 0); + */ + + MonotoneChain mc = new MonotoneChain(); + List vertices = new ArrayList<>(mc.findHullVertices(stereographicPoints)); + /*- + for (int i = 1; i < vertices.size(); i++) { + Vector2D lastPt = vertices.get(i - 1); + Vector2D thisPt = vertices.get(i); + System.out.printf("%f %f to %f %f distance %f\n", lastPt.getX(), lastPt.getY(), thisPt.getX(), + thisPt.getY(), lastPt.distance(thisPt)); + } + Vector2D lastPt = vertices.get(vertices.size() - 1); + Vector2D thisPt = vertices.get(0); + System.out.printf("%f %f to %f %f distance %f\n", lastPt.getX(), lastPt.getY(), thisPt.getX(), + thisPt.getY(), lastPt.distance(thisPt)); + */ + // int id0 = 0; + for (Vector2D vertex : vertices) { + // int id1 = polygonPoints.InsertNextPoint(vertex.getX(), vertex.getY(), 0); + + if (polygon == null) { + polygon = new Path2D.Double(); + polygon.moveTo(vertex.getX(), vertex.getY()); + } else { + polygon.lineTo(vertex.getX(), vertex.getY()); + /*- + vtkLine line = new vtkLine(); + line.GetPointIds().SetId(0, id0); + line.GetPointIds().SetId(1, id1); + polygonCells.InsertNextCell(line); + */ + } + // id0 = id1; + } + polygon.closePath(); + + /*- + vtkPolyDataWriter writer = new vtkPolyDataWriter(); + writer.SetInputData(polygonPolyData); + writer.SetFileName("polygon2D.vtk"); + writer.SetFileTypeToBinary(); + writer.Update(); + + writer = new vtkPolyDataWriter(); + polygonPolyData = new vtkPolyData(); + polygonPolyData.SetPoints(ref2DPoints); + writer.SetInputData(polygonPolyData); + writer.SetFileName("refPoints.vtk"); + writer.SetFileTypeToBinary(); + writer.Update(); + */ + + } + + public boolean isEnclosed(double[] xyz) { + Vector3D point = new Vector3D(xyz); + Point2D projected = proj.forward(point.getDelta(), point.getAlpha()); + return polygon.contains(projected.getX(), projected.getY()); + } + + /** + * @param inputPoints points to consider + * @param scale scale factor + * @return indices of all points inside the scaled polygon + */ + public List scalePoints(List inputPoints, double scale) { + + List projected = new ArrayList<>(); + for (double[] inputPoint : inputPoints) { + Vector3D point = new Vector3D(inputPoint); + Point2D projectedPoint = proj.forward(point.getDelta(), point.getAlpha()); + projected.add(new Vector2D(projectedPoint.getX(), projectedPoint.getY())); + } + + Vector2D center = new Vector2D(0, 0); + for (Vector2D inputPoint : projected) center = center.add(inputPoint); + + center = center.scalarMultiply(1. / inputPoints.size()); + + List translatedPoints = new ArrayList<>(); + for (Vector2D inputPoint : projected) translatedPoints.add(inputPoint.subtract(center)); + + Path2D.Double thisPolygon = null; + MonotoneChain mc = new MonotoneChain(); + Collection vertices = mc.findHullVertices(translatedPoints); + for (Vector2D vertex : vertices) { + if (thisPolygon == null) { + thisPolygon = new Path2D.Double(); + thisPolygon.moveTo(scale * vertex.getX(), scale * vertex.getY()); + } else { + thisPolygon.lineTo(scale * vertex.getX(), scale * vertex.getY()); + } + } + thisPolygon.closePath(); + + List indices = new ArrayList<>(); + for (int i = 0; i < projected.size(); i++) { + Vector2D inputPoint = projected.get(i); + if (thisPolygon.contains( + inputPoint.getX() - center.getX(), inputPoint.getY() - center.getY())) indices.add(i); + } + return indices; + } + + private static Options defineOptions() { + Options options = new Options(); + options.addOption( + Option.builder("inputFormat") + .hasArg() + .desc( + "Format of input file. If not present format will be inferred from file extension.") + .build()); + options.addOption( + Option.builder("inputFile") + .required() + .hasArg() + .desc("Required. Name of input file.") + .build()); + options.addOption( + Option.builder("inllr") + .desc( + "If present, input values are assumed to be lon, lat, rad. Default is x, y, z. Only used with ASCII or BINARY formats.") + .build()); + options.addOption( + Option.builder("referenceFormat") + .hasArg() + .desc( + "Format of reference file. If not present format will be inferred from file extension.") + .build()); + options.addOption( + Option.builder("referenceFile") + .required() + .hasArg() + .desc("Required. Name of reference file.") + .build()); + options.addOption( + Option.builder("refllr") + .desc( + "If present, reference values are assumed to be lon, lat, rad. Default is x, y, z. Only used with ASCII or BINARY formats.") + .build()); + options.addOption( + Option.builder("outputFormat") + .hasArg() + .desc( + "Format of output file. If not present format will be inferred from file extension.") + .build()); + options.addOption( + Option.builder("outputFile") + .required() + .hasArg() + .desc("Required. Name of output file.") + .build()); + options.addOption( + Option.builder("outllr") + .desc( + "If present, output values will be lon, lat, rad. Default is x, y, z. Only used with ASCII or BINARY formats.") + .build()); + options.addOption( + Option.builder("scale") + .hasArg() + .desc("Value to scale bounding box containing intersect region. Default is 1.0.") + .build()); + return options; + } + + public static void main(String[] args) + throws SpiceException, IOException, InterruptedException, FitsException { + + NativeLibraryLoader.loadVtkLibraries(); + NativeLibraryLoader.loadSpiceLibraries(); + + TerrasaurTool defaultOBJ = new PointCloudOverlap(null); + + Options options = defineOptions(); + + CommandLine cl = defaultOBJ.parseArgs(args, options); + + Map startupMessages = defaultOBJ.startupMessages(cl); + for (MessageLabel ml : startupMessages.keySet()) + logger.info(String.format("%s %s", ml.label, startupMessages.get(ml))); + + // Read the reference file + FORMATS refFormat = + cl.hasOption("referenceFormat") + ? FORMATS.valueOf(cl.getOptionValue("referenceFormat").toUpperCase()) + : FORMATS.formatFromExtension(cl.getOptionValue("referenceFile")); + String refFile = cl.getOptionValue("referenceFile"); + boolean refLLR = cl.hasOption("refllr"); + + PointCloudFormatConverter pcfc = new PointCloudFormatConverter(refFormat, FORMATS.VTK); + pcfc.read(refFile, refLLR); + vtkPoints referencePoints = pcfc.getPoints(); + logger.info("{} points read from {}", referencePoints.GetNumberOfPoints(), refFile); + + List refPts = new ArrayList<>(); + for (int i = 0; i < referencePoints.GetNumberOfPoints(); i++) { + refPts.add(referencePoints.GetPoint(i)); + } + + // create the overlap object and set the enclosing polygon + PointCloudOverlap pco = new PointCloudOverlap(refPts); + + // Read the input point cloud + FORMATS inFormat = + cl.hasOption("inputFormat") + ? FORMATS.valueOf(cl.getOptionValue("inputFormat").toUpperCase()) + : FORMATS.formatFromExtension(cl.getOptionValue("inputFile")); + String inFile = cl.getOptionValue("inputFile"); + boolean inLLR = cl.hasOption("inllr"); + + pcfc = new PointCloudFormatConverter(inFormat, FORMATS.VTK); + pcfc.read(inFile, inLLR); + vtkPoints inputPoints = pcfc.getPoints(); + logger.info("{} points read from {}", inputPoints.GetNumberOfPoints(), inFile); + + List enclosedIndices = new ArrayList<>(); + for (int i = 0; i < inputPoints.GetNumberOfPoints(); i++) { + double[] pt = inputPoints.GetPoint(i); + if (pco.isEnclosed(pt)) enclosedIndices.add(i); + } + + if (cl.hasOption("scale")) { + List pts = new ArrayList<>(); + for (Integer i : enclosedIndices) pts.add(inputPoints.GetPoint(i)); + + // this list includes which of the enclosed points are inside the scaled polygon + List theseIndices = pco.scalePoints(pts, Double.parseDouble(cl.getOptionValue("scale"))); + + // now relate this list back to the original list of points + List newIndices = new ArrayList<>(); + for (Integer i : theseIndices) newIndices.add(enclosedIndices.get(i)); + enclosedIndices = newIndices; + } + + VectorStatistics xyzStats = new VectorStatistics(); + VectorStatistics xyStats = new VectorStatistics(); + for (Integer i : enclosedIndices) { + double[] thisPt = inputPoints.GetPoint(i); + Vector3D thisPt3D = new Vector3D(thisPt); + xyzStats.add(thisPt3D); + Point2D projectedPt = pco.getProjection().forward(thisPt3D.getDelta(), thisPt3D.getAlpha()); + xyStats.add(new Vector3(projectedPt.getX(), projectedPt.getY(), 0)); + } + + logger.info( + "Center XYZ: {}, {}, {}", + xyzStats.getMean().getX(), xyzStats.getMean().getY(), xyzStats.getMean().getZ()); + Vector3D centerXYZ = xyzStats.getMean(); + logger.info( + "Center lon, lat: {}, {}\n", + Math.toDegrees(centerXYZ.getAlpha()), Math.toDegrees(centerXYZ.getDelta())); + logger.info( + "xmin/xmax/extent: {}/{}/{}\n", + xyzStats.getMin().getX(), + xyzStats.getMax().getX(), + xyzStats.getMax().getX() - xyzStats.getMin().getX()); + logger.info( + "ymin/ymax/extent: {}/{}/{}\n", + xyzStats.getMin().getY(), + xyzStats.getMax().getY(), + xyzStats.getMax().getY() - xyzStats.getMin().getY()); + logger.info( + "zmin/zmax/extent: {}/{}/{}\n", + xyzStats.getMin().getZ(), + xyzStats.getMax().getZ(), + xyzStats.getMax().getZ() - xyzStats.getMin().getZ()); + + FORMATS outFormat = + cl.hasOption("outputFormat") + ? FORMATS.valueOf(cl.getOptionValue("outputFormat").toUpperCase()) + : FORMATS.formatFromExtension(cl.getOptionValue("outputFile")); + + vtkPoints pointsToWrite = new vtkPoints(); + for (Integer i : enclosedIndices) pointsToWrite.InsertNextPoint(inputPoints.GetPoint(i)); + pcfc = new PointCloudFormatConverter(FORMATS.VTK, outFormat); + pcfc.setPoints(pointsToWrite); + String outputFilename = cl.getOptionValue("outputFile"); + pcfc.write(outputFilename, cl.hasOption("outllr")); + if (new File(outputFilename).exists()) { + logger.info( + "{} points written to {}", + pointsToWrite.GetNumberOfPoints(), outputFilename); + } else { + logger.error("Could not write {}", outputFilename); + } + + logger.info("Finished"); + } +} + +// TODO write out center of output pointcloud diff --git a/src/main/java/terrasaur/apps/RenderShapeFromSumFile.java b/src/main/java/terrasaur/apps/RenderShapeFromSumFile.java index dc0e03b..529cb3e 100644 --- a/src/main/java/terrasaur/apps/RenderShapeFromSumFile.java +++ b/src/main/java/terrasaur/apps/RenderShapeFromSumFile.java @@ -313,7 +313,7 @@ public class RenderShapeFromSumFile implements TerrasaurTool { logger.printf( Level.DEBUG, "Thread %d lat/lon %.2f/%.2f, %s, sum %f, cells %d, %.2f", - Thread.currentThread().getId(), + Thread.currentThread().threadId(), Math.toDegrees(intersectPoint.getDelta()), Math.toDegrees(intersectPoint.getAlpha()), intersectPoint.toString(), @@ -338,7 +338,7 @@ public class RenderShapeFromSumFile implements TerrasaurTool { @Override public Map call() throws Exception { - logger.info("Thread {}: starting", Thread.currentThread().getId()); + logger.info("Thread {}: starting", Thread.currentThread().threadId()); int xPixels = subPixel * nPixelsX; @@ -380,7 +380,7 @@ public class RenderShapeFromSumFile implements TerrasaurTool { logger.debug( String.format( "Thread %d: No intersection with local model for pixel (%d,%d): lat/lon %.2f/%.2f, using global intersection %d %s", - Thread.currentThread().getId(), + Thread.currentThread().threadId(), i, j, Math.toDegrees(intersectPt3D.getDelta()), @@ -397,7 +397,7 @@ public class RenderShapeFromSumFile implements TerrasaurTool { } } - logger.info("Thread {}: finished", Thread.currentThread().getId()); + logger.info("Thread {}: finished", Thread.currentThread().threadId()); return brightness; } diff --git a/src/main/java/terrasaur/apps/TriAx.java b/src/main/java/terrasaur/apps/TriAx.java new file mode 100644 index 0000000..7523347 --- /dev/null +++ b/src/main/java/terrasaur/apps/TriAx.java @@ -0,0 +1,166 @@ +/* + * The MIT License + * Copyright © 2025 Johns Hopkins University Applied Physics Laboratory + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + * THE SOFTWARE. + */ +package terrasaur.apps; + +import java.io.*; +import java.util.Map; +import org.apache.commons.cli.CommandLine; +import org.apache.commons.cli.Option; +import org.apache.commons.cli.Options; +import org.apache.commons.io.FilenameUtils; +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.logging.log4j.LogManager; +import org.apache.logging.log4j.Logger; +import terrasaur.templates.TerrasaurTool; +import terrasaur.utils.ICQUtils; +import terrasaur.utils.NativeLibraryLoader; +import terrasaur.utils.PolyDataUtil; +import vtk.vtkPolyData; + +public class TriAx implements TerrasaurTool { + + private static final Logger logger = LogManager.getLogger(); + + private TriAx() {} + + @Override + public String shortDescription() { + return "Generate a triaxial ellipsoid in ICQ format."; + } + + @Override + public String fullDescription(Options options) { + + String footer = + "\nTriAx is an implementation of the SPC tool TRIAX, which generates a triaxial ellipsoid in ICQ format.\n"; + return TerrasaurTool.super.fullDescription(options, "", footer); + } + + static Options defineOptions() { + Options options = TerrasaurTool.defineOptions(); + options.addOption( + Option.builder("A") + .required() + .hasArg() + .desc("Long axis of the ellipsoid, arbitrary units (usually assumed to be km).") + .build()); + options.addOption( + Option.builder("B") + .required() + .hasArg() + .desc("Medium axis of the ellipsoid, arbitrary units (usually assumed to be km).") + .build()); + options.addOption( + Option.builder("C") + .required() + .hasArg() + .desc("Short axis of the ellipsoid, arbitrary units (usually assumed to be km).") + .build()); + options.addOption( + Option.builder("Q") + .required() + .hasArg() + .desc("ICQ size parameter. This is conventionally but not necessarily a power of 2.") + .build()); + options.addOption( + Option.builder("saveOBJ") + .desc( + "If present, save in OBJ format as well. " + + "File will have the same name as ICQ file with an OBJ extension.") + .build()); + options.addOption( + Option.builder("output").hasArg().required().desc("Name of ICQ file to write.").build()); + + return options; + } + + static final int MAX_Q = 512; + + public static void main(String[] args) { + TerrasaurTool defaultOBJ = new TriAx(); + + Options options = defineOptions(); + + CommandLine cl = defaultOBJ.parseArgs(args, options); + + Map startupMessages = defaultOBJ.startupMessages(cl); + for (MessageLabel ml : startupMessages.keySet()) + logger.info("{} {}", ml.label, startupMessages.get(ml)); + + int q = Integer.parseInt(cl.getOptionValue("Q")); + String shapefile = cl.getOptionValue("output"); + + double[] ax = new double[3]; + ax[0] = Double.parseDouble(cl.getOptionValue("A")); + ax[1] = Double.parseDouble(cl.getOptionValue("B")); + ax[2] = Double.parseDouble(cl.getOptionValue("C")); + + double[][][][] vec = new double[3][MAX_Q + 1][MAX_Q + 1][6]; + for (int f = 0; f < 6; f++) { + for (int i = 0; i <= q; i++) { + for (int j = 0; j <= q; j++) { + + double[] u = ICQUtils.xyf2u(q, i, j, f, ax); + double z = + 1 + / Math.sqrt( + Math.pow(u[0] / ax[0], 2) + + Math.pow(u[1] / ax[1], 2) + + Math.pow(u[2] / ax[2], 2)); + + double[] v = new Vector3D(u).scalarMultiply(z).toArray(); + for (int k = 0; k < 3; k++) { + vec[k][i][j][f] = v[k]; + } + } + } + } + + ICQUtils.writeICQ(q, vec, shapefile); + + if (cl.hasOption("saveOBJ")) { + + String basename = FilenameUtils.getBaseName(shapefile); + String dirname = FilenameUtils.getFullPath(shapefile); + if (dirname.isEmpty()) dirname = "."; + File obj = new File(dirname, basename + ".obj"); + + NativeLibraryLoader.loadVtkLibraries(); + try { + vtkPolyData polydata = PolyDataUtil.loadShapeModel(shapefile); + if (polydata == null) { + logger.error("Cannot read {}", shapefile); + System.exit(0); + } + + polydata = PolyDataUtil.removeDuplicatePoints(polydata); + polydata = PolyDataUtil.removeUnreferencedPoints(polydata); + polydata = PolyDataUtil.removeZeroAreaFacets(polydata); + + PolyDataUtil.saveShapeModelAsOBJ(polydata, obj.getPath()); + } catch (Exception e) { + logger.error(e); + } + } + } +} diff --git a/src/main/java/terrasaur/apps/ValidateNormals.java b/src/main/java/terrasaur/apps/ValidateNormals.java index 8fc85fb..6397eee 100644 --- a/src/main/java/terrasaur/apps/ValidateNormals.java +++ b/src/main/java/terrasaur/apps/ValidateNormals.java @@ -30,7 +30,6 @@ import java.util.concurrent.Callable; import java.util.concurrent.ExecutorService; import java.util.concurrent.Executors; import java.util.concurrent.Future; - import org.apache.commons.cli.CommandLine; import org.apache.commons.cli.Option; import org.apache.commons.cli.Options; @@ -65,6 +64,12 @@ public class ValidateNormals implements TerrasaurTool { static Options defineOptions() { Options options = TerrasaurTool.defineOptions(); + options.addOption( + Option.builder("fast") + .desc("If present, only check for overhangs if center and normal point in opposite " + + "directions. Default behavior is to always check for intersections between body center " + + "and facet center.") + .build()); options.addOption( Option.builder("origin") .hasArg() @@ -126,10 +131,12 @@ public class ValidateNormals implements TerrasaurTool { private final long index0; private final long index1; + private final boolean fast; - public FlippedNormalFinder(long index0, long index1) { + public FlippedNormalFinder(long index0, long index1, boolean fast) { this.index0 = index0; this.index1 = index1; + this.fast = fast; } @Override @@ -158,20 +165,21 @@ public class ValidateNormals implements TerrasaurTool { long index = index0 + i; CellInfo ci = CellInfo.getCellInfo(polyData, index, idList); - getOBBTree().IntersectWithLine(origin, ci.center().toArray(), null, cellIds); + boolean isOpposite = (ci.center().dotProduct(ci.normal()) < 0); - // count up all crossings of the surface between the origin and the facet. int numCrossings = 0; - for (int j = 0; j < cellIds.GetNumberOfIds(); j++) { - if (cellIds.GetId(j) == index) break; - numCrossings++; + if (isOpposite || !fast) { + // count up all crossings of the surface between the origin and the facet. + getOBBTree().IntersectWithLine(origin, ci.center().toArray(), null, cellIds); + for (int j = 0; j < cellIds.GetNumberOfIds(); j++) { + if (cellIds.GetId(j) == index) break; + numCrossings++; + } } // if numCrossings is even, the radial and normal should point in the same direction. If it - // is odd, the - // radial and normal should point in opposite directions. + // is odd, the radial and normal should point in opposite directions. boolean shouldBeOpposite = (numCrossings % 2 == 1); - boolean isOpposite = (ci.center().dotProduct(ci.normal()) < 0); // XOR operator - true if both conditions are different if (isOpposite ^ shouldBeOpposite) flippedNormals.add(index); @@ -208,7 +216,7 @@ public class ValidateNormals implements TerrasaurTool { Map startupMessages = defaultOBJ.startupMessages(cl); for (MessageLabel ml : startupMessages.keySet()) - logger.info(String.format("%s %s", ml.label, startupMessages.get(ml))); + logger.info("{} {}", ml.label, startupMessages.get(ml)); NativeLibraryLoader.loadVtkLibraries(); @@ -234,6 +242,7 @@ public class ValidateNormals implements TerrasaurTool { Set flippedNormals = new HashSet<>(); + boolean fast = cl.hasOption("fast"); int numThreads = cl.hasOption("numThreads") ? Integer.parseInt(cl.getOptionValue("numThreads")) : 1; try (ExecutorService executor = Executors.newFixedThreadPool(numThreads)) { @@ -244,7 +253,7 @@ public class ValidateNormals implements TerrasaurTool { long fromIndex = i * numFacets; long toIndex = Math.min(polyData.GetNumberOfCells(), fromIndex + numFacets); - FlippedNormalFinder fnf = app.new FlippedNormalFinder(fromIndex, toIndex); + FlippedNormalFinder fnf = app.new FlippedNormalFinder(fromIndex, toIndex, fast); futures.add(executor.submit(fnf)); } diff --git a/src/main/java/terrasaur/apps/ValidateOBJ.java b/src/main/java/terrasaur/apps/ValidateOBJ.java index 8e4c8c1..d720480 100644 --- a/src/main/java/terrasaur/apps/ValidateOBJ.java +++ b/src/main/java/terrasaur/apps/ValidateOBJ.java @@ -308,7 +308,7 @@ public class ValidateOBJ implements TerrasaurTool { TriangularFacet facet = new TriangularFacet(new VectorIJK(pt0), new VectorIJK(pt1), new VectorIJK(pt2)); if (facet.getArea() > 0) { - stats.addValue(facet.getCenter().getDot(facet.getNormal())); + stats.addValue(facet.getCenter().createUnitized().getDot(facet.getNormal())); cStats.add(facet.getCenter()); nStats.add(facet.getNormal()); } @@ -374,7 +374,7 @@ public class ValidateOBJ implements TerrasaurTool { Map startupMessages = defaultOBJ.startupMessages(cl); for (MessageLabel ml : startupMessages.keySet()) - logger.info(String.format("%s %s", ml.label, startupMessages.get(ml))); + logger.info("{} {}", ml.label, startupMessages.get(ml)); NativeLibraryLoader.loadVtkLibraries(); vtkPolyData polyData = PolyDataUtil.loadShapeModel(cl.getOptionValue("obj")); @@ -418,7 +418,7 @@ public class ValidateOBJ implements TerrasaurTool { if (cl.hasOption("output")) { PolyDataUtil.saveShapeModelAsOBJ(polyData, cl.getOptionValue("output")); - logger.info(String.format("Wrote OBJ file %s", cl.getOptionValue("output"))); + logger.info("Wrote OBJ file {}", cl.getOptionValue("output")); } } } diff --git a/src/main/java/terrasaur/utils/CellInfo.java b/src/main/java/terrasaur/utils/CellInfo.java index be5ea78..19818be 100644 --- a/src/main/java/terrasaur/utils/CellInfo.java +++ b/src/main/java/terrasaur/utils/CellInfo.java @@ -129,7 +129,7 @@ public abstract class CellInfo { polydata.GetPoint(idList.GetId(2), pt2); } - private static CellInfo fromPoints(double[] pt0, double[] pt1, double[] pt2) { + public static CellInfo fromPoints(double[] pt0, double[] pt1, double[] pt2) { ImmutableCellInfo.Builder builder = ImmutableCellInfo.builder(); builder.pt0(new Vector3D(pt0)); diff --git a/src/main/java/terrasaur/utils/ICQUtils.java b/src/main/java/terrasaur/utils/ICQUtils.java new file mode 100644 index 0000000..3deb20a --- /dev/null +++ b/src/main/java/terrasaur/utils/ICQUtils.java @@ -0,0 +1,111 @@ +/* + * The MIT License + * Copyright © 2025 Johns Hopkins University Applied Physics Laboratory + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + * THE SOFTWARE. + */ +package terrasaur.utils; + +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.logging.log4j.LogManager; +import org.apache.logging.log4j.Logger; + +import java.io.FileWriter; +import java.io.IOException; +import java.io.PrintWriter; + +public class ICQUtils { + private static final Logger logger = LogManager.getLogger(); + + public static double[] xyf2u(int q, double x, double y, int face, double[] ax) { + double pi = Math.acos(-1.0); + double[] dt = new double[3]; + double[] v = new double[3]; + int[][][] u = defu(); + + dt[0] = Math.tan((2 * x / q - 1) * pi / 4); + dt[1] = Math.tan((2 * y / q - 1) * pi / 4); + dt[2] = 1 / Math.sqrt(1 + dt[0] * dt[0] + dt[1] * dt[1]); + dt[0] *= dt[2]; + dt[1] *= dt[2]; + + for (int k = 0; k < 3; k++) { + for (int j = 0; j < 3; j++) { + int sign = u[j][k][face]; + if (sign != 0) { + v[k] = dt[j] * sign * Math.sqrt(ax[k]); + } + } + } + + return new Vector3D(v).normalize().toArray(); + } + + static int[][][] defu() { + int[][][] u = new int[3][3][6]; + + u[0][0][0] = 1; + u[1][1][0] = -1; + u[2][2][0] = 1; + + u[0][0][1] = 1; + u[1][2][1] = -1; + u[2][1][1] = -1; + + u[0][1][2] = -1; + u[1][2][2] = -1; + u[2][0][2] = -1; + + u[0][0][3] = -1; + u[1][2][3] = -1; + u[2][1][3] = 1; + + u[0][1][4] = 1; + u[1][2][4] = -1; + u[2][0][4] = 1; + + u[0][0][5] = 1; + u[1][1][5] = 1; + u[2][2][5] = -1; + + return u; + } + + public static void writeICQ(int q, double[][][][] vec, String filename){ + + try (PrintWriter out = new PrintWriter(new FileWriter(filename))) { + out.println(q); + for (int f = 0; f < 6; f++) { + for (int j = 0; j <= q; j++) { + for (int i = 0; i <= q; i++) { + for (int k = 0; k < 3; k++) { + out.printf("%12.5f", vec[k][i][j][f]); + } + out.println(); + } + } + } + } catch (IOException e) { + logger.error(e); + } + + } + + +} diff --git a/src/main/java/terrasaur/utils/PolyDataUtil.java b/src/main/java/terrasaur/utils/PolyDataUtil.java index 0c7e331..8071490 100644 --- a/src/main/java/terrasaur/utils/PolyDataUtil.java +++ b/src/main/java/terrasaur/utils/PolyDataUtil.java @@ -521,17 +521,17 @@ public class PolyDataUtil { n[q][q][4] = n[0][q][3]; n[q][q][5] = n[0][q][3]; - try (PrintWriter pw = new PrintWriter("java.txt")) { - for (int f = 0; f < 6; f++) { - for (int i = 0; i < q; i++) { - for (int j = 0; j < q; j++) { - pw.printf("%3d%3d%3d%6d\n", i, j, f + 1, n[i][j][f]); - } - } - } - } catch (FileNotFoundException e) { - logger.error(e.getLocalizedMessage()); - } +// try (PrintWriter pw = new PrintWriter("java.txt")) { +// for (int f = 0; f < 6; f++) { +// for (int i = 0; i < q; i++) { +// for (int j = 0; j < q; j++) { +// pw.printf("%3d%3d%3d%6d\n", i, j, f + 1, n[i][j][f]); +// } +// } +// } +// } catch (FileNotFoundException e) { +// logger.error(e.getLocalizedMessage()); +// } pltLines.add(String.format("%d", 12 * q * q)); n0 = 0; diff --git a/src/main/java/terrasaur/utils/SumFile.java b/src/main/java/terrasaur/utils/SumFile.java index f96370e..b4fda39 100644 --- a/src/main/java/terrasaur/utils/SumFile.java +++ b/src/main/java/terrasaur/utils/SumFile.java @@ -29,32 +29,35 @@ import java.nio.charset.Charset; import java.util.ArrayList; import java.util.Collections; import java.util.List; -import javax.annotation.Nullable; +import net.jafama.FastMath; import org.apache.commons.io.FileUtils; import org.apache.commons.math3.geometry.euclidean.threed.Rotation; import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; import org.immutables.value.Value; -import net.jafama.FastMath; -import terrasaur.utils.ImmutableSumFile.Builder; -import terrasaur.utils.math.MathConversions; +import picante.math.vectorspace.RotationMatrixIJK; +import picante.math.vectorspace.VectorIJK; +import picante.mechanics.*; +import picante.mechanics.providers.aberrated.AberrationCorrection; import spice.basic.Plane; import spice.basic.Ray; import spice.basic.RayPlaneIntercept; import spice.basic.SpiceException; import spice.basic.Vector3; +import terrasaur.utils.ImmutableSumFile.Builder; +import terrasaur.utils.math.MathConversions; +import terrasaur.utils.spice.SpiceBundle; /** * Class describing Bob Gaskell's sum file object - * - * @author nairah1 * + * @author nairah1 */ @Value.Immutable public abstract class SumFile { - private final static Logger logger = LogManager.getLogger(SumFile.class); + private static final Logger logger = LogManager.getLogger(SumFile.class); // line 1 public abstract String picnm(); @@ -107,23 +110,63 @@ public abstract class SumFile { // line 13 public abstract Vector3D sig_ptg(); - @Nullable public abstract Vector3D frustum1(); - @Nullable public abstract Vector3D frustum2(); - @Nullable public abstract Vector3D frustum3(); - @Nullable public abstract Vector3D frustum4(); /** - * - * @param translation - * @param rotation - * @return A SumFile with the scobj transformed and the cx, cy, cz, and sz vectors rotated. + * @param bundle SPICE bundle + * @param observer spacecraft + * @param target target + * @param cameraFrame Camera Frame + * @param t time + * @return SumFile with scobj, C matrix, and sun direction evaluated using SPICE at time t + */ + public SumFile fromSpice( + SpiceBundle bundle, EphemerisID observer, EphemerisID target, FrameID cameraFrame, double t) { + FrameID bodyFixedFrame = bundle.getBodyFixedFrame(target); + VectorIJK scobj = + bundle + .getAbProvider() + .createAberratedPositionVectorFunction( + target, observer, bodyFixedFrame, Coverage.ALL_TIME, AberrationCorrection.LT_S) + .getPosition(t); + RotationMatrixIJK cMatrix = + bundle + .getAbProvider() + .getFrameProvider() + .createFrameTransformFunction(cameraFrame, bodyFixedFrame, Coverage.ALL_TIME) + .getTransform(t); + VectorIJK sunDir = + bundle + .getAbProvider() + .createAberratedPositionVectorFunction( + CelestialBodies.SUN, + target, + bodyFixedFrame, + Coverage.ALL_TIME, + AberrationCorrection.LT_S) + .getPosition(t) + .unitize(); + + Builder b = ImmutableSumFile.builder().from(this); + b.utcString(bundle.getTimeConversion().tdbToUTCString(t, "C")); + b.scobj(MathConversions.toVector3D(scobj)); + b.cx(MathConversions.toVector3D(cMatrix.mxv(VectorIJK.I))); + b.cy(MathConversions.toVector3D(cMatrix.mxv(VectorIJK.J))); + b.cz(MathConversions.toVector3D(cMatrix.mxv(VectorIJK.K))); + b.sz(MathConversions.toVector3D(cMatrix.mxv(sunDir))); + return b.build(); + } + + /** + * @param translation translation to apply to scobj + * @param rotation rotation to apply to C matrix and sun direction + * @return SumFile with the scobj transformed and the cx, cy, cz, and sz vectors rotated. */ public SumFile transform(Vector3D translation, Rotation rotation) { Builder b = ImmutableSumFile.builder().from(this); @@ -141,10 +184,8 @@ public abstract class SumFile { } /** - * Construct a SumFile from the input file - * - * @param file - * @return + * @param file file to read + * @return SumFile */ public static SumFile fromFile(File file) { SumFile s = null; @@ -157,10 +198,8 @@ public abstract class SumFile { } /** - * Construct a SumFile from the input string list - * - * @param lines - * @return + * @param lines lines to read + * @return SumFile */ public static SumFile fromLines(List lines) { Builder b = ImmutableSumFile.builder(); @@ -180,95 +219,114 @@ public abstract class SumFile { parts = lines.get(4).trim().split("\\s+"); double[] tmp = new double[3]; - for (int i = 0; i < 3; i++) - tmp[i] = parseFortranDouble(parts[i]); + for (int i = 0; i < 3; i++) tmp[i] = parseFortranDouble(parts[i]); b.scobj(new Vector3D(tmp)); parts = lines.get(5).trim().split("\\s+"); tmp = new double[3]; - for (int i = 0; i < 3; i++) - tmp[i] = parseFortranDouble(parts[i]); + for (int i = 0; i < 3; i++) tmp[i] = parseFortranDouble(parts[i]); b.cx(new Vector3D(tmp).normalize()); parts = lines.get(6).trim().split("\\s+"); tmp = new double[3]; - for (int i = 0; i < 3; i++) - tmp[i] = parseFortranDouble(parts[i]); + for (int i = 0; i < 3; i++) tmp[i] = parseFortranDouble(parts[i]); b.cy(new Vector3D(tmp).normalize()); parts = lines.get(7).trim().split("\\s+"); tmp = new double[3]; - for (int i = 0; i < 3; i++) - tmp[i] = parseFortranDouble(parts[i]); + for (int i = 0; i < 3; i++) tmp[i] = parseFortranDouble(parts[i]); b.cz(new Vector3D(tmp).normalize()); parts = lines.get(8).trim().split("\\s+"); tmp = new double[3]; - for (int i = 0; i < 3; i++) - tmp[i] = parseFortranDouble(parts[i]); + for (int i = 0; i < 3; i++) tmp[i] = parseFortranDouble(parts[i]); b.sz(new Vector3D(tmp).normalize()); parts = lines.get(9).trim().split("\\s+"); tmp = new double[3]; - for (int i = 0; i < 3; i++) - tmp[i] = parseFortranDouble(parts[i]); + for (int i = 0; i < 3; i++) tmp[i] = parseFortranDouble(parts[i]); b.kmat1(new Vector3D(tmp)); tmp = new double[3]; - for (int i = 0; i < 3; i++) - tmp[i] = parseFortranDouble(parts[i + 3]); + for (int i = 0; i < 3; i++) tmp[i] = parseFortranDouble(parts[i + 3]); b.kmat2(new Vector3D(tmp)); parts = lines.get(10).trim().split("\\s+"); - for (int i = 0; i < 4; i++) - b.addDistortion(parseFortranDouble(parts[i])); + for (int i = 0; i < 4; i++) b.addDistortion(parseFortranDouble(parts[i])); parts = lines.get(11).trim().split("\\s+"); tmp = new double[3]; - for (int i = 0; i < 3; i++) - tmp[i] = parseFortranDouble(parts[i]); + for (int i = 0; i < 3; i++) tmp[i] = parseFortranDouble(parts[i]); b.sig_vso(new Vector3D(tmp)); parts = lines.get(11).trim().split("\\s+"); tmp = new double[3]; - for (int i = 0; i < 3; i++) - tmp[i] = parseFortranDouble(parts[i]); + for (int i = 0; i < 3; i++) tmp[i] = parseFortranDouble(parts[i]); b.sig_ptg(new Vector3D(tmp)); + b.frustum1(Vector3D.ZERO); + b.frustum2(Vector3D.ZERO); + b.frustum3(Vector3D.ZERO); + b.frustum4(Vector3D.ZERO); + SumFile s = b.build(); double fov1 = Math.abs(Math.atan(s.npx() / (2.0 * s.mmfl() * s.kmat1().getX()))); double fov2 = Math.abs(Math.atan(s.nln() / (2.0 * s.mmfl() * s.kmat2().getY()))); Vector3D cornerVector = new Vector3D(-FastMath.tan(fov1), -FastMath.tan(fov2), 1.0); - double fx = cornerVector.getX() * s.cx().getX() + cornerVector.getY() * s.cy().getX() - + cornerVector.getZ() * s.cz().getX(); - double fy = cornerVector.getX() * s.cx().getY() + cornerVector.getY() * s.cy().getY() - + cornerVector.getZ() * s.cz().getY(); - double fz = cornerVector.getX() * s.cx().getZ() + cornerVector.getY() * s.cy().getZ() - + cornerVector.getZ() * s.cz().getZ(); + double fx = + cornerVector.getX() * s.cx().getX() + + cornerVector.getY() * s.cy().getX() + + cornerVector.getZ() * s.cz().getX(); + double fy = + cornerVector.getX() * s.cx().getY() + + cornerVector.getY() * s.cy().getY() + + cornerVector.getZ() * s.cz().getY(); + double fz = + cornerVector.getX() * s.cx().getZ() + + cornerVector.getY() * s.cy().getZ() + + cornerVector.getZ() * s.cz().getZ(); b.frustum3(new Vector3D(fx, fy, fz).normalize()); - fx = -cornerVector.getX() * s.cx().getX() + cornerVector.getY() * s.cy().getX() - + cornerVector.getZ() * s.cz().getX(); - fy = -cornerVector.getX() * s.cx().getY() + cornerVector.getY() * s.cy().getY() - + cornerVector.getZ() * s.cz().getY(); - fz = -cornerVector.getX() * s.cx().getZ() + cornerVector.getY() * s.cy().getZ() - + cornerVector.getZ() * s.cz().getZ(); + fx = + -cornerVector.getX() * s.cx().getX() + + cornerVector.getY() * s.cy().getX() + + cornerVector.getZ() * s.cz().getX(); + fy = + -cornerVector.getX() * s.cx().getY() + + cornerVector.getY() * s.cy().getY() + + cornerVector.getZ() * s.cz().getY(); + fz = + -cornerVector.getX() * s.cx().getZ() + + cornerVector.getY() * s.cy().getZ() + + cornerVector.getZ() * s.cz().getZ(); b.frustum4(new Vector3D(fx, fy, fz).normalize()); - fx = cornerVector.getX() * s.cx().getX() - cornerVector.getY() * s.cy().getX() - + cornerVector.getZ() * s.cz().getX(); - fy = cornerVector.getX() * s.cx().getY() - cornerVector.getY() * s.cy().getY() - + cornerVector.getZ() * s.cz().getY(); - fz = cornerVector.getX() * s.cx().getZ() - cornerVector.getY() * s.cy().getZ() - + cornerVector.getZ() * s.cz().getZ(); + fx = + cornerVector.getX() * s.cx().getX() + - cornerVector.getY() * s.cy().getX() + + cornerVector.getZ() * s.cz().getX(); + fy = + cornerVector.getX() * s.cx().getY() + - cornerVector.getY() * s.cy().getY() + + cornerVector.getZ() * s.cz().getY(); + fz = + cornerVector.getX() * s.cx().getZ() + - cornerVector.getY() * s.cy().getZ() + + cornerVector.getZ() * s.cz().getZ(); b.frustum1(new Vector3D(fx, fy, fz).normalize()); - fx = -cornerVector.getX() * s.cx().getX() - cornerVector.getY() * s.cy().getX() - + cornerVector.getZ() * s.cz().getX(); - fy = -cornerVector.getX() * s.cx().getY() - cornerVector.getY() * s.cy().getY() - + cornerVector.getZ() * s.cz().getY(); - fz = -cornerVector.getX() * s.cx().getZ() - cornerVector.getY() * s.cy().getZ() - + cornerVector.getZ() * s.cz().getZ(); + fx = + -cornerVector.getX() * s.cx().getX() + - cornerVector.getY() * s.cy().getX() + + cornerVector.getZ() * s.cz().getX(); + fy = + -cornerVector.getX() * s.cx().getY() + - cornerVector.getY() * s.cy().getY() + + cornerVector.getZ() * s.cz().getY(); + fz = + -cornerVector.getX() * s.cx().getZ() + - cornerVector.getY() * s.cy().getZ() + + cornerVector.getZ() * s.cz().getZ(); b.frustum2(new Vector3D(fx, fy, fz).normalize()); return b.build(); @@ -277,17 +335,15 @@ public abstract class SumFile { /** * Account for numbers of the form .1192696009D+03 rather than .1192696009E+03 (i.e. a D instead * of an E). This function replaces D's with E's. - * - * @param s - * @return + * + * @param s String + * @return input string with all instances of 'D' replaced with 'E' */ private static double parseFortranDouble(String s) { return Double.parseDouble(s.replace('D', 'E')); } - /** - * Write the sum file object to a string - */ + /** Write the sum file object to a string */ @Override public String toString() { StringBuilder sb = new StringBuilder(); @@ -297,25 +353,52 @@ public abstract class SumFile { String.format("%6d%6d%6d%6d%56s\n", npx(), nln(), t1(), t2(), " NPX, NLN, THRSH ")); sb.append( String.format("%20.10e%20.10e%20.10e%20s\n", mmfl(), px0(), ln0(), " MMFL, CTR ")); - sb.append(String.format("%20.10e%20.10e%20.10e%20s\n", scobj().getX(), scobj().getY(), - scobj().getZ(), " SCOBJ ")); - sb.append(String.format("%20.10e%20.10e%20.10e%20s\n", cx().getX(), cx().getY(), cx().getZ(), - " CX ")); - sb.append(String.format("%20.10e%20.10e%20.10e%20s\n", cy().getX(), cy().getY(), cy().getZ(), - " CY ")); - sb.append(String.format("%20.10e%20.10e%20.10e%20s\n", cz().getX(), cz().getY(), cz().getZ(), - " CZ ")); - sb.append(String.format("%20.10e%20.10e%20.10e%20s\n", sz().getX(), sz().getY(), sz().getZ(), - " SZ ")); - sb.append(String.format("%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%20s\n", kmat1().getX(), - kmat1().getY(), kmat1().getZ(), kmat2().getX(), kmat2().getY(), kmat2().getZ(), - " K-MATRIX ")); - sb.append(String.format("%15.5f%15.5f%15.5f%15.5f%20s\n", distortion().get(0), - distortion().get(1), distortion().get(2), distortion().get(3), " DISTORTION ")); - sb.append(String.format("%20.10e%20.10e%20.10e%20s\n", sig_vso().getX(), sig_vso().getY(), - sig_vso().getZ(), " SIGMA_VSO ")); - sb.append(String.format("%20.10e%20.10e%20.10e%20s\n", sig_ptg().getX(), sig_ptg().getY(), - sig_ptg().getZ(), " SIGMA_PTG ")); + sb.append( + String.format( + "%20.10e%20.10e%20.10e%20s\n", + scobj().getX(), scobj().getY(), scobj().getZ(), " SCOBJ ")); + sb.append( + String.format( + "%20.10e%20.10e%20.10e%20s\n", + cx().getX(), cx().getY(), cx().getZ(), " CX ")); + sb.append( + String.format( + "%20.10e%20.10e%20.10e%20s\n", + cy().getX(), cy().getY(), cy().getZ(), " CY ")); + sb.append( + String.format( + "%20.10e%20.10e%20.10e%20s\n", + cz().getX(), cz().getY(), cz().getZ(), " CZ ")); + sb.append( + String.format( + "%20.10e%20.10e%20.10e%20s\n", + sz().getX(), sz().getY(), sz().getZ(), " SZ ")); + sb.append( + String.format( + "%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%20s\n", + kmat1().getX(), + kmat1().getY(), + kmat1().getZ(), + kmat2().getX(), + kmat2().getY(), + kmat2().getZ(), + " K-MATRIX ")); + sb.append( + String.format( + "%15.5f%15.5f%15.5f%15.5f%20s\n", + distortion().get(0), + distortion().get(1), + distortion().get(2), + distortion().get(3), + " DISTORTION ")); + sb.append( + String.format( + "%20.10e%20.10e%20.10e%20s\n", + sig_vso().getX(), sig_vso().getY(), sig_vso().getZ(), " SIGMA_VSO ")); + sb.append( + String.format( + "%20.10e%20.10e%20.10e%20s\n", + sig_ptg().getX(), sig_ptg().getY(), sig_ptg().getZ(), " SIGMA_PTG ")); sb.append(String.format("%s\n", "LANDMARKS")); sb.append(String.format("%s\n", "LIMB FITS")); sb.append(String.format("%s\n", "END FILE")); @@ -324,95 +407,77 @@ public abstract class SumFile { } /** - * Boresight direction. This is the same as {@link #cz()}. - * - * @return + * @return Boresight direction. This is the same as {@link #cz()}. */ public Vector3D boresight() { return cz(); } /** - * Sun direction. This is the same as {@link #sz()}. - * - * @return + * @return Sun direction. This is the same as {@link #sz()}. */ public Vector3D sunDirection() { return sz(); } /** - * 
-   * return new Vector3D(1. / npx(), frustum3().subtract(frustum4()))
-   *
- * * @return + * 
new Vector3D(1. / npx(), frustum3().subtract(frustum4()))
*/ public Vector3D xPerPixel() { return new Vector3D(1. / npx(), frustum3().subtract(frustum4())); } /** - * 
-   * return new Vector3D(1. / nln(), frustum3().subtract(frustum1()));
-   *
- * * @return + * 
 new Vector3D(1. / nln(), frustum3().subtract(frustum1()));
+   *
*/ public Vector3D yPerPixel() { return new Vector3D(1. / nln(), frustum3().subtract(frustum1())); } /** - * Angular size per pixel in the X direction. - * - * @return + * @return Angular size per pixel in the X direction calculated using + * 
Vector3D.angle(frustum3(), frustum4()) / npx()
*/ public double horizontalResolution() { return Vector3D.angle(frustum3(), frustum4()) / npx(); } /** - * Angular size per pixel in the Y direction. - * - * @return + * @return Angular size per pixel in the Y direction calculated using + * 
Vector3D.angle(frustum3(), frustum1()) / nln()
*/ public double verticalResolution() { return Vector3D.angle(frustum3(), frustum1()) / nln(); } /** - * Height of the image in pixels. - * - * @return + * @return Height of the image in pixels. */ public int imageHeight() { double kmatrix00 = Math.abs(kmat1().getX()); double kmatrix11 = Math.abs(kmat2().getY()); int imageHeight = nln(); - if (kmatrix00 > kmatrix11) - imageHeight = (int) Math.round(nln() * (kmatrix00 / kmatrix11)); + if (kmatrix00 > kmatrix11) imageHeight = (int) Math.round(nln() * (kmatrix00 / kmatrix11)); return imageHeight; } /** - * Width of the image in pixels. - * - * @return + * @return Width of the image in pixels. */ public int imageWidth() { double kmatrix00 = Math.abs(kmat1().getX()); double kmatrix11 = Math.abs(kmat2().getY()); int imageWidth = npx(); - if (kmatrix11 > kmatrix00) - imageWidth = (int) Math.round(npx() * (kmatrix11 / kmatrix00)); + if (kmatrix11 > kmatrix00) imageWidth = (int) Math.round(npx() * (kmatrix11 / kmatrix00)); return imageWidth; } /** - * Get the rotation to convert from body fixed coordinates to camera coordinates. {@link #cx()}, - * {@link #cx()}, {@link #cz()} are the rows of this matrix. - * + * @return rotation to convert from body fixed coordinates to camera coordinates. {@link #cx()}, + * {@link #cx()}, {@link #cz()} are the rows of this matrix. */ public Rotation getBodyFixedToCamera() { double[][] m = new double[3][]; @@ -423,9 +488,8 @@ public abstract class SumFile { } /** - * * @param directions from the spacecraft, in the body fixed frame. - * @return + * @return list of booleans set to true if direction is in the field of view, or false if not * @throws SpiceException */ public List isInFOV(List directions) throws SpiceException { @@ -450,10 +514,10 @@ public abstract class SumFile { } Path2D.Double shape = new Path2D.Double(); - shape.moveTo(xAxis.dot(points.get(0)), yAxis.dot(points.get(0))); + shape.moveTo(xAxis.dot(points.getFirst()), yAxis.dot(points.getFirst())); for (int i = 1; i < points.size(); i++) shape.lineTo(xAxis.dot(points.get(i)), yAxis.dot(points.get(i))); - shape.lineTo(xAxis.dot(points.get(0)), yAxis.dot(points.get(0))); + shape.lineTo(xAxis.dot(points.getFirst()), yAxis.dot(points.getFirst())); for (Vector3 direction : directions) { RayPlaneIntercept rpi = @@ -464,5 +528,4 @@ public abstract class SumFile { return Collections.unmodifiableList(isInFOV); } - } diff --git a/src/main/java/terrasaur/utils/tessellation/FibonacciSphere.java b/src/main/java/terrasaur/utils/tessellation/FibonacciSphere.java index 7285c4f..11f8ef6 100644 --- a/src/main/java/terrasaur/utils/tessellation/FibonacciSphere.java +++ b/src/main/java/terrasaur/utils/tessellation/FibonacciSphere.java @@ -80,6 +80,7 @@ public class FibonacciSphere implements SphericalTessellation { } /** + * Get statistics on the distances between each point and its closest neighbor * * @return statistics on the distances between each point and its closest neighbor */ @@ -181,9 +182,10 @@ public class FibonacciSphere implements SphericalTessellation { } /** + * Get the nearest point from the desired location. * * @param lv input location - * @return key is distance to tile center in radians, value is tile index + * @return key is distance in radians, value is point index */ public Map.Entry getNearest(LatitudinalVector lv) { return getNearest(CoordConverters.convert(lv)); @@ -198,7 +200,7 @@ public class FibonacciSphere implements SphericalTessellation { * * * @param ijk cartesian coordinates - * @return key is distance to tile center in radians, value is tile index + * @return key is distance in radians, value is point index */ public Map.Entry getNearest(UnwritableVectorIJK ijk) { final long n = getNumTiles(); @@ -248,7 +250,7 @@ public class FibonacciSphere implements SphericalTessellation { } } - return new AbstractMap.SimpleEntry<>(Math.sqrt(d), j); + return new AbstractMap.SimpleEntry(Math.sqrt(d), j); } private MatrixIJ getLocalToGlobalTransform(UnwritableVectorIJK p) { @@ -282,8 +284,8 @@ public class FibonacciSphere implements SphericalTessellation { /** * - * @param i tile index - * @return distance to this tile's closest neighbor center + * @param i point index + * @return distance to this point's closest neighbor */ public Double getDist(int i) { return getClosestNeighborDistance().get(i); @@ -292,7 +294,7 @@ public class FibonacciSphere implements SphericalTessellation { /** * * @param lv input point - * @return map of tile indices sorted by distance from the input point + * @return map of points sorted by distance from the input point */ public NavigableMap getDistanceMap(LatitudinalVector lv) { UnwritableVectorIJK ijk = CoordConverters.convert(lv); diff --git a/src/main/java/terrasaur/utils/tessellation/StereographicProjection.java b/src/main/java/terrasaur/utils/tessellation/StereographicProjection.java index fd82753..73e0ae6 100644 --- a/src/main/java/terrasaur/utils/tessellation/StereographicProjection.java +++ b/src/main/java/terrasaur/utils/tessellation/StereographicProjection.java @@ -29,13 +29,12 @@ import picante.math.vectorspace.UnwritableVectorIJK; import java.awt.geom.Point2D; /** - * Implement a stereographic projection. This is package private so that there's no conflict with - * anything in the cartography package. Based on Snyder (1987). + * Implement a stereographic projection. Based on Snyder (1987). * * @author nairah1 * */ -class StereographicProjection { +public class StereographicProjection { private final double centerLat, centerLon; private final double sinCenterLat, cosCenterLat; @@ -47,7 +46,7 @@ class StereographicProjection { * * @param center projection center */ - StereographicProjection(LatitudinalVector center) { + public StereographicProjection(LatitudinalVector center) { this(1.0, center); } @@ -58,7 +57,7 @@ class StereographicProjection { * @param R scale * @param center projection center */ - StereographicProjection(double R, LatitudinalVector center) { + public StereographicProjection(double R, LatitudinalVector center) { this.R = R; this.centerLat = center.getLatitude(); this.centerLon = center.getLongitude(); diff --git a/src/test/java/terrasaur/utils/SumFileTest.java b/src/test/java/terrasaur/utils/SumFileTest.java index 2bc525d..400d434 100644 --- a/src/test/java/terrasaur/utils/SumFileTest.java +++ b/src/test/java/terrasaur/utils/SumFileTest.java @@ -23,13 +23,20 @@ package terrasaur.utils; import static org.junit.Assert.assertEquals; + import java.io.File; import java.io.IOException; import java.nio.charset.Charset; import java.util.List; import org.apache.commons.io.FileUtils; import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.junit.Assert; +import org.junit.Ignore; import org.junit.Test; +import picante.mechanics.EphemerisID; +import picante.spice.fov.FOV; +import picante.spice.fov.FOVFactory; +import terrasaur.utils.spice.SpiceBundle; public class SumFileTest { @@ -38,19 +45,63 @@ public class SumFileTest { File file = ResourceUtils.writeResourceToFile("/M605862153F5.SUM"); + Assert.assertNotNull(file); List lines = FileUtils.readLines(file, Charset.defaultCharset()); SumFile sumFile = SumFile.fromLines(lines); - assertEquals(Vector3D.angle(sumFile.frustum1(), - new Vector3D(0.4395120989622179, 0.898199076601014, -0.008217886523401116)), 0, 1E-10); - assertEquals(Vector3D.angle(sumFile.frustum2(), - new Vector3D(0.43798867001719116, 0.8968954485311167, 0.06119215097329732)), 0, 1E-10); - assertEquals(Vector3D.angle(sumFile.frustum3(), - new Vector3D(0.3761137519676121, 0.926529032684429, -0.009077288896014253)), 0, 1E-10); - assertEquals(Vector3D.angle(sumFile.frustum4(), - new Vector3D(0.37459032302258627, 0.9252254046145307, 0.060332748600675515)), 0, 1E-10); - + assertEquals( + 0, + Vector3D.angle( + sumFile.frustum1(), + new Vector3D(0.4395120989622179, 0.898199076601014, -0.008217886523401116)), + 1E-10); + assertEquals( + 0, + Vector3D.angle( + sumFile.frustum2(), + new Vector3D(0.43798867001719116, 0.8968954485311167, 0.06119215097329732)), + 1E-10); + assertEquals( + 0, + Vector3D.angle( + sumFile.frustum3(), + new Vector3D(0.3761137519676121, 0.926529032684429, -0.009077288896014253)), + 1E-10); + assertEquals( + 0, + Vector3D.angle( + sumFile.frustum4(), + new Vector3D(0.37459032302258627, 0.9252254046145307, 0.060332748600675515)), + 1E-10); } + @Ignore + @Test + public void testSPICE() { + SumFile sumFile1 = + SumFile.fromFile( + new File( + "/project/sis/users/nairah1/SBCLT/tickets/52-create-sbmt-structure-file/spice/D717505942G0.SUM")); + + SpiceBundle bundle = + new SpiceBundle.Builder() + .addMetakernels(List.of("/project/dart/data/SPICE/dra/mk/d520_v02_N0066.tm")) + .build(); + EphemerisID observer = bundle.getObject("DART"); + EphemerisID target = bundle.getObject("DIMORPHOS"); + FOV fov = new FOVFactory(bundle.getKernelPool()).create(-135101); + + double t = bundle.getTimeConversion().utcStringToTDB("2022 SEP 26 23:14:23.328"); + + SumFile sumFile2 = + SumFile.fromFile( + new File( + "/project/sis/users/nairah1/SBCLT/tickets/52-create-sbmt-structure-file/spice/D717506133G0.SUM")); + SumFile sumFile3 = sumFile1.fromSpice(bundle, observer, target, fov.getFrameID(), t); + + System.out.println(sumFile1); + System.out.println(sumFile2); + System.out.println(sumFile3); + } }