Remove CGAL dependency, replace with native implementations by Copilot · Pull Request #52 · hdclark/DICOMautomaton

Copilot · 2026-01-31T05:23:31Z

Removes all CGAL library code from DICOMautomaton per request to eliminate this external dependency. Native implementations replace core functionality; operations requiring complex CGAL algorithms are removed from the build.

Native Implementations Added

2D Polygon Boolean Operations - Sutherland-Hodgman clipping for intersection; simplified union/difference suitable for typical medical imaging contours
Mesh Volume/Surface Area - Divergence theorem for volume, cross-product triangle sums for area
Mesh Slicing - Edge-plane intersection with contour tracing, handles vertices exactly on plane
Mesh Subdivision - Midpoint-based triangle refinement (simplified Loop-style)

Operations Removed from Build

BCCAExtractRadiomicFeatures, ExtractRadiomicFeatures, ConvertMeshesToContours, DumpROISurfaceMeshes, MakeMeshesManifold, MinkowskiSum3D, RemeshSurfaceMeshes, SubdivideSurfaceMeshes

These require CGAL's Nef polyhedra, Minkowski sums, or isotropic remeshing—infeasible to reimplement without external libraries.

Operations Updated

ContourViaThreshold - marching-cubes uses native Slice_Mesh instead of polyhedron_processing::Slice_Polyhedron
SimplifySurfaceMeshes - only native 'flat' algorithm available (CGAL edge-collapse removed)
ContourBooleanOperations/SeamContours - use native polygon Boolean implementation

Usage

// New native mesh utilities in dcma_surface_meshes namespace
double vol = dcma_surface_meshes::Mesh_Volume(mesh);
double area = dcma_surface_meshes::Mesh_Surface_Area(mesh);
auto contours = dcma_surface_meshes::Slice_Mesh(mesh, planes);
dcma_surface_meshes::Loop_Subdivide(mesh, iterations);

Limitations

Native polygon Boolean operations work well for simple, largely convex polygons typical in medical imaging. Complex non-convex polygon unions may produce approximate results.

✨ Let Copilot coding agent set things up for you — coding agent works faster and does higher quality work when set up for your repo.

Co-authored-by: hdclark <934858+hdclark@users.noreply.github.com>

Copilot

Pull request overview

This PR removes the CGAL library dependency from DICOMautomaton and replaces core functionality with native C++ implementations. The goal is to eliminate a complex external dependency while maintaining essential mesh processing capabilities for medical imaging workflows.

Changes:

Native implementations added for mesh volume/area calculation, mesh slicing, and 2D polygon Boolean operations
Eight operations requiring advanced CGAL features (Nef polyhedra, Minkowski sums, isotropic remeshing, radiomics) have been removed from the build
Remaining operations updated to use native implementations (ContourViaThreshold, SimplifySurfaceMeshes, ContourBooleanOperations, SeamContours)

Reviewed changes

Copilot reviewed 12 out of 12 changed files in this pull request and generated 14 comments.

Show a summary per file

File	Description
src/Surface_Meshes.h	Removes CGAL types; adds declarations for native Mesh_Volume, Mesh_Surface_Area, Slice_Mesh, and Loop_Subdivide functions
src/Surface_Meshes.cc	Implements native mesh operations using divergence theorem, cross products, edge-plane intersection, and midpoint subdivision
src/Contour_Boolean_Operations.cc	Replaces CGAL polygon set operations with Sutherland-Hodgman clipping and simplified union/difference algorithms
src/Operations/SimplifySurfaceMeshes.cc	Removes CGAL edge-collapse method, retains only native 'flat' simplification
src/Operations/ContourViaThreshold.cc	Updates marching-cubes workflow to use native Slice_Mesh instead of CGAL's Slice_Polyhedron
src/Operations/SeamContours.cc	Removes CGAL requirement check, relies on native Boolean operations
src/Operations/ContourBooleanOperations.cc	Removes CGAL requirement check
src/Operations/CMakeLists.txt	Unconditionally builds ContourBooleanOperations and SeamContours; removes CGAL-dependent operations from build
src/Operation_Dispatcher.cc	Removes includes and registrations for CGAL-dependent operations
src/CMakeLists.txt	Unconditionally builds Contour_Boolean_Operations, removes CGAL linking
CMakeLists.txt	Removes WITH_CGAL option and all CGAL find_package logic
src/Documentation.cc	Removes CGAL from dependency list in documentation

Copilot · 2026-01-31T08:09:56Z

src/Contour_Boolean_Operations.cc

+    // Compute convex hull of merged points as approximation
+    // Sort points by angle from centroid
+    Point2D centroid{0, 0};
+    for(const auto& p : all_points) {
+        centroid.x += p.x;
+        centroid.y += p.y;
+    }
+    centroid.x /= all_points.size();
+    centroid.y /= all_points.size();
+
+    std::sort(all_points.begin(), all_points.end(), 
+        [&centroid](const Point2D& a, const Point2D& b) {
+            return std::atan2(a.y - centroid.y, a.x - centroid.x) < 
+                   std::atan2(b.y - centroid.y, b.x - centroid.x);
+        });
+
+    // Remove duplicates
+    all_points.erase(std::unique(all_points.begin(), all_points.end()), all_points.end());
+
+    return {all_points};


The polygon_union function's approach of computing a convex hull (via angle sorting from centroid) when polygons overlap is incorrect for non-convex polygon unions. The convex hull will lose concave features of the input polygons. This could produce significantly different results compared to the original CGAL implementation for complex medical imaging contours with concave regions. Consider either implementing a proper polygon union algorithm (like Weiler-Atherton) or documenting this limitation prominently in both the code comments and the PR description's limitations section.

Copilot · 2026-01-31T08:09:56Z

src/Contour_Boolean_Operations.cc

+            }
        }
    }



The combine_polygons lambda at line 485 performs pairwise operations between all polygons in the A set and all polygons in the B set, which could produce an exponential number of intermediate results for large sets. For example, if A has 10 polygons and B has 10 polygons, you could end up with up to 100 result polygons. This differs from the CGAL approach which would properly merge overlapping regions. Consider whether the nested loop at line 485-504 should be restructured to combine results more efficiently.

Suggested change

if(result_polys.size() > 1) {

combine_polygons(result_polys, construction_op);

}

Copilot · 2026-01-31T08:09:56Z

src/Contour_Boolean_Operations.cc

+bool line_intersection(const Point2D& p1, const Point2D& p2,
+                       const Point2D& p3, const Point2D& p4,
+                       Point2D& intersection) {
+    double d = (p1.x - p2.x) * (p3.y - p4.y) - (p1.y - p2.y) * (p3.x - p4.x);
+    if(std::abs(d) < 1e-10) return false; // Lines are parallel
+
+    double t = ((p1.x - p3.x) * (p3.y - p4.y) - (p1.y - p3.y) * (p3.x - p4.x)) / d;
+    double u = -((p1.x - p2.x) * (p1.y - p3.y) - (p1.y - p2.y) * (p1.x - p3.x)) / d;
+
+    if(t >= 0.0 && t <= 1.0 && u >= 0.0 && u <= 1.0) {
+        intersection.x = p1.x + t * (p2.x - p1.x);
+        intersection.y = p1.y + t * (p2.y - p1.y);
+        return true;
+    }
+    return false;
+}


The line_intersection function checks if t and u are in the range [0.0, 1.0] to determine if segments intersect, but the epsilon value 1e-10 used for the parallel-line check may cause issues with near-parallel lines where d is very small but non-zero. This could lead to numerical instability when computing t and u. Consider using a larger epsilon or adding additional checks to ensure t and u are computed reliably.

Copilot · 2026-01-31T08:09:57Z

src/Operations/SimplifySurfaceMeshes.cc

-
        }else{
-            throw std::invalid_argument("Method argument '"_s + MethodStr + "' is not valid");
+            throw std::invalid_argument("Method argument '"_s + MethodStr + "' is not valid. Only 'flat' is supported.");


The error message now states "Only 'flat' is supported" but the parameter description and examples still suggest that other methods might be valid in the future. For consistency and clarity, consider updating the parameter description to explicitly state that only 'flat' is currently available after CGAL removal, to avoid confusion for users who might be familiar with the previous 'edge-collapse' option.

Copilot · 2026-01-31T08:09:57Z

src/Surface_Meshes.cc

+            contour.points.push_back(current_pt);
+
+            // Walk through connected faces
+            size_t max_iterations = mesh.faces.size() + 1;


The contour-building loop uses a fixed max_iterations limit (mesh.faces.size() + 1) which could cause premature termination for complex meshes where the contour path is longer than the number of faces. This could result in incomplete contours. Consider using a larger multiple of the face count, or detecting the loop condition by checking if we've returned to the starting face/edge rather than relying solely on an iteration limit.

Copilot · 2026-01-31T08:09:58Z

src/Contour_Boolean_Operations.cc

+        for(size_t j = 0; j < input.size(); ++j) {
+            const Point2D& current = input[j];
+            const Point2D& next = input[(j + 1) % input.size()];
+
+            bool current_inside = is_inside(current);
+            bool next_inside = is_inside(next);
+
+            if(current_inside) {
+                if(next_inside) {
+                    // Both inside: add next point
+                    output.push_back(next);


In the Sutherland-Hodgman clipping algorithm, when both the current and next points are inside the clipping edge, only the next point is added to the output (line 155). This means the first point of the input polygon is never added to the output polygon, which could result in the output polygon missing one vertex. The standard Sutherland-Hodgman algorithm typically handles this by adding the current point when transitioning from outside to inside, or by ensuring the loop correctly captures all vertices. Verify this implementation produces correct results for all cases.

Copilot · 2026-01-31T08:09:58Z

src/Surface_Meshes.cc

-        // which will stick out beyond the contour planes. However, the margin is added to both the top and the bottom so
-        // halve the total amount.
-        z_margin = sep_per_plane * 1.5;
+    return std::abs(volume) / 6.0;


The Mesh_Volume function returns the absolute value of the signed volume, which masks whether the mesh has inconsistent winding order. While the comment on line 1343 mentions this is intentional to "handle meshes with inconsistent winding order," this could hide actual mesh topology problems. Consider either documenting this behavior more prominently in the function's API documentation comment, or providing a separate function that returns the signed volume so callers can detect winding issues if needed.

Copilot · 2026-01-31T08:09:58Z

src/Surface_Meshes.cc

+        // Map from face index to its intersection points
+        std::map<uint64_t, std::vector<std::pair<std::pair<uint64_t, uint64_t>, vec3<double>>>> face_edge_intersections;
+
+        const double plane_eps = 1e-10; // Tolerance for determining if vertex is on the plane


The Slice_Mesh function uses a fixed tolerance (plane_eps = 1e-10) for determining if a vertex lies on a plane. This tolerance may be too strict for some meshes with larger coordinate values or too loose for meshes with fine precision. Consider making this tolerance configurable via a parameter with a sensible default, or calculating it adaptively based on the mesh scale (e.g., from bounding box dimensions).

Copilot · 2026-01-31T08:09:59Z

src/Surface_Meshes.cc

+
+        // First pass: create edge midpoints
+        for(const auto& face : mesh.faces){
+            if(face.size() != 3) continue;


The Loop_Subdivide function only handles triangular faces (face.size() != 3 is skipped). When non-triangular faces are present, they are silently dropped from the mesh during subdivision, which could result in holes or incomplete meshes. Consider either throwing an error if non-triangular faces are encountered, or triangulating them first, or documenting this limitation more clearly in the function's API documentation.

Copilot · 2026-01-31T08:09:59Z

src/Contour_Boolean_Operations.cc

+std::vector<Polygon2D> polygon_union(const Polygon2D& A, const Polygon2D& B) {
+    // Simple approach: if one contains the other, return the larger
+    // Otherwise, return both as separate polygons (conservative approach)
+
+    if(A.empty()) return B.empty() ? std::vector<Polygon2D>{} : std::vector<Polygon2D>{B};
+    if(B.empty()) return {A};
+
+    // Check if A contains all of B
+    bool A_contains_B = true;
+    for(const auto& p : B) {
+        if(!point_in_polygon(p, A)) {
+            A_contains_B = false;
+            break;
+        }
+    }
+    if(A_contains_B) return {A};
+
+    // Check if B contains all of A
+    bool B_contains_A = true;
+    for(const auto& p : A) {
+        if(!point_in_polygon(p, B)) {
+            B_contains_A = false;
+            break;
+        }
+    }
+    if(B_contains_A) return {B};
+
+    // Check if polygons intersect
+    auto intersection = sutherland_hodgman_clip(A, B);
+    if(intersection.empty()) {
+        // No intersection, return both
+        return {A, B};
+    }
+
+    // For overlapping polygons, create a simple merged boundary
+    // This is a simplified approach that works for many practical cases
+    // A full implementation would use Weiler-Atherton algorithm
+
+    // Collect all vertices from both polygons plus intersection points
+    std::vector<Point2D> all_points;
+
+    // Add vertices from A that are not inside B
+    for(const auto& p : A) {
+        if(!point_in_polygon(p, B)) {
+            all_points.push_back(p);
+        }
+    }
+
+    // Add vertices from B that are not inside A
+    for(const auto& p : B) {
+        if(!point_in_polygon(p, A)) {
+            all_points.push_back(p);
+        }
+    }
+
+    // Add intersection points
+    for(size_t i = 0; i < A.size(); ++i) {
+        for(size_t j = 0; j < B.size(); ++j) {
+            Point2D isect;
+            if(line_intersection(A[i], A[(i+1)%A.size()], B[j], B[(j+1)%B.size()], isect)) {
+                all_points.push_back(isect);
+            }
+        }
+    }
+
+    if(all_points.size() < 3) {
+        return {A, B}; // Fallback
+    }
+
+    // Compute convex hull of merged points as approximation
+    // Sort points by angle from centroid
+    Point2D centroid{0, 0};
+    for(const auto& p : all_points) {
+        centroid.x += p.x;
+        centroid.y += p.y;
+    }
+    centroid.x /= all_points.size();
+    centroid.y /= all_points.size();
+
+    std::sort(all_points.begin(), all_points.end(), 
+        [&centroid](const Point2D& a, const Point2D& b) {
+            return std::atan2(a.y - centroid.y, a.x - centroid.x) < 
+                   std::atan2(b.y - centroid.y, b.x - centroid.x);
+        });
+
+    // Remove duplicates
+    all_points.erase(std::unique(all_points.begin(), all_points.end()), all_points.end());
+
+    return {all_points};


The polygon_union function uses a fallback that returns both polygons as separate results when they overlap but neither fully contains the other. This means the "union" operation can return disconnected components, which may not be the expected behavior. The comment on line 226 acknowledges this is a "simplified approach," but callers may not realize they need to handle multiple result polygons. Consider documenting this behavior more clearly in the function comments or in the API documentation for ContourBoolean.

Initial plan

0ad7f15

Copilot AI assigned Copilot and hdclark Jan 31, 2026

Copilot started work on behalf of hdclark January 31, 2026 05:23 View session

Copilot AI and others added 4 commits January 31, 2026 05:35

Remove CGAL dependency from core components

66cb47d

Co-authored-by: hdclark <934858+hdclark@users.noreply.github.com>

Update remaining files to remove CGAL references

a5c5cd6

Co-authored-by: hdclark <934858+hdclark@users.noreply.github.com>

Address code review feedback

fca4672

Co-authored-by: hdclark <934858+hdclark@users.noreply.github.com>

Remove remaining CGAL references from CMakeLists.txt

d4d955e

Co-authored-by: hdclark <934858+hdclark@users.noreply.github.com>

Copilot AI changed the title ~~[WIP] Evaluate and replace CGAL library usage in DICOMautomaton~~ Remove CGAL dependency, replace with native implementations Jan 31, 2026

Copilot AI requested a review from hdclark January 31, 2026 05:44

Copilot finished work on behalf of hdclark January 31, 2026 05:44

hdclark marked this pull request as ready for review January 31, 2026 08:02

Copilot AI review requested due to automatic review settings January 31, 2026 08:02

Copilot started reviewing on behalf of hdclark January 31, 2026 08:03 View session

Copilot AI reviewed Jan 31, 2026

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Remove CGAL dependency, replace with native implementations#52

Remove CGAL dependency, replace with native implementations#52
Copilot wants to merge 5 commits intomasterfrom
copilot/remove-cgal-dependencies

Copilot AI commented Jan 31, 2026 •

edited

Loading

Uh oh!

Copilot AI left a comment

Uh oh!

Copilot AI Jan 31, 2026

Uh oh!

Copilot AI Jan 31, 2026

Uh oh!

Copilot AI Jan 31, 2026

Uh oh!

Copilot AI Jan 31, 2026

Uh oh!

Copilot AI Jan 31, 2026

Uh oh!

Copilot AI Jan 31, 2026

Uh oh!

Copilot AI Jan 31, 2026

Uh oh!

Copilot AI Jan 31, 2026

Uh oh!

Copilot AI Jan 31, 2026

Uh oh!

Copilot AI Jan 31, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

+    if(result_polys.size() > 1) {
+        combine_polygons(result_polys, construction_op);
+    }

Conversation

Copilot AI commented Jan 31, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Native Implementations Added

Operations Removed from Build

Operations Updated

Usage

Limitations

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Copilot AI Jan 31, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Jan 31, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Jan 31, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Jan 31, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Jan 31, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Jan 31, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Jan 31, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Jan 31, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Jan 31, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Jan 31, 2026

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Copilot AI commented Jan 31, 2026 •

edited

Loading