Refactor meshgen and CleanPSLG for Improved Breakline Processing and Mesh Quality

@geodata/private/shoelace.m

@@ -1,9 +1,34 @@
-function area = shoelace(x,y)
-n = length(x);
-xp = [x; x(1)];
-yp = [y; y(1)];
-area = 0;
-for i = 1:n
-    area = area + det([xp(i), xp(i+1); yp(i), yp(i+1)]);
+function area = shoelace(x, y)
+    % Computes the area of a polygon using the Shoelace formula.
+    % Handles NaNs by computing areas of individual sub-polygons separately.
+
+    % Ensure x and y are column vectors
+    x = x(:);
+    y = y(:);
+
+    % Identify NaN indices to split polygons
+    nanIndices = isnan(x) | isnan(y);
+
+    % Find segment start and end indices
+    segmentStart = find([true; nanIndices(1:end-1)] & ~nanIndices);
+    segmentEnd = find(~nanIndices & [nanIndices(2:end); true]);
+
+    % Initialize total area
+    area = 0;
+
+    % Process each sub-polygon separately
+    for i = 1:length(segmentStart)
+        xi = x(segmentStart(i):segmentEnd(i));
+        yi = y(segmentStart(i):segmentEnd(i));
+
+        % Ensure the polygon is closed (first and last points should match)
+        if ~isequal(xi(1), xi(end)) || ~isequal(yi(1), yi(end))
+            xi = [xi; xi(1)];
+            yi = [yi; yi(1)];
+        end
+
+        % Apply Shoelace formula
+        Ai = 0.5 * abs(sum(xi(1:end-1) .* yi(2:end) - xi(2:end) .* yi(1:end-1)));
+        area = area + Ai; % Accumulate area from all sub-polygons
+    end
 end
-area = 1/2*abs(area);

@meshgen/private/Get_line_edges.m

@@ -1,17 +1,22 @@
 function edges = Get_line_edges(nodes)
-%GET_LINE_EDGES Generate edges from a sequence of nodes in a closed loop.
+%GET_LINE_EDGES Generate edges from a sequence of nodes.
 %   edges = GET_LINE_EDGES(nodes) takes a list of nodes (vertices) as input
-%   and generates a list of edges where each node is connected to the next,
-%   forming a closed loop. The output 'edges' is an Mx2 matrix, where M is
-%   the number of edges, and each row represents an edge defined by the
-%   indices of its two endpoints.
+%   and generates a list of edges where each node is connected to the next.
+%   If the nodes form a closed loop, the last node is connected back to the first.
+%   The output 'edges' is an Mx2 matrix, where M is the number of edges,
+%   and each row represents an edge defined by the indices of its two endpoints.
+%
+%   Assumes 'nodes' is an Nx2 array (or Nx3 for 3D), where each row is a coordinate.
 
     % Number of nodes
-    nNodes = length(nodes);
+    nNodes = size(nodes, 1);
 
     % Generate edges connecting each node to the next
     edges = [(1:nNodes-1)', (2:nNodes)'];
 
-    % Add the edge connecting the last node back to the first
-    edges(end+1, :) = [nNodes, 1];
+    % Check if the shape is a closed loop
+    if all(nodes(1, :) == nodes(end, :))
+        % Add the edge connecting the last node back to the first
+        edges(end+1, :) = [nNodes, 1];
+    end
 end

@meshgen/private/filter_polygon_constraints.m

@@ -1,46 +1,47 @@
 function [pfix, egfix] = filter_polygon_constraints(pfix, egfix, ibouboxes, box_number)
-%FILTER_CONSTRAINTS Removes edge constraints based on bounding box criteria.
-%   This function filters out edge constraints (egfix) where at least one
-%   endpoint (from pfix) is outside the specified bounding box (ibouboxes)
-%   for the given box_number. It also adjusts the edges based on nested
-%   bounding boxes, if applicable.
-% remove bars if one point is outside
-node1=pfix(egfix(:,1),:);
-node2=pfix(egfix(:,2),:);
-iboubox = ibouboxes{box_number};
-% to enlarge or shrink the box, you must make sure bbox is equi
-% spaced
-[ty,tx]=my_interpm(iboubox(:,2),iboubox(:,1),100/111e3);
-iboubox = [tx,ty]; 
-buffer_size = 1.0;
-iboubox(:,1) = buffer_size*iboubox(:,1)+(1-buffer_size)*mean(iboubox(1:end-1,1));
-iboubox(:,2) = buffer_size*iboubox(:,2)+(1-buffer_size)*mean(iboubox(1:end-1,2));
-inside_node1 = inpoly(node1,iboubox(1:end-1,:)) ;
-inside_node2 = inpoly(node2,iboubox(1:end-1,:)) ;
-inside = inside_node1 .* inside_node2;
-% Get all points inside inner boxes and consider these outside for
-% all the nested boxes.
-for bn = box_number+1:length(ibouboxes)
-    % enlarge nest
-    iboubox = ibouboxes{bn}(1:end-1,:);
-    [ty,tx]=my_interpm(iboubox(:,2),iboubox(:,1),100/111e3);
-    iboubox = [tx,ty];
-    iboubox(:,1) = 1.25*iboubox(:,1)+(1-1.25)*mean(iboubox(1:end-1,1));
-    iboubox(:,2) = 1.25*iboubox(:,2)+(1-1.25)*mean(iboubox(1:end-1,2));
-    inside_node1 = inpoly(node1,iboubox);
-    inside_node2 = inpoly(node2,iboubox);
-    inside2 = inside_node1 .* inside_node2;
-    inside(find(inside2)) = false;
-end
-egfix(~inside,:) = [];
-tegfix=egfix';
-uid = unique(tegfix(:));
-tuid = length(uid);
-% remove nfix outside iboubox
-if tuid > 0
-    % remove pfix if outside domain
-    pfix = pfix(uid,:);
-end
-egfix = renumberEdges(egfix);
-
-end
+%FILTER_POLYGON_CONSTRAINTS Removes edges and points outside a buffered bounding polygon.
+
+    % Suppress polyshape inconsistency warning
+    warning('off', 'MATLAB:polyshape:repairedBySimplify');
+
+    % Extract bounding polygon and apply buffer
+    iboubox = ibouboxes{box_number};
+    buffered_poly = polybuffer(polyshape(iboubox), 100 / 111000); % ~100m buffer in degrees
+    [tx, ty] = boundary(buffered_poly);
+    iboubox = [tx, ty];
+
+    % Extract node positions from edges
+    node1 = pfix(egfix(:,1), :);
+    node2 = pfix(egfix(:,2), :);
+
+    % Identify edges where at least one endpoint is inside
+    inside = inpolygon(node1(:,1), node1(:,2), iboubox(:,1), iboubox(:,2)) | ...
+             inpolygon(node2(:,1), node2(:,2), iboubox(:,1), iboubox(:,2));
+
+    % Process nested bounding boxes (if applicable)
+    for bn = box_number+1:length(ibouboxes)
+        nested_iboubox = ibouboxes{bn};
+        nested_buffered_poly = polybuffer(polyshape(nested_iboubox), 125 / 111000); % ~125m buffer
+        [tx, ty] = boundary(nested_buffered_poly);
+        nested_iboubox = [tx, ty];
+
+        % Identify edges where both endpoints are inside the nested box (remove these)
+        inside_nested = inpolygon(node1(:,1), node1(:,2), nested_iboubox(:,1), nested_iboubox(:,2)) & ...
+                        inpolygon(node2(:,1), node2(:,2), nested_iboubox(:,1), nested_iboubox(:,2));
+        inside(inside_nested) = false;
+    end
+
+    % Retain only edges that pass filtering
+    egfix = egfix(inside, :);
+
+    % Retain only necessary points
+    if ~isempty(egfix)
+        uid = unique(egfix(:)); % Unique point indices used in edges
+        pfix = pfix(uid, :); % Retain only those points
+        egfix = renumberEdges(egfix); % Renumber edges to reflect reduced point set
+    end
+
+    % Restore warnings
+    warning('on', 'MATLAB:polyshape:repairedBySimplify');
+
+end

utilities/drawedge2.m

@@ -19,6 +19,3 @@ function drawedge2(pp,e2,color)
 
 
 end
-
-
-