From 3b389bb7f77b8846edd137787d664030d70f44b5 Mon Sep 17 00:00:00 2001 From: tiffanymatthe <49164666+tiffanymatthe@users.noreply.github.com> Date: Sat, 25 Feb 2023 16:48:23 -0800 Subject: [PATCH 01/20] CompareMeshes: Added operation to compare meshes Adds functionality to compare meshes based on Hausdorff distance and surface area --- src/Operation_Dispatcher.cc | 2 + src/Operations/CMakeLists.txt | 1 + src/Operations/CompareMeshes.cc | 120 ++++++++++++++++++++++++++++++++ src/Operations/CompareMeshes.h | 16 +++++ 4 files changed, 139 insertions(+) create mode 100644 src/Operations/CompareMeshes.cc create mode 100644 src/Operations/CompareMeshes.h diff --git a/src/Operation_Dispatcher.cc b/src/Operation_Dispatcher.cc index 2d37c02f..7790ad4b 100644 --- a/src/Operation_Dispatcher.cc +++ b/src/Operation_Dispatcher.cc @@ -37,6 +37,7 @@ #include "Operations/BuildLexiconInteractively.h" #include "Operations/ClusterDBSCAN.h" #include "Operations/CombineMeshes.h" +#include "Operations/CompareMeshes.h" #include "Operations/ComparePixels.h" #include "Operations/ContourBasedRayCastDoseAccumulate.h" #include "Operations/ContourSimilarity.h" @@ -280,6 +281,7 @@ std::map Known_Operations(){ out["CellularAutomata"] = std::make_pair(OpArgDocCellularAutomata, CellularAutomata); out["ClusterDBSCAN"] = std::make_pair(OpArgDocClusterDBSCAN, ClusterDBSCAN); out["CombineMeshes"] = std::make_pair(OpArgDocCombineMeshes, CombineMeshes); + out["CompareMeshes"] = std::make_pair(OpArgDocCompareMeshes, CompareMeshes); out["ComparePixels"] = std::make_pair(OpArgDocComparePixels, ComparePixels); out["ContourBasedRayCastDoseAccumulate"] = std::make_pair(OpArgDocContourBasedRayCastDoseAccumulate, ContourBasedRayCastDoseAccumulate); out["ContouringAides"] = std::make_pair(OpArgDocContouringAides, ContouringAides); diff --git a/src/Operations/CMakeLists.txt b/src/Operations/CMakeLists.txt index 3492e102..8f86992b 100644 --- a/src/Operations/CMakeLists.txt +++ b/src/Operations/CMakeLists.txt @@ -16,6 +16,7 @@ add_library(Operations_objs OBJECT BuildLexiconInteractively.cc ClusterDBSCAN.cc CombineMeshes.cc + CompareMeshes.cc ComparePixels.cc ContourBasedRayCastDoseAccumulate.cc ContouringAides.cc diff --git a/src/Operations/CompareMeshes.cc b/src/Operations/CompareMeshes.cc new file mode 100644 index 00000000..5ea59abb --- /dev/null +++ b/src/Operations/CompareMeshes.cc @@ -0,0 +1,120 @@ +#include +#include +#include //Needed for exit() calls. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "YgorMath.h" //Needed for vec3 class. +#include "YgorMathIOOBJ.h" +#include "YgorMisc.h" //Needed for FUNCINFO, FUNCWARN, FUNCERR macros. + +#include "Explicator.h" + +#include "../Contour_Boolean_Operations.h" +#include "../Structs.h" +#include "../Regex_Selectors.h" + +#include "../Simple_Meshing.h" +#include "../Surface_Meshes.h" + +#include "CompareMeshes.h" + +OperationDoc OpArgDocCompareMeshes(){ + OperationDoc out; + out.name = "CompareMeshes"; + + out.desc = + "This routine calculates various metrics of difference between two meshes and prints it to the terminal output." + ; + + + out.args.emplace_back(); + out.args.back() = SMWhitelistOpArgDoc(); + out.args.back().name = "MeshSelection1"; + out.args.back().default_val = "#-0"; + out.args.emplace_back(); + out.args.back() = SMWhitelistOpArgDoc(); + out.args.back().name = "MeshSelection2"; + out.args.back().default_val = "#-1"; + + return out; +} + +bool CompareMeshes(Drover &DICOM_data, + const OperationArgPkg& OptArgs, + std::map& /*InvocationMetadata*/, + const std::string& FilenameLex){ + + Explicator X(FilenameLex); + + //---------------------------------------------- User Parameters -------------------------------------------------- + const auto MeshSelection1Str = OptArgs.getValueStr("MeshSelection1").value(); + const auto MeshSelection2Str = OptArgs.getValueStr("MeshSelection2").value(); + + //----------------------------------------------------------------------------------------------------------------- + auto SMs_all = All_SMs( DICOM_data ); + auto SMs1 = Whitelist( SMs_all, MeshSelection1Str ); + auto SMs2 = Whitelist( SMs_all, MeshSelection2Str ); + + + if (SMs1.size() < 1 || SMs2.size() < 1) { + FUNCERR("Must select at least two meshes."); + } + if (SMs1.size() > 1 || SMs2.size() > 1) { + FUNCWARN("Can only calculate the metrics of two meshes, only looking at last element of each selected"); + } + + std::shared_ptr mesh1 = *SMs1.front(); + std::shared_ptr mesh2 = *SMs2.front(); + + double max_distance = 0; + + FUNCINFO("Iterating through " << size(mesh1->meshes.vertices) << " and " << size(mesh2->meshes.vertices) << " vertices."); + + // O(n^2) time, might be way to speed this up if this is ever a bottleneck + for (auto & vertex : mesh1->meshes.vertices) { + // find closest vertex on mesh2 that corresponds to vertex + double min_distance = -1; + for (auto & vertex2 : mesh2->meshes.vertices) { + double distance = vertex.distance(vertex2); + if (min_distance == -1) { + min_distance = distance; + } else { + min_distance = std::min(distance, min_distance); + } + } + max_distance = std::max(min_distance, max_distance); + } + + double second_max_distance = 0; + + for (auto & vertex : mesh2->meshes.vertices) { + // find closest vertex on mesh2 that corresponds to vertex + double min_distance = -1; + for (auto & vertex2 : mesh1->meshes.vertices) { + double distance = vertex.distance(vertex2); + if (min_distance == -1) { + min_distance = distance; + } else { + min_distance = std::min(distance, min_distance); + } + } + second_max_distance = std::max(min_distance, second_max_distance); + } + + FUNCINFO("HAUSDORFF DISTANCE: " << max_distance << " or " << second_max_distance); + + // print out surface areas + FUNCINFO("SURFACE AREA: First mesh = " << mesh1->meshes.surface_area() << ", second mesh = " << mesh2->meshes.surface_area()); + FUNCINFO("SURFACE AREA difference: " << mesh1->meshes.surface_area() - mesh2->meshes.surface_area()); + + return true; +} diff --git a/src/Operations/CompareMeshes.h b/src/Operations/CompareMeshes.h new file mode 100644 index 00000000..b5df74f1 --- /dev/null +++ b/src/Operations/CompareMeshes.h @@ -0,0 +1,16 @@ +// CompareMeshes.h. + +#pragma once + +#include +#include + +#include "../Structs.h" + + +OperationDoc OpArgDocCompareMeshes(); + +bool CompareMeshes(Drover &DICOM_data, + const OperationArgPkg& /*OptArgs*/, + std::map& /*InvocationMetadata*/, + const std::string& /*FilenameLex*/); From 5e7aa2c8d840a9848f9b24929f0187d163d6469a Mon Sep 17 00:00:00 2001 From: misprit7 Date: Sun, 5 Mar 2023 14:24:22 -0800 Subject: [PATCH 02/20] Added script to compare meshes with reconstructed versions --- scripts/compare_all_meshes.sh | 103 ++++++++++++++++++++++++++++++++++ 1 file changed, 103 insertions(+) create mode 100755 scripts/compare_all_meshes.sh diff --git a/scripts/compare_all_meshes.sh b/scripts/compare_all_meshes.sh new file mode 100755 index 00000000..a7d0a2b3 --- /dev/null +++ b/scripts/compare_all_meshes.sh @@ -0,0 +1,103 @@ + +POSITIONAL_ARGS=() +OPERATIONS="" +VERBOSE=0 + +while [[ $# -gt 0 ]]; do + case $1 in + -s|--show) + OPERATIONS+="-o SDL_Viewer" + shift # past argument + ;; + -v|--verbose) + VERBOSE=1 + shift # past argument + ;; + -*|--*) + echo "Unknown option $1" + exit 1 + ;; + *) + POSITIONAL_ARGS+=("$1") # save positional arg + shift # past argument + ;; + esac +done + + +declare -a shapes=("Sphere(10)" "Sphere(5)" "Sphere(7)") +# declare -a shapes=("Sphere(10)" ") + +printf "| %-15s | %-13s | %-13s | %-15s | %-15s | %-10s |\n" Shape "Hausdorff 1" "Hausdorff 2" "Surface Area 1" "Aurface Area 2" "Area Diff" +for shape in "${shapes[@]}" +do + # echo "Computing differences for $shape" + OUTPUT="$(dicomautomaton_dispatcher \ + -v \ + -o GenerateMeshes \ + -p resolution="0.25, 0.25, 0.25" \ + -p MeshLabel="original" \ + -p Objects="$shape" \ + -o GenerateSyntheticImages \ + -p NumberOfImages='128' \ + -p NumberOfRows='64', \ + -p NumberOfColumns='64' \ + -p NumberOfChannels='1' \ + -p SliceThickness='1' \ + -p SpacingBetweenSlices='1' \ + -p VoxelWidth='1.0' \ + -p VoxelHeight='1.0' \ + -p ImageAnchor='-14.0, -14.0, -14.0' \ + -p ImagePosition='0.0, 0.0, 0.0' \ + -p ImageOrientationColumn='1.0, 0.0, 0.0' \ + -p ImageOrientationRow='0.0, 1.0, 0.0' \ + -p InstanceNumber='1' \ + -p AcquisitionNumber='1' \ + -p VoxelValue='0.0' \ + -p StipleValue='nan' \ + -o ConvertMeshesToContours \ + -p ROILabel="original contours" \ + -p MeshSelection="last" \ + -p ImageSelection="last" \ + -o ConvertContoursToMeshes \ + -p NormalizedROILabelRegex=".*" \ + -p ROILabelRegex=".*" \ + -p MeshLabel="reconstructed" \ + -p Method="direct" \ + -o CompareMeshes \ + -p MeshSelection1="#-0" \ + -p MeshSelection2="#-1" \ + $OPERATIONS 2>&1)" + if [ "$VERBOSE" == 1 ]; then + echo "$OUTPUT" + fi + HD_LINE=$(grep "HAUSDORFF DISTANCE" <<< "$OUTPUT") + SA_LINE=$(grep "SURFACE AREA: First mesh" <<< "$OUTPUT") + SA_DIFF_LINE=$(grep "SURFACE AREA difference" <<< "$OUTPUT") + + HD1=$(sed -r "s/.*([0-9]+\.[0-9]+) or.*/\1/g" <<< $HD_LINE) + HD2=$(sed -r "s/.*or ([0-9]+\.[0-9]+).*/\1/g" <<< $HD_LINE) + + SA1=$(sed -r "s/.*First mesh = ([0-9]+\.[0-9]+).*/\1/g" <<< $SA_LINE) + SA2=$(sed -r "s/.*second mesh = ([0-9]+\.[0-9]+).*/\1/g" <<< $SA_LINE) + + SA_DIFF=$(sed -r "s/.*difference: (-?[0-9]+\.[0-9]+).*/\1/g" <<< $SA_DIFF_LINE) + + # echo "$HD1" + # echo "$HD2" + # echo "$SA1" + # echo "$SA2" + # echo "$SA_DIFF" + printf "| %-15s | %-13.3f | %-13.3f | %-15.3f | %-15.3f | %-10.3f |\n" $shape $HD1 $HD2 $SA1 $SA2 $SA_DIFF + + # echo "Finished computing differences for $shape" +done + + + + + + + + + From 5bf3f07dd456aa64ddb861ee2b40d3af26d6b658 Mon Sep 17 00:00:00 2001 From: Renu Date: Sun, 5 Mar 2023 16:18:18 -0800 Subject: [PATCH 03/20] Added volume comparison --- src/Operations/CompareMeshes.cc | 38 +++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/src/Operations/CompareMeshes.cc b/src/Operations/CompareMeshes.cc index 5ea59abb..5fbd5705 100644 --- a/src/Operations/CompareMeshes.cc +++ b/src/Operations/CompareMeshes.cc @@ -109,6 +109,44 @@ bool CompareMeshes(Drover &DICOM_data, } second_max_distance = std::max(min_distance, second_max_distance); } + + mesh1->meshes.convert_to_triangles(); + mesh2->meshes.convert_to_triangles(); + + + double volume1 = 0; + double volume2 = 0; + + decltype(mesh1->meshes.faces)* faces = &(mesh1->meshes.faces); + for (auto & fv : *faces){ + + const auto P_A = mesh1->meshes.vertices.at( fv[0] ); + const auto P_B = mesh1->meshes.vertices.at( fv[1] ); + const auto P_C = mesh1->meshes.vertices.at( fv[2] ); + + volume1 += (-P_C.x*P_B.y*P_A.z + P_B.x*P_C.y*P_A.z + + P_C.x*P_A.y*P_B.z - P_A.x*P_C.y*P_B.z - + P_B.x*P_A.y*P_C.z + P_A.x*P_B.y*P_C.z)/(6.0); + + } + + FUNCINFO("vol1 = "<meshes.faces); + for (auto & fv : *faces){ + + const auto P_A = mesh2->meshes.vertices.at( fv[0] ); + const auto P_B = mesh2->meshes.vertices.at( fv[1] ); + const auto P_C = mesh2->meshes.vertices.at( fv[2] ); + + volume2 += (-P_C.x*P_B.y*P_A.z + P_B.x*P_C.y*P_A.z + + P_C.x*P_A.y*P_B.z - P_A.x*P_C.y*P_B.z - + P_B.x*P_A.y*P_C.z + P_A.x*P_B.y*P_C.z)/(6.0); + } + + FUNCINFO("vol2 = "< Date: Sun, 12 Mar 2023 14:27:58 -0700 Subject: [PATCH 04/20] Add more shapes to bash script (#5) * add conic box issue * add correspondence issue * add more shapes --- scripts/compare_all_meshes.sh | 45 ++++++++++++++++++++++++++++------- 1 file changed, 37 insertions(+), 8 deletions(-) diff --git a/scripts/compare_all_meshes.sh b/scripts/compare_all_meshes.sh index a7d0a2b3..4baf6f20 100755 --- a/scripts/compare_all_meshes.sh +++ b/scripts/compare_all_meshes.sh @@ -24,18 +24,47 @@ while [[ $# -gt 0 ]]; do esac done - -declare -a shapes=("Sphere(10)" "Sphere(5)" "Sphere(7)") -# declare -a shapes=("Sphere(10)" ") - -printf "| %-15s | %-13s | %-13s | %-15s | %-15s | %-10s |\n" Shape "Hausdorff 1" "Hausdorff 2" "Surface Area 1" "Aurface Area 2" "Area Diff" -for shape in "${shapes[@]}" +tri_force_shape=" +subtract(){ + join(){ + plane(0,0,0, 1,1,0, 10); + plane(0,0,0, 0,1,1, 10); + plane(0,0,0, 1,0,1, 10); + + plane(10,10,10, -1,-1,-0, 10); + plane(10,10,10, -0,-1,-1, 10); + plane(10,10,10, -1,-0,-1, 10); + }; + + join(){ + plane(0,0,0, 1,2,0, 10); + plane(0,0,0, 0,1,2, 10); + plane(0,0,0, 2,0,1, 10); + + plane(10,10,10, -2,-1,-0, 10); + plane(10,10,10, -0,-2,-1, 10); + plane(10,10,10, -1,-0,-2, 10); + }; +}; +" + +default_res="0.25,0.25,0.25" + +declare -a shape_labels=("Sphere(10)" "aa_box(1.0,2.0,4.0)" "tri-force" "chamfered(box-sphere)") +declare -a shapes=("Sphere(10)" "aa_box(1.0,2.0,4.0)" "$tri_force_shape" "chamfer_subtract(2.0){aa_box(8.0,20.0,20.0);sphere(12.0);}") +declare -a resolutions=("$default_res" "0.25,0.25,0.75" "0.25,0.25,0.75" "$default_res") + +printf "| %-20s | %-13s | %-13s | %-15s | %-15s | %-10s |\n" Shape "Hausdorff 1" "Hausdorff 2" "Surface Area 1" "Surface Area 2" "Area Diff" +for i in ${!shapes[@]} do + shape=${shapes[$i]} + res=${resolutions[$i]} + shape_label=${shape_labels[$i]} # echo "Computing differences for $shape" OUTPUT="$(dicomautomaton_dispatcher \ -v \ -o GenerateMeshes \ - -p resolution="0.25, 0.25, 0.25" \ + -p resolution=$res \ -p MeshLabel="original" \ -p Objects="$shape" \ -o GenerateSyntheticImages \ @@ -88,7 +117,7 @@ do # echo "$SA1" # echo "$SA2" # echo "$SA_DIFF" - printf "| %-15s | %-13.3f | %-13.3f | %-15.3f | %-15.3f | %-10.3f |\n" $shape $HD1 $HD2 $SA1 $SA2 $SA_DIFF + printf "| %-20s | %-13.3f | %-13.3f | %-15.3f | %-15.3f | %-10.3f |\n" $shape_label $HD1 $HD2 $SA1 $SA2 $SA_DIFF # echo "Finished computing differences for $shape" done From 0470ee6b74fc23cc7f54cc9d6a8ffb0785f2af39 Mon Sep 17 00:00:00 2001 From: tiffanymatthe <49164666+tiffanymatthe@users.noreply.github.com> Date: Sun, 12 Mar 2023 14:41:25 -0700 Subject: [PATCH 05/20] Add box mesh with conic issue in compare mesh script (#4) --- .github/workflows/main.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index d62c8779..c6058cb7 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -27,6 +27,7 @@ jobs: rsync -avP /out/ ./Windows/ strip ./Windows/usr/bin/*exe strip ./Windows/usr/lib/*a + git config --global --add safe.directory '*' || true - name: test run: | export DEBIAN_FRONTEND='noninteractive' @@ -47,7 +48,7 @@ jobs: ./integration_tests/Run.sh -v - name: create_zip run: | - export short_sha="$(git rev-parse --short HEAD)" + export short_sha="$(git rev-parse --short HEAD || printf 'unknown')" export DEBIAN_FRONTEND="noninteractive" apt-get update --yes apt-get install --yes --no-install-recommends coreutils binutils findutils rsync openssh-client patchelf zip @@ -86,6 +87,7 @@ jobs: run: | ( ( while true ; do sleep 225 ; printf '\n\n' ; df -h ; printf '\n\n' ; free || true ; printf '\n\n' ; done ) &) ./docker/builders/ci/implementation_ci_debian_oldstable.sh + git config --global --add safe.directory '*' || true - name: test run: | ./integration_tests/Run.sh -v From c30cb50a11bdd0393684a17941a29668487a8f78 Mon Sep 17 00:00:00 2001 From: Renu Rajamagesh <60906270+Renu-R@users.noreply.github.com> Date: Thu, 23 Mar 2023 13:47:50 -0700 Subject: [PATCH 06/20] added vol comparison to script (#6) * added vol comparison to script * added calculations of centroids and centroid shift * fixed volmue calc err * addressed comments on PR * fixed typo * another typo fixed, spacing fixed * corrected variable name * added some spacing * added spacing * typo fix --- scripts/compare_all_meshes.sh | 22 ++++++++----- src/Operations/CompareMeshes.cc | 56 +++++++++++++++++++++++---------- 2 files changed, 55 insertions(+), 23 deletions(-) diff --git a/scripts/compare_all_meshes.sh b/scripts/compare_all_meshes.sh index 4baf6f20..91990064 100755 --- a/scripts/compare_all_meshes.sh +++ b/scripts/compare_all_meshes.sh @@ -1,4 +1,4 @@ - +#!/bin/bash POSITIONAL_ARGS=() OPERATIONS="" VERBOSE=0 @@ -54,7 +54,8 @@ declare -a shape_labels=("Sphere(10)" "aa_box(1.0,2.0,4.0)" "tri-force" "chamfer declare -a shapes=("Sphere(10)" "aa_box(1.0,2.0,4.0)" "$tri_force_shape" "chamfer_subtract(2.0){aa_box(8.0,20.0,20.0);sphere(12.0);}") declare -a resolutions=("$default_res" "0.25,0.25,0.75" "0.25,0.25,0.75" "$default_res") -printf "| %-20s | %-13s | %-13s | %-15s | %-15s | %-10s |\n" Shape "Hausdorff 1" "Hausdorff 2" "Surface Area 1" "Surface Area 2" "Area Diff" +printf "| %-20s | %-13s | %-13s | %-15s | %-15s | %-10s | %-15s | %-15s | %-10s| \n\n" Shape "Hausdorff 1" "Hausdorff 2" "Surface Area 1" "Surface Area 2" "Area Diff" "Volume 1" "Volume 2" "Volume Diff" + for i in ${!shapes[@]} do shape=${shapes[$i]} @@ -102,7 +103,10 @@ do fi HD_LINE=$(grep "HAUSDORFF DISTANCE" <<< "$OUTPUT") SA_LINE=$(grep "SURFACE AREA: First mesh" <<< "$OUTPUT") - SA_DIFF_LINE=$(grep "SURFACE AREA difference" <<< "$OUTPUT") + SA_DIFF_LINE=$(grep "SURFACE AREA (%) difference" <<< "$OUTPUT") + VOL_LINE=$(grep "VOLUME: First mesh" <<< "$OUTPUT") + VOL_DIFF_LINE=$(grep "VOLUME (%) difference" <<< "$OUTPUT") + HD1=$(sed -r "s/.*([0-9]+\.[0-9]+) or.*/\1/g" <<< $HD_LINE) HD2=$(sed -r "s/.*or ([0-9]+\.[0-9]+).*/\1/g" <<< $HD_LINE) @@ -112,12 +116,16 @@ do SA_DIFF=$(sed -r "s/.*difference: (-?[0-9]+\.[0-9]+).*/\1/g" <<< $SA_DIFF_LINE) + V1=$(sed -r "s/.*First mesh = ([0-9]+\.[0-9]+).*/\1/g" <<< $VOL_LINE) + V2=$(sed -r "s/.*second mesh = ([0-9]+\.[0-9]+).*/\1/g" <<< $VOL_LINE) + + V_DIFF=$(sed -r "s/.*difference: (-?[0-9]+\.[0-9]+).*/\1/g" <<< $VOL_DIFF_LINE) + # echo "$HD1" # echo "$HD2" - # echo "$SA1" - # echo "$SA2" - # echo "$SA_DIFF" - printf "| %-20s | %-13.3f | %-13.3f | %-15.3f | %-15.3f | %-10.3f |\n" $shape_label $HD1 $HD2 $SA1 $SA2 $SA_DIFF + + printf "| %-20s | %-13.3f | %-13.3f | %-15.3f | %-15.3f | %-10.3f | %-15.3f | %-15.3f | %-10.3f |\n\n" $shape_label $HD1 $HD2 $SA1 $SA2 $SA_DIFF $V1 $V2 $V_DIFF + # echo "Finished computing differences for $shape" done diff --git a/src/Operations/CompareMeshes.cc b/src/Operations/CompareMeshes.cc index 5fbd5705..83d3ef83 100644 --- a/src/Operations/CompareMeshes.cc +++ b/src/Operations/CompareMeshes.cc @@ -80,10 +80,10 @@ bool CompareMeshes(Drover &DICOM_data, FUNCINFO("Iterating through " << size(mesh1->meshes.vertices) << " and " << size(mesh2->meshes.vertices) << " vertices."); // O(n^2) time, might be way to speed this up if this is ever a bottleneck - for (auto & vertex : mesh1->meshes.vertices) { + for (auto& vertex : mesh1->meshes.vertices) { // find closest vertex on mesh2 that corresponds to vertex double min_distance = -1; - for (auto & vertex2 : mesh2->meshes.vertices) { + for (auto& vertex2 : mesh2->meshes.vertices) { double distance = vertex.distance(vertex2); if (min_distance == -1) { min_distance = distance; @@ -96,10 +96,10 @@ bool CompareMeshes(Drover &DICOM_data, double second_max_distance = 0; - for (auto & vertex : mesh2->meshes.vertices) { + for (auto& vertex : mesh2->meshes.vertices) { // find closest vertex on mesh2 that corresponds to vertex double min_distance = -1; - for (auto & vertex2 : mesh1->meshes.vertices) { + for (auto& vertex2 : mesh1->meshes.vertices) { double distance = vertex.distance(vertex2); if (min_distance == -1) { min_distance = distance; @@ -109,7 +109,31 @@ bool CompareMeshes(Drover &DICOM_data, } second_max_distance = std::max(min_distance, second_max_distance); } + + // Calculate centroids by finding the average of all the vertices in the surface mesh + //Assumes that contours vertices are distributed evenly on the surface mesh which is not strictly true. + double sumx = 0, sumy = 0, sumz = 0; + for (auto& vertex : mesh1->meshes.vertices) { + sumx += vertex.x; + sumy += vertex.y; + sumz += vertex.z; + } + vec3 centroid1 = vec3(sumx/size(mesh1->meshes.vertices) , sumy/size(mesh1->meshes.vertices) , sumz/size(mesh1->meshes.vertices)); + + sumx = 0, sumy= 0, sumz = 0; + for (auto& vertex : mesh2->meshes.vertices) { + sumx += vertex.x; + sumy += vertex.y; + sumz += vertex.z; + } + vec3 centroid2 = vec3(sumx/size(mesh1->meshes.vertices) , sumy/size(mesh1->meshes.vertices) , sumz/size(mesh1->meshes.vertices)); + + const double centroid_shift = sqrt(pow((centroid2.x-centroid1.x),2) + pow((centroid2.y-centroid1.y),2) + pow((centroid2.y-centroid1.y),2)); + //Converting to a triangular mesh to ensure that each face is made up of 3 vertices for volume calculation + //Total volume is calculated by summing the signed volme of the triganular prism made by each face and the origin as the apex. + //This method returns a finite volume even if the mesh is open, so watertightness needs to be checked separately. + mesh1->meshes.convert_to_triangles(); mesh2->meshes.convert_to_triangles(); @@ -118,7 +142,7 @@ bool CompareMeshes(Drover &DICOM_data, double volume2 = 0; decltype(mesh1->meshes.faces)* faces = &(mesh1->meshes.faces); - for (auto & fv : *faces){ + for (auto& fv : *faces){ const auto P_A = mesh1->meshes.vertices.at( fv[0] ); const auto P_B = mesh1->meshes.vertices.at( fv[1] ); @@ -128,31 +152,31 @@ bool CompareMeshes(Drover &DICOM_data, P_C.x*P_A.y*P_B.z - P_A.x*P_C.y*P_B.z - P_B.x*P_A.y*P_C.z + P_A.x*P_B.y*P_C.z)/(6.0); - } - - FUNCINFO("vol1 = "<meshes.faces); - for (auto & fv : *faces){ + for (auto& fv : *faces){ const auto P_A = mesh2->meshes.vertices.at( fv[0] ); const auto P_B = mesh2->meshes.vertices.at( fv[1] ); const auto P_C = mesh2->meshes.vertices.at( fv[2] ); - volume2 += (-P_C.x*P_B.y*P_A.z + P_B.x*P_C.y*P_A.z + + volume2 += (-P_C.x*P_B.y*P_A.z + P_B.x*P_C.y*P_A.z + P_C.x*P_A.y*P_B.z - P_A.x*P_C.y*P_B.z - P_B.x*P_A.y*P_C.z + P_A.x*P_B.y*P_C.z)/(6.0); - } + } - FUNCINFO("vol2 = "<meshes.surface_area() - mesh2->meshes.surface_area()); + FUNCINFO("SURFACE AREA (%) difference: " << (mesh1->meshes.surface_area() - mesh2->meshes.surface_area())*100/mesh1->meshes.surface_area()); + FUNCINFO("VOLUME: First mesh = " <get().Get_Signed_Area() ); const auto m2_area = std::abs( m2_cop_it->get().Get_Signed_Area() ); - YLOGWARN("Found overlapping upper-plane contours, trimmed smallest-area contour"); + YLOGWARN("Found intersecting upper-plane contours, trimmed smallest-area contour"); if(m1_area < m2_area){ m1_cop_it = m_cops.erase(m1_cop_it); m2_cop_it = std::next(m1_cop_it); @@ -322,14 +354,14 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, } for(auto l1_cop_it = std::begin(l_cops); l1_cop_it != std::end(l_cops); ){ for(auto l2_cop_it = std::next(l1_cop_it); l2_cop_it != std::end(l_cops); ){ - if(projected_contours_overlap(*l_cp_it, *l1_cop_it, - *l_cp_it, *l2_cop_it)){ + if(projected_contours_intersect(*l_cp_it, *l1_cop_it, + *l_cp_it, *l2_cop_it)){ // Cull the smaller contour. const auto l1_area = std::abs( l1_cop_it->get().Get_Signed_Area() ); const auto l2_area = std::abs( l2_cop_it->get().Get_Signed_Area() ); - YLOGWARN("Found overlapping lower-plane contours, trimmed smallest-area contour"); + YLOGWARN("Found intersecting lower-plane contours, trimmed smallest-area contour"); if(l1_area < l2_area){ l1_cop_it = l_cops.erase(l1_cop_it); l2_cop_it = std::next(l1_cop_it); From e7733c6407eece09b2a165b6b241c8f3de1c17b3 Mon Sep 17 00:00:00 2001 From: DominicKlukas Date: Thu, 23 Mar 2023 19:10:48 -0700 Subject: [PATCH 09/20] Created a new operation which lets us create custom contours, for testing --- src/Operation_Dispatcher.cc | 2 + src/Operations/CreateCustomContour.cc | 64 +++++++++++++++++++++++++++ src/Operations/CreateCustomContour.h | 16 +++++++ src/Simple_Meshing.cc | 62 ++++++++++++++++++++++++++ 4 files changed, 144 insertions(+) create mode 100644 src/Operations/CreateCustomContour.cc create mode 100644 src/Operations/CreateCustomContour.h diff --git a/src/Operation_Dispatcher.cc b/src/Operation_Dispatcher.cc index 7790ad4b..7bfd9207 100644 --- a/src/Operation_Dispatcher.cc +++ b/src/Operation_Dispatcher.cc @@ -67,6 +67,7 @@ #include "Operations/CopyTables.h" #include "Operations/CopyPoints.h" #include "Operations/CountVoxels.h" +#include "Operations/CreateCustomContour.h" #include "Operations/CropImageDoseToROIs.h" #include "Operations/CropImages.h" #include "Operations/CropROIDose.h" @@ -310,6 +311,7 @@ std::map Known_Operations(){ out["CopyTables"] = std::make_pair(OpArgDocCopyTables, CopyTables); out["CopyPoints"] = std::make_pair(OpArgDocCopyPoints, CopyPoints); out["CountVoxels"] = std::make_pair(OpArgDocCountVoxels, CountVoxels); + out["CreateCustomContour"] = std::make_pair(OpArgDocCreateCustomContour, CreateCustomContour); out["CropImageDoseToROIs"] = std::make_pair(OpArgDocCropImageDoseToROIs, CropImageDoseToROIs); out["CropImages"] = std::make_pair(OpArgDocCropImages, CropImages); out["CropROIDose"] = std::make_pair(OpArgDocCropROIDose, CropROIDose); diff --git a/src/Operations/CreateCustomContour.cc b/src/Operations/CreateCustomContour.cc new file mode 100644 index 00000000..d8db29d5 --- /dev/null +++ b/src/Operations/CreateCustomContour.cc @@ -0,0 +1,64 @@ + +//CopyContours.cc - A part of DICOMautomaton 2021. Written by hal clark. + +#include +#include +#include + +#include "../Structs.h" +#include "YgorLog.h" + +#include "CreateCustomContour.h" + +OperationDoc OpArgDocCreateCustomContour(){ + OperationDoc out; + out.name = "CreateCustomContour"; + + out.desc = + "This operation creates a new contour from a list of x, y, and z coordinates."; + + out.args.emplace_back(); + out.args.back().name = "ROILabel"; + out.args.back().desc = "A label to attach to the copied ROI contours."; + out.args.back().default_val = "unspecified"; + out.args.back().expected = true; + + out.args.emplace_back(); + out.args.back().name = "XValues"; + out.args.back().desc = "List of X coordinates of the points in the contour"; + out.args.back().default_val = "0"; + out.args.back().expected = true; + + out.args.emplace_back(); + out.args.back().name = "YValues"; + out.args.back().desc = "List of Y coordinates of the points in the contour"; + out.args.back().default_val = "0"; + out.args.back().expected = true; + + out.args.emplace_back(); + out.args.back().name = "ZValues"; + out.args.back().desc = "List of Z coordinates of the points in the contour"; + out.args.back().default_val = "0"; + out.args.back().expected = true; + return out; +} + +bool CreateCustomContour(Drover &DICOM_data, + const OperationArgPkg& OptArgs, + std::map& /*InvocationMetadata*/, + const std::string& FilenameLex){ + + //---------------------------------------------- User Parameters -------------------------------------------------- + const auto ROILabel = OptArgs.getValueStr("ROILabel").value(); + const auto XValues = OptArgs.getValueStr("XValues").value(); + const auto YValues = OptArgs.getValueStr("YValues").value(); + const auto ZValues = OptArgs.getValueStr("ZValues").value(); + //----------------------------------------------------------------------------------------------------------------- + + //Stuff references to all contours into a list. Remember that you can still address specific contours through + // the original holding containers (which are not modified here). + YLOGINFO("The X values are " << XValues); + YLOGINFO("The Y values are " << YValues); + YLOGINFO("The Z values are " << YValues); + return true; +} diff --git a/src/Operations/CreateCustomContour.h b/src/Operations/CreateCustomContour.h new file mode 100644 index 00000000..dcbf0aba --- /dev/null +++ b/src/Operations/CreateCustomContour.h @@ -0,0 +1,16 @@ +// CreateCustomContour.h. + +#pragma once + +#include +#include + +#include "../Structs.h" + + +OperationDoc OpArgDocCreateCustomContour(); + +bool CreateCustomContour(Drover &DICOM_data, + const OperationArgPkg& /*OptArgs*/, + std::map& /*InvocationMetadata*/, + const std::string& FilenameLex); diff --git a/src/Simple_Meshing.cc b/src/Simple_Meshing.cc index defc7c67..1f5e2694 100644 --- a/src/Simple_Meshing.cc +++ b/src/Simple_Meshing.cc @@ -745,3 +745,65 @@ Minimally_Amalgamate_Contours( return amal; } +contour_of_points +Minimally_Amalgamate_Contours_N_to_N( + const vec3 &ortho_unit, + const vec3 &pseudo_vert_offset, + std::list>> B ){ + if(B.empty()) { + throw std::invalid_argument("No contours supplied in B. Cannot continue."); + } + + const auto has_consistent_orientation = [&]( std::reference_wrapper> cop_refw ) -> bool { + // Determine if the contour has a consistent orientation as that provided by the user. + auto l_ortho_unit = ortho_unit; + try{ + l_ortho_unit = cop_refw.get().Estimate_Planar_Normal(); + }catch(const std::exception &){} + return (0.0 < ortho_unit.Dot(l_ortho_unit)); + }; + + // Confirm all contour orientations are consistent. + std::list< contour_of_points > reversed_cops_storage; + for(auto crw_it = std::begin(B); crw_it != std::end(B); ){ + if(!has_consistent_orientation(*crw_it)){ + YLOGWARN("Found contour with inconsistent orientation. Making a reversed copy. This may discard information"); + reversed_cops_storage.emplace_back( crw_it->get() ); + reversed_cops_storage.back().points.reverse(); + (*crw_it) = std::ref(reversed_cops_storage.back()); + } + ++crw_it; + } + for(const auto &cop_refw : B){ + if(cop_refw.get().closed == false){ + throw std::invalid_argument("Found open contour. Refusing to continue."); + } + } + contour_of_points amal = B.front().get(); + B.pop_front(); + if(amal.points.size() < 3) { + throw std::invalid_argument("Seed contour in B contains insufficient vertices. Cannot continue."); + } + + std::list pseudo_verts; + const auto is_a_pseudo_vert = [&]( decltype(std::begin(amal.points)) it ) -> bool { + for(const auto &pv : pseudo_verts){ + if( it == pv ) return true; + } + return false; + }; + + while(true) { + double shortest_criteria = std::numeric_limits::infinity(); + auto closest_cop = std::begin(B); + decltype(std::begin(amal.points)) closest_a_v_it; + decltype(std::begin(B.front().get().points)) closest_b_it; + + for(auto b_it = std::begin(B); b_it != std::end(B); ++b_it){ + auto a_v_it = std::begin(amal.points); + if(is_a_pseudo_vert(a_v_it)) { + + } + } + } +} \ No newline at end of file From 4e1d55c956ef05d4298f02eecc97fb4408b1bd80 Mon Sep 17 00:00:00 2001 From: DominicKlukas Date: Thu, 23 Mar 2023 20:35:01 -0700 Subject: [PATCH 10/20] Implemented first version of CreateCustomContour --- src/Operations/CMakeLists.txt | 1 + src/Operations/CreateCustomContour.cc | 54 ++++++++++++++++++++++++++- src/Operations/CreateCustomContour.h | 4 +- src/Simple_Meshing.cc | 4 +- 4 files changed, 59 insertions(+), 4 deletions(-) diff --git a/src/Operations/CMakeLists.txt b/src/Operations/CMakeLists.txt index 8f86992b..d9445905 100644 --- a/src/Operations/CMakeLists.txt +++ b/src/Operations/CMakeLists.txt @@ -46,6 +46,7 @@ add_library(Operations_objs OBJECT CopyTables.cc CopyPoints.cc CountVoxels.cc + CreateCustomContour.cc CropImageDoseToROIs.cc CropImages.cc CropROIDose.cc diff --git a/src/Operations/CreateCustomContour.cc b/src/Operations/CreateCustomContour.cc index d8db29d5..acf5f1d1 100644 --- a/src/Operations/CreateCustomContour.cc +++ b/src/Operations/CreateCustomContour.cc @@ -59,6 +59,58 @@ bool CreateCustomContour(Drover &DICOM_data, // the original holding containers (which are not modified here). YLOGINFO("The X values are " << XValues); YLOGINFO("The Y values are " << YValues); - YLOGINFO("The Z values are " << YValues); + YLOGINFO("The Z values are " << ZValues); + + contour_of_points cop; + contour_collection cc; + + cop.closed = true; + + std::vector x_values = ValueStringToDoubleList(XValues); + std::vector y_values = ValueStringToDoubleList(YValues); + std::vector z_values = ValueStringToDoubleList(ZValues); + + if(x_values.size()!=y_values.size()||y_values.size()!=z_values.size()) { + YLOGWARN("Each list must have the same number of numbers"); + return true; + } + + if(x_values.size() < 3) { + YLOGWARN("We require more than 3 points to be in a contour"); + return true; + } + + for(int i = 0; i < x_values.size(); i++) { + const vec3 p(static_cast(x_values[i]), + static_cast(y_values[i]), + static_cast(z_values[i])); + cop.points.emplace_back(p); + } + cop.metadata["ROILabel"] = ROILabel; + + cc.contours.emplace_back(cop); + + DICOM_data.Ensure_Contour_Data_Allocated(); + DICOM_data.contour_data->ccs.emplace_back(cc); return true; } + +std::vector ValueStringToDoubleList(std::string input) { + std::vector values; + size_t pos = 0; + std::string delimiter = " "; + std::string token; + while ((pos = input.find(delimiter)) != std::string::npos) { + token = input.substr(0, pos); + try { + values.emplace_back(stod(token)); + std::string error_message = "Invalid Input"; + } catch(std::string error_message) { + YLOGWARN(error_message); + return std::vector(); + } + input.erase(0, pos + delimiter.length()); + } + return values; +} + diff --git a/src/Operations/CreateCustomContour.h b/src/Operations/CreateCustomContour.h index dcbf0aba..74b06a02 100644 --- a/src/Operations/CreateCustomContour.h +++ b/src/Operations/CreateCustomContour.h @@ -2,7 +2,7 @@ #pragma once -#include +#include #include #include "../Structs.h" @@ -14,3 +14,5 @@ bool CreateCustomContour(Drover &DICOM_data, const OperationArgPkg& /*OptArgs*/, std::map& /*InvocationMetadata*/, const std::string& FilenameLex); + +std::vector ValueStringToDoubleList(std::string input); diff --git a/src/Simple_Meshing.cc b/src/Simple_Meshing.cc index 1f5e2694..e5732168 100644 --- a/src/Simple_Meshing.cc +++ b/src/Simple_Meshing.cc @@ -745,7 +745,7 @@ Minimally_Amalgamate_Contours( return amal; } -contour_of_points +/*contour_of_points Minimally_Amalgamate_Contours_N_to_N( const vec3 &ortho_unit, const vec3 &pseudo_vert_offset, @@ -806,4 +806,4 @@ Minimally_Amalgamate_Contours_N_to_N( } } } -} \ No newline at end of file +}*/ \ No newline at end of file From f502f44d9cfa8d263678af2c2f14845c1fe08868 Mon Sep 17 00:00:00 2001 From: tiffanymatthe <49164666+tiffanymatthe@users.noreply.github.com> Date: Wed, 5 Apr 2023 16:23:42 -0700 Subject: [PATCH 11/20] Fix capping issue (#9) * CI: GitHub: work around dir ownership. For some reason the ownership of the root repo dire is not owned by the current user. (Maybe to make it de-facto read-only? Or maybe an artifact of containerization? Something else?) This commit silences git warnings about ownership that would otherwise result in irrelevant fatal errors. * CI: GitHub: tweak dir ownership workaround. Not sure which build stages / enviornments actually need this workaround. Either way, defer until git package has been explicitly installed. * initial changes to add higher plane * switch from pairings to pairs * add set removal * initial try * final before new simple try * almost there * done * clean up * simplified a lot * final final clean * add boolean --------- Co-authored-by: Haley Clark Co-authored-by: Hal Clark --- src/Operations/ConvertContoursToMeshes.cc | 49 ++++++++++++++--------- 1 file changed, 31 insertions(+), 18 deletions(-) diff --git a/src/Operations/ConvertContoursToMeshes.cc b/src/Operations/ConvertContoursToMeshes.cc index 96ca1dfa..c4353ea2 100644 --- a/src/Operations/ConvertContoursToMeshes.cc +++ b/src/Operations/ConvertContoursToMeshes.cc @@ -311,7 +311,9 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, // Identify whether there are adjacent planes within the contour spacing on either side. auto l_cp_it = std::prev(m_cp_it); - //auto h_cp_it = std::next(m_cp_it); + auto h_cp_it = std::next(m_cp_it); + + bool cap_roof_of_m_cops = false; if(l_cp_it == std::cend(ucps)){ continue; @@ -319,10 +321,13 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, const auto l_cp_dist = std::abs(m_cp_it->Get_Signed_Distance_To_Point(l_cp_it->R_0)); if((1.5 * contour_sep) < l_cp_dist) l_cp_it = std::cend(ucps); } - //if(h_cp_it != std::cend(ucps)){ - // const auto h_cp_dist = std::abs(m_cp_it->Get_Signed_Distance_To_Point(h_cp_it->R_0)); - // if((1.5 * contour_sep) < h_cp_dist) h_cp_it = std::cend(ucps); - //} + if(h_cp_it != std::cend(ucps)){ + const auto h_cp_dist = std::abs(m_cp_it->Get_Signed_Distance_To_Point(h_cp_it->R_0)); + if((1.5 * contour_sep) < h_cp_dist) { + h_cp_it = std::cend(ucps); + cap_roof_of_m_cops = true; + } + } auto l_cops = locate_contours_on_plane(*l_cp_it); if( (l_cops.size() == 0) && (m_cops.size() == 0) ){ @@ -543,22 +548,24 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, //YLOGINFO("Performing N-to-N meshing.."); auto ofst_upper = m_cp_it->N_0 * contour_sep * -0.49; auto ofst_lower = m_cp_it->N_0 * contour_sep * 0.49; + + // will modify contours in m_cops and l_cops, ok if only processing in one direction auto amal_upper = Minimally_Amalgamate_Contours(m_cp_it->N_0, ofst_upper, pcs.upper); auto amal_lower = Minimally_Amalgamate_Contours(m_cp_it->N_0, ofst_lower, pcs.lower); - /* - // Leaving this here for future debugging, for which it will no-doubt be needed... - { - const auto amal_cop_str = amal_upper.write_to_string(); - const auto fname = Get_Unique_Sequential_Filename("/tmp/amal_upper_", 6, ".txt"); - OverwriteStringToFile(amal_cop_str, fname); - } - { - const auto amal_cop_str = amal_lower.write_to_string(); - const auto fname = Get_Unique_Sequential_Filename("/tmp/amal_lower_", 6, ".txt"); - OverwriteStringToFile(amal_cop_str, fname); - } - */ + /* + // Leaving this here for future debugging, for which it will no-doubt be needed... + { + const auto amal_cop_str = amal_upper.write_to_string(); + const auto fname = Get_Unique_Sequential_Filename("/tmp/amal_upper_", 6, ".txt"); + OverwriteStringToFile(amal_cop_str, fname); + } + { + const auto amal_cop_str = amal_lower.write_to_string(); + const auto fname = Get_Unique_Sequential_Filename("/tmp/amal_lower_", 6, ".txt"); + OverwriteStringToFile(amal_cop_str, fname); + } + */ auto new_faces = Estimate_Contour_Correspondence(std::ref(amal_upper), std::ref(amal_lower)); const auto old_face_count = amesh.vertices.size(); @@ -570,6 +577,12 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, const auto f_C = static_cast(fs[2] + old_face_count); amesh.faces.emplace_back( std::vector{{f_A, f_B, f_C}} ); } + } + } + + if (cap_roof_of_m_cops) { + for (auto &cop : m_cops) { + close_hole_in_roof(cop); } } } From fbea4dfd903c8bfac9424fc54bb311a499a535bd Mon Sep 17 00:00:00 2001 From: DominicKlukas <69833981+DominicKlukas@users.noreply.github.com> Date: Sat, 8 Apr 2023 17:00:43 -0700 Subject: [PATCH 12/20] Upgraded CreateCustomContour, Fixed an N-to-N bug which results in open meshes in pathalogical cases (#12) * Upgraded CreateCustomContour * Fixed N-to-N bug which results open meshes in some pathalogical cases * Changed Namespace * Updated Docs --- src/Operations/CreateCustomContour.cc | 91 ++++++++++++++++++--------- src/Operations/CreateCustomContour.h | 2 +- src/Simple_Meshing.cc | 13 ++-- 3 files changed, 70 insertions(+), 36 deletions(-) diff --git a/src/Operations/CreateCustomContour.cc b/src/Operations/CreateCustomContour.cc index acf5f1d1..7137e680 100644 --- a/src/Operations/CreateCustomContour.cc +++ b/src/Operations/CreateCustomContour.cc @@ -3,7 +3,6 @@ #include #include -#include #include "../Structs.h" #include "YgorLog.h" @@ -25,19 +24,19 @@ OperationDoc OpArgDocCreateCustomContour(){ out.args.emplace_back(); out.args.back().name = "XValues"; - out.args.back().desc = "List of X coordinates of the points in the contour"; + out.args.back().desc = "List of contours separated by | character, where each countour is a list of doubles separated by spaces. Contains x-values"; out.args.back().default_val = "0"; out.args.back().expected = true; out.args.emplace_back(); out.args.back().name = "YValues"; - out.args.back().desc = "List of Y coordinates of the points in the contour"; + out.args.back().desc = "List of contours separated by | character, where each countour is a list of doubles separated by spaces. Contains y-values"; out.args.back().default_val = "0"; out.args.back().expected = true; out.args.emplace_back(); out.args.back().name = "ZValues"; - out.args.back().desc = "List of Z coordinates of the points in the contour"; + out.args.back().desc = "List of contours separated by | character, where each countour is a list of doubles separated by spaces. Contains z-values"; out.args.back().default_val = "0"; out.args.back().expected = true; return out; @@ -61,56 +60,88 @@ bool CreateCustomContour(Drover &DICOM_data, YLOGINFO("The Y values are " << YValues); YLOGINFO("The Z values are " << ZValues); - contour_of_points cop; contour_collection cc; - cop.closed = true; + std::vector x_contours = SplitString(XValues, "|"); + std::vector y_contours = SplitString(YValues, "|"); + std::vector z_contours = SplitString(ZValues, "|"); + - std::vector x_values = ValueStringToDoubleList(XValues); - std::vector y_values = ValueStringToDoubleList(YValues); - std::vector z_values = ValueStringToDoubleList(ZValues); - - if(x_values.size()!=y_values.size()||y_values.size()!=z_values.size()) { - YLOGWARN("Each list must have the same number of numbers"); + if(x_contours.size()!=y_contours.size()||y_contours.size()!=z_contours.size()) { + YLOGWARN("Ensure same number of contours for each dimension"); return true; } - if(x_values.size() < 3) { - YLOGWARN("We require more than 3 points to be in a contour"); + if(x_contours.size() < 1) { + YLOGWARN("We require that there is at least one contour."); return true; } - for(int i = 0; i < x_values.size(); i++) { - const vec3 p(static_cast(x_values[i]), - static_cast(y_values[i]), - static_cast(z_values[i])); - cop.points.emplace_back(p); + DICOM_data.Ensure_Contour_Data_Allocated(); + + std::vector> x_cops; + std::vector> y_cops; + std::vector> z_cops; + + for(int i = 0; i < x_contours.size(); i++) { + std::vector x_cop_string = SplitString(x_contours[i], " "); + std::vector y_cop_string = SplitString(y_contours[i], " "); + std::vector z_cop_string = SplitString(z_contours[i], " "); + if(x_cop_string.size()!=y_cop_string.size()||y_cop_string.size()!=z_cop_string.size()) { + YLOGWARN("Ensure each contour has the same number of points for each dimension"); + } + if(x_cop_string.size() < 3) { + YLOGWARN("Each contour must have at least 3 points"); + } + std::vector x_cop; + std::vector y_cop; + std::vector z_cop; + for(int j = 0; j < x_cop_string.size(); j++) { + try { + x_cop.emplace_back(stod(x_cop_string[j])); + y_cop.emplace_back(stod(y_cop_string[j])); + z_cop.emplace_back(stod(z_cop_string[j])); + std::string error_message = "Invalid Input"; + } catch(std::string error_message) { + return true; + } + } + x_cops.emplace_back(x_cop); + y_cops.emplace_back(y_cop); + z_cops.emplace_back(z_cop); } - cop.metadata["ROILabel"] = ROILabel; - cc.contours.emplace_back(cop); + for(int i = 0; i < x_cops.size(); i++) { + contour_of_points cop; + for(int j = 0; j < x_cops[i].size(); j++) { + const vec3 p(static_cast(x_cops[i][j]), + static_cast(y_cops[i][j]), + static_cast(z_cops[i][j])); + cop.points.emplace_back(p); + } + cop.closed = true; + cop.metadata["ROILabel"] = ROILabel; + cc.contours.emplace_back(cop); + } - DICOM_data.Ensure_Contour_Data_Allocated(); DICOM_data.contour_data->ccs.emplace_back(cc); return true; } -std::vector ValueStringToDoubleList(std::string input) { - std::vector values; +std::vector SplitString(std::string input, std::string delimiter) { + std::vector values; size_t pos = 0; - std::string delimiter = " "; std::string token; while ((pos = input.find(delimiter)) != std::string::npos) { token = input.substr(0, pos); try { - values.emplace_back(stod(token)); + values.emplace_back(token); std::string error_message = "Invalid Input"; } catch(std::string error_message) { - YLOGWARN(error_message); - return std::vector(); + return std::vector(); } input.erase(0, pos + delimiter.length()); } + values.emplace_back(input); return values; -} - +} \ No newline at end of file diff --git a/src/Operations/CreateCustomContour.h b/src/Operations/CreateCustomContour.h index 74b06a02..7425d2f8 100644 --- a/src/Operations/CreateCustomContour.h +++ b/src/Operations/CreateCustomContour.h @@ -15,4 +15,4 @@ bool CreateCustomContour(Drover &DICOM_data, std::map& /*InvocationMetadata*/, const std::string& FilenameLex); -std::vector ValueStringToDoubleList(std::string input); +std::vector SplitString(std::string input, std::string delimiter); diff --git a/src/Simple_Meshing.cc b/src/Simple_Meshing.cc index e5732168..7cf2e636 100644 --- a/src/Simple_Meshing.cc +++ b/src/Simple_Meshing.cc @@ -646,14 +646,17 @@ Minimally_Amalgamate_Contours( auto a_v1_it = std::prev(std::end(amal.points)); auto a_v2_it = std::begin(amal.points); - // Disregard this edge if any of the verts are ficticious. - if( is_a_pseudo_vert(a_v1_it) - || is_a_pseudo_vert(a_v2_it) ){ - continue; - } const auto a_end = std::end(amal.points); while(a_v2_it != a_end){ + // Disregard this edge if any of the verts are ficticious. + if( is_a_pseudo_vert(a_v1_it) + || is_a_pseudo_vert(a_v2_it) ){ + YLOGWARN("We are here again!"); + a_v1_it = a_v2_it; + ++a_v2_it; + continue; + } auto b_v1_it = std::prev(std::end(b_it->get().points)); auto b_v2_it = std::begin(b_it->get().points); const auto b_end = std::end(b_it->get().points); From ea3aa5458b046bec2319bede22b81dcf0d2601f3 Mon Sep 17 00:00:00 2001 From: Renu Rajamagesh <60906270+Renu-R@users.noreply.github.com> Date: Sun, 9 Apr 2023 12:37:27 -0700 Subject: [PATCH 13/20] Nested contour handling * check areas * attempting capping and tiling * bug fix * handled one-to-two case * two-to-two case * compiles * commented outer tiling for debugging * tried some test shapes * changed capping of pipe like structures * bug fix * removed redundant code * removed log statements * removed more log statements * addressed comments from PR * fixed spacing * fixed spacing * more changes from PR --- src/Operations/ConvertContoursToMeshes.cc | 161 ++++++++++++++++++---- 1 file changed, 138 insertions(+), 23 deletions(-) diff --git a/src/Operations/ConvertContoursToMeshes.cc b/src/Operations/ConvertContoursToMeshes.cc index c4353ea2..298297db 100644 --- a/src/Operations/ConvertContoursToMeshes.cc +++ b/src/Operations/ConvertContoursToMeshes.cc @@ -396,6 +396,9 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, std::list pairs; const auto set_union_is_empty = [](const std::set &A, const std::set &B) -> bool { + //if both sets are empty their union is not empty + if (A.empty() && B.empty()) return false; + for(const auto &a : A){ if(B.count(a) != 0){ return false; @@ -468,7 +471,6 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, // Convert from integer numbering to direct pairing info. for(const auto &p : pairs){ - pairings.emplace_back(); for(const auto &u : p.upper){ pairings.back().upper.emplace_back( *std::next( std::begin(m_cops), u ) ); @@ -520,6 +522,22 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, return; }; + const auto contours_are_enclosed = [&](const plane &pln, cop_refw_t A, cop_refw_t B) -> bool { + return (projected_contours_overlap(pln, A, pln, B) && !projected_contours_intersect(pln, A, pln, B)); + }; + + const auto add_faces_to_mesh = [&](cop_refw_t cop_refw_A, cop_refw_t cop_refw_B,std::vector> new_faces) -> void{ + const auto old_face_count = amesh.vertices.size(); + for(const auto &p : cop_refw_A.get().points) amesh.vertices.emplace_back(p); + for(const auto &p : cop_refw_B.get().points) amesh.vertices.emplace_back(p); + for(const auto &fs : new_faces){ + const auto f_A = static_cast(fs[0] + old_face_count); + const auto f_B = static_cast(fs[1] + old_face_count); + const auto f_C = static_cast(fs[2] + old_face_count); + amesh.faces.emplace_back( std::vector{{f_A, f_B, f_C}} ); + } + }; + // Estimate connectivity and append triangles. for(auto &pcs : pairings){ const auto N_upper = pcs.upper.size(); @@ -527,25 +545,130 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, //YLOGINFO("Processing contour map from " << N_upper << " to " << N_lower); if( (N_upper != 0) && (N_lower == 0) ){ - for(const auto &cop_refw : pcs.upper) close_hole_in_floor(cop_refw); + //If the upper plane contains 2 contours and one is enclosed in the other, + //tile the contours together instead of closing the floor + //this routine is for pipe like structures. + if(N_upper == 2){ + if(contours_are_enclosed(*m_cp_it, pcs.upper.front(), pcs.upper.back())){ + auto new_faces = Estimate_Contour_Correspondence(pcs.upper.front(), pcs.upper.back()); + add_faces_to_mesh(pcs.upper.front(), pcs.upper.back(), new_faces); + } + }else{ + for(const auto &cop_refw : pcs.upper) close_hole_in_floor(cop_refw); + } - }else if( (N_upper == 0) && (N_lower == 1) ){ - for(const auto &cop_refw : pcs.lower) close_hole_in_roof(cop_refw); + }else if( (N_upper == 0) && (N_lower != 0) ){ + //if the upper plane contains 2 contours and one is enclosed in the other, + //tile the contours together instead of closing the roof + //this routine is for pipe like structures. + if (N_lower == 2){ + if(contours_are_enclosed(*l_cp_it, pcs.lower.front(), pcs.lower.back())){ + auto new_faces = Estimate_Contour_Correspondence(pcs.lower.front(), pcs.lower.back()); + add_faces_to_mesh(pcs.lower.front(), pcs.lower.back(), new_faces); + } + }else{ + for(const auto &cop_refw : pcs.lower) close_hole_in_roof(cop_refw); + } }else if( (N_upper == 1) && (N_lower == 1) ){ auto new_faces = Estimate_Contour_Correspondence(pcs.upper.front(), pcs.lower.front()); + add_faces_to_mesh(pcs.upper.front(), pcs.lower.front(), new_faces); + + }else if((N_upper == 2) && (N_lower == 1) ){ + //check if the upper plane contains enclosed contours + if(contours_are_enclosed(*m_cp_it, pcs.upper.front(), pcs.upper.back())){ + //get contour areas to determine which contour is enclosed by the other + auto contour = std::begin(pcs.upper); + const auto c1_area = std::abs(contour->get().Get_Signed_Area()); + ++contour; + const auto c2_area = std::abs(contour->get().Get_Signed_Area()); + + //cap the smaller contour and tile the larger contour with the lower plane contour + if (c1_area < c2_area){ + close_hole_in_floor(pcs.upper.front()); + auto new_faces = Estimate_Contour_Correspondence(pcs.upper.back(), pcs.lower.front()); + add_faces_to_mesh(pcs.upper.back(), pcs.lower.front(), new_faces); + }else{ + close_hole_in_floor(pcs.upper.back()); + auto new_faces = Estimate_Contour_Correspondence(pcs.upper.front(), pcs.lower.front()); + add_faces_to_mesh(pcs.upper.front(), pcs.lower.front(), new_faces); + } + }else{ + //do two-to-one branching + } - const auto old_face_count = amesh.vertices.size(); - for(const auto &p : pcs.upper.front().get().points) amesh.vertices.emplace_back(p); - for(const auto &p : pcs.lower.front().get().points) amesh.vertices.emplace_back(p); - for(const auto &fs : new_faces){ - const auto f_A = static_cast(fs[0] + old_face_count); - const auto f_B = static_cast(fs[1] + old_face_count); - const auto f_C = static_cast(fs[2] + old_face_count); - amesh.faces.emplace_back( std::vector{{f_A, f_B, f_C}} ); + }else if( (N_upper == 1) && (N_lower == 2) ){ + //check if the lower plane contains enclosed contours + if(contours_are_enclosed(*l_cp_it, pcs.lower.front(), pcs.lower.back())){ + //get contour areas to determine which contour is enclosed by the other + auto contour = std::begin(pcs.lower); + const auto c1_area = std::abs(contour->get().Get_Signed_Area()); + ++contour; + const auto c2_area = std::abs(contour->get().Get_Signed_Area()); + + //cap the smaller contour and tile the larger contour with the lower plane contour + if (c1_area < c2_area){ + close_hole_in_roof(pcs.lower.front()); + auto new_faces = Estimate_Contour_Correspondence(pcs.lower.back(), pcs.upper.front()); + add_faces_to_mesh(pcs.lower.back(), pcs.upper.front(), new_faces); + }else{ + close_hole_in_roof(pcs.lower.back()); + auto new_faces = Estimate_Contour_Correspondence(pcs.lower.front(), pcs.upper.front()); + add_faces_to_mesh(pcs.lower.front(), pcs.upper.front(), new_faces); + } + }else{ + //do one-to-two branching } + }else{ //YLOGINFO("Performing N-to-N meshing.."); + + //routine for hollow structures with an inner contour and an outer contour on both planes + if((N_upper == 2) && (N_lower == 2)){ + //check if both planes have enclosed contours + if (contours_are_enclosed(*m_cp_it, pcs.upper.front(), pcs.upper.back()) + && contours_are_enclosed(*l_cp_it, pcs.lower.front(), pcs.lower.back())){ + //check areas to determine which the inner and outer contours are on both planes + + //assume 2nd contour is the inner to begin with + auto upper_inner = pcs.upper.back(); + auto upper_outer = pcs.upper.front(); + + auto contour = std::begin(pcs.upper); + auto c1_area = std::abs(contour->get().Get_Signed_Area()); + ++contour; + auto c2_area = std::abs(contour->get().Get_Signed_Area()); + + if (c1_area < c2_area){ + upper_inner = pcs.upper.front(); + upper_outer = pcs.upper.back(); + } + + //do the same for the lower plane + auto lower_inner = pcs.lower.back(); + auto lower_outer = pcs.lower.front(); + + contour = std::begin(pcs.lower); + c1_area = std::abs(contour->get().Get_Signed_Area()); + ++contour; + c2_area = std::abs(contour->get().Get_Signed_Area()); + + if (c1_area < c2_area){ + lower_inner = pcs.lower.front(); + lower_outer = pcs.lower.back(); + } + + //connect inner contours together and outer contours together + auto new_faces = Estimate_Contour_Correspondence(lower_inner, upper_inner); + add_faces_to_mesh(lower_inner, upper_inner, new_faces); + + new_faces = Estimate_Contour_Correspondence(lower_outer, upper_outer); + add_faces_to_mesh(lower_outer, upper_outer, new_faces); + //move to next iteration of for loop since we have tiled it + continue; + } + } + auto ofst_upper = m_cp_it->N_0 * contour_sep * -0.49; auto ofst_lower = m_cp_it->N_0 * contour_sep * 0.49; @@ -567,19 +690,11 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, } */ auto new_faces = Estimate_Contour_Correspondence(std::ref(amal_upper), std::ref(amal_lower)); + add_faces_to_mesh(std::ref(amal_upper), std::ref(amal_lower), new_faces); + } - const auto old_face_count = amesh.vertices.size(); - for(const auto &p : amal_upper.points) amesh.vertices.emplace_back(p); - for(const auto &p : amal_lower.points) amesh.vertices.emplace_back(p); - for(const auto &fs : new_faces){ - const auto f_A = static_cast(fs[0] + old_face_count); - const auto f_B = static_cast(fs[1] + old_face_count); - const auto f_C = static_cast(fs[2] + old_face_count); - amesh.faces.emplace_back( std::vector{{f_A, f_B, f_C}} ); - } - } } - + //caps contours that have no corresponding contours on the lower plane if (cap_roof_of_m_cops) { for (auto &cop : m_cops) { close_hole_in_roof(cop); From cb12c4600c7e983ef3b4bfcb6ee2a9e368569358 Mon Sep 17 00:00:00 2001 From: tiffanymatthe <49164666+tiffanymatthe@users.noreply.github.com> Date: Mon, 10 Apr 2023 07:15:44 -0700 Subject: [PATCH 14/20] Manifold metric (#15) --- scripts/compare_all_meshes.sh | 10 +-- src/Operations/CompareMeshes.cc | 153 ++++++++++++++++++++++++++++++-- 2 files changed, 153 insertions(+), 10 deletions(-) diff --git a/scripts/compare_all_meshes.sh b/scripts/compare_all_meshes.sh index 91990064..8e63ddd6 100755 --- a/scripts/compare_all_meshes.sh +++ b/scripts/compare_all_meshes.sh @@ -54,7 +54,7 @@ declare -a shape_labels=("Sphere(10)" "aa_box(1.0,2.0,4.0)" "tri-force" "chamfer declare -a shapes=("Sphere(10)" "aa_box(1.0,2.0,4.0)" "$tri_force_shape" "chamfer_subtract(2.0){aa_box(8.0,20.0,20.0);sphere(12.0);}") declare -a resolutions=("$default_res" "0.25,0.25,0.75" "0.25,0.25,0.75" "$default_res") -printf "| %-20s | %-13s | %-13s | %-15s | %-15s | %-10s | %-15s | %-15s | %-10s| \n\n" Shape "Hausdorff 1" "Hausdorff 2" "Surface Area 1" "Surface Area 2" "Area Diff" "Volume 1" "Volume 2" "Volume Diff" +printf "| %-20s | %-13s | %-13s | %-15s | %-15s | %-10s | %-10s | %-10s | %-10s| %-8s | %-8s | \n\n" Shape "Hausdorff 1" "Hausdorff 2" "Surface Area 1" "Surface Area 2" "Area Diff" "Vol 1" "Vol 2" "Vol Diff" "Manifold 1" "Manifold 2" for i in ${!shapes[@]} do @@ -106,7 +106,7 @@ do SA_DIFF_LINE=$(grep "SURFACE AREA (%) difference" <<< "$OUTPUT") VOL_LINE=$(grep "VOLUME: First mesh" <<< "$OUTPUT") VOL_DIFF_LINE=$(grep "VOLUME (%) difference" <<< "$OUTPUT") - + MA_LINE=$(grep "MANIFOLD: First mesh" <<< "$OUTPUT") HD1=$(sed -r "s/.*([0-9]+\.[0-9]+) or.*/\1/g" <<< $HD_LINE) HD2=$(sed -r "s/.*or ([0-9]+\.[0-9]+).*/\1/g" <<< $HD_LINE) @@ -121,10 +121,10 @@ do V_DIFF=$(sed -r "s/.*difference: (-?[0-9]+\.[0-9]+).*/\1/g" <<< $VOL_DIFF_LINE) - # echo "$HD1" - # echo "$HD2" + MA1=$(sed -r "s/.*First mesh = ([0-9]+).*/\1/g" <<< $MA_LINE) + MA2=$(sed -r "s/.*second mesh = ([0-9]+).*/\1/g" <<< $MA_LINE) - printf "| %-20s | %-13.3f | %-13.3f | %-15.3f | %-15.3f | %-10.3f | %-15.3f | %-15.3f | %-10.3f |\n\n" $shape_label $HD1 $HD2 $SA1 $SA2 $SA_DIFF $V1 $V2 $V_DIFF + printf "| %-20s | %-13.3f | %-13.3f | %-15.3f | %-15.3f | %-10.3f | %-10.3f | %-10.3f | %-10.3f | %-8s | %-8s |\n\n" $shape_label $HD1 $HD2 $SA1 $SA2 $SA_DIFF $V1 $V2 $V_DIFF $MA1 $MA2 # echo "Finished computing differences for $shape" diff --git a/src/Operations/CompareMeshes.cc b/src/Operations/CompareMeshes.cc index 83d3ef83..04aa3004 100644 --- a/src/Operations/CompareMeshes.cc +++ b/src/Operations/CompareMeshes.cc @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -48,6 +49,140 @@ OperationDoc OpArgDocCompareMeshes(){ return out; } +std::vector> get_face_edges(const std::vector &face) { + std::vector> edges({ + std::set({face[0], face[1]}), + std::set({face[1], face[2]}), + std::set({face[2], face[0]}) + }); + + return edges; +} + +// will consolidate different indices that point to same vertex +// will not modify mesh +std::vector> GetCleanFaces(const std::shared_ptr &mesh) { + std::map, std::vector> vertex_to_indices; + int i = 0; + auto &vertices = mesh.get()->meshes.vertices; + for (auto &vertex : vertices) { + vertex_to_indices[vertex].emplace_back(i); + if (vertex_to_indices[vertex].size() > 1) YLOGINFO("Vertex " << vertex << " has " << vertex_to_indices[vertex].size()); + i++; + } + + std::vector> new_faces; + + // loop through faces and access vertex for each one. pick first index + for (auto &face : mesh.get()->meshes.faces) { + std::vector new_face; + for (auto &vertex_index : face) { + new_face.emplace_back(vertex_to_indices[vertices[vertex_index]].front()); + } + new_faces.emplace_back(new_face); + } + + return new_faces; +} + +// returns true if mesh is edge manifold +// it is edge manifold when every edge is connected to 2 faces +// https://www.mathworks.com/help/lidar/ref/surfacemesh.isedgemanifold.html +bool IsEdgeManifold(const std::shared_ptr &mesh) { + std::map, int> edge_counts; + auto mesh_faces = GetCleanFaces(mesh); + int face_count = 0; + // YLOGINFO("Printing vertices for mesh"); + // for (auto &face : mesh.get()->meshes.faces) { + // YLOGINFO(face_count << " face vertex indices " << face[0] << ", " << face[1] << ", " << face[2] << " with size " << face.size()); + // face_count++; + // } + + face_count = 0; + for (auto &face : mesh_faces) { + // assumes each face has 3 vertices + auto edges = get_face_edges(face); + + for (auto &edge : edges) { + edge_counts[edge] += 1; + // auto itt = edge.begin(); + // YLOGINFO("Edge count = " << edge_counts[edge] << " for edge " << *itt << ", " << *(std::next(itt))); + if (edge_counts[edge] > 2) { + // auto it = edge.begin(); + // YLOGINFO("edge count > 2 for edge " << *it << ", " << *(std::next(it)) << " for face #" << face_count); + return false; + } + } + face_count++; + } + + for (auto &key_value_pair : edge_counts) { + if (key_value_pair.second != 2) { + // YLOGINFO("edge count is " << key_value_pair.second); + return false; + } + } + + return true; +} + +// returns true if mesh if vertex manifold +// it is vertex manifold when each vertex's faces form an open or closed fan +// https://www.mathworks.com/help/lidar/ref/surfacemesh.isvertexmanifold.html +bool IsVertexManifold(const std::shared_ptr&mesh) { + // for each face, find faces with same edges and add to search + // keep searching until you hit a face you look for before (closed) or ended + // is there still faces unsearched? + + // store vertex to face hashmap + // for each vertex, store edge to list of faces + // start with a face + // make sure two edges have adjacent faces (and add to queue) + // search adjacent faces (only one edge now) + // have a list of visited edges to avoid re-searching + // keep searching until no more faces (make sure number of faces searched is equal to total faces) + + using face_type = std::vector; + auto mesh_faces = GetCleanFaces(mesh); + std::map> vertex_to_faces; + + for (auto &mesh_face : mesh_faces) { + vertex_to_faces[mesh_face[0]].emplace_back(mesh_face); + vertex_to_faces[mesh_face[1]].emplace_back(mesh_face); + vertex_to_faces[mesh_face[2]].emplace_back(mesh_face); + } + + for (auto &vertex_faces_pair : vertex_to_faces) { + auto vertex = vertex_faces_pair.first; + auto faces = vertex_faces_pair.second; + + std::map, std::vector> edge_to_faces; + for (auto &face : faces) { + auto edges = get_face_edges(face); + for (auto &edge : edges) edge_to_faces[edge].emplace_back(face); + } + + std::set visited_faces; + std::queue face_queue; + face_queue.push(faces[0]); + + while (face_queue.size() != 0) { + auto face = face_queue.front(); + face_queue.pop(); + if (visited_faces.find(face) != visited_faces.end()) continue; // already visited + visited_faces.emplace(face); + auto edges = get_face_edges(face); + for (auto &edge : edges) { + auto adjacent_faces = edge_to_faces[edge]; + for (auto &adj_face : adjacent_faces) face_queue.push(adj_face); + } + } + + if (visited_faces.size() != faces.size()) return false; + } + return true; +} + bool CompareMeshes(Drover &DICOM_data, const OperationArgPkg& OptArgs, std::map& /*InvocationMetadata*/, @@ -167,16 +302,24 @@ bool CompareMeshes(Drover &DICOM_data, P_B.x*P_A.y*P_C.z + P_A.x*P_B.y*P_C.z)/(6.0); } - + bool v_manifold_1 = IsVertexManifold(mesh1); + bool v_manifold_2 = IsVertexManifold(mesh2); + bool e_manifold_1 = IsEdgeManifold(mesh1); + bool e_manifold_2 = IsEdgeManifold(mesh2); + FUNCINFO("Vertex manifoldness: " << v_manifold_1 << " and " << v_manifold_2); + FUNCINFO("Edge manifoldness: " << e_manifold_1 << " and " << e_manifold_2); + + bool manifold1 = v_manifold_1 && e_manifold_1; + bool manifold2 = v_manifold_2 && e_manifold_2; FUNCINFO("HAUSDORFF DISTANCE: " << max_distance << " or " << second_max_distance); FUNCINFO("SURFACE AREA: First mesh = " << mesh1->meshes.surface_area() << ", second mesh = " << mesh2->meshes.surface_area()); FUNCINFO("SURFACE AREA (%) difference: " << (mesh1->meshes.surface_area() - mesh2->meshes.surface_area())*100/mesh1->meshes.surface_area()); FUNCINFO("VOLUME: First mesh = " < &cop, const vec3 &cont_normal) { + //The rotation matrix R that rotates the contour normal to the direction of + //the z-axis can be applied to all the contour points to move them to the XY plane + auto unitNormContour = cont_normal.unit(); + + vec3 unitNormDesired = vec3(0, 0, 1); + if((unitNormDesired.Cross(unitNormContour)).length() == 0) { + return; + }; + + //obtaining R matrix as described here https://math.stackexchange.com/a/2672702 + num_array I = num_array().identity(3); + auto k = unitNormContour + unitNormDesired; + auto K = k.to_num_array(); + auto scale = 2 / (K.transpose() * K).read_coeff(0, 0); + auto R = K * K.transpose() *= (scale); + R = R - I; + //apply rotation to all points + for(auto &point : cop.points) { + point = (R * point.to_num_array()).to_vec3(); + } +} + +// retrieve element in stack underneath top element +int +Next_To_Top(std::stack &s) { + auto p = s.top(); + s.pop(); + auto res = s.top(); + s.push(p); + return res; +} + +// returns 0 if collinear, 1 if CW, 2 if CCW +// CW if looking from p, we go CW from q to r +// CCW if looking from p, we go CCW from q to r +int +orientation(const vec3 p, const vec3 q, const vec3 r) { + auto val = (q.y - p.y) * (r.x - q.x) - (q.x - p.x) * (r.y - q.y); + if(val == 0) + return 0; + return (val > 0) ? 1 : 2; +} + +// MAIN METHODS + +// returns contour of points of convex hull by following https://www.geeksforgeeks.org/convex-hull-using-graham-scan/ +contour_of_points +Contour_From_Convex_Hull_2(std::vector> &cops, const vec3 &cont_normal) { + if(cops.size() == 0) { + YLOGINFO("Returning empty contour of points."); + return contour_of_points(); + } + + contour_of_points merged_cop; + + for(auto &cop : cops) { + merged_cop.points.insert(merged_cop.points.end(), cop.points.begin(), cop.points.end()); + } + + contour_of_points rotated_cop(merged_cop); + Rotate_Contour_Onto_XY(rotated_cop, cont_normal); + + // creating vector from list for quick access + std::vector> points(rotated_cop.points.begin(), rotated_cop.points.end()); + + // create vector of associated indices to manipulate when finding convex hull + // useful for recreating contour of points for convex hull with same order as input + std::vector indices(points.size()); + std::iota(std::begin(indices), std::end(indices), 0); + + // returns point based on index from indices[i] + const auto get_point = [&](int i) -> vec3 { + return points[indices[i]]; + }; + + // find leftest and lowest point + double ymin = get_point(0).y; + int min_i = 0; + + for(int i = 1; i < indices.size(); ++i) { + int y = get_point(i).y; + if(y < ymin || (ymin == y && get_point(i).x < get_point(min_i).x)) { + ymin = y; + min_i = i; + } + } + std::iter_swap(indices.begin() + 0, indices.begin() + min_i); + + // sort rest of points to form a simple closed path CCW + const auto p0 = get_point(0); + + // returns true if point A is closer to p0 than point B + const auto compare_points = [&](const int index_A, const int index_B) -> bool { + const auto point_A = points[index_A]; + const auto point_B = points[index_B]; + const int o = orientation(p0, point_A, point_B); + + if(o == 0) + return (point_B.sq_dist(p0) > point_A.sq_dist(p0)) ? true : false; + + return (o == 2) ? true : false; + }; + + std::sort(indices.begin() + 1, indices.end(), [&](const int A, const int B) -> bool { + return compare_points(A, B); + }); + + int m = 1; // size of array for concave hull + + // only keep furthest point if multiple points form the same angle with p0 + for(int i = 1; i < points.size(); ++i) { + while(i < points.size() - 1 && orientation(p0, get_point(i), get_point(i + 1)) == 0) { + i++; + } + + indices[m] = indices[i]; + m++; + } + + if(m < 3) { + throw std::runtime_error("Convex hull is not possible for less than 3 points"); + } + + // remove concave points + std::stack s; + s.push(0); + s.push(1); + s.push(2); + + for(int i = 3; i < m; ++i) { + while(s.size() > 1 && orientation(get_point(Next_To_Top(s)), get_point(s.top()), get_point(i)) != 2) { + s.pop(); + } + s.push(i); + } + + // stack contains output points + std::vector kept_indices; + while(!s.empty()) { + kept_indices.emplace_back(indices[s.top()]); + s.pop(); + } + + contour_of_points convex_hull_contour; + + std::vector> merged_points(merged_cop.points.begin(), merged_cop.points.end()); + + for(auto i : kept_indices) { + convex_hull_contour.points.emplace_back(merged_points[i]); + } + + convex_hull_contour.closed = true; + + return convex_hull_contour; +} + +// assumes there are only two contours in cops +// modifies convex hull by adding midpoints and including more original contour points +// returns closed modified convex hull, two non-convex hull contours (one for each original contour), and list of midpoints +// the non-convex hull contours start and end with closest point on convex hull +std::tuple, contour_of_points, contour_of_points, std::vector>> +Modify_Convex_Hull(const contour_of_points &convex_hull, const std::vector> cops, + const vec3 &pseudo_vert_offset) { + vec3 start_A; + vec3 end_A; + vec3 start_B; + vec3 end_B; + vec3 midpoint_for_end_A; + vec3 midpoint_for_start_A; + + bool found_first = false; + + int prev_contour_index = -1; + vec3 prev_point; + int first_point_contour_index = 0; + + for(auto it = convex_hull.points.begin(); it != convex_hull.points.end(); it++) { + int curr_contour_index = 0; + for(auto &cop : cops) { + if(std::find(cop.points.begin(), cop.points.end(), *it) != cop.points.end()) { + if(prev_contour_index != -1 && prev_contour_index != curr_contour_index) { + vec3 midpoint = (*it + prev_point) / 2 + pseudo_vert_offset; + if(!found_first) { + end_A = prev_point; + start_B = *it; + midpoint_for_end_A = midpoint; + found_first = true; + } else { + end_B = prev_point; + start_A = *it; + midpoint_for_start_A = midpoint; + } + } + break; + } + curr_contour_index += 1; + } + if(it == convex_hull.points.begin()) { + first_point_contour_index = curr_contour_index; + } + prev_contour_index = curr_contour_index; + prev_point = *it; + + // check last point and first point + if(it == std::prev(convex_hull.points.end()) && first_point_contour_index != curr_contour_index) { + end_B = *it; + start_A = *(convex_hull.points.begin()); + vec3 midpoint = (*it + *(convex_hull.points.begin())) / 2 + pseudo_vert_offset; + midpoint_for_start_A = midpoint; + break; + } + } + + // swap points depending on orientation + bool convex_hull_ccw = convex_hull.Is_Counter_Clockwise(); + if(cops[first_point_contour_index].Is_Counter_Clockwise() != convex_hull_ccw) { + std::swap(start_A, end_A); + std::swap(midpoint_for_start_A, midpoint_for_end_A); + } + if(cops[1 - first_point_contour_index].Is_Counter_Clockwise() != convex_hull_ccw) { + std::swap(start_B, end_B); + } + + auto d1 = start_A.sq_dist(midpoint_for_start_A); + auto d2 = end_A.sq_dist(midpoint_for_end_A); + if (std::max(d1,d2) > 2 * std::min(d1,d2)) throw std::runtime_error("Complex 2 to 1 meshing is not suitable for this contour. Aborted"); + + const auto get_convex_and_non_convex_points = + [&](const int contour_index, vec3 start_point, + vec3 end_point) -> std::tuple>, std::list>> { + std::list> split_section_1; + std::list> consecutive_section_2; + std::list> split_section_3; + + bool is_consecutive_section_convex_hull = false; + bool found_start_point = false; + bool found_end_point = false; + bool filling_split_section = true; + bool on_split_section_3 = false; + + for(auto point : cops[contour_index].points) { + if(point == start_point) { + if(found_end_point) { // already found end point before, so doing 3rd portion + on_split_section_3 = true; + filling_split_section = true; + } else { + filling_split_section = false; + } + found_start_point = true; + } + + if(point == end_point) { + if(found_start_point) { // already found start point, so doing 3rd portion + is_consecutive_section_convex_hull = true; + filling_split_section = true; + on_split_section_3 = true; + consecutive_section_2.emplace_back(point); // don't want the end point to be in the third portion + } else { + filling_split_section = false; + split_section_1.emplace_back(point); + } + found_end_point = true; + continue; // since we placed the point in the previous section already + } + + if(filling_split_section) { + if(on_split_section_3) { + split_section_3.emplace_back(point); + } else { + split_section_1.emplace_back(point); + } + } else { + consecutive_section_2.emplace_back(point); + } + } + + split_section_3.insert(split_section_3.end(), split_section_1.begin(), split_section_1.end()); + + std::list> &convex_hull_points = split_section_3; + std::list> &other_points = consecutive_section_2; + + if(is_consecutive_section_convex_hull) { + std::swap(convex_hull_points, other_points); + } + + // add start and end points to non convex hull points to ensure proper meshing + other_points.emplace_front(convex_hull_points.back()); + other_points.emplace_back(convex_hull_points.front()); + + return std::make_tuple(convex_hull_points, other_points); + }; + + // create better convex hull by chaining original points + auto [convex_hull_A, other_A] = get_convex_and_non_convex_points(first_point_contour_index, start_A, end_A); + auto [convex_hull_B, other_B] = get_convex_and_non_convex_points(1 - first_point_contour_index, start_B, end_B); + + // flip contour B if end point is closer to midpoint for end A + if(end_B.sq_dist(midpoint_for_end_A) < start_B.sq_dist(midpoint_for_end_A)) { + convex_hull_B.reverse(); + other_B.reverse(); + } + + // add midpoints in contour + contour_of_points modified_convex_cop; + modified_convex_cop.points.insert(modified_convex_cop.points.end(), convex_hull_A.begin(), convex_hull_A.end()); + modified_convex_cop.points.emplace_back(midpoint_for_end_A); + modified_convex_cop.points.insert(modified_convex_cop.points.end(), convex_hull_B.begin(), convex_hull_B.end()); + modified_convex_cop.points.emplace_back(midpoint_for_start_A); + modified_convex_cop.closed = true; + + std::vector> midpoints; + midpoints.emplace_back(midpoint_for_end_A); + midpoints.emplace_back(midpoint_for_start_A); + + return std::make_tuple(modified_convex_cop, other_A, other_B, midpoints); +} + +// creates faces by connecting non convex hull points to the midpoints based on distance +// returns faces and associated ordered points +std::tuple>, std::vector>> +Mesh_Inner_Points_With_Midpoints(const std::list> &left_points, + const std::list> &right_points, + const std::vector> &midpoints, const vec3 &pseudo_vert_offset) { + if(midpoints.size() != 2) { + YLOGWARN("Unable to handle " << midpoints.size() << " midpoints at this time."); + throw std::runtime_error("Unable to handle !=2 midpoints at this time."); + } + + // compare distances by using midpoints on the same plane + auto midpoint1_flattened = midpoints[0] - pseudo_vert_offset; + auto midpoint2_flattened = midpoints[1] - pseudo_vert_offset; + + int midpoint1_position = left_points.size() + right_points.size(); + int midpoint2_position = midpoint1_position + 1; + + std::vector> faces; + + // adds face + // iterators must be of same vector + const auto add_points_with_midpoint = + [&](std::list>::iterator point1_it, std::list>::iterator point2_it, + std::list>::iterator beg_it, const int offset, const int midpoint_index) -> void { + const auto v_a = static_cast(offset + std::distance(beg_it, point1_it)); + const auto v_b = static_cast(offset + std::distance(beg_it, point2_it)); + const auto v_c = static_cast(midpoint_index); + faces.emplace_back(std::array { { v_a, v_b, v_c } }); + }; + + const auto add_point_with_two_midpoints = [&](std::list>::iterator point1_it, + std::list>::iterator beg_it, const int offset, + const int midpoint1_index, const int midpoint2_index) -> void { + const auto v_a = static_cast(offset + std::distance(beg_it, point1_it)); + const auto v_b = static_cast(midpoint1_index); + const auto v_c = static_cast(midpoint2_index); + faces.emplace_back(std::array { { v_a, v_b, v_c } }); + }; + + const auto create_faces_with_points = [&](std::list> points, const int offset) -> void { + bool switched = false; + auto closer_midpoint = midpoint1_flattened; + auto closer_position = midpoint1_position; + auto further_midpoint = midpoint2_flattened; + auto further_position = midpoint2_position; + + if(points.begin()->sq_dist(closer_midpoint) > points.begin()->sq_dist(further_midpoint)) { + std::swap(closer_midpoint, further_midpoint); + std::swap(closer_position, further_position); + } + + for(auto it = std::next(points.begin()); it != points.end(); it++) { + // assumes once switched, will not go back + auto point = *it; + if(point.sq_dist(closer_midpoint) < point.sq_dist(further_midpoint) && !switched) { + // create face with prev point and closer_midpoint + add_points_with_midpoint(it, std::prev(it), points.begin(), offset, closer_position); + } else { + if(!switched) { + // connect to both points (two faces) + add_points_with_midpoint(it, std::prev(it), points.begin(), offset, closer_position); + add_point_with_two_midpoints(it, points.begin(), offset, closer_position, further_position); + switched = true; + } else { + // connect to prev point and further_midpoint + add_points_with_midpoint(it, std::prev(it), points.begin(), offset, further_position); + } + } + } + }; + + create_faces_with_points(left_points, 0); + create_faces_with_points(right_points, left_points.size()); + + // make vector of all points (order is important for faces to be parsed correctly) + std::vector> all_points(left_points.begin(), left_points.end()); + all_points.insert(all_points.end(), right_points.begin(), right_points.end()); + all_points.emplace_back(midpoints[0]); + all_points.emplace_back(midpoints[1]); + + return std::make_tuple(faces, all_points); +} + +// meshes 2 to 1 (will need to mesh convex hull contour with other contour outside of this routine) +// returns non-convex hull faces, non-convex hull points, and convex hull contour (with midpoints) +std::tuple>, std::vector>, contour_of_points> +Mesh_With_Convex_Hull_2(const std::list>> &A, + const vec3 &ortho_unit_A, const vec3 &pseudo_vert_offset) { + if(A.size() != 2) { + throw std::runtime_error("Convex hull is currently only possible for 2 contours. Aborted."); + } + std::vector> list_of_cops; + for(auto &cop : A) { + auto new_cop = cop.get(); + if (new_cop.points.front() == new_cop.points.back()) new_cop.points.pop_back(); + list_of_cops.emplace_back(new_cop); + } + auto convex_hull_cop = Contour_From_Convex_Hull_2(list_of_cops, ortho_unit_A); + auto [modified_convex_cop, left_points, right_points, midpoints] = + Modify_Convex_Hull(convex_hull_cop, list_of_cops, pseudo_vert_offset); + auto [faces_from_non_convex, ordered_non_convex_points] = + Mesh_Inner_Points_With_Midpoints(left_points.points, right_points.points, midpoints, pseudo_vert_offset); + + return std::make_tuple(faces_from_non_convex, ordered_non_convex_points, modified_convex_cop); +} \ No newline at end of file diff --git a/src/Complex_Branching_Meshing.h b/src/Complex_Branching_Meshing.h new file mode 100644 index 00000000..0e464591 --- /dev/null +++ b/src/Complex_Branching_Meshing.h @@ -0,0 +1,61 @@ +//Complex_Branching_Meshing.h - A part of DICOMautomaton 2020. Written by hal clark. + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include //Needed for exit() calls. +#include //Needed for std::pair. +#include +#include + +#include "YgorMath.h" //Needed for vec3 class. + +// returns contour of points of convex hull by following +// https://www.geeksforgeeks.org/convex-hull-using-graham-scan/ +contour_of_points +Contour_From_Convex_Hull_2( + std::vector> &cops, + const vec3 &normal +); + +// assumes there are only two contours in cops +// modifies convex hull by adding midpoints and including more original contour points +// returns closed modified convex hull, two non-convex hull contours (one for each original contour), and list of midpoints +// the non-convex hull contours start and end with closest point on convex hull +std::tuple,contour_of_points,contour_of_points, std::vector>> +Modify_Convex_Hull( + const contour_of_points &convex_hull, + const std::vector> cops, + const vec3 &pseudo_vert_offset +); + +// creates faces by connecting non convex hull points to the midpoints based on distance +// returns faces and associated ordered points +std::tuple>, std::vector>> +Mesh_Inner_Points_With_Midpoints( + const std::list> &left_points, + const std::list> &right_points, + const std::vector> &midpoints, + const vec3 &pseudo_vert_offset +); + +// meshes 2 to 1 (will need to mesh convex hull contour with other contour outside of this routine) +// returns non-convex hull faces, non-convex hull points, and convex hull contour (with midpoints) +std::tuple>, std::vector>, contour_of_points> +Mesh_With_Convex_Hull_2( + const std::list>> &A, + const vec3 &ortho_unit_A, + const vec3 &pseudo_vert_offset +); \ No newline at end of file diff --git a/src/Operations/ConvertContoursToMeshes.cc b/src/Operations/ConvertContoursToMeshes.cc index 298297db..8cf42691 100644 --- a/src/Operations/ConvertContoursToMeshes.cc +++ b/src/Operations/ConvertContoursToMeshes.cc @@ -27,6 +27,7 @@ #include "../Simple_Meshing.h" #include "../Surface_Meshes.h" +#include "../Complex_Branching_Meshing.h" #include "ConvertContoursToMeshes.h" @@ -149,7 +150,6 @@ OperationDoc OpArgDocConvertContoursToMeshes(){ } - bool ConvertContoursToMeshes(Drover &DICOM_data, const OperationArgPkg& OptArgs, std::map& /*InvocationMetadata*/, @@ -537,12 +537,31 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, amesh.faces.emplace_back( std::vector{{f_A, f_B, f_C}} ); } }; + + const auto add_faces_and_vertices = [&]( + std::vector> &new_faces, + std::vector> &points + ) -> void { + const auto old_face_count = amesh.vertices.size(); + // for(const auto &p : points) amesh.vertices.emplace_back(p); + amesh.vertices.insert(amesh.vertices.end(), points.begin(), points.end()); + for(const auto &fs : new_faces){ + const auto f_A = static_cast(fs[0] + old_face_count); + const auto f_B = static_cast(fs[1] + old_face_count); + const auto f_C = static_cast(fs[2] + old_face_count); + amesh.faces.emplace_back( std::vector{{f_A, f_B, f_C}} ); + } + }; // Estimate connectivity and append triangles. for(auto &pcs : pairings){ const auto N_upper = pcs.upper.size(); const auto N_lower = pcs.lower.size(); - //YLOGINFO("Processing contour map from " << N_upper << " to " << N_lower); + + // useful for addition of midpoints in complex branching + auto ofst_upper = m_cp_it->N_0 * contour_sep * -0.49; + auto ofst_lower = m_cp_it->N_0 * contour_sep * 0.49; + // YLOGINFO("Processing contour map from " << N_upper << " to " << N_lower); if( (N_upper != 0) && (N_lower == 0) ){ //If the upper plane contains 2 contours and one is enclosed in the other, @@ -573,7 +592,6 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, }else if( (N_upper == 1) && (N_lower == 1) ){ auto new_faces = Estimate_Contour_Correspondence(pcs.upper.front(), pcs.lower.front()); add_faces_to_mesh(pcs.upper.front(), pcs.lower.front(), new_faces); - }else if((N_upper == 2) && (N_lower == 1) ){ //check if the upper plane contains enclosed contours if(contours_are_enclosed(*m_cp_it, pcs.upper.front(), pcs.upper.back())){ @@ -594,9 +612,16 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, add_faces_to_mesh(pcs.upper.front(), pcs.lower.front(), new_faces); } }else{ - //do two-to-one branching + try { + auto [faces, points, amal_upper] = Mesh_With_Convex_Hull_2(pcs.upper, m_cp_it->N_0, ofst_upper); + add_faces_and_vertices(faces, points); + auto amal_lower = pcs.lower.begin()->get(); + auto new_faces = Estimate_Contour_Correspondence(std::ref(amal_upper), std::ref(amal_lower)); + add_faces_to_mesh(std::ref(amal_upper), std::ref(amal_lower), new_faces); + } catch (const std::runtime_error& error) { + goto generic_n_to_n_meshing; + } } - }else if( (N_upper == 1) && (N_lower == 2) ){ //check if the lower plane contains enclosed contours if(contours_are_enclosed(*l_cp_it, pcs.lower.front(), pcs.lower.back())){ @@ -617,9 +642,16 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, add_faces_to_mesh(pcs.lower.front(), pcs.upper.front(), new_faces); } }else{ - //do one-to-two branching + try { + auto [faces, points, amal_lower] = Mesh_With_Convex_Hull_2(pcs.lower, l_cp_it->N_0, ofst_lower); + add_faces_and_vertices(faces, points); + auto amal_upper = pcs.upper.begin()->get(); + auto new_faces = Estimate_Contour_Correspondence(std::ref(amal_upper), std::ref(amal_lower)); + add_faces_to_mesh(std::ref(amal_upper), std::ref(amal_lower), new_faces); + } catch (const std::runtime_error& error) { + goto generic_n_to_n_meshing; + } } - }else{ //YLOGINFO("Performing N-to-N meshing.."); @@ -667,32 +699,27 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, //move to next iteration of for loop since we have tiled it continue; } + } else { + generic_n_to_n_meshing: + auto amal_upper = Minimally_Amalgamate_Contours(m_cp_it->N_0, ofst_upper, pcs.upper); + auto amal_lower = Minimally_Amalgamate_Contours(m_cp_it->N_0, ofst_lower, pcs.lower); + /* + // Leaving this here for future debugging, for which it will no-doubt be needed... + { + const auto amal_cop_str = amal_upper.write_to_string(); + const auto fname = Get_Unique_Sequential_Filename("/tmp/amal_upper_", 6, ".txt"); + OverwriteStringToFile(amal_cop_str, fname); + } + { + const auto amal_cop_str = amal_lower.write_to_string(); + const auto fname = Get_Unique_Sequential_Filename("/tmp/amal_lower_", 6, ".txt"); + OverwriteStringToFile(amal_cop_str, fname); + } + */ + auto new_faces = Estimate_Contour_Correspondence(std::ref(amal_upper), std::ref(amal_lower)); + add_faces_to_mesh(std::ref(amal_upper), std::ref(amal_lower), new_faces); } - - auto ofst_upper = m_cp_it->N_0 * contour_sep * -0.49; - auto ofst_lower = m_cp_it->N_0 * contour_sep * 0.49; - - // will modify contours in m_cops and l_cops, ok if only processing in one direction - auto amal_upper = Minimally_Amalgamate_Contours(m_cp_it->N_0, ofst_upper, pcs.upper); - auto amal_lower = Minimally_Amalgamate_Contours(m_cp_it->N_0, ofst_lower, pcs.lower); - - /* - // Leaving this here for future debugging, for which it will no-doubt be needed... - { - const auto amal_cop_str = amal_upper.write_to_string(); - const auto fname = Get_Unique_Sequential_Filename("/tmp/amal_upper_", 6, ".txt"); - OverwriteStringToFile(amal_cop_str, fname); - } - { - const auto amal_cop_str = amal_lower.write_to_string(); - const auto fname = Get_Unique_Sequential_Filename("/tmp/amal_lower_", 6, ".txt"); - OverwriteStringToFile(amal_cop_str, fname); - } - */ - auto new_faces = Estimate_Contour_Correspondence(std::ref(amal_upper), std::ref(amal_lower)); - add_faces_to_mesh(std::ref(amal_upper), std::ref(amal_lower), new_faces); } - } //caps contours that have no corresponding contours on the lower plane if (cap_roof_of_m_cops) { diff --git a/src/Simple_Meshing.cc b/src/Simple_Meshing.cc index 7cf2e636..4af379ee 100644 --- a/src/Simple_Meshing.cc +++ b/src/Simple_Meshing.cc @@ -651,8 +651,7 @@ Minimally_Amalgamate_Contours( while(a_v2_it != a_end){ // Disregard this edge if any of the verts are ficticious. if( is_a_pseudo_vert(a_v1_it) - || is_a_pseudo_vert(a_v2_it) ){ - YLOGWARN("We are here again!"); + || is_a_pseudo_vert(a_v2_it) ){ a_v1_it = a_v2_it; ++a_v2_it; continue; From 9ccbbce302ebbb9f0f51984d6a577291a3b6851f Mon Sep 17 00:00:00 2001 From: Xander Naumenko <33139843+misprit7@users.noreply.github.com> Date: Mon, 10 Apr 2023 17:13:39 -0700 Subject: [PATCH 16/20] Graph tiling (#16) * Started graph approach * Added graph construction and conversion back to index list * Added recursive graph search * Finished tenetative graph implementation * Fixed graph search specific bugs * Fixed indexing issues * Got tiling working for big contours * Got tiling working on large shapes * Cleaned up code for PR * Addressed PR comments * Fixed indexing issues --- compile_and_install.sh | 5 +- src/Operations/ConvertContoursToMeshes.cc | 19 +- src/Simple_Meshing.cc | 393 +++++++++++++++++++++- src/Simple_Meshing.h | 4 + 4 files changed, 409 insertions(+), 12 deletions(-) diff --git a/compile_and_install.sh b/compile_and_install.sh index 22a84dea..18aed003 100755 --- a/compile_and_install.sh +++ b/compile_and_install.sh @@ -16,8 +16,8 @@ shopt -s nocasematch # The temporary location in which to build. #BUILDROOT="/home/hal/Builds/DICOMautomaton/" -BUILDROOT="/tmp/dcma_build" -#BUILDROOT="build" +# BUILDROOT="/tmp/dcma_build" +BUILDROOT="build" DISTRIBUTION="auto" # used to tailor the build process for specific distributions/environments. @@ -190,6 +190,7 @@ elif [[ "${DISTRIBUTION}" =~ .*debian.* ]] ; then -DWITH_ASAN=OFF \ -DWITH_TSAN=OFF \ -DWITH_MSAN=OFF \ + -DCMAKE_EXPORT_COMPILE_COMMANDS=1 \ ../ fi JOBS=$(nproc) diff --git a/src/Operations/ConvertContoursToMeshes.cc b/src/Operations/ConvertContoursToMeshes.cc index 8cf42691..d46a164c 100644 --- a/src/Operations/ConvertContoursToMeshes.cc +++ b/src/Operations/ConvertContoursToMeshes.cc @@ -569,7 +569,7 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, //this routine is for pipe like structures. if(N_upper == 2){ if(contours_are_enclosed(*m_cp_it, pcs.upper.front(), pcs.upper.back())){ - auto new_faces = Estimate_Contour_Correspondence(pcs.upper.front(), pcs.upper.back()); + auto new_faces = Tile_Contours(pcs.upper.front(), pcs.upper.back()); add_faces_to_mesh(pcs.upper.front(), pcs.upper.back(), new_faces); } }else{ @@ -582,7 +582,7 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, //this routine is for pipe like structures. if (N_lower == 2){ if(contours_are_enclosed(*l_cp_it, pcs.lower.front(), pcs.lower.back())){ - auto new_faces = Estimate_Contour_Correspondence(pcs.lower.front(), pcs.lower.back()); + auto new_faces = Tile_Contours(pcs.lower.front(), pcs.lower.back()); add_faces_to_mesh(pcs.lower.front(), pcs.lower.back(), new_faces); } }else{ @@ -590,7 +590,8 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, } }else if( (N_upper == 1) && (N_lower == 1) ){ - auto new_faces = Estimate_Contour_Correspondence(pcs.upper.front(), pcs.lower.front()); + /* auto new_faces = Estimate_Contour_Correspondence(pcs.upper.front(), pcs.lower.front()); */ + auto new_faces = Tile_Contours(pcs.upper.front(), pcs.lower.front()); add_faces_to_mesh(pcs.upper.front(), pcs.lower.front(), new_faces); }else if((N_upper == 2) && (N_lower == 1) ){ //check if the upper plane contains enclosed contours @@ -604,11 +605,11 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, //cap the smaller contour and tile the larger contour with the lower plane contour if (c1_area < c2_area){ close_hole_in_floor(pcs.upper.front()); - auto new_faces = Estimate_Contour_Correspondence(pcs.upper.back(), pcs.lower.front()); + auto new_faces = Tile_Contours(pcs.upper.back(), pcs.lower.front()); add_faces_to_mesh(pcs.upper.back(), pcs.lower.front(), new_faces); }else{ close_hole_in_floor(pcs.upper.back()); - auto new_faces = Estimate_Contour_Correspondence(pcs.upper.front(), pcs.lower.front()); + auto new_faces = Tile_Contours(pcs.upper.front(), pcs.lower.front()); add_faces_to_mesh(pcs.upper.front(), pcs.lower.front(), new_faces); } }else{ @@ -634,11 +635,11 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, //cap the smaller contour and tile the larger contour with the lower plane contour if (c1_area < c2_area){ close_hole_in_roof(pcs.lower.front()); - auto new_faces = Estimate_Contour_Correspondence(pcs.lower.back(), pcs.upper.front()); + auto new_faces = Tile_Contours(pcs.lower.back(), pcs.upper.front()); add_faces_to_mesh(pcs.lower.back(), pcs.upper.front(), new_faces); }else{ close_hole_in_roof(pcs.lower.back()); - auto new_faces = Estimate_Contour_Correspondence(pcs.lower.front(), pcs.upper.front()); + auto new_faces = Tile_Contours(pcs.lower.front(), pcs.upper.front()); add_faces_to_mesh(pcs.lower.front(), pcs.upper.front(), new_faces); } }else{ @@ -691,10 +692,10 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, } //connect inner contours together and outer contours together - auto new_faces = Estimate_Contour_Correspondence(lower_inner, upper_inner); + auto new_faces = Tile_Contours(lower_inner, upper_inner); add_faces_to_mesh(lower_inner, upper_inner, new_faces); - new_faces = Estimate_Contour_Correspondence(lower_outer, upper_outer); + new_faces = Tile_Contours(lower_outer, upper_outer); add_faces_to_mesh(lower_outer, upper_outer, new_faces); //move to next iteration of for loop since we have tiled it continue; diff --git a/src/Simple_Meshing.cc b/src/Simple_Meshing.cc index 4af379ee..ee1f6245 100644 --- a/src/Simple_Meshing.cc +++ b/src/Simple_Meshing.cc @@ -16,6 +16,8 @@ #include #include #include +#include //Needed for double max +#include #include //Needed for exit() calls. #include //Needed for std::pair. @@ -34,6 +36,395 @@ #include "Simple_Meshing.h" +// This can be useful when debuggin with gdb to prevent +// Obviously don't enable this when tiling large objects +/* #pragma GCC push_options */ +/* #pragma GCC optimize ("O0") */ + +/* + * Class to represent a single node in a vertex + */ +class node { + public: + // Indices of A and B that this node represents + size_t iA, iB; + // Weight of choosing the next vertex to be on A/B + // In our case this corresponds to the area of (iA, iB, iA+1) + // or (iA, iB, iB+1) respectively + double wA, wB; + + node(size_t ia, size_t ib, + vec3& point_A, + vec3& point_B, + vec3& point_A_next, + vec3& point_B_next){ + iA = ia; iB = ib; + const vec3 ab = point_A - point_B; + // Weights are surface area of potential next faces + /* wA = 0.5 * ab.Cross(point_B - point_A_next).length(); */ + /* wB = 0.5 * ab.Cross(point_A - point_B_next).length(); */ + + // Perimeter + wA = ab.length() + + (point_A - point_A_next).length() + + (point_B - point_A_next).length(); + + wB = ab.length() + + (point_A - point_B_next).length() + + (point_B - point_B_next).length(); + } +}; + +/* + * Convert top and bottom paths to not be overlapping by swapping sections + * + * path_top: top path to modify + * path_bottom: bottom path to modify + */ +void convert_nonoverlapping( + std::vector& path_top, std::vector& path_bottom){ + for(size_t k = 0; k < path_top.size(); ++k){ + if(path_top[k] > path_bottom[k]){ + std::swap(path_top[k], path_bottom[k]); + } + } +} + +/* + * Calculates lowest cost path between assuming a specific start vertex + * + * nodes: graph of nodes + * k: index of row to start at + * path_top: top bounding path to search below + * path_bottom: bottom bounding path to search above + * + * Assumes bottom_path[i] > top_path[i] for all i + * + * returns (minimal path, cost of minimal path) + */ +std::pair, double> +single_path(std::vector>& nodes, size_t k, + std::vector& path_top, std::vector& path_bottom){ + size_t m = nodes.size()/2, n = nodes[0].size(); + + convert_nonoverlapping(path_top, path_bottom); + + // i, j, distance + std::vector> distances(2*m, + std::vector(n, std::numeric_limits::max())); + + std::set> visited = {}; + + // Current minimum distance among nodes + typedef std::pair> node_dist; + std::priority_queue, std::greater> min_dist; + distances[k][0] = 0; + min_dist.push({0, {k, 0}}); + + + // Loop to get distances of each relevant node + while(!min_dist.empty() && (min_dist.top().second.first != m + k || min_dist.top().second.second != n-1)){ + size_t i = min_dist.top().second.first; + size_t j = min_dist.top().second.second; + min_dist.pop(); + if(i == m+k && j == n-1) break; + + // Bottom of the current path + // Since we only directly get the top, have to compute this + size_t bottom = j == n-1 ? m+k : path_bottom[j+1]; + // Traverse next nodes and set their distances based on current + // Try down + // k instead of n because we can't go past endpoint + if(i+1 <= m + k && i+1 >= path_top[j] && i+1 <= bottom && !visited.count({i+1, j})){ + distances[i+1][j] = std::min(distances[i][j] + nodes[i][j].wA, distances[i+1][j]); + min_dist.push({distances[i+1][j], {i+1, j}}); + visited.insert({i+1, j}); + } + // Try right + if(j+1 < n && i >= path_top[j+1] && i <= bottom && !visited.count({i, j+1})){ + distances[i][j+1] = std::min(distances[i][j] + nodes[i][j].wB, distances[i][j+1]); + min_dist.push({distances[i][j+1], {i, j+1}}); + visited.insert({i, j+1}); + } + } + + // Reconstruct best path + std::vector ret(n); + ret[0] = k; + + for(size_t j = n-1, i = m + k; j > 0; --j){ + while(i > 0 && j > 0 && distances[i-1][j] + nodes[i-1][j].wA <= distances[i][j-1] + nodes[i][j-1].wB) --i; + ret[j] = i; + } + // Cost of path is length to end node + return {ret, distances[m + k][n-1]}; +} + +/* + * Calculates optimal path between two indices + * + * nodes: graph of nodes + * i: Upper index + * j: lower index, assumes i, double> +paths_between(std::vector>& nodes, size_t i, size_t j, + std::vector& path_i, std::vector& path_j){ + // If we're at something of length 2 or less, stop recursion + if(j-i < 2) return {std::vector(), std::numeric_limits::max()}; + size_t k = (i + j) / 2; + + // SINGLEPATH(k, G'(i, j)) + auto [path_k, cost_k] = single_path(nodes, k, path_i, path_j); + + // PATHSBETWEEN(i, k) + // PATHSBETWEEN(k, j) + auto [path_before, cost_before] = paths_between(nodes, i, k, path_i, path_k); + auto [path_after, cost_after] = paths_between(nodes, k, j, path_k, path_j); + + // Return minimum of paths calculated + if(cost_before < cost_k && cost_before < cost_after){ + return {path_before, cost_before}; + } + else if(cost_after < cost_k && cost_after < cost_before){ + return {path_after, cost_after}; + } + else { + return {path_k, cost_k}; + } +} + +/* + * Calculates optimal path through graph representation of tiling + * + * nodes: graph of nodes, size 2m x n + * + * Returns minimal path + */ +std::vector all_paths(std::vector>& nodes){ + // Start with G' as best bounds + std::vector top_path(nodes[0].size(), 0); + std::vector bottom_path(nodes[0].size(), nodes.size()); + + // nodes.size() should always be even + // There's some weirdness with the indexing relative to the paper since + // they label their vertices for 0,1,...,m, so have m+1 vertices + size_t m = nodes.size()/2; + + // SINGLEPATH(0,G') + // SINGLEPATH(m,G') + auto [path_start, cost_start] = single_path(nodes, 0, top_path, bottom_path); + auto [path_end, cost_end] = single_path(nodes, m-1, top_path, bottom_path); + + // PATHSBETWEEN(0, m) + auto [path_between, cost_between] = paths_between(nodes, 0, m-1, path_start, path_end); + + // Return minimum of paths calculated + if(cost_start < cost_end && cost_start < cost_between){ + return path_start; + } + else if(cost_end < cost_start && cost_end < cost_between){ + return path_end; + } + else { + return path_between; + } +} + +/* + * Quick and dirty way of visualizing graphs with a colored path, this is purely for + * debugging and definitey isn't the most elegantly written, but the output looks pretty + */ +void print_path(std::vector> nodes, std::vector path){ + printf("Path: "); + for(size_t p : path){ + printf("%zu ", p); + } + printf("\n"); + + path.push_back(path[0] + nodes.size()/2); + const char* white = "\033[0;37m"; + const char* green = "\033[0;32m"; + + for(size_t i = 0; i < nodes.size(); ++i){ + for(size_t j = 0; j < nodes[0].size(); ++j){ + const char* c_node; + if(i >= path[j] && i <= path[j+1]){ + c_node = green; + } else c_node = white; + const char* c_edge; + if(i >= path[j] && i == path[j+1] && j < nodes[0].size()-1){ + c_edge = green; + } else c_edge = white; + printf("-- %so%s --%-6.2lf%s", c_node, c_edge, nodes[i][j].wB, white); + } + printf("\n"); + for(size_t j = 0; j < nodes[0].size(); ++j){ + const char* c; + if(i >= path[j] && i < path[j+1] && path[j] != path[j+1]){ + c = green; + } else c = white; + printf(" %s|%s ", c, white); + } + printf("\n"); + for(size_t j = 0; j < nodes[0].size(); ++j){ + const char* c; + if(i >= path[j] && i < path[j+1] && path[j] != path[j+1]){ + c = green; + } else c = white; + printf(" %s%-6.2lf%s ", c, nodes[i][j].wA, white); + } + printf("\n"); + for(size_t j = 0; j < nodes[0].size(); ++j){ + const char* c; + if(i >= path[j] && i < path[j+1] && path[j] != path[j+1]){ + c = green; + } else c = white; + printf(" %s|%s ", c, white); + } + printf("\n"); + } +} + +// Method taken from here: +// https://www.cs.jhu.edu/~misha/Fall13b/Papers/Fuchs77.pdf +std::vector< std::array > Tile_Contours( + std::reference_wrapper> A, + std::reference_wrapper> B ){ + + contour_of_points contour_A = A.get(); + contour_of_points contour_B = B.get(); + + const auto N_A = contour_A.points.size(); + const auto N_B = contour_B.points.size(); + if( N_A == 0 ){ + throw std::invalid_argument("Contour A contains no vertices. Cannot continue."); + } + if( N_B == 0 ){ + throw std::invalid_argument("Contour B contains no vertices. Cannot continue."); + } + + // Determine contour orientations. Single-vertex contours can take any orientation, so use a reasonable default. + auto ortho_unit_A = vec3(0.0, 0.0, 1.0); + auto ortho_unit_B = vec3(0.0, 0.0, 1.0); + try{ + ortho_unit_A = contour_A.Estimate_Planar_Normal(); + }catch(const std::exception &){} + try{ + ortho_unit_B = contour_B.Estimate_Planar_Normal(); + }catch(const std::exception &){} + + // Ensure the contours have the same orientation. + if(ortho_unit_A.Dot(ortho_unit_B) <= 0.0){ + // Handle special cases + if( (N_A == 1) && (N_B != 1) ){ + ortho_unit_A = ortho_unit_B; + }else if( (N_A != 1) && (N_B == 1) ){ + ortho_unit_B = ortho_unit_A; + + // Otherwise, flip one of the contours, recurse, and adjust the face labels to point to the actual vertex + // layout. + // + // Note: This effectively ignores contour orientation altogether. + }else{ + YLOGWARN("Ignoring adjacent contours with opposite orientations. Recursing.."); + std::reverse( std::begin(contour_B.points), std::end(contour_B.points) ); + auto faces = Tile_Contours(A, std::ref(contour_B) ); + + for(auto &face_arr: faces){ + for(auto &v_i : face_arr){ + if(N_A <= v_i) v_i = N_A + (N_B - 1) - (v_i - N_A); + } + } + return faces; + } + } + + /*************************************************************************** + * Graph theory search + **************************************************************************/ + + std::vector> nodes; + + + // Keep track of where we are looping over points + auto iter_A = std::begin(contour_A.points); + + // const convenience iterators + const auto begin_A = std::begin(contour_A.points); + const auto begin_B = std::begin(contour_B.points); + const auto end_A = std::end(contour_A.points); + const auto end_B = std::end(contour_B.points); + + auto p_i_next = *iter_A; + iter_A = std::next(iter_A); + + size_t m = *begin_A == contour_A.points.back() ? N_A - 1 : N_A; + size_t n = *begin_B == contour_B.points.back() ? N_A - 1 : N_A; + + // We repeat over A twice, since to search the graph we need to search repeats + // Don't loop over all points since last point is repeat of first + for(size_t i = 0; i < 2 * m; ++i, ++iter_A){ + + auto p_i = p_i_next; + // If at end of iterator loop back to beginning + if(std::next(iter_A) == end_A) iter_A = begin_A; + p_i_next = *iter_A; + + nodes.push_back(std::vector()); + + auto iter_B = std::begin(contour_B.points); + auto p_j_next = *iter_B; + iter_B = std::next(iter_B); + for(size_t j = 0; j < n; ++j, ++iter_B){ + + auto p_j = p_j_next; + // If at end of iterator loop back to beginning + p_j_next = iter_B == end_B ? *begin_B : *iter_B; + + nodes[i].push_back(node(i % m, j, p_i, p_j, p_i_next, p_j_next)); + } + nodes[i].push_back(node(i % m, 0, p_i, *begin_B, p_i_next, *std::next(begin_B))); + } + + // Useful to check progress of tiling if its going too slow + printf("Finding Optimal path, m=%zu, n=%zu\n", m, n); + + // Execute shortest path graph search + std::vector optimal_path = all_paths(nodes); + + + // Even for debugging, don't use this for dense graphs as it will print a lot to stdout + if (m < 15 && n < 15) + print_path(nodes, optimal_path); + + std::vector> ret; + + // Reconstruct paths + for(size_t j = 0; j <= n; ++j){ + size_t end_i = (j == n ? optimal_path[0] + m : optimal_path[j+1]); + for(size_t i = optimal_path[j]; i <= end_i; ++i){ + // If end then the next node is to the right, otherwise down + // On last column we can't go right any more so skip + size_t next_index; + + if(i < end_i) next_index = nodes[i+1][j].iA; + else if(j == n) continue; + else next_index = nodes[i][j+1].iB + N_A; + + ret.push_back({nodes[i][j].iA, nodes[i][j].iB + N_A, next_index}); + } + } + + return ret; + +} + +/* #pragma GCC pop_options */ // Low-level routine that joins the vertices of two contours. // Returns a list of faces where the vertex indices refer to A followed by B. @@ -808,4 +1199,4 @@ Minimally_Amalgamate_Contours_N_to_N( } } } -}*/ \ No newline at end of file +}*/ diff --git a/src/Simple_Meshing.h b/src/Simple_Meshing.h index e5b858b3..a55c769e 100644 --- a/src/Simple_Meshing.h +++ b/src/Simple_Meshing.h @@ -20,6 +20,10 @@ #include #include +std::vector< std::array > Tile_Contours( + std::reference_wrapper> A, + std::reference_wrapper> B ); + // Low-level routine that joins the vertices of two contours. // Returns a list of faces where the vertex indices refer to A followed by B. std::vector< std::array > From f4062cfd410732df8119dca4cd8be04147998c72 Mon Sep 17 00:00:00 2001 From: Renu Rajamagesh <60906270+Renu-R@users.noreply.github.com> Date: Mon, 10 Apr 2023 17:35:02 -0700 Subject: [PATCH 17/20] replaced tiling methods (#17) --- src/Operations/ConvertContoursToMeshes.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Operations/ConvertContoursToMeshes.cc b/src/Operations/ConvertContoursToMeshes.cc index d46a164c..b2a1cf64 100644 --- a/src/Operations/ConvertContoursToMeshes.cc +++ b/src/Operations/ConvertContoursToMeshes.cc @@ -617,7 +617,7 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, auto [faces, points, amal_upper] = Mesh_With_Convex_Hull_2(pcs.upper, m_cp_it->N_0, ofst_upper); add_faces_and_vertices(faces, points); auto amal_lower = pcs.lower.begin()->get(); - auto new_faces = Estimate_Contour_Correspondence(std::ref(amal_upper), std::ref(amal_lower)); + auto new_faces = Tile_Contours(std::ref(amal_upper), std::ref(amal_lower)); add_faces_to_mesh(std::ref(amal_upper), std::ref(amal_lower), new_faces); } catch (const std::runtime_error& error) { goto generic_n_to_n_meshing; @@ -647,7 +647,7 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, auto [faces, points, amal_lower] = Mesh_With_Convex_Hull_2(pcs.lower, l_cp_it->N_0, ofst_lower); add_faces_and_vertices(faces, points); auto amal_upper = pcs.upper.begin()->get(); - auto new_faces = Estimate_Contour_Correspondence(std::ref(amal_upper), std::ref(amal_lower)); + auto new_faces = Tile_Contours(std::ref(amal_upper), std::ref(amal_lower)); add_faces_to_mesh(std::ref(amal_upper), std::ref(amal_lower), new_faces); } catch (const std::runtime_error& error) { goto generic_n_to_n_meshing; @@ -717,7 +717,7 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, OverwriteStringToFile(amal_cop_str, fname); } */ - auto new_faces = Estimate_Contour_Correspondence(std::ref(amal_upper), std::ref(amal_lower)); + auto new_faces = Tile_Contours(std::ref(amal_upper), std::ref(amal_lower)); add_faces_to_mesh(std::ref(amal_upper), std::ref(amal_lower), new_faces); } } From 0dc5276f9a3ded8a0854dc9a56e2a71fe941dbdc Mon Sep 17 00:00:00 2001 From: Xander Naumenko <33139843+misprit7@users.noreply.github.com> Date: Tue, 11 Apr 2023 07:44:58 -0700 Subject: [PATCH 18/20] Fixed limites for graph construction (#18) --- src/Simple_Meshing.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Simple_Meshing.cc b/src/Simple_Meshing.cc index ee1f6245..ed50ec35 100644 --- a/src/Simple_Meshing.cc +++ b/src/Simple_Meshing.cc @@ -364,7 +364,7 @@ std::vector< std::array > Tile_Contours( iter_A = std::next(iter_A); size_t m = *begin_A == contour_A.points.back() ? N_A - 1 : N_A; - size_t n = *begin_B == contour_B.points.back() ? N_A - 1 : N_A; + size_t n = *begin_B == contour_B.points.back() ? N_B - 1 : N_B; // We repeat over A twice, since to search the graph we need to search repeats // Don't loop over all points since last point is repeat of first From cb58a4ae941b4116584e208cb541e5c7eb0b2fd2 Mon Sep 17 00:00:00 2001 From: tiffanymatthe <49164666+tiffanymatthe@users.noreply.github.com> Date: Sun, 23 Apr 2023 19:48:00 -0700 Subject: [PATCH 19/20] clean code (#19) --- scripts/compare_all_meshes.sh | 4 +-- src/Operations/CompareMeshes.cc | 32 ++++++++--------------- src/Operations/ConvertContoursToMeshes.cc | 10 ++++--- 3 files changed, 18 insertions(+), 28 deletions(-) diff --git a/scripts/compare_all_meshes.sh b/scripts/compare_all_meshes.sh index 8e63ddd6..a08bae9d 100755 --- a/scripts/compare_all_meshes.sh +++ b/scripts/compare_all_meshes.sh @@ -61,7 +61,7 @@ do shape=${shapes[$i]} res=${resolutions[$i]} shape_label=${shape_labels[$i]} - # echo "Computing differences for $shape" + OUTPUT="$(dicomautomaton_dispatcher \ -v \ -o GenerateMeshes \ @@ -126,8 +126,6 @@ do printf "| %-20s | %-13.3f | %-13.3f | %-15.3f | %-15.3f | %-10.3f | %-10.3f | %-10.3f | %-10.3f | %-8s | %-8s |\n\n" $shape_label $HD1 $HD2 $SA1 $SA2 $SA_DIFF $V1 $V2 $V_DIFF $MA1 $MA2 - - # echo "Finished computing differences for $shape" done diff --git a/src/Operations/CompareMeshes.cc b/src/Operations/CompareMeshes.cc index 04aa3004..e655d921 100644 --- a/src/Operations/CompareMeshes.cc +++ b/src/Operations/CompareMeshes.cc @@ -92,11 +92,6 @@ bool IsEdgeManifold(const std::shared_ptr &mesh) { std::map, int> edge_counts; auto mesh_faces = GetCleanFaces(mesh); int face_count = 0; - // YLOGINFO("Printing vertices for mesh"); - // for (auto &face : mesh.get()->meshes.faces) { - // YLOGINFO(face_count << " face vertex indices " << face[0] << ", " << face[1] << ", " << face[2] << " with size " << face.size()); - // face_count++; - // } face_count = 0; for (auto &face : mesh_faces) { @@ -105,11 +100,7 @@ bool IsEdgeManifold(const std::shared_ptr &mesh) { for (auto &edge : edges) { edge_counts[edge] += 1; - // auto itt = edge.begin(); - // YLOGINFO("Edge count = " << edge_counts[edge] << " for edge " << *itt << ", " << *(std::next(itt))); if (edge_counts[edge] > 2) { - // auto it = edge.begin(); - // YLOGINFO("edge count > 2 for edge " << *it << ", " << *(std::next(it)) << " for face #" << face_count); return false; } } @@ -118,7 +109,6 @@ bool IsEdgeManifold(const std::shared_ptr &mesh) { for (auto &key_value_pair : edge_counts) { if (key_value_pair.second != 2) { - // YLOGINFO("edge count is " << key_value_pair.second); return false; } } @@ -130,10 +120,6 @@ bool IsEdgeManifold(const std::shared_ptr &mesh) { // it is vertex manifold when each vertex's faces form an open or closed fan // https://www.mathworks.com/help/lidar/ref/surfacemesh.isvertexmanifold.html bool IsVertexManifold(const std::shared_ptr&mesh) { - // for each face, find faces with same edges and add to search - // keep searching until you hit a face you look for before (closed) or ended - // is there still faces unsearched? - // store vertex to face hashmap // for each vertex, store edge to list of faces // start with a face @@ -246,14 +232,18 @@ bool CompareMeshes(Drover &DICOM_data, } // Calculate centroids by finding the average of all the vertices in the surface mesh - //Assumes that contours vertices are distributed evenly on the surface mesh which is not strictly true. + // Assumes that contours vertices are distributed evenly on the surface mesh which is not strictly true. double sumx = 0, sumy = 0, sumz = 0; for (auto& vertex : mesh1->meshes.vertices) { sumx += vertex.x; sumy += vertex.y; sumz += vertex.z; } - vec3 centroid1 = vec3(sumx/size(mesh1->meshes.vertices) , sumy/size(mesh1->meshes.vertices) , sumz/size(mesh1->meshes.vertices)); + vec3 centroid1 = vec3( + sumx/size(mesh1->meshes.vertices), + sumy/size(mesh1->meshes.vertices), + sumz/size(mesh1->meshes.vertices) + ); sumx = 0, sumy= 0, sumz = 0; for (auto& vertex : mesh2->meshes.vertices) { @@ -265,9 +255,9 @@ bool CompareMeshes(Drover &DICOM_data, const double centroid_shift = sqrt(pow((centroid2.x-centroid1.x),2) + pow((centroid2.y-centroid1.y),2) + pow((centroid2.y-centroid1.y),2)); - //Converting to a triangular mesh to ensure that each face is made up of 3 vertices for volume calculation - //Total volume is calculated by summing the signed volme of the triganular prism made by each face and the origin as the apex. - //This method returns a finite volume even if the mesh is open, so watertightness needs to be checked separately. + // Converting to a triangular mesh to ensure that each face is made up of 3 vertices for volume calculation + // Total volume is calculated by summing the signed volme of the triganular prism made by each face and the origin as the apex. + // This method returns a finite volume even if the mesh is open, so watertightness needs to be checked separately. mesh1->meshes.convert_to_triangles(); mesh2->meshes.convert_to_triangles(); @@ -306,8 +296,8 @@ bool CompareMeshes(Drover &DICOM_data, bool v_manifold_2 = IsVertexManifold(mesh2); bool e_manifold_1 = IsEdgeManifold(mesh1); bool e_manifold_2 = IsEdgeManifold(mesh2); - FUNCINFO("Vertex manifoldness: " << v_manifold_1 << " and " << v_manifold_2); - FUNCINFO("Edge manifoldness: " << e_manifold_1 << " and " << e_manifold_2); + FUNCINFO("Vertex manifoldness (first vs. second): " << v_manifold_1 << " and " << v_manifold_2); + FUNCINFO("Edge manifoldness (first vs. second): " << e_manifold_1 << " and " << e_manifold_2); bool manifold1 = v_manifold_1 && e_manifold_1; bool manifold2 = v_manifold_2 && e_manifold_2; diff --git a/src/Operations/ConvertContoursToMeshes.cc b/src/Operations/ConvertContoursToMeshes.cc index b2a1cf64..4ec3e0ef 100644 --- a/src/Operations/ConvertContoursToMeshes.cc +++ b/src/Operations/ConvertContoursToMeshes.cc @@ -150,6 +150,7 @@ OperationDoc OpArgDocConvertContoursToMeshes(){ } + bool ConvertContoursToMeshes(Drover &DICOM_data, const OperationArgPkg& OptArgs, std::map& /*InvocationMetadata*/, @@ -471,6 +472,7 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, // Convert from integer numbering to direct pairing info. for(const auto &p : pairs){ + pairings.emplace_back(); for(const auto &u : p.upper){ pairings.back().upper.emplace_back( *std::next( std::begin(m_cops), u ) ); @@ -543,7 +545,6 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, std::vector> &points ) -> void { const auto old_face_count = amesh.vertices.size(); - // for(const auto &p : points) amesh.vertices.emplace_back(p); amesh.vertices.insert(amesh.vertices.end(), points.begin(), points.end()); for(const auto &fs : new_faces){ const auto f_A = static_cast(fs[0] + old_face_count); @@ -593,7 +594,7 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, /* auto new_faces = Estimate_Contour_Correspondence(pcs.upper.front(), pcs.lower.front()); */ auto new_faces = Tile_Contours(pcs.upper.front(), pcs.lower.front()); add_faces_to_mesh(pcs.upper.front(), pcs.lower.front(), new_faces); - }else if((N_upper == 2) && (N_lower == 1) ){ + }else if( (N_upper == 2) && (N_lower == 1) ){ //check if the upper plane contains enclosed contours if(contours_are_enclosed(*m_cp_it, pcs.upper.front(), pcs.upper.back())){ //get contour areas to determine which contour is enclosed by the other @@ -633,7 +634,7 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, const auto c2_area = std::abs(contour->get().Get_Signed_Area()); //cap the smaller contour and tile the larger contour with the lower plane contour - if (c1_area < c2_area){ + if(c1_area < c2_area){ close_hole_in_roof(pcs.lower.front()); auto new_faces = Tile_Contours(pcs.lower.back(), pcs.upper.front()); add_faces_to_mesh(pcs.lower.back(), pcs.upper.front(), new_faces); @@ -657,7 +658,7 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, //YLOGINFO("Performing N-to-N meshing.."); //routine for hollow structures with an inner contour and an outer contour on both planes - if((N_upper == 2) && (N_lower == 2)){ + if( (N_upper == 2) && (N_lower == 2) ){ //check if both planes have enclosed contours if (contours_are_enclosed(*m_cp_it, pcs.upper.front(), pcs.upper.back()) && contours_are_enclosed(*l_cp_it, pcs.lower.front(), pcs.lower.back())){ @@ -704,6 +705,7 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, generic_n_to_n_meshing: auto amal_upper = Minimally_Amalgamate_Contours(m_cp_it->N_0, ofst_upper, pcs.upper); auto amal_lower = Minimally_Amalgamate_Contours(m_cp_it->N_0, ofst_lower, pcs.lower); + /* // Leaving this here for future debugging, for which it will no-doubt be needed... { From 54c8b3fbcc047f386c7ed9596eebf6da592da6cd Mon Sep 17 00:00:00 2001 From: Xander Naumenko <33139843+misprit7@users.noreply.github.com> Date: Sun, 23 Apr 2023 20:12:52 -0700 Subject: [PATCH 20/20] Added downsampling (#20) --- src/Operations/ConvertContoursToMeshes.cc | 12 ++++++++++++ src/Simple_Meshing.cc | 4 ++-- 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/src/Operations/ConvertContoursToMeshes.cc b/src/Operations/ConvertContoursToMeshes.cc index 4ec3e0ef..643b86cc 100644 --- a/src/Operations/ConvertContoursToMeshes.cc +++ b/src/Operations/ConvertContoursToMeshes.cc @@ -146,6 +146,14 @@ OperationDoc OpArgDocConvertContoursToMeshes(){ "contours" }; out.args.back().samples = OpArgSamples::Exhaustive; + out.args.emplace_back(); + out.args.back().name = "Downsample"; + out.args.back().desc = "Whether to downsample before meshing. Tiling approach is quite computationally" + "expensive so often useful. Each contour will have at most this many points. If" + "0 then no downsampling."; + out.args.back().default_val = "0"; + out.args.back().expected = true; + return out; } @@ -163,6 +171,7 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, const auto ROILabelRegex = OptArgs.getValueStr("ROILabelRegex").value(); const auto MeshLabel = OptArgs.getValueStr("MeshLabel").value(); const auto MethodStr = OptArgs.getValueStr("Method").value(); + const auto Downsample = std::stoi(OptArgs.getValueStr("Downsample").value()); //----------------------------------------------------------------------------------------------------------------- const auto NormalizedMeshLabel = X(MeshLabel); @@ -191,6 +200,9 @@ bool ConvertContoursToMeshes(Drover &DICOM_data, // Contours. for(auto &c : cc_refw.get().contours){ if(c.points.empty()) continue; + if(Downsample != 0){ + c = c.Resample_LTE_Evenly_Along_Perimeter(Downsample); + } auto cop_base_ptr = reinterpret_cast *>(&c); cops.push_back( std::ref(*cop_base_ptr) ); } diff --git a/src/Simple_Meshing.cc b/src/Simple_Meshing.cc index ed50ec35..607a05e9 100644 --- a/src/Simple_Meshing.cc +++ b/src/Simple_Meshing.cc @@ -399,8 +399,8 @@ std::vector< std::array > Tile_Contours( // Even for debugging, don't use this for dense graphs as it will print a lot to stdout - if (m < 15 && n < 15) - print_path(nodes, optimal_path); + /* if (m < 15 && n < 15) */ + /* print_path(nodes, optimal_path); */ std::vector> ret;