diff --git a/@CleanPSLG/CleanPSLG.m b/@CleanPSLG/CleanPSLG.m new file mode 100644 index 00000000..6760fa9b --- /dev/null +++ b/@CleanPSLG/CleanPSLG.m @@ -0,0 +1,628 @@ +classdef CleanPSLG + properties + Vertices + Segments + Tolerance = 0.001; + AngleThreshold = 10; + DistanceThreshold = 0.001; % Distance threshold for nearby edges + PreservedSegments = []; % Logical array marking segments to preserve + end + + methods + function obj = CleanPSLG(vertices, segments, tol, angle_thresh, dist_thresh) + obj.Vertices = vertices; + obj.Segments = segments; + if nargin > 2, obj.Tolerance = tol; end + if nargin > 3, obj.AngleThreshold = angle_thresh; end + if nargin > 4, obj.DistanceThreshold = dist_thresh; + else obj.DistanceThreshold = obj.Tolerance; end + + % Initialize preserved segments array + obj.PreservedSegments = false(size(segments, 1), 1); + end + + function obj = mergeVertices(obj) + [obj.Vertices, obj.Segments] = obj.fix_geo(obj.Vertices, obj.Segments, obj.Tolerance); + end + + function obj = dropIntersectingEdges(obj) + obj.Segments = obj.drop_intersecting_edges(obj.Vertices, obj.Segments); + end + + function obj = pruneEncroachingEdges(obj) + obj.Segments = obj.drop_encroaching_edges(obj.Vertices, obj.Segments, obj.AngleThreshold, obj.Tolerance); + end + + function obj = dropNearbyEdges(obj) + obj.Segments = obj.drop_nearby_edges(obj.Vertices, obj.Segments, obj.DistanceThreshold, obj.PreservedSegments); + end + + function obj = preserveLongChains(obj, minChainLength) + % Preserve contiguous chains of segments longer than minChainLength + % These chains will be protected during dropNearbyEdges + + % Identify chains of connected segments + chains = obj.identifyContiguousChains(); + + % Calculate the total length of each chain and mark for preservation + numChains = length(chains); + chainLengths = zeros(numChains, 1); + preservedChains = false(numChains, 1); + + for i = 1:numChains + chainLengths(i) = obj.calculateChainLength(chains{i}); + % Convert length to kilometers (assuming coordinates are in degrees) + chainLengthKm = chainLengths(i) * 111; % Approximate conversion from degrees to km + if chainLengthKm >= minChainLength + preservedChains(i) = true; + end + end + + % Mark segments in long chains as preserved + obj.PreservedSegments = false(size(obj.Segments, 1), 1); + for i = 1:numChains + if preservedChains(i) + obj.PreservedSegments(chains{i}) = true; + end + end + + fprintf('Identified %d chains, preserving %d chains longer than %.1f km\n', ... + numChains, sum(preservedChains), minChainLength); + fprintf('Preserving %d segments out of %d total segments\n', ... + sum(obj.PreservedSegments), length(obj.PreservedSegments)); + end + + function obj = deleteShortChains(obj, maxChainLength) + % Delete chains of segments shorter than maxChainLength (in kilometers) + % Identify chains of connected segments + chains = obj.identifyContiguousChains(); + + % Calculate the total length of each chain and identify short ones to delete + numChains = length(chains); + chainLengths = zeros(numChains, 1); + shortChains = false(numChains, 1); + segmentsToDelete = []; + + for i = 1:numChains + chainSegs = chains{i}; + chainLengths(i) = obj.calculateChainLength(chainSegs); + % Convert length to kilometers (assuming coordinates are in degrees) + chainLengthKm = chainLengths(i) * 111; % Approximate conversion from degrees to km + + % Mark short chains for deletion + if chainLengthKm < maxChainLength + shortChains(i) = true; + segmentsToDelete = [segmentsToDelete; chainSegs]; + end + end + + % Remove the short chain segments + if ~isempty(segmentsToDelete) + keepMask = true(size(obj.Segments, 1), 1); + keepMask(segmentsToDelete) = false; + obj.Segments = obj.Segments(keepMask, :); + + % Also update the preserved segments array if it exists + if ~isempty(obj.PreservedSegments) + obj.PreservedSegments = obj.PreservedSegments(keepMask); + end + end + + fprintf('Identified %d chains, deleted %d chains shorter than %.1f km\n', ... + numChains, sum(shortChains), maxChainLength); + fprintf('Deleted %d segments out of %d original segments (%.1f%%)\n', ... + length(segmentsToDelete), length(segmentsToDelete) + size(obj.Segments, 1), ... + 100 * length(segmentsToDelete) / (length(segmentsToDelete) + size(obj.Segments, 1))); + end + + function chains = identifyContiguousChains(obj) + % Identify all contiguous chains of segments + segments = obj.Segments; + vertices = obj.Vertices; + + % Create vertex adjacency lists + n = size(vertices, 1); + vertexToSegments = cell(n, 1); + + for i = 1:size(segments, 1) + v1 = segments(i, 1); + v2 = segments(i, 2); + vertexToSegments{v1} = [vertexToSegments{v1}, i]; + vertexToSegments{v2} = [vertexToSegments{v2}, i]; + end + + % Find connected components using depth-first search + visited = false(size(segments, 1), 1); + chains = {}; + + for i = 1:size(segments, 1) + if ~visited(i) + % Start a new chain + chain = []; + stack = i; + + % DFS to find all connected segments + while ~isempty(stack) + current = stack(end); + stack(end) = []; + + if ~visited(current) + visited(current) = true; + chain = [chain; current]; + + % Find connected segments through shared vertices + v1 = segments(current, 1); + v2 = segments(current, 2); + + % Add connected segments to stack + for connectedSegIdx = [vertexToSegments{v1}, vertexToSegments{v2}] + if ~visited(connectedSegIdx) + stack = [stack, connectedSegIdx]; + end + end + end + end + + chains{end+1} = chain; + end + end + end + + function length = calculateChainLength(obj, chainIndices) + % Calculate the total length of a chain of segments + length = 0; + for i = 1:numel(chainIndices) + segIdx = chainIndices(i); + v1 = obj.Vertices(obj.Segments(segIdx, 1), :); + v2 = obj.Vertices(obj.Segments(segIdx, 2), :); + length = length + norm(v2 - v1); + end + end + + function plotPSLG(obj) + figure; + subplot(1,2,1); + scatter(obj.Vertices(:,1), obj.Vertices(:,2), 10, 'b', 'filled'); hold on; + for i = 1:size(obj.Segments,1) + pts = obj.Vertices(obj.Segments(i,:),:); + plot(pts(:,1), pts(:,2), 'k-'); + end + title('Processed PSLG'); xlabel('Longitude'); ylabel('Latitude'); axis equal; + end + + function obj = cleanUnusedVertices(obj) + % Remove vertices not referenced by any segments and renumber segments accordingly + % Find all vertex indices that are used in segments + usedIndices = unique(obj.Segments(:)); + numOriginalVertices = size(obj.Vertices, 1); + + % Create a mapping from old indices to new indices + maxIndex = max(max(obj.Segments)); + vertexMap = zeros(maxIndex, 1); + vertexMap(usedIndices) = 1:length(usedIndices); + + % Update the segments with new vertex indices + obj.Segments = vertexMap(obj.Segments); + + % Keep only used vertices + obj.Vertices = obj.Vertices(usedIndices, :); + + fprintf('Removed %d unused vertices (%.1f%% reduction)\n', ... + numOriginalVertices - length(usedIndices), ... + 100 * (1 - length(usedIndices)/numOriginalVertices)); + end + + end + + methods (Static) + function [newVertices, newSegments] = fix_geo(vertices, segments, tol) + % Merge close vertices + n = size(vertices,1); + parent = (1:n)'; + function r = find_parent(i) + while parent(i) ~= i, parent(i) = parent(parent(i)); i = parent(i); end + r = i; + end + function union(i, j) + ri = find_parent(i); + rj = find_parent(j); + if ri ~= rj, parent(rj) = ri; end + end + for i = 1:n-1 + for j = i+1:n + if norm(vertices(i,:) - vertices(j,:)) < tol + union(i,j); + end + end + end + mapping = arrayfun(@find_parent, (1:n)'); + [uniqueMappings, ~, newMapping] = unique(mapping); + newVertices = vertices(uniqueMappings, :); + newSegments = arrayfun(@(x) newMapping(x), segments); + newSegments = unique(sort(newSegments, 2), 'rows'); + end + + function newSegments = drop_intersecting_edges(vertices, segments) + % Drop intersecting edges + nseg = size(segments,1); + drop = false(nseg,1); + for i = 1:nseg + for j = i+1:nseg + if any(segments(i,:) == segments(j,:)) + continue; + end + p = vertices(segments(i,1),:); + r = vertices(segments(i,2),:); + q = vertices(segments(j,1),:); + s = vertices(segments(j,2),:); + if CleanPSLG.segments_intersect(p, r, q, s) + drop(i) = true; drop(j) = true; + end + end + end + newSegments = segments(~drop,:); + end + + function flag = segments_intersect(p, r, q, s) + o1 = CleanPSLG.orientation(p, r, q); + o2 = CleanPSLG.orientation(p, r, s); + o3 = CleanPSLG.orientation(q, s, p); + o4 = CleanPSLG.orientation(q, s, r); + flag = (o1 ~= o2) && (o3 ~= o4); + end + + function o = orientation(p, q, r) + val = (q(2)-p(2))*(r(1)-q(1)) - (q(1)-p(1))*(r(2)-q(2)); + o = (val > 0) - (val < 0); + end + + function newSegments = drop_encroaching_edges(vertices, segments, angle_thresh_deg, tol) + nseg = size(segments,1); + drop = false(nseg,1); + cos_thresh = cosd(angle_thresh_deg); + for i = 1:nseg + if drop(i), continue; end + for j = i+1:nseg + if drop(j), continue; end + if any(segments(i,:) == segments(j,:)) + continue; + end + v1 = vertices(segments(i,:),:); + v2 = vertices(segments(j,:),:); + unit1 = (v1(2,:) - v1(1,:)) / norm(v1(2,:) - v1(1,:)); + unit2 = (v2(2,:) - v2(1,:)) / norm(v2(2,:) - v2(1,:)); + if abs(dot(unit1, unit2)) > cos_thresh + drop(i) = true; + end + end + end + newSegments = segments(~drop,:); + end + + function newSegments = drop_nearby_edges(vertices, segments, dist_thresh, preservedSegments) + % Drop edges that are close to other edges - optimized version + % Only considers perpendicular distances between edges + % Respects segments marked for preservation + nseg = size(segments, 1); + if nseg < 2 + newSegments = segments; + return; + end + + % If preservedSegments not provided, assume no segments are preserved + if nargin < 4 || isempty(preservedSegments) + preservedSegments = false(nseg, 1); + end + + % Calculate segment midpoints and bounding boxes for spatial binning + midpoints = zeros(nseg, 2); + bbox = zeros(nseg, 4); % [minX, minY, maxX, maxY] + + for i = 1:nseg + v1 = vertices(segments(i,1), :); + v2 = vertices(segments(i,2), :); + midpoints(i, :) = (v1 + v2) / 2; + bbox(i, :) = [min(v1(1), v2(1)), min(v1(2), v2(2)), ... + max(v1(1), v2(1)), max(v1(2), v2(2))]; + end + + % Create spatial grid for efficient lookup + % Determine grid dimensions + domain = [min(bbox(:,1)), min(bbox(:,2)), max(bbox(:,3)), max(bbox(:,4))]; + gridSize = max(1, ceil(dist_thresh * 10)); % Adjust based on dist_thresh + + numCellsX = max(1, ceil((domain(3) - domain(1)) / gridSize)); + numCellsY = max(1, ceil((domain(4) - domain(2)) / gridSize)); + + % Assign segments to grid cells + cellAssignment = cell(numCellsX, numCellsY); + for i = 1:nseg + % Determine which cells this segment overlaps + minCellX = max(1, min(numCellsX, floor((bbox(i,1) - domain(1)) / gridSize) + 1)); + minCellY = max(1, min(numCellsY, floor((bbox(i,2) - domain(2)) / gridSize) + 1)); + maxCellX = max(1, min(numCellsX, ceil((bbox(i,3) - domain(1)) / gridSize))); + maxCellY = max(1, min(numCellsY, ceil((bbox(i,4) - domain(2)) / gridSize))); + + % Add segment to all overlapping cells + for cx = minCellX:maxCellX + for cy = minCellY:maxCellY + cellAssignment{cx, cy} = [cellAssignment{cx, cy}; i]; + end + end + end + + % Process segments + drop = false(nseg, 1); + + % Fast check for potential nearby segments using grid + for cx = 1:numCellsX + for cy = 1:numCellsY + cellSegments = cellAssignment{cx, cy}; + numCellSegs = length(cellSegments); + + if numCellSegs < 2 + continue; % Skip cells with 0 or 1 segment + end + + % Check all pairs of segments in this cell + for ii = 1:numCellSegs-1 + i = cellSegments(ii); + if drop(i) || preservedSegments(i) + continue; % Skip already dropped or preserved segments + end + + for jj = ii+1:numCellSegs + j = cellSegments(jj); + if drop(j) || preservedSegments(j) || any(segments(i,:) == segments(j,:)) + continue; % Skip already dropped, preserved, or connected segments + end + + % Quick bounding box check + if bbox(i,1) > bbox(j,3) + dist_thresh || ... + bbox(j,1) > bbox(i,3) + dist_thresh || ... + bbox(i,2) > bbox(j,4) + dist_thresh || ... + bbox(j,2) > bbox(i,4) + dist_thresh + continue; % Bounding boxes too far apart + end + + % Detailed perpendicular distance check + v1 = vertices(segments(i,:), :); + v2 = vertices(segments(j,:), :); + + % Get perpendicular distances + [perp_dist, is_perp] = CleanPSLG.perpendicular_segment_distance(v1, v2); + + % Only consider perpendicular distances + if is_perp && perp_dist < dist_thresh + % Drop the longer segment (unless preserved) + len1 = norm(v1(2,:) - v1(1,:)); + len2 = norm(v2(2,:) - v2(1,:)); + if len1 > len2 + drop(i) = true; + break; % Break inner loop, segment i is dropped + else + drop(j) = true; + end + end + end + end + end + end + + newSegments = segments(~drop,:); + end + + function [dist, is_perpendicular] = perpendicular_segment_distance(v1, v2) + % Calculate perpendicular distance between two line segments + % Returns the minimum perpendicular distance and whether it actually occurs + % v1 and v2 are 2x2 matrices containing start and end points of segments + + % Get segment vectors + p1 = v1(1,:); + p2 = v1(2,:); + p3 = v2(1,:); + p4 = v2(2,:); + + % Direction vectors of the segments + dir1 = p2 - p1; + dir2 = p4 - p3; + + % Initialize values + dist = Inf; + is_perpendicular = false; + + % Normalize directions + len1 = norm(dir1); + len2 = norm(dir2); + + if len1 < 1e-10 || len2 < 1e-10 + % One of the segments is almost a point + return; + end + + unit_dir1 = dir1 / len1; + unit_dir2 = dir2 / len2; + + % Check if segments are nearly parallel (not perpendicular) + dot_product = abs(dot(unit_dir1, unit_dir2)); + if dot_product > 0.9 % Approximately parallel if dot product > 0.9 (< 26 degrees difference) + % Segments are nearly parallel - not considering them perpendicular + return; + end + + % For the first segment, calculate perpendicular distance to the second segment + d1 = CleanPSLG.point_to_segment_perp_distance(p1, p3, p4); + d2 = CleanPSLG.point_to_segment_perp_distance(p2, p3, p4); + + % For the second segment, calculate perpendicular distance to the first segment + d3 = CleanPSLG.point_to_segment_perp_distance(p3, p1, p2); + d4 = CleanPSLG.point_to_segment_perp_distance(p4, p1, p2); + + % Get the minimum perpendicular distance + min_perp_dist = min([d1, d2, d3, d4]); + + % Only consider it a perpendicular case if the minimum distance is finite + if isfinite(min_perp_dist) + dist = min_perp_dist; + is_perpendicular = true; + end + end + + function dist = point_to_segment_perp_distance(p, s1, s2) + % Calculate perpendicular distance from point p to segment (s1,s2) + % Returns Inf if projection is outside the segment + + v = s2 - s1; % Segment vector + w = p - s1; % Vector from segment start to point + + % Check if we can project onto the segment + c1 = dot(w, v); + if c1 <= 0 + % p projects outside the segment, beyond s1 + dist = Inf; + return; + end + + c2 = dot(v, v); + if c2 <= c1 + % p projects outside the segment, beyond s2 + dist = Inf; + return; + end + + % Calculate projection parameter and perpendicular distance + b = c1 / c2; % Projection parameter (0 to 1) + pb = s1 + b * v; % Projected point on segment + dist = norm(p - pb); + end + + function d = fast_segment_distance(v1, v2) + % Optimized calculation of minimum distance between two line segments + + % First do quick endpoint checks + d_endpoints = min([ + norm(v1(1,:) - v2(1,:)), + norm(v1(1,:) - v2(2,:)), + norm(v1(2,:) - v2(1,:)), + norm(v1(2,:) - v2(2,:)) + ]); + + % Early exit if endpoints are close enough + if d_endpoints < 1e-10 + d = 0; + return; + end + + % Get line segment vectors and cross-segment vector + p1 = v1(1,:); + p2 = v1(2,:); + p3 = v2(1,:); + p4 = v2(2,:); + + v13 = p1 - p3; + v43 = p4 - p3; + v21 = p2 - p1; + + % Check if either segment is actually a point + d43 = sum(v43.^2); + d21 = sum(v21.^2); + + % Handle special cases + if d43 < 1e-10 && d21 < 1e-10 % Both segments are points + d = norm(p1 - p3); + return; + end + + if d21 < 1e-10 % First segment is a point + d = CleanPSLG.point_line_distance(p1, p3, p4); + return; + end + + if d43 < 1e-10 % Second segment is a point + d = CleanPSLG.point_line_distance(p3, p1, p2); + return; + end + + % Compute the parameters of the closest points + d1343 = sum(v13 .* v43); + d4321 = sum(v43 .* v21); + d1321 = sum(v13 .* v21); + d4343 = d43; + d2121 = d21; + + denom = d2121 * d4343 - d4321 * d4321; + + if abs(denom) < 1e-10 % Lines are parallel + % Use perpendicular distance between parallel lines + d_p1l2 = CleanPSLG.point_line_distance(p1, p3, p4); + d = d_p1l2; + return; + end + + numer = d1343 * d4321 - d1321 * d4343; + + mua = numer / denom; + mub = (d1343 + d4321 * mua) / d4343; + + % Clamp parameters to segment bounds + mua = max(0, min(1, mua)); + mub = max(0, min(1, mub)); + + % Compute the closest points on each segment + pa = p1 + mua * v21; + pb = p3 + mub * v43; + + % Return the distance between these points + d = norm(pa - pb); + end + + function d = segment_distance(v1, v2) + % Calculate minimum distance between two line segments + % v1 and v2 are 2x2 matrices containing start and end points of segments + + % Vector directions of the segments + d1 = v1(2,:) - v1(1,:); + d2 = v2(2,:) - v2(1,:); + + % Point-to-point distances + d_p1p1 = norm(v1(1,:) - v2(1,:)); + d_p1p2 = norm(v1(1,:) - v2(2,:)); + d_p2p1 = norm(v1(2,:) - v2(1,:)); + d_p2p2 = norm(v1(2,:) - v2(2,:)); + + % Point-to-line distances + d_p1l2 = CleanPSLG.point_line_distance(v1(1,:), v2(1,:), v2(2,:)); + d_p2l2 = CleanPSLG.point_line_distance(v1(2,:), v2(1,:), v2(2,:)); + d_p3l1 = CleanPSLG.point_line_distance(v2(1,:), v1(1,:), v1(2,:)); + d_p4l1 = CleanPSLG.point_line_distance(v2(2,:), v1(1,:), v1(2,:)); + + % Minimum distance is the smallest of all distances + d = min([d_p1p1, d_p1p2, d_p2p1, d_p2p2, d_p1l2, d_p2l2, d_p3l1, d_p4l1]); + end + + function d = point_line_distance(p, l1, l2) + % Calculate perpendicular distance from point p to line through l1 and l2 + % Check if projection falls within line segment + v = l2 - l1; + w = p - l1; + + c1 = dot(w, v); + if c1 <= 0 + d = norm(p - l1); + return; + end + + c2 = dot(v, v); + if c2 <= c1 + d = norm(p - l2); + return; + end + + b = c1 / c2; + pb = l1 + b * v; + d = norm(p - pb); + end + + + end +end \ No newline at end of file diff --git a/@geodata/geodata.m b/@geodata/geodata.m index c20bb9d3..4b2d1bcc 100644 --- a/@geodata/geodata.m +++ b/@geodata/geodata.m @@ -41,6 +41,7 @@ innerb % height of inner mainlandb_type % type of mainland (lake or river or ocean) innerb_type % type of inner (lake or river or ocean) + linestrings % unclosed vector features (e.g., thalwegs) weirs % weir crestlines weirPfix % boundaries of weir weirEgfix % edges of weir @@ -62,9 +63,35 @@ high_fidelity % performs a 1D mesh generation step to form pfix and egfix prior to 2D meshing end - + + %% **✅ Static Method: Handles NaN-Separated Segments** + methods (Static) + function plotStartEnd(vector, startColor, endColor) + if isempty(vector) || size(vector,2) < 2 + return; % Skip if invalid + end + + % Find NaN delimiters + nanIndices = [find(isnan(vector(:,1))); size(vector,1) + 1]; + startIndices = [1; nanIndices(1:end-1) + 1]; + endIndices = nanIndices - 1; + + % Remove invalid indices (e.g., NaNs at start or end) + validMask = (startIndices <= endIndices) & (startIndices > 0) & (endIndices > 0); + startIndices = startIndices(validMask); + endIndices = endIndices(validMask); + + % Plot markers at start and end of each segment + if ~isempty(startIndices) && ~isempty(endIndices) + m_plot(vector(startIndices,1), vector(startIndices,2), 'go', 'MarkerFaceColor', startColor, 'MarkerSize', 6); + m_plot(vector(endIndices,1), vector(endIndices,2), 'rx', 'MarkerFaceColor', endColor, 'MarkerSize', 6); + end + end + end + + methods - + function obj = geodata(varargin) % Class constructor to parse NetCDF DEM data, NaN-delimited vector, % or shapefile that defines polygonal boundary of meshing @@ -79,7 +106,7 @@ %addOptional(p,'weirs',defval); %addOptional(p,'pslg',defval); %addOptional(p,'boubox',defval); - + % Check for m_map dir M_MAP_EXISTS=0 ; if exist('m_proj','file')==2 @@ -88,7 +115,7 @@ if M_MAP_EXISTS~=1 error('Where''s m_map? Please read the user guide') end - + % Check for utilties dir UTIL_DIR_EXISTS=0 ; if exist('inpoly.m','file') @@ -97,7 +124,7 @@ if UTIL_DIR_EXISTS~=1 error('Where''s the utilities directory? Please read the user guide') end - + % Check for dataset dir DATASET_DIR_EXISTS=0 ; if exist('datasets','dir')==7 @@ -106,10 +133,10 @@ if DATASET_DIR_EXISTS~=1 warning('We suggest you to place your files in a directory called "datasets". Please read the user guide') end - - + + p = inputParser; - + defval = 0; % placeholder value if arg is not passed. % add name/value pairs addOptional(p,'bbox',defval); @@ -125,7 +152,7 @@ addOptional(p,'shapefile_3d',defval); addOptional(p,'high_fidelity',defval); - + % parse the inputs parse(p,varargin{:}); % store the inputs as a struct @@ -209,7 +236,7 @@ obj.window = 5; end case('shapefile_3d') - obj.shapefile_3d = inp.(fields{i}) ; + obj.shapefile_3d = inp.(fields{i}) ; case('high_fidelity') obj.high_fidelity = inp.(fields{i}); case('weirs') @@ -271,7 +298,7 @@ if size(obj.bbox,1) == 1 && isempty(obj.pslg) obj = ParseDEM(obj,'bbox'); end - + if obj.bbox == 0 error('No bbox supplied. If you are using the pslg option then you need to supply a bbox.') elseif size(obj.bbox,1) == 2 @@ -287,9 +314,9 @@ obj.bbox = [min(obj.boubox(:,1)) max(obj.boubox(:,1)) min(obj.boubox(:,2)) max(obj.boubox(:,2))] ; end - + obj = ParseShoreline(obj) ; - + % kjr Add the weir faux islands to the inner geometry if ~isempty(obj.weirPfix) idx = [0; cumsum(weirLength)']+1 ; @@ -302,22 +329,22 @@ end obj.inner = [obj.inner ; NaN NaN ; tmp ] ; end - + % Ensure inpoly flip is correct obj = check_connectedness_inpoly(obj); - + disp(['Read in meshing boundary: ',obj.contourfile]); - + if obj.BACKUPdemfile~=0 obj = ParseDEM(obj,'BackUp') ; end - + obj = ParseDEM(obj) ; - - - + + + end - + function obj = ParseShoreline(obj) % Read in the geometric meshing boundary information from a % ESRI-shapefile @@ -325,10 +352,10 @@ if ~iscell(obj.contourfile) obj.contourfile = {obj.contourfile}; end - + polygon_struct = Read_shapefile( obj.contourfile, [], ... obj.bbox, obj.gridspace, obj.boubox, 0, obj.shapefile_3d); - + % Unpack data from function Read_Shapefile()s obj.outer = polygon_struct.outer; obj.mainland = polygon_struct.mainland; @@ -337,15 +364,17 @@ obj.innerb = polygon_struct.innerb; obj.mainlandb_type = polygon_struct.mainlandb_type; obj.innerb_type = polygon_struct.innerb_type; - + + obj.linestrings = polygon_struct.linestrings; + % Read in the geometric meshing boundary information from a % NaN-delimited vector. elseif ~isempty(obj.pslg) - + % Handle the case for user defined mesh boundary information polygon_struct = Read_shapefile( [], obj.pslg, ... obj.bbox, obj.gridspace, obj.boubox, 0, obj.shapefile_3d); - + % Unpack data from function Read_Shapefile()s obj.outer = polygon_struct.outer; obj.mainland = polygon_struct.mainland; @@ -354,14 +383,14 @@ obj.innerb = polygon_struct.innerb; obj.mainlandb_type = polygon_struct.mainlandb_type; obj.innerb_type = polygon_struct.innerb_type; - + else % set outer to the boubox obj.outer = obj.boubox; end obj = ClassifyShoreline(obj) ; end - + function obj = ClassifyShoreline(obj) % Helper function to... % 1) Read the data from supplied file @@ -369,34 +398,34 @@ % 3) Ensure point spacing is adequate to support mesh sizes % and ensure computational efficiency when calculating the % signed distance function. - + % Check if no mainland segments, set outer % to the boubox. if isempty(obj.mainland) obj.outer = [ ]; obj.outer = obj.boubox; end - + % Make sure the shoreline components have spacing of % gridspace/spacing [la,lo] = my_interpm(obj.outer(:,2),obj.outer(:,1),... obj.gridspace/obj.spacing); obj.outer = []; obj.outer(:,1) = lo; obj.outer(:,2) = la; outerbox = obj.outer(1:find(isnan(obj.outer(:,1)),1,'first'),:); - + if ~isempty(obj.mainland) [la,lo] = my_interpm(obj.mainland(:,2),obj.mainland(:,1),... obj.gridspace/obj.spacing); obj.mainland = []; obj.mainland(:,1) = lo; obj.mainland(:,2) = la; end - + if ~isempty(obj.inner) [la,lo]=my_interpm(obj.inner(:,2),obj.inner(:,1),... obj.gridspace/obj.spacing); obj.inner = []; obj.inner(:,1) = lo; obj.inner(:,2) = la; end clearvars lo la - + % Smooth the coastline (apply moving average filter). if obj.window > 1 disp(['Smoothing coastline with ' ... @@ -413,32 +442,32 @@ else disp('No smoothing of coastline enabled') end - + % KJR: Coarsen portions of outer, mainland % and inner outside bbox. iboubox = bbox_to_bou(obj.bbox); iboubox(:,1) = 1.10*iboubox(:,1)+(1-1.10)*mean(iboubox(1:end-1,1)); iboubox(:,2) = 1.10*iboubox(:,2)+(1-1.10)*mean(iboubox(1:end-1,2)); - + % Coarsen outer obj.outer = coarsen_polygon(obj.outer,iboubox); - + % Coarsen inner and move parts that overlap with bounding box % to mainland if ~isempty(obj.inner) obj.inner = coarsen_polygon(obj.inner,iboubox); end - + % Coarsen mainland and remove parts that overlap with bounding % box (note this doesn't change the polygon used for inpoly, % only changes the distance function used for edgefx) if ~isempty(obj.mainland) obj.mainland = coarsen_polygon(obj.mainland,iboubox); end - - + + end - + function obj = ParseDEM(obj,varargin) % obj = ParseDEM(obj,varargin) % set 'BackUp' for reading backupdem @@ -448,23 +477,23 @@ backup = 1; fname = obj.BACKUPdemfile ; end - + % Process the DEM for the meshing region. if ~isempty(fname) - + % Find name of x, y (one that has 1 dimensions) % z values (use one that has 2 dimensions) [xvn, yvn, zvn] = getdemvarnames(fname); - + % Read x and y x = double(ncread(fname,xvn)); y = double(ncread(fname,yvn)); - + if any(strcmp(varargin,'bbox')) obj.bbox = [min(x) max(x); min(y) max(y)]; return; end - + modbox = [0 0]; if obj.bbox(1,2) > 180 && obj.bbox(1,1) < 180 && min(x) < 0 % bbox straddles 180/-180 line @@ -494,7 +523,7 @@ end It = find(x >= bboxt(1,1) & x <= bboxt(1,2)); I = [I; It]; - + % At this point, detect the memory footprint of the % DEM subset and compute the required stride necessary % to satisfy the memory requirements @@ -517,7 +546,7 @@ 'of ' num2str(STRIDE) ' to match h0']) end end - + end % grab only the portion that was requested if STRIDE > 1 && STRIDE_MEM > STRIDE_H0 @@ -557,7 +586,7 @@ % Determine bottom left corner of DEM % (after possible flipping of DEM packing) obj.x0y0 = [x(1),y(1)]; - + % check for any invalid values bad = isnan(demz); if sum(bad(:)) > 0 && ~backup @@ -571,7 +600,7 @@ demz(bad) = NaN ; end end - + % creating back up interpolant if backup obj.Fb2 = griddedInterpolant({x,y},demz,... @@ -590,15 +619,15 @@ % clear data from memory clear x y demz end - + disp(['Read in demfile ',fname]); end - + % Handle the case of no dem if isempty(obj.x0y0) obj.x0y0 = [obj.bbox(1,1), obj.bbox(2,1)]; end - + function [xvn, yvn, zvn] = getdemvarnames(fname) % Define well-known variables for longitude and latitude % coordinates in Digital Elevation Model NetCDF file (CF @@ -634,16 +663,16 @@ error('Could not locate z coordinate in DEM') ; end end - + end - - - + + + function obj = check_connectedness_inpoly(obj) % Check for connected polygons and if not connected, % whether to flip the inpoly result to ensure meshing % on the coastal side of the polygon. - + obj.inpoly_flip = 0; % return if outer polygon is connected shpEnd = find(isnan(obj.outer(:,1))); @@ -653,11 +682,11 @@ sum(obj.outer(shpEnd(loc+1)-1,:))) > eps disp('Warning: Shapefile is unconnected... continuing anyway') end - + % check for inpoly goodness read the GSHHS checker ps = Read_shapefile( {'GSHHS_l_L1'}, [], ... obj.bbox, obj.gridspace, obj.boubox, 0, 0 ); - + % make a "fake" tester grid x = linspace(obj.bbox(1,1),obj.bbox(1,2),100); y = linspace(obj.bbox(2,1),obj.bbox(2,2),100); @@ -673,20 +702,20 @@ disp(['Shapefile inpoly is inconsistent ' ... 'with GHSSS test file, flipping the inpoly test']) end - + % if flooplaind meshing, flip the inpoly test if obj.fp obj.inpoly_flip = mod(1,obj.inpoly_flip); end end - + function obj = close(obj) % Clips the mainland segment with the boubox. % Performs a breadth-first search given a seed position % of the piecewise-straight line graph (PSLG) that is used to define the meshing boundary. % This returns back an updated geodata class instance with the outer boundary clipped with the boubox. % kjr,und,chl 2018 - + if isempty(obj.Fb) warning('ALERT: Meshing boundary is not a polygon!') warning('ALERT: DEM is required to clip line segment with polygon') @@ -709,27 +738,27 @@ cands2=cands(in,:) ; seed = cands2(50,:) ; end - + geom = [obj.mainland; obj.boubox] ; - + [NODE,PSLG]=getnan2(geom) ; - + [NODE2,PSLG2,PART2] = bfsgeo2(NODE,PSLG,seed) ; - + POLY = extdom_polygon(PSLG2(PART2{1},:),NODE2,-1) ; - + new_outer = cell2mat(POLY') ; - + [la,lo] = my_interpm(new_outer(:,2),new_outer(:,1),((obj.h0/2)/111e3)) ; - + new_outer = [lo la] ; - + obj.outer = new_outer ; - + % reset this to default obj.inpoly_flip = 0 ; end - + function obj = extractContour(obj,ilev) % Extract a geometric contour from the DEM at elevation ilev. % obj = extractContour(obj,ilev) @@ -737,33 +766,32 @@ % Can use to get the mean sea level contour, e.g.; % gdat = geodata('pslg',0,'h0',min_el,'dem',dem); % make the dummy gdat for the dem extents; % lmsl = extractContour(gdat,0); %using the dummy gdat with dem info to get the 'lmsl' gdat with the 0-m contour. - + [node,edge] = ... getiso2( obj.Fb.GridVectors{1},obj.Fb.GridVectors{2},... double(obj.Fb.Values'),ilev) ; - + polyline = cell2mat(extdom_polygon(edge,node,-1,1,10)') ; - + obj = geodata('pslg',polyline,'bbox',obj.bbox,... 'h0',obj.h0,'dem',obj.demfile) ; - + end - - function plot(obj,type,projection,holdon) - % plot(obj,type,projection,holdon) - % Plot geodata class info + + function plot(obj, type, projection, holdon) + % plot(obj, type, projection, holdon) + % Plot geodata class info with start/end markers % % Inputs: % obj : geodata class object [required input] % - % optional inputs =>.... - % type : i) 'shp' [default] - plots the shapelines (shoreline) only - % ii) 'dem' - plots the dem bathy in addition to shapelines - % iii) 'omega' - hatches the meshing domain in addition to plotting shapelines - % projection : choose the projection type from m_map options - % (default is Mercator) - % holdon : plot on a new (= 0 [default]) or existing (= 1) figure? - + % Optional inputs: + % type : 'shp' [default] - plots the shapelines (shoreline) only, + % 'dem' - plots the DEM bathy in addition to shapelines, + % 'omega' - hatches the meshing domain in addition to plotting shapelines. + % projection : choose the projection type from m_map options (default is Mercator) + % holdon : plot on a new (0 [default]) or existing (1) figure? + if nargin == 1 || isempty(type) type = 'shp'; end @@ -773,41 +801,37 @@ function plot(obj,type,projection,holdon) if nargin < 4 || isempty(holdon) holdon = 0; end - - % plotting on new or existing figure? + + % Set up projection and figure if ~holdon - % setup the projection bufx = 0.2*(obj.bbox(1,2) - obj.bbox(1,1)); bufy = 0.2*(obj.bbox(2,2) - obj.bbox(2,1)); - if ~isempty(regexp(projection,'ste')) - m_proj(projection,'lat',min(obj.bbox(2,:)),... - 'long',mean(obj.bbox(1,:)),'radius',... - min(179.9,1.20*max(diff(obj.bbox(2,:))))); + if ~isempty(regexp(projection, 'ste')) + m_proj(projection, 'lat', min(obj.bbox(2,:)), 'long', mean(obj.bbox(1,:)), 'radius', ... + min(179.9, 1.20*max(diff(obj.bbox(2,:))))); else lmin = -180; lmax = +180; - if obj.bbox(1,2) > 180; lmax = 360; lmin = 0; end - lon1 = max(lmin,obj.bbox(1,1) - bufx); - lon2 = min(lmax,obj.bbox(1,2) + bufx); - lat1 = max(- 90,obj.bbox(2,1) - bufy); - lat2 = min(+ 90,obj.bbox(2,2) + bufy); - m_proj(projection,... - 'long',[lon1, lon2],'lat',[lat1, lat2]); + if obj.bbox(1,2) > 180 + lmax = 360; lmin = 0; + end + lon1 = max(lmin, obj.bbox(1,1) - bufx); + lon2 = min(lmax, obj.bbox(1,2) + bufx); + lat1 = max(-90, obj.bbox(2,1) - bufy); + lat2 = min(+90, obj.bbox(2,2) + bufy); + m_proj(projection, 'long', [lon1, lon2], 'lat', [lat1, lat2]); end - % plot on new figure figure; - colori = 'g-'; % set island color - colorm = 'r-'; % set mainland color + colori = 'g-'; % island color + colorm = 'r-'; % mainland color else - colori = 'b-'; % set island color - colorm = 'm-'; % set mainland color + colori = 'b-'; % island color for holdon + colorm = 'm-'; % mainland color for holdon end hold on - % select optional types + + % Select additional plotting options switch type - case('dem') - % interpolate DEM's bathy linearly onto our - % edgefunction grid (or a coarsened version of it for - % memory considerations) + case 'dem' mem = inf; stride = obj.h0/111e3; while mem > 1 xx = obj.x0y0(1):stride:obj.bbox(1,2); @@ -816,55 +840,44 @@ function plot(obj,type,projection,holdon) mem = xs.bytes*ys.bytes/1e9; stride = stride*2; end - [demx,demy] = ndgrid(xx,yy); - demz = obj.Fb(demx,demy); - m_fastscatter(demx(:),demy(:),demz(:)); - cb = colorbar; ylabel(cb,'topo-bathy depth [m]') - case('omega') - % hatch the meshing domain, Omega - [demx,demy] = ndgrid(obj.x0y0(1):obj.h0/111e3:obj.bbox(1,2), ... + [demx, demy] = ndgrid(xx, yy); + demz = obj.Fb(demx, demy); + m_fastscatter(demx(:), demy(:), demz(:)); + cb = colorbar; ylabel(cb, 'topo-bathy depth [m]'); + case 'omega' + [demx, demy] = ndgrid(obj.x0y0(1):obj.h0/111e3:obj.bbox(1,2), ... obj.x0y0(2):obj.h0/111e3:obj.bbox(2,2)); - edges = Get_poly_edges( [obj.outer; obj.inner] ); - in = inpoly([demx(:),demy(:)],[obj.outer; obj.inner], edges); + edges = Get_poly_edges([obj.outer; obj.inner]); + in = inpoly([demx(:), demy(:)], [obj.outer; obj.inner], edges); long = demx(~in); lati = demy(~in); - m_hatch(obj.boubox(1:end-1,1),... - obj.boubox(1:end-1,2),'cross',45,0.05); - m_plot(long,lati,'.','Color','white') + m_hatch(obj.boubox(1:end-1,1), obj.boubox(1:end-1,2), 'cross', 45, 0.05); + m_plot(long, lati, '.', 'Color', 'white'); end + + % Plot mainland, inner boundaries, and (if available) line strings + h1 = []; h2 = []; h3 = []; h4 = []; if ~isempty(obj.mainland) && obj.mainland(1) ~= 0 - h1 = m_plot(obj.mainland(:,1),obj.mainland(:,2),... - colorm,'linewi',1); hold on; + h1 = m_plot(obj.mainland(:,1), obj.mainland(:,2), colorm, 'linewi', 1); + geodata.plotStartEnd(obj.mainland, 'g', 'r'); % Mark start (green) and end (red) end if ~isempty(obj.inner) && obj.inner(1) ~= 0 - h2 = m_plot(obj.inner(:,1),obj.inner(:,2),... - colori,'linewi',1); hold on; - end - if ~isempty(obj.weirs) - if isstruct(obj.weirs) ~= 0 - for ii =1 : length(obj.weirs) - h3 = m_plot(obj.weirs.X,obj.weirs.Y,... - 'm-','linewi',1); hold on; - end - else - for ii =1 : length(obj.weirs) - h3 = m_plot(obj.weirs{ii}(:,1),obj.weirs{ii}(:,2),... - 'm-','linewi',1); hold on; - end - end + h2 = m_plot(obj.inner(:,1), obj.inner(:,2), colori, 'linewi', 1); + geodata.plotStartEnd(obj.inner, 'g', 'r'); end - [la,lo] = my_interpm(obj.boubox(:,2),obj.boubox(:,1),... - 0.5*obj.h0/111e3); - m_plot(lo,la,'k--','linewi',2); - if ~holdon - m_grid('xtick',10,'tickdir','out','yaxislocation','left','fontsize',10); + if ~isempty(obj.linestrings) + h4 = m_plot(obj.linestrings(:,1), obj.linestrings(:,2), 'c--', 'linewi', 1); + geodata.plotStartEnd(obj.linestrings, 'g', 'r'); % Mark start (green) and end (red) end - if exist('h1','var') && exist('h2','var') && exist('h3','var') - legend([h1 h2,h3],{'mainland' 'inner' 'weirs'},'Location','NorthWest') - elseif exist('h1','var') && exist('h2','var') - legend([h1 h2],{'mainland' 'inner'},'Location','NorthWest') + + % Plot the outer boundary hatch (dashed line) + [la, lo] = my_interpm(obj.boubox(:,2), obj.boubox(:,1), 0.5*obj.h0/111e3); + m_plot(lo, la, 'k--', 'linewi', 2); + + if ~holdon + m_grid('xtick', 10, 'tickdir', 'out', 'yaxislocation', 'left', 'fontsize', 10); end end - + end - + end diff --git a/@geodata/private/Read_shapefile.m b/@geodata/private/Read_shapefile.m index 4102c6f6..7ba12a12 100644 --- a/@geodata/private/Read_shapefile.m +++ b/@geodata/private/Read_shapefile.m @@ -1,297 +1,137 @@ -function polygon_struct = Read_shapefile( finputname, polygon, bbox, ... - h0, boubox, plot_on, shapefile_3d) -% Read_shapefile: Reads a shapefile or a NaN-delimited vector -% containing polygons and/or segments in the the desired region -% of interest. Classifies the vector data as either a -% mainland, island, or outer boundary. The program will automatically -% trim islands that are have an area smaller than 4*h0^2 and will in gaps -% in the vector that are larger than h0/2. This is necessary for the -% use of boundary rejection method in dpoly. +function polygon_struct = Read_shapefile(finputname, ~, bbox, h0, boubox, plot_on, ~) +%% ============================ +% Initialize Data Structures +% ============================ +SG = []; +loop = 1; minus = 0; +tolerance = 1e-5; +min_area_inner = 4 * h0^2; +min_area_mainland = 100 * h0^2; -% INPUTS: -% finputname : file name(s) of the shapefile listed in a cell array -% polygon : a NaN-delimited vector of polygons and/or segments. -% bbox : the bounding box that we want to extract out -% h0 : minimum edge length in region of interest. -% plot_on : plot the final polygon or not (=1 or =0) -% -% OUTPUTS: -% polygon_struct : a structure containing the vector features identified as -% either islands, mainland, or outer. -% Written by William Pringle and Keith Roberts, CHL,UND, 2017 -% Edits by Keith Roberts, July 2018. -%% Loop over all the filenames and get the shapefile within bbox -SG = []; -if bbox(1,2) > 180 && bbox(1,1) < 180 - % bbox straddles 180/-180 meridian - loop = 2; minus = 0; -elseif all(bbox(1,:) > 180) - % beyond 180 in 0 to 360 format - loop = 1; minus = 1; -else - loop = 1; minus = 0; -end -if (size(finputname,1)~=0) - for fname = finputname - for nn = 1:loop - bboxt = bbox'; - if loop == 2 - if nn == 1 - bboxt(2,1) = 180; - else - bboxt(1,1) = -180; bboxt(2,1) = bboxt(2,1) - 360; - end - end - if minus - bboxt(:,1) = bboxt(:,1) - 360; - end - % Read the structure - try - % The shaperead is faster if it is available - S = shaperead(fname{1},'BoundingBox',bboxt); - % Get rid of unwanted components; - D = struct2cell(S); - S = cell2struct(D(3:4,:)',{'X','Y'},2); - disp('Read shapefile with shaperead') - sr = 1; - catch - % If only m_shaperead is available or if some error occured - % with shaperead (e.g., 3D shapefile, m_shaperead may work) - disp('Reading shapefile with m_shaperead') - % This uses m_map (slower but free) - S = m_shaperead(fname{1},reshape(bboxt',4,1)); - % Let's just keep the x-y data - D = S.ncst; - if isfield(S,'dbf') || isfield(S,'dbfdata') - code = S.dbfdata(:,1); - S = cell2struct([D code]',{'points' 'type'},1); - else - S = cell2struct(D','points',1); - end - sr = 0; - end - if ~isempty(S) - % Keep the following polygons - SG = [SG; S]; - end - end - end -else - sr = 1 ; - % convert NaN-delimited vector to struct - count = 1; j=1; - for i = 1 : length(polygon) - % the end of the segment - - if(isnan(polygon(i,1))==1) - % put NaN at end - SG(count,1).X(:,j) =NaN; - SG(count,1).Y(:,j) =NaN; - % reset - j = 1 ; count = count + 1; - - continue - else - % keep going - SG(count,1).X(:,j) = polygon(i,1); - SG(count,1).Y(:,j) = polygon(i,2); - j=j+1; +if bbox(1,2) > 180 && bbox(1,1) < 180, loop = 2; end +if all(bbox(1,:) > 180), minus = 1; end + +%% ============================ +% Read Shapefile Data +% ============================ +for fname = finputname + for nn = 1:loop + bboxt = bbox'; + if loop == 2, bboxt(2-((nn-1)*1),1) = 180 - (nn-1)*360; end + if minus, bboxt(:,1) = bboxt(:,1) - 360; end + + try + S = shaperead(fname{1}, 'BoundingBox', bboxt); + S = rmfield(S, setdiff(fieldnames(S), {'X', 'Y'})); + catch + S = m_shaperead(fname{1}, reshape(bboxt', 4, 1)); + S = cell2struct(S.ncst', 'points', 1); end + SG = [SG; S]; end end -% If we don't have an outer polygon already then make it by bbox -polygon_struct.outer = boubox; -% Densify the outer polygon (fills gaps larger than half min edgelength). -[latout,lonout] = my_interpm(polygon_struct.outer(:,2),... - polygon_struct.outer(:,1),h0/2); -polygon_struct.outer = []; -polygon_struct.outer(:,1) = lonout; -polygon_struct.outer(:,2) = latout; -%% Find whether the polygons are wholey inside the bbox.. -%% Set as islands or mainland -disp('Partitioning the boundary into islands, mainland, ocean') -polygon_struct.inner = []; -polygon_struct.mainland = []; -polygon_struct.innerb = []; -polygon_struct.mainlandb = []; -polygon_struct.innerb_type = []; -polygon_struct.mainlandb_type = []; -edges = Get_poly_edges( polygon_struct.outer ); +%% ============================ +% Handle Empty Data +% ============================ +if isempty(SG) + polygon_struct = struct('outer', boubox, 'inner', [], 'mainland', [], ... + 'linestrings', [], 'mainlandb', [], 'innerb', [], ... + 'mainlandb_type', [], 'innerb_type', []); + return; +end -if isempty(SG); return; end +polygon_struct = struct('outer', boubox, 'inner', [], 'mainland', [], ... + 'linestrings', [], 'mainlandb', [], 'innerb', [], ... + 'mainlandb_type', [], 'innerb_type', []); +edges = Get_poly_edges(polygon_struct.outer); -if sr - tmpM = [[SG.X]',[SG.Y]'] ; % MAT - if bbox(1,2) > 180 - tmpM(tmpM(:,1) < 0,1) = tmpM(tmpM(:,1) < 0,1) + 360; - end - for i = 1 : length(SG) - dims(i) = length(SG(i).X) ; - end - tmpC = mat2cell(tmpM,dims); % TO CELL -else - tmpC = struct2cell(SG)'; - - % Remove empty cells from shapefile - if any(cellfun(@isempty,tmpC(:,1))) - warning('Empty objects in shapefile have been removed') - tmpC = tmpC(~cellfun(@isempty,tmpC(:,1)),:); - end +%% ============================ +% Extract & Classify Shapes +% ============================ +tmpC = struct2cell(SG)'; +tmpC = tmpC(~cellfun(@isempty, tmpC(:,1)),:); +tmpC = cellfun(@(row) row(:,1:2), tmpC, 'UniformOutput', false); - tmpCC = []; nn = 0; - for ii = 1:size(tmpC,1) - % may have NaNs inside - isnan1 = find(isnan(tmpC{ii,1}(:,1))); - if isempty(isnan1) - isnan1 = length(tmpC{ii,1})+1; - elseif isnan1(end) ~= length(tmpC{ii,1}) - isnan1(end+1) = length(tmpC{ii,1})+1; - end - is = 1; - for jj = 1:length(isnan1) - nn = nn + 1; - ie = isnan1(jj)-1; - tmpCC{nn,1} = tmpC{ii,1}(is:ie,:); - tmpCC{nn,2} = tmpC{ii,2}; - is = isnan1(jj)+1; - end - end - if ~isempty(tmpCC) - tmpC = tmpCC; - end - tmpM = cell2mat(tmpC(:,1)); - if size(tmpM,2) == 3 - tmpM = tmpM(:,1:2); - end -end -% Get current polygon -% Check proportion of polygon that is within bbox -tmpIn = inpoly(tmpM,polygon_struct.outer, edges); -tmpInC = mat2cell(tmpIn,cellfun(@length,tmpC(:,1))); +valid_polygons = {}; +line_strings = {}; -j = 0 ; k = 0 ; height = []; new_islandb = []; new_mainb = []; -new_islandb_type = []; new_mainb_type = []; -for i = 1 : size(tmpC,1) - if sr - % using shaperead - points = tmpC{i,1}(1:end-1,:) ; - In = tmpInC{i,1}(1:end-1) ; - else - % using m_shaperead - points = tmpC{i,1}(1:end,:) ; - if size(points,2) == 3 - points = points(:,1:2); - end - if shapefile_3d - % if 3-D shapefile - height = points(:,3); - points = points(:,1:2); - type = tmpC{i,2}; - if strcmp(type,'BA040') - type = 'ocean'; - elseif strcmp(type,'BH080') - type = 'lake'; - elseif strcmp(type,'BH140') - type = 'river'; - end - else - height = []; - end - In = tmpInC{i,1}(1:end) ; - end - if bbox(1,2) > 180 - lond = abs(diff(points(:,1))); - if any(lond > 350) - points(points(:,1) > 180,1) = 0; - end - end - % lets calculate the area of the - % feature using the shoelace algorithm and decided whether to keep or - % not based on area. - if all(points(1,:) == points(end,:)) - area = shoelace(points(:,1),points(:,2)); - else - area = 999; % not a polygon - end - if length(find(In == 1)) == length(points) - % Wholey inside box - if area < 4*h0^2 % too small, then don't consider it. - continue; - end - % Set as island (with NaN delimiter) - k = k + 1 ; - new_island{k} = [points; NaN NaN]; - if ~isempty(height) - new_islandb{k} = [points height; NaN NaN NaN]; - new_islandb_type{k} = type; - end +for i = 1:size(tmpC,1) + points = tmpC{i,1}; + + % **Check for Polygon Closure Using Tolerance** + first_point = points(1, :); + distances = sqrt(sum((points(2:end, :) - first_point).^2, 2)); + + if any(distances < tolerance) + valid_polygons{end+1} = points; else - %Partially inside box - if area < 100*h0^2 % too small, then don't consider it. - continue; - end - % Set as mainland - j = j + 1 ; - new_main{j} = [points; NaN NaN]; - if ~isempty(height) - new_mainb{j} = [points height; NaN NaN NaN]; - new_mainb_type{j} = type; - end + line_strings{end+1} = [points; NaN NaN]; end end -if k > 0 - polygon_struct.inner = [polygon_struct.inner; cell2mat(new_island')]; - polygon_struct.innerb = cell2mat(new_islandb'); - polygon_struct.innerb_type = new_islandb_type; + +if ~isempty(line_strings) + polygon_struct.linestrings = cell2mat(line_strings'); end -if j > 0 - polygon_struct.mainland = [polygon_struct.mainland; cell2mat(new_main')]; - polygon_struct.mainlandb = cell2mat(new_mainb'); - polygon_struct.mainlandb_type = new_mainb_type; + +%% ============================ +% Compute Polygon Areas and Assign to Inner or Mainland +% ============================ +inner_polys = {}; +mainland_polys = {}; +innerb_types = {}; +mainlandb_types = {}; + +% **Progress Bar for Large Datasets** +total_polygons = numel(valid_polygons); +if total_polygons > 1000 + f = waitbar(0, 'Processing polygons...'); end -% Merge overlapping mainland and inner -if exist('polybool','file') || exist('polyshape','file') - if ~isempty(new_mainb) && k > 0 - polym = polygon_struct.mainland; - mergei = false(k,1); - for kk = 1:k - polyi = new_island{kk}; - IA = find(ismembertol(polym,polyi,1e-5,'ByRows',true)); - if length(IA) > 2 - if exist('polyshape','file') - polyout = union(polyshape(polym),polyshape(polyi)); - polym = polyout.Vertices; - else - [x,y] = polybool('union',polym(:,1),polym(:,2),... - polyi(:,1),polyi(:,2)); - polym = [x,y]; - end - mergei(kk) = 1; - end - end - if ~isnan(polym(end,1)); polym(end+1,:) = NaN; end - polygon_struct.mainland = polym; - new_island(mergei) = []; - polygon_struct.inner = cell2mat(new_island'); + +for i = 1:total_polygons + points = valid_polygons{i}; + area = shoelace(points(:,1), points(:,2)); + inside_bbox = all(inpoly(points, polygon_struct.outer, edges)); + + if inside_bbox && abs(area) >= min_area_inner + inner_polys{end+1} = points; + innerb_types{end+1} = 'inner'; + elseif abs(area) >= min_area_mainland + mainland_polys{end+1} = points; + mainlandb_types{end+1} = 'mainland'; end -else - warning(['no polyshape or polybool available to merge possible ' ... - 'overlapping of mainland and inner']) -end -% Remove parts of inner and mainland overlapping with outer -polygon_struct.outer = [polygon_struct.outer; polygon_struct.mainland]; -%% Plot the map -if plot_on >= 1 && ~isempty(polygon_struct) - figure(1); - hold on - plot(polygon_struct.outer(:,1),polygon_struct.outer(:,2)) - if ~isempty(polygon_struct.inner) - plot(polygon_struct.inner(:,1),polygon_struct.inner(:,2)) + + % **Update Progress Bar** + if exist('f', 'var') && mod(i, 500) == 0 + waitbar(i / total_polygons, f); end - if ~isempty(polygon_struct.mainland) - plot(polygon_struct.mainland(:,1),polygon_struct.mainland(:,2)) +end + +if exist('f', 'var'), close(f); end % Close progress bar if it exists + +polygon_struct.inner = cell2mat(cellfun(@(x) [x; NaN NaN], inner_polys, 'UniformOutput', false)'); +polygon_struct.mainland = cell2mat(cellfun(@(x) [x; NaN NaN], mainland_polys, 'UniformOutput', false)'); + +%% ============================ +% Initialize `mainlandb` and `innerb` Correctly +% ============================ +polygon_struct.mainlandb = polygon_struct.mainland; +polygon_struct.innerb = polygon_struct.inner; +polygon_struct.mainlandb_type = mainlandb_types; +polygon_struct.innerb_type = innerb_types; + +%% ============================ +% Plot Results (Optional) +% ============================ +if plot_on >= 1 + figure(1); hold on; + plot(polygon_struct.outer(:,1), polygon_struct.outer(:,2), 'k'); + plot(polygon_struct.inner(:,1), polygon_struct.inner(:,2), 'b'); + plot(polygon_struct.mainland(:,1), polygon_struct.mainland(:,2), 'r'); + if ~isempty(polygon_struct.linestrings) + plot(polygon_struct.linestrings(:,1), polygon_struct.linestrings(:,2), 'c--'); end + axis equal end -%EOF + end diff --git a/@geodata/private/shoelace.m b/@geodata/private/shoelace.m index 8b1ae8d5..ba73bc39 100644 --- a/@geodata/private/shoelace.m +++ b/@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); \ No newline at end of file diff --git a/@meshgen/meshgen.m b/@meshgen/meshgen.m index dba9459b..4ce9428e 100644 --- a/@meshgen/meshgen.m +++ b/@meshgen/meshgen.m @@ -16,156 +16,51 @@ % % You should have received a copy of the GNU General Public License % along with this program. If not, see . - % - % Available options: - % ef % edgefx class - % bou % geodata class - % h0 % minimum edge length (optional if bou exists) - % bbox % bounding box [xmin,ymin; xmax,ymax] (manual specification, no bou) - % proj % structure containing the m_map projection info - % plot_on % flag to plot (def: 1) or not (0) - % nscreen % how many it to plot and write temp files (def: 5) - % itmax % maximum number of iterations. - % pfix % fixed node positions (nfix x 2 ) - % egfix % edge constraints - % outer % meshing boundary (manual specification, no bou) - % inner % island boundaries (manual specification, no bou) - % mainland % the shoreline boundary (manual specification, no bou) - % fixboxes % a flag that indicates which boxes will use fixed constraints - % memory_gb % memory in GB allowed to use for initial rejector - % cleanup % logical flag or string to trigger cleaning of topology (default is on). - % direc_smooth % logical flag to trigger direct smoothing of mesh in the cleanup - % dj_cutoff % the cutoff area fraction for disjoint portions to delete - % qual_tol % tolerance for the accepted negligible change in quality - % enforceWeirs % whether or not to enforce weirs in meshgen - % enforceMin % whether or not to enfore minimum edgelength for all edgefxs - % delaunay_elim_on_exit % whether or not to run delaunay_elim on exit of meshgen - % improve_with_reduced_quality % whether or not to allow mesh improvements with decreases in mesh quality - % improve_boundary % run a gradient desc. to improve the boundary conformity - % high_fidelity % flag to form pfix and egfix for this domain - % + properties - fd % handle to distance function - fh % handle to edge function - h0 % minimum edge length - edgefx % edgefx class - bbox % bounding box [xmin,ymin; xmax,ymax] - pfix % fixed node positions (nfix x 2 ) - egfix % edge constraints - fixboxes % a flag that indicates which boxes will use fixed constraints - plot_on % flag to plot (def: 1) or not (0) - nscreen % how many it to plot and write temp files (def: 5) - bou % geodata class - ef % edgefx class - itmax % maximum number of iterations. - outer % meshing boundary (manual specification, no bou) - inner % island boundaries (manual specification, no bou) - mainland % the shoreline boundary (manual specification, no bou) - boubox % the bbox as a polygon 2-tuple - inpoly_flip % used to flip the inpoly test to determine the signed distance. - memory_gb % memory in GB allowed to use for initial rejector - cleanup % logical flag or string to trigger cleaning of topology (default is on). - direc_smooth % logical flag to trigger direct smoothing of mesh in the cleanup - dj_cutoff % the cutoff area fraction for disjoint portions to delete - grd = msh(); % create empty mesh class to return p and t in. - qual % mean, lower 3rd sigma, and the minimum element quality. - qual_tol % tolerance for the accepted negligible change in quality - proj % structure containing the m_map projection info - anno % Approx. Nearest Neighbor search object. - annData % datat contained with KD-tree in anno - Fb % bathymetry data interpolant - enforceWeirs % whether or not to enforce weirs in meshgen - enforceMin % whether or not to enfore minimum edgelength for all edgefxs - improve_boundary % improve the boundary representation - high_fidelity % flag to form pfix and egfix for this domain - delaunay_elim_on_exit % whether or not to run delaunay_elim on exit of meshgen - improve_with_reduced_quality % whether or not to allow mesh improvements with decreases in mesh quality + % User-defined inputs and internal parameters + fd % Handle to distance function + fh % Handle to edge function(s) + h0 % Minimum edge length (scalar or array) + bbox % Bounding box [xmin,ymin; xmax,ymax] (or cell array) + pfix % Fixed node positions (nfix x 2) + egfix % Fixed edge constraints (indices into pfix) + fixboxes % Flag indicating which boxes use fixed constraints + plot_on % Plotting flag (1 = on, 0 = off) + nscreen % Frequency of plotting and writing temp files + bou % Geodata object(s) + ef % Edgefx class instance(s) + itmax % Maximum number of iterations for mesh improvement + outer % Outer boundary (cell array if multiple boxes) + inner % Island boundaries (cell array) + mainland % Shoreline boundary (cell array) + boubox % Bounding box as a polygon (cell array) + inpoly_flip % Flag to flip the inpoly test for signed distances + memory_gb % Memory in GB for initial rejection method + cleanup % Flag/string to trigger mesh cleanup (default 'default') + direc_smooth % Flag to trigger direct smoothing during cleanup + dj_cutoff % Cutoff area fraction for deleting disjoint portions + grd = msh(); % Mesh container (of type msh) for nodes p and triangles t + qual % Quality metrics for mesh elements + qual_tol % Tolerance for negligible mesh quality change + proj % Projection structure for m_map + anno % Approximate Nearest Neighbor search object(s) + annData % Data contained in the KD-tree(s) + Fb % Bathymetry data interpolant(s) + enforceWeirs % Flag to enforce weirs in mesh generation + enforceMin % Flag to enforce minimum edge length for all ef's + improve_boundary % Flag to improve the boundary representation + high_fidelity % Cell array or scalar flag for high-fidelity mesh + delaunay_elim_on_exit % Flag to run delaunay_elim on exit + improve_with_reduced_quality % Allow improvements with reduced quality end - - - - + methods - - - function obj = plot(obj) - if ~isempty(obj.pfix) && ~isempty(obj.egfix) - figure; hold on; % Initialize figure and hold for multiple plots - - % Plot the bounding boxes for visual aid - for box_number = 1:length(obj.boubox) - iboubox = obj.boubox{box_number}; - plot(iboubox(:,1),iboubox(:,2), 'g-', 'LineWidth', 2, 'DisplayName', 'Bounding Box'); - - % Plot outer boundaries - touter = obj.outer(box_number); - if ~isempty(touter) - tedges = Get_poly_edges(touter{1}); % Assuming Get_poly_edges is defined elsewhere - [touter, ~] = filter_polygon_constraints(touter{1}, tedges, obj.boubox, box_number); - plot(touter(:,1), touter(:,2), 'bx', 'DisplayName', 'Outer Boundary'); - end - - % Plot inner boundaries - tinner = obj.inner(box_number); - if ~isempty(tinner) && ~isempty(tinner{1}) - tedges = Get_poly_edges(tinner{1}); % Assuming Get_poly_edges is defined elsewhere - [tinner, ~] = filter_polygon_constraints(tinner{1}, tedges, obj.boubox, box_number); - plot(tinner(:,1), tinner(:,2), 'rx', 'DisplayName', 'Inner Boundary'); - end - end - - % Plot the fixed constraints as squares and edges in black - if exist('drawedge2', 'file') == 2 % Check if drawedge2 exists - drawedge2(obj.pfix, obj.egfix, [0,0,0]); % Assuming drawedge2 is a custom function - else - % Alternative plotting if drawedge2 is not available - plot(obj.pfix(:,1), obj.pfix(:,2), 'ks', 'MarkerSize', 8, 'MarkerFaceColor', 'k', 'DisplayName', 'Fixed Points'); - for i = 1:size(obj.egfix, 1) - plot(obj.pfix(obj.egfix(i,:), 1), obj.pfix(obj.egfix(i,:), 2), 'k-', 'LineWidth', 2, 'DisplayName', 'Fixed Edges'); - end - end - - axis equal; - title('Mesh Generation Constraints are black lines'); - xlabel('Longitude'); - ylabel('Latitude'); - legend('show', 'Location', 'bestoutside'); - else - disp('No constraints to plot!'); - end - end - - - - - % class constructor/default grd generation options + %% Constructor function obj = meshgen(varargin) - % Check for m_map dir - M_MAP_EXISTS=0; - if exist('m_proj','file')==2 - M_MAP_EXISTS=1 ; - end - if M_MAP_EXISTS~=1 - error('Where''s m_map? Chief, you need to read the user guide') - end - - - % Check for utilties dir - UTIL_DIR_EXISTS=0 ; - if exist('inpoly.m','file') - UTIL_DIR_EXISTS=1; - end - if UTIL_DIR_EXISTS~=1 - error('Where''s the utilities directory? Chief, you need to read the user guide') - end - - + % Parse input arguments and set defaults. p = inputParser; - % unpack options and set default ones, catch errors. - - - defval = 0; % placeholder value if arg is not passed. - % add name/value pairs + defval = 0; % placeholder for unspecified args addOptional(p,'h0',defval); addOptional(p,'bbox',defval); addOptional(p,'fh',defval); @@ -187,392 +82,384 @@ addOptional(p,'big_mesh',defval); addOptional(p,'proj',defval); addOptional(p,'qual_tol',defval); - addOptional(p,'enforceWeirs',1); + addOptional(p,'enforceWeirs',0); addOptional(p,'enforceMin',1); addOptional(p,'delaunay_elim_on_exit',1); addOptional(p,'improve_with_reduced_quality',0); - - - % parse the inputs parse(p,varargin{:}); - - - %if isempty(varargin); return; end - % store the inputs as a struct - inp = p.Results; - - - % kjr...order these argument so they are processed in a predictable - % manner. Process the general opts first, then the OceanMesh - % classes...then basic non-critical options. - inp = orderfields(inp,{'h0','bbox','enforceWeirs','enforceMin',... + inp = orderfields(p.Results,{'h0','bbox','enforceWeirs','enforceMin',... 'delaunay_elim_on_exit','improve_with_reduced_quality',... - 'fh',... - 'inner','outer','mainland',... - 'bou','ef',... %<--OceanMesh classes come after - 'egfix','pfix','fixboxes',... + 'fh','inner','outer','mainland',... + 'bou','ef','egfix','pfix','fixboxes',... 'plot_on','nscreen','itmax',... 'memory_gb','qual_tol','cleanup',... 'direc_smooth','dj_cutoff',... 'big_mesh','proj'}); - % get the fieldnames of the edge functions + + % Loop through options and assign to object properties. fields = fieldnames(inp); - % loop through and determine which args were passed. - % also, assign reasonable default values if some options were - % not assigned. - for i = 1 : numel(fields) - type = fields{i}; - switch type - % parse aux options first - case('h0') - % Provide in meters - obj.h0 = inp.(fields{i}); - case('fh') - if isa(inp.(fields{i}),'function_handle') - obj.fh = inp.(fields{i}); + for i = 1:numel(fields) + switch fields{i} + case 'h0' + obj.h0 = inp.h0; + case 'fh' + if isa(inp.fh, 'function_handle') + obj.fh = inp.fh; end - % can't check for errors here yet. - case('bbox') - obj.bbox= inp.(fields{i}); + case 'bbox' + obj.bbox = inp.bbox; if iscell(obj.bbox) - % checking bbox extents ob_min = obj.bbox{1}(:,1); ob_max = obj.bbox{1}(:,2); for ii = 2:length(obj.bbox) - if any(obj.bbox{ii}(:,1) < ob_min) || ... - any(obj.bbox{ii}(:,2) > ob_max) - error(['Outer bbox must contain all ' ... - 'inner bboxes: inner box #' ... - num2str(ii) ' violates this']) + if any(obj.bbox{ii}(:,1) < ob_min) || any(obj.bbox{ii}(:,2) > ob_max) + error(['Outer bbox must contain all inner bboxes: inner box #' num2str(ii) ' violates this']) end end end - - - % if user didn't pass anything explicitly for - % bounding box make it empty so it can be populated - % from ef as a cell-array - if obj.bbox(1)==0 - obj.bbox = []; - end - case('pfix') - obj.pfix= inp.(fields{i}); - if obj.pfix(1)~=0 - obj.pfix(:,:) = inp.(fields{i}); - else - obj.pfix = []; - end + if obj.bbox(1)==0, obj.bbox = []; end + case 'pfix' + obj.pfix = inp.pfix; + if ~isempty(obj.pfix) && obj.pfix(1) == 0, obj.pfix = []; end if obj.enforceWeirs - for j = 1 : length(obj.bou) - if ~isempty(obj.bou{j}.weirPfix) - obj.pfix = [obj.pfix ; obj.bou{j}.weirPfix]; + for j = 1:length(obj.bou) + if ~isempty(obj.bou{j}.weirPfix) + obj.pfix = [obj.pfix; obj.bou{j}.weirPfix]; end end end - case('egfix') - obj.egfix= inp.(fields{i}); - if ~isempty(obj.egfix) && obj.egfix(1)~=0 - obj.egfix = inp.(fields{i}); - else - obj.egfix = []; - end + case 'egfix' + obj.egfix = inp.egfix; + if ~isempty(obj.egfix) && obj.egfix(1)==0, obj.egfix = []; end if obj.enforceWeirs - for j = 1 : length(obj.bou) + for j = 1:length(obj.bou) if ~isempty(obj.bou{j}.weirEgfix) && ~isempty(obj.egfix) - obj.egfix = [obj.egfix ; obj.bou{j}.weirEgfix+max(obj.egfix(:))]; - elseif isempty(obj.egfix) - obj.egfix = obj.bou{j}.weirEgfix; + obj.egfix = [obj.egfix; obj.bou{j}.weirEgfix + max(obj.egfix(:))]; + else + obj.egfix = obj.bou{j}.weirEgfix; end end end obj.egfix = renumberEdges(obj.egfix); - case('fixboxes') - obj.fixboxes= inp.(fields{i}); - case('bou') - % got it from user arg - if obj.outer~=0, continue; end - - obj.outer = {} ; - obj.inner = {} ; - obj.mainland = {} ; - - obj.bou = inp.(fields{i}); - - % handle when not a cell + case 'fixboxes' + obj.fixboxes = inp.fixboxes; + case 'bou' + if obj.outer ~= 0, continue; end + obj.outer = {}; obj.inner = {}; obj.mainland = {}; + obj.bou = inp.bou; if ~iscell(obj.bou) - boutemp = obj.bou; - obj.bou = cell(1); - obj.bou{1} = boutemp; + obj.bou = {obj.bou}; end - - % then the geodata class was provide, unpack for ee = 1:length(obj.bou) - try - arg = obj.bou{ee} ; - catch - arg = obj.bou; - end - if isa(arg,'geodata') - + arg = obj.bou{ee}; + if isa(arg, 'geodata') obj.high_fidelity{ee} = obj.bou{ee}.high_fidelity; - obj.outer{ee} = obj.bou{ee}.outer; obj.inner{ee} = obj.bou{ee}.inner; - - % save bathy interpolant to meshgen if ~isempty(obj.bou{ee}.Fb) - obj.Fb{ee} = obj.bou{ee}.Fb ; + obj.Fb{ee} = obj.bou{ee}.Fb; end - - if ~isempty(obj.inner{ee}) && ... - obj.inner{ee}(1)~= 0 - obj.outer{ee} = [obj.outer{ee}; - obj.inner{ee}]; + if ~isempty(obj.inner{ee}) && obj.inner{ee}(1) ~= 0 + obj.outer{ee} = [obj.outer{ee}; obj.inner{ee}]; end obj.mainland{ee} = obj.bou{ee}.mainland; obj.boubox{ee} = obj.bou{ee}.boubox; obj.inpoly_flip{ee} = obj.bou{ee}.inpoly_flip; - end end - - - case('ef') - tmp = inp.(fields{i}); + case 'ef' + tmp = inp.ef; if isa(tmp, 'function_handle') - error('Please specify your edge function handle through the name/value pair fh'); + error('Please specify your edge function handle via the name/value pair fh'); end obj.ef = tmp; - - - % handle when not a cell if ~iscell(obj.ef) - eftemp = obj.ef; - obj.ef = cell(1); - obj.ef{1} = eftemp; + obj.ef = {obj.ef}; end - - - % Gather boxes from ef class. - for ee = 1 : length(obj.ef) - if isa(obj.ef{ee},'edgefx') + for ee = 1:length(obj.ef) + if isa(obj.ef{ee}, 'edgefx') obj.bbox{ee} = obj.ef{ee}.bbox; end end - - - % checking bbox extents if iscell(obj.bbox) ob_min = obj.bbox{1}(:,1); ob_max = obj.bbox{1}(:,2); for ii = 2:length(obj.bbox) - if any(obj.bbox{ii}(:,1) < ob_min) || ... - any(obj.bbox{ii}(:,2) > ob_max) - error(['Outer bbox must contain all ' ... - 'inner bboxes: inner box #' ... - num2str(ii) ' violates this']) + if any(obj.bbox{ii}(:,1) < ob_min) || any(obj.bbox{ii}(:,2) > ob_max) + error(['Outer bbox must contain all inner bboxes: inner box #' num2str(ii) ' violates this']) end end end - - - % kjr 2018 June: get h0 from edge functions for ee = 1:length(obj.ef) - if isa(obj.ef{ee},'edgefx') + if isa(obj.ef{ee}, 'edgefx') obj.h0(ee) = obj.ef{ee}.h0; end end - - - % kjr 2018 smooth the outer automatically - if length(obj.ef) > 1 - % kjr 2020, ensure the min. sizing func is - % used - if obj.enforceMin - obj.ef = enforce_min_ef(obj.ef); - end - obj.ef = smooth_outer(obj.ef,obj.Fb); + if length(obj.ef) > 1 && obj.enforceMin + obj.ef = enforce_min_ef(obj.ef); end - - - % Save the ef interpolants into the edgefx + obj.ef = smooth_outer(obj.ef, obj.Fb); for ee = 1:length(obj.ef) - if isa(obj.ef{ee},'edgefx') - obj.fh{ee} = @(p)obj.ef{ee}.F(p); + if isa(obj.ef{ee}, 'edgefx') + obj.fh{ee} = @(p) obj.ef{ee}.F(p); end end - - - case('plot_on') - obj.plot_on= inp.(fields{i}); - case('nscreen') - obj.nscreen= inp.(fields{i}); - if obj.nscreen ~=0 - obj.nscreen = inp.(fields{i}); + case 'plot_on' + obj.plot_on = inp.plot_on; + case 'nscreen' + obj.nscreen = inp.nscreen; + if obj.nscreen ~= 0 obj.plot_on = 1; else - obj.nscreen = 5; % default + obj.nscreen = 5; end - case('itmax') - obj.itmax= inp.(fields{i}); - if obj.itmax ~=0 - obj.itmax = inp.(fields{i}); - else + case 'itmax' + obj.itmax = inp.itmax; + if obj.itmax == 0 obj.itmax = 100; - warning('No itmax specified, itmax set to 100'); + warning('No itmax specified; defaulting to 100'); end - case('qual_tol') - obj.qual_tol = inp.(fields{i}); - if obj.qual_tol ~=0 - obj.qual_tol = inp.(fields{i}); - else - obj.qual_tol = 0.01; - end - case('inner') - if ~isa(obj.bou,'geodata') - obj.inner = inp.(fields{i}); + case 'qual_tol' + obj.qual_tol = inp.qual_tol; + if obj.qual_tol == 0, obj.qual_tol = 0.01; end + case 'inner' + if ~isa(obj.bou, 'geodata') + obj.inner = inp.inner; end - case('outer') - if ~isa(obj.bou,'geodata') - obj.outer = inp.(fields{i}); - if obj.inner(1)~=0 + case 'outer' + if ~isa(obj.bou, 'geodata') + obj.outer = inp.outer; + if obj.inner(1) ~= 0 obj.outer = [obj.outer; obj.inner]; end end - case('mainland') - if ~isa(obj.bou,'geodata') - obj.mainland = inp.(fields{i}); + case 'mainland' + if ~isa(obj.bou, 'geodata') + obj.mainland = inp.mainland; end - case('memory_gb') - if ~isa(obj.bou,'memory_gb') - obj.memory_gb = inp.(fields{i}); - end - case('cleanup') - obj.cleanup = inp.(fields{i}); + case 'memory_gb' + obj.memory_gb = inp.memory_gb; + case 'cleanup' + obj.cleanup = inp.cleanup; if isempty(obj.cleanup) || obj.cleanup == 0 obj.cleanup = 'none'; elseif obj.cleanup == 1 obj.cleanup = 'default'; end - case('dj_cutoff') - obj.dj_cutoff = inp.(fields{i}); - case('direc_smooth') - obj.direc_smooth = inp.(fields{i}); - case('proj') - obj.proj = inp.(fields{i}); - % default CPP - if obj.proj == 0; obj.proj = 'equi'; end + case 'dj_cutoff' + obj.dj_cutoff = inp.dj_cutoff; + case 'direc_smooth' + obj.direc_smooth = inp.direc_smooth; + case 'proj' + obj.proj = inp.proj; + if obj.proj == 0, obj.proj = 'equi'; end if ~isempty(obj.bbox) - % kjr Oct 2018 use outer coarsest box for - % multiscale meshing - lon_mi = obj.bbox{1}(1,1)-obj.h0(1)/1110; - lon_ma = obj.bbox{1}(1,2)+obj.h0(1)/1110; - lat_mi = obj.bbox{1}(2,1)-obj.h0(1)/1110; - lat_ma = obj.bbox{1}(2,2)+obj.h0(1)/1110; + lon_mi = obj.bbox{1}(1,1) - obj.h0(1)/1110; + lon_ma = obj.bbox{1}(1,2) + obj.h0(1)/1110; + lat_mi = obj.bbox{1}(2,1) - obj.h0(1)/1110; + lat_ma = obj.bbox{1}(2,2) + obj.h0(1)/1110; else - lon_mi = -180; lon_ma = 180; - lat_mi = -90; lat_ma = 90; + lon_mi = -180; lon_ma = 180; lat_mi = -90; lat_ma = 90; end - % Set up projected space - dmy = msh() ; + dmy = msh(); dmy.p(:,1) = [lon_mi; lon_ma]; dmy.p(:,2) = [lat_mi; lat_ma]; - del = setProj(dmy,1,obj.proj) ; - case('enforceWeirs') - obj.enforceWeirs = inp.(fields{i}); - case('enforceMin') - obj.enforceMin = inp.(fields{i}); - case('delaunay_elim_on_exit') - obj.delaunay_elim_on_exit = inp.(fields{i}); - case('improve_with_reduced_quality') - obj.improve_with_reduced_quality = inp.(fields{i}); + setProj(dmy,1,obj.proj); + case 'enforceWeirs' + obj.enforceWeirs = inp.enforceWeirs; + case 'enforceMin' + obj.enforceMin = inp.enforceMin; + case 'delaunay_elim_on_exit' + obj.delaunay_elim_on_exit = inp.delaunay_elim_on_exit; + case 'improve_with_reduced_quality' + obj.improve_with_reduced_quality = inp.improve_with_reduced_quality; end end - - - if isempty(varargin); return; end - - - % error checking - if isempty(obj.boubox) && ~isempty(obj.bbox) - % Make the bounding box 5 x 2 matrix in clockwise order if - % it isn't present. This case must be when the user is - % manually specifying the PSLG. - obj.boubox{1} = [obj.bbox(1,1) obj.bbox(2,1); - obj.bbox(1,1) obj.bbox(2,2); ... - obj.bbox(1,2) obj.bbox(2,2); - obj.bbox(1,2) obj.bbox(2,1); ... - obj.bbox(1,1) obj.bbox(2,1); NaN NaN]; + + % Essential input checks. + if any(obj.h0 == 0) + error('h0 (minimum edge length) was not correctly specified!'); end - if any(obj.h0==0), error('h0 was not correctly specified!'), end - if isempty(obj.outer), error('no outer boundary specified!'), end - if isempty(obj.bbox), error('no bounding box specified!'), end - obj.fd = @dpoly; % <-default distance fx accepts p and pv (outer polygon). - % kjr build ANN object into meshgen - obj = createANN(obj) ; - - + if isempty(obj.outer) + error('No outer boundary specified!'); + end + if isempty(obj.bbox) + error('No bounding box specified!'); + end + + % Set default distance function and build ANN. + obj.fd = @dpoly; + obj = createANN(obj); global MAP_PROJECTION MAP_COORDS MAP_VAR_LIST - obj.grd.proj = MAP_PROJECTION ; - obj.grd.coord = MAP_COORDS ; - obj.grd.mapvar = MAP_VAR_LIST ; - - % Check and update geometry for high fidelity option + obj.grd.proj = MAP_PROJECTION; + obj.grd.coord = MAP_COORDS; + obj.grd.mapvar = MAP_VAR_LIST; + + % Generate breakline constraints from shoreline and island boundaries. + [tpfix, tegfix] = obj.generateBreaklineConstraints(); + + % Check if any entry in high_fidelity is 3 + if any(cellfun(@(x) isscalar(x) && x == 3, obj.high_fidelity)) + disp(' Pruning breakline connectivity due to high-fidelity mode 3...'); + + % Set tolerance and angle threshold + tol = (obj.h0(end)) / 111e3 ; + + maxChainLengthKm = 2.5; + + % Apply CleanPSLG processing with timing + disp(' Starting PSLG cleaning process...'); + + % Apply the filtering pipeline + pslg = CleanPSLG(tpfix, tegfix, tol, 0, tol); + pslg = pslg.dropNearbyEdges(); + pslg = pslg.deleteShortChains(maxChainLengthKm); + pslg = pslg.cleanUnusedVertices(); + + + % Update points and edges + tpfix = pslg.Vertices; + tegfix = pslg.Segments; + + [tpfix, tegfix] = fixgeo2(tpfix, tegfix); + end + + % Check for duplicate fixed points. + if ~isempty(tpfix) + checkfixp = setdiff(tpfix, fixmesh(tpfix), 'rows'); + if ~isempty(checkfixp) + error('Duplicate fixed points detected, cannot proceed'); + end + end + + % Update object with user-supplied and generated constraints. + obj.pfix = [obj.pfix; tpfix]; + if isempty(obj.egfix) + obj.egfix = tegfix; + else + obj.egfix = [obj.egfix; tegfix + max(obj.egfix(:))]; + end + + if ~isempty(obj.pfix) + disp(['Using ', num2str(size(obj.pfix,1)), ' fixed points.']); + end + if ~isempty(obj.egfix) + if max(obj.egfix(:)) > size(obj.pfix,1) + error('FATAL: egfix indices exceed the number of fixed points.'); + end + disp(['Using ', num2str(size(obj.egfix,1)), ' fixed edges.']); + warning('Please verify fixed constraints using plot.'); + end + end + + + + %% Private method: generateBreaklineConstraints + function [tpfix, tegfix] = generateBreaklineConstraints(obj) + % Generates breakline constraints (fixed points and edges) from the + % provided mainland, inner (island) boundaries, and line strings. + % + % Outputs: + % tpfix - Accumulated fixed points from polygon boundaries and line strings. + % tegfix - Corresponding fixed edge constraints. + tpfix = []; tegfix = []; + for box_num = 1:length(obj.h0) + % Collect polygon data from mainland and inner boundaries. + polys = {}; + if ~isempty(obj.mainland{box_num}) + polys{end+1} = obj.mainland{box_num}; + end + if ~isempty(obj.inner{box_num}) && obj.inner{box_num}(1) ~= 0 + polys{end+1} = obj.inner{box_num}; + end - - % High fidelity - formation of point & edge constraints - if obj.high_fidelity{box_num} + % Collect NaN-separated line strings + lineStrings = {}; + if ~isempty(obj.bou{box_num}.linestrings) + lineStrings{end+1} = obj.bou{box_num}.linestrings; % Store line strings + end - % Define polygons based on available geometry data - ml = obj.mainland{box_num}; - il = obj.inner{box_num}; - polys = {}; - if ~isempty(il), polys{end+1} = il; end - if ~isempty(ml), polys{end+1} = ml; end - - if obj.cleanup == 1 - warning('Setting cleanup to 0 since high_fidelity mode is on'); + % Process only if high-fidelity is enabled for this box. + if obj.high_fidelity{box_num} + if isequal(obj.cleanup, 1) + warning('Disabling cleanup since high-fidelity mode is active.'); obj.cleanup = 0; end - disp(['Redistributing vertices for box #' num2str(box_num)]); - - % Convert polygon data, removing NaNs and duplicates - poly = cell2mat(polys'); - D = nandelim_to_cell(poly); % Custom function to handle NaN-delimited cells - D = D(~cellfun(@(p) all(isnan(p(:))), D)); % Remove cells with all NaNs - areas = cellfun(@(d) polyarea(d(:, 1), d(:, 2)), D); - [~, uniqueIdx] = uniquetol(areas, 1e-12); % Remove duplicate polygons + disp(['Redistributing vertices for box #', num2str(box_num)]); + + % Process Polygons: Concatenate and remove NaNs + poly_all = cell2mat(polys'); + D = nandelim_to_cell(poly_all); % Split polygons at NaNs. + D = D(~cellfun(@(p) all(isnan(p(:))), D)); % Remove empty segments + areas = cellfun(@(d) polyarea(d(:,1), d(:,2)), D); + [~, uniqueIdx] = uniquetol(areas, 1e-12); D = D(uniqueIdx); - - % Process each polygon + + % Process each polygon separately for i = 1:length(D) - my_poly = D{i}; - if size(my_poly, 1) > 2 - tedges = Get_line_edges(my_poly); - [pts, bnde] = filter_polygon_constraints(my_poly, tedges, obj.boubox, box_num); + current_poly = D{i}; + if size(current_poly,1) > 2 + polyEdges = Get_line_edges(current_poly); + [pts, bnde] = filter_polygon_constraints(current_poly, polyEdges, obj.boubox, box_num); if isempty(bnde), continue; end - poly_split = extdom_polygon(bnde, pts, 0, 1); - - for ii = 1:length(poly_split) - points = unique(poly_split{ii}, 'rows', 'stable'); - % Option 2 simply constrains the vector as - % it is in the file. + + % Split polygon if necessary. + polySplits = extdom_polygon(bnde, pts, 0, 1); + for j = 1:length(polySplits) + points = unique(polySplits{j}, 'rows', 'stable'); + + % Generate fixed constraints based on high_fidelity option. if obj.high_fidelity{box_num} == 2 tmp_pfix = points; - tmp_egfix = Get_poly_edges([points; NaN,NaN]); - % Option 1: resamples based on the local - % mesh resolution - elseif obj.high_fidelity{box_num} == 1 - [tmp_pfix, tmp_egfix] = mesh1d(points, obj.fh, obj.h0./111e3,... - [], obj.boubox, box_num, []); + tmp_egfix = Get_poly_edges([points; NaN, NaN]); + elseif obj.high_fidelity{box_num} == 1 || obj.high_fidelity{box_num} == 3 + [tmp_pfix, tmp_egfix] = mesh1d(points, obj.fh, obj.h0./111e3, [], obj.boubox, box_num, []); + else + continue; end - - if size(tmp_pfix, 1) > 2 + + if size(tmp_pfix,1) > 2 [tmp_pfix, tmp_egfix] = fixgeo2(tmp_pfix, tmp_egfix); - if max(tmp_egfix(:)) ~= size(tmp_pfix, 1), continue; end + if max(tmp_egfix(:)) ~= size(tmp_pfix,1), continue; end + tpfix = [tpfix; tmp_pfix]; + if isempty(tegfix) + tegfix = tmp_egfix; + else + tegfix = [tegfix; tmp_egfix + max(tegfix(:))]; + end + end + end + end + end + + % Process NaN-separated Line Strings + for i = 1:length(lineStrings) + current_lines = lineStrings{i}; + if size(current_lines,1) > 1 + % Convert NaN-separated line strings into individual segments + lineSegments = nandelim_to_cell(current_lines); % Split at NaNs + lineSegments = lineSegments(~cellfun(@isempty, lineSegments)); % Remove empty segments + + for j = 1:length(lineSegments) + lineSeg = unique(lineSegments{j}, 'rows', 'stable'); + if size(lineSeg,1) < 2, continue; end % Ensure it's valid + + % Generate fixed constraints + if obj.high_fidelity{box_num} == 2 + tmp_pfix = lineSeg; + tmp_egfix = Get_poly_edges([lineSeg; NaN, NaN]); + elseif obj.high_fidelity{box_num} == 1 || obj.high_fidelity{box_num} == 3 + [tmp_pfix, tmp_egfix] = mesh1d(lineSeg, obj.fh, (3*obj.h0)./111e3, [], obj.boubox, box_num, []); + else + continue; + end + + if size(tmp_pfix,1) > 1 + [tmp_pfix, tmp_egfix] = fixgeo2(tmp_pfix, tmp_egfix); + if max(tmp_egfix(:)) ~= size(tmp_pfix,1), continue; end tpfix = [tpfix; tmp_pfix]; if isempty(tegfix) tegfix = tmp_egfix; @@ -585,37 +472,80 @@ end end end - - % Final geometry check and update - [tpfix, tegfix] = fixgeo2(tpfix, tegfix); - checkfixp = setdiff(tpfix, fixmesh(tpfix), 'rows'); - if ~isempty(checkfixp), error('Duplicate fixed points detected, cannot proceed'); end - % Update object with user-passed fixed points - obj.pfix = [obj.pfix; tpfix]; - - if isempty(obj.egfix) - obj.egfix = tegfix; - else - obj.egfix = [obj.egfix; tegfix + max(obj.egfix(:))]; - end - - if ~isempty(obj.pfix) - disp(['Using ', num2str(length(obj.pfix)), ' fixed points.']); - end - - if ~isempty(obj.egfix) - if max(obj.egfix(:)) > length(obj.pfix) - error('FATAL: egfix does not index correctly into pfix.'); + % Final adjustment: re-index and remove duplicate fixed constraints. + [tpfix, tegfix] = fixgeo2(tpfix, tegfix); + end + + + function obj = plot(obj) + if ~isempty(obj.pfix) && ~isempty(obj.egfix) + figure; hold on; grid on; + + % Set colors and styles + bboxColor = [0, 0.7, 0]; % Green for bounding boxes + outerColor = [0, 0, 1]; % Blue for outer boundaries + innerColor = [1, 0, 0]; % Red for inner boundaries + edgeColor = [0, 0, 0]; % Black for fixed edges + fixedPointColor = [0, 0, 0]; % Black for fixed points + + % Loop over bounding boxes and plot constraints + for box_number = 1:length(obj.boubox) + iboubox = obj.boubox{box_number}; + + % Plot bounding box with slight transparency + plot(iboubox(:,1), iboubox(:,2), '-', 'Color', bboxColor, 'LineWidth', 2, ... + 'DisplayName', 'Bounding Box', 'HandleVisibility', 'off'); + + % Plot outer boundary (if available) + touter = obj.outer(box_number); + if ~isempty(touter) + tedges = Get_poly_edges(touter{1}); + [touter, ~] = filter_polygon_constraints(touter{1}, tedges, obj.boubox, box_number); + plot(touter(:,1), touter(:,2), 'o', 'Color', outerColor, 'MarkerSize', 5, ... + 'DisplayName', 'Outer Boundary'); + end + + % Plot inner boundary (if available) + tinner = obj.inner(box_number); + if ~isempty(tinner) && ~isempty(tinner{1}) + tedges = Get_poly_edges(tinner{1}); + [tinner, ~] = filter_polygon_constraints(tinner{1}, tedges, obj.boubox, box_number); + plot(tinner(:,1), tinner(:,2), 'x', 'Color', innerColor, 'MarkerSize', 5, ... + 'DisplayName', 'Inner Boundary'); + end + end + + % Plot fixed edges and points + if exist('drawedge2', 'file') == 2 + drawedge2(obj.pfix, obj.egfix, edgeColor); + else + % Plot fixed points + scatter(obj.pfix(:,1), obj.pfix(:,2), 40, fixedPointColor, 'filled', ... + 'DisplayName', 'Fixed Points'); + + % Plot fixed edges with a thicker line + for i = 1:size(obj.egfix, 1) + plot(obj.pfix(obj.egfix(i,:), 1), obj.pfix(obj.egfix(i,:), 2), '-', ... + 'Color', edgeColor, 'LineWidth', 2.5, 'DisplayName', 'Fixed Edges'); + end end - disp(['Using ', num2str(length(obj.egfix)), ' fixed edges.']); - warning('Please check fixed constraints: plot(mshopts)'); + + % Improve plot formatting + axis equal; + title('Constrained Breaklines and Mesh Constraints', 'FontWeight', 'bold'); + xlabel('Longitude', 'FontSize', 12); + ylabel('Latitude', 'FontSize', 12); + + % Remove duplicate legend entries + legendEntries = findobj(gca, '-property', 'DisplayName'); + else + disp('No constraints to plot!'); end end - - - % Creates Approximate nearest neighbor objects on start-up - function obj = createANN(obj) + + + function obj = createANN(obj) box_vec = 1:length(obj.bbox); for box_num = box_vec if ~iscell(obj.outer) @@ -626,8 +556,6 @@ dataset(isnan(obj.outer{box_num}(:,1)),:) = []; end if all(abs(obj.bbox{box_num}(1,:)) == 180) - % This line removes the line that can appear in the - % center for a global mesh dataset(abs(dataset(:,1)) > 180-1e-6,:) = []; dataset(abs(dataset(:,1)) < 1e-6,:) = []; end @@ -635,11 +563,60 @@ dataset(isnan(dataset(:,1)),:) = []; dmy = ann(dataset'); obj.anno{box_num} = dmy; - obj.annData{box_num}=dataset; + obj.annData{box_num} = dataset; end end - - + + function mesh_out = collapse_thin_triangles(obj, aspect_ratio_threshold) + % Identify and collapse thin triangles in the mesh + % aspect_ratio_threshold: Defines what is considered "thin" + + tri = obj.t; % Get triangle connectivity + nodes = obj.p; % Get node coordinates + + num_tri = size(tri, 1); + + for i = 1:num_tri + % Get triangle node indices + n1 = tri(i, 1); + n2 = tri(i, 2); + n3 = tri(i, 3); + + % Compute edge lengths + e1 = norm(nodes(n2, :) - nodes(n1, :)); % Edge 1-2 + e2 = norm(nodes(n3, :) - nodes(n2, :)); % Edge 2-3 + e3 = norm(nodes(n1, :) - nodes(n3, :)); % Edge 3-1 + + % Find longest edge + [longest_edge, idx] = max([e1, e2, e3]); + + % Compute aspect ratio (shortest to longest) + shortest_edge = min([e1, e2, e3]); + aspect_ratio = shortest_edge / longest_edge; + + % If the triangle is too thin, collapse it + if aspect_ratio < aspect_ratio_threshold + % Find the node opposite the longest edge + switch idx + case 1, opposite_node = n3; edge_nodes = [n1, n2]; + case 2, opposite_node = n1; edge_nodes = [n2, n3]; + case 3, opposite_node = n2; edge_nodes = [n3, n1]; + end + + % Compute midpoint of longest edge + midpoint = mean(nodes(edge_nodes, :), 1); + + % Move the opposite node to the midpoint (collapsing the triangle) + nodes(opposite_node, :) = midpoint; + end + end + + % Update mesh with modified nodes + mesh_out = mesh; + mesh_out.p = nodes; + end + + function obj = build(obj) % 2-D Mesh Generator using Distance Functions. % Checking existence of major inputs @@ -655,7 +632,7 @@ delIT = 0 ; delImp = 2; imp = 10; % number of iterations to do mesh improvements (delete/add) EXIT_QUALITY = 0.30; % minimum quality necessary to terminate if iter < itmax - + % unpack initial points. p = obj.grd.p; if isempty(p) @@ -702,7 +679,7 @@ [ys;ys])/(2/sqrt(3)*h0_l)) + ... floor(1e3*m_lldist([0;ed],... [ys;ys])/(2/sqrt(3)*h0_l)); - + else nx = floor(1e3*m_lldist([st;ed],... [ys;ys])/(2/sqrt(3)*h0_l)); @@ -723,14 +700,14 @@ y(ns:ne) = ys; end ns = ne+1; ys = ys + dy; - + end - + st = ed; ed = st + blklen; p1 = [x(:) y(:)]; clear x y - - + + %% 2. Remove points outside the region, apply the rejection method p1 = p1(feval(obj.fd,p1,obj,box_num) < geps,:); % Keep only d<0 points r0 = 1./feval(fh_l,p1).^2; % Probability to keep point @@ -754,7 +731,7 @@ obj.grd.b = []; h0_l = obj.h0(end); % finest h0 (in case of a restart of meshgen.build). end - + nfix = length(obj.pfix); negfix = length(obj.egfix); if ~isempty(obj.pfix); p = [obj.pfix; p]; end % kjr July 2023, set to these values for better convg. @@ -763,15 +740,15 @@ deltat=0.10; end - % Check if any boxes are set to high-fidelity - % If so turn off heal_fixed_edges + % Check if any boxes are set to high-fidelity + % If so turn off heal_fixed_edges HIGH_FIDELITY_MODE = 0; for i = 1 : length(obj.h0) if obj.high_fidelity{i} HIGH_FIDELITY_MODE = 1; end end - + N = size(p,1); % Number of points N disp(['Number of initial points after rejection is ',num2str(N)]); %% Iterate @@ -818,7 +795,7 @@ if ~mod(it,imp+1) && ((obj.qual(it,1) - obj.qual(it-1,1) < -0.10) || ... (~obj.improve_with_reduced_quality && ... (N - length(p_before_improve))/length(p_before_improve) < -0.10)) - + disp('Mesh improvement was unsuccessful...rewinding...'); p = p_before_improve; N = size(p,1); % Number of points changed @@ -921,8 +898,8 @@ [ideal_bars(:,1),ideal_bars(:,2)] = ... % needs to be in non-projected m_xy2ll(ideal_bars(:,1),ideal_bars(:,2)); % coordinates hbars = 0*ideal_bars(:,1); - - + + for box_num = 1:length(obj.h0) % For each bbox, find the bars that are in and calculate if ~iscell(obj.fh) % their ideal lengths. fh_l = obj.fh; @@ -939,12 +916,12 @@ end hbars(inside) = feval(fh_l,ideal_bars(inside,:)); % Ideal lengths end - - + + L0 = hbars*Fscale*median(L)/median(hbars); % L0 = Desired lengths using ratio of medians scale factor LN = L./L0; % LN = Normalized bar lengths - - + + % Mesh improvements (deleting and addition) p_before_improve = p; if ~mod(it,imp) % @@ -955,8 +932,8 @@ % Remove elements with small connectivity nn = get_small_connectivity(p,t); disp(['Deleting ' num2str(length(nn)) ' due to small connectivity']) - - + + % Remove points that are too close (< LN = 0.5) if any(LN < 0.5) % Do not delete pfix too close. @@ -964,8 +941,8 @@ disp(['Deleting ' num2str(length(nn1)) ' points too close together']) nn = unique([nn; nn1]); end - - + + % Split long edges however many times to % better lead to LN of 1 if any(LN > 2) @@ -998,22 +975,22 @@ continue; end end - - + + F = (1-LN.^4).*exp(-LN.^4)./LN; % Bessens-Heckbert edge force - F(isinf(F)) = 0; + F(isinf(F)) = 0; Fvec = F*[1,1].*barvec; - - + + Ftot = full(sparse(bars(:,[1,1,2,2]),ones(size(F))*[1,2,1,2],[Fvec,-Fvec],N,2)); Ftot(1:nfix,:) = 0; % Force = 0 at fixed points - + pt = pt + deltat*Ftot; % Update node positions - - + + [p(:,1),p(:,2)] = m_xy2ll(pt(:,1),pt(:,2)); - - + + %7. Bring outside points back to the boundary d = feval(obj.fd,p,obj,[],1); ix = d > 0; % Find points outside (d>0) ix(1:nfix) = 0; @@ -1032,12 +1009,12 @@ end alpha = alpha / 0.5; end - - + + % 8. Termination criterion: Exceed itmax it = it + 1 ; - - + + if ( it > obj.itmax ) % Do the final deletion of small connectivity if obj.delaunay_elim_on_exit @@ -1053,12 +1030,12 @@ %% disp('Finished iterating...'); fprintf(1,' ------------------------------------------------------->\n') ; - - + + %% Doing the final cleaning and fixing to the mesh... % Always save the mesh! save('Precleaned_grid.mat','it','p','t'); - + % Clean up the mesh if specified if ~strcmp(obj.cleanup,'none') % Put the mesh class into the grd part of meshgen and clean @@ -1080,13 +1057,13 @@ obj.grd.pfix = obj.pfix ; obj.grd.egfix= obj.egfix ; end - - + + % Check element order, important for the global meshes crossing % -180/180 boundary obj.grd = CheckElementOrder(obj.grd); - - + + if obj.plot_on figure; plot(obj.qual,'linewi',2); hold on @@ -1103,8 +1080,8 @@ %%%%%%%%%%%%%%%%%%%%%%%%%% % Auxiliary subfunctions % %%%%%%%%%%%%%%%%%%%%%%%%%% - - + + function [t,p] = delaunay_elim(p,fd,geps,final) % Removing mean to reduce the magnitude of the points to % help the convex calc @@ -1137,8 +1114,8 @@ [p(:,1),p(:,2)] = m_xy2ll(pt1(:,1),pt1(:,2)); end end - - + + function nn = get_small_connectivity(p,t) % Get node connectivity (look for 4) [~, enum] = VertToEle(t); @@ -1150,7 +1127,7 @@ return; end - function del = heal_fixed_edges(p,t,egfix) + function del = heal_fixed_edges(p,t,egfix) % kjr april2019 % if there's a triangle with a low geometric quality that % contains a fixed edge, remove the non-fixed vertex @@ -1166,12 +1143,11 @@ badtria = t(dmy,:); del = badtria(badtria > nfix) ; end - + end % end mesh generator - - + + end % end methods - - -end % end class + +end % end class \ No newline at end of file diff --git a/@meshgen/private/Get_line_edges.m b/@meshgen/private/Get_line_edges.m index bb666b2f..a5cc43f4 100644 --- a/@meshgen/private/Get_line_edges.m +++ b/@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 diff --git a/@meshgen/private/filter_polygon_constraints.m b/@meshgen/private/filter_polygon_constraints.m index b1c50d86..4fcac7c6 100644 --- a/@meshgen/private/filter_polygon_constraints.m +++ b/@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 \ No newline at end of file diff --git a/utilities/drawedge2.m b/utilities/drawedge2.m index ff72d721..01dade9a 100644 --- a/utilities/drawedge2.m +++ b/utilities/drawedge2.m @@ -19,6 +19,3 @@ function drawedge2(pp,e2,color) end - - -