module fortplot_contour_regions !! Contour polygon decomposition for extracting filled regions !! !! This module provides functions for decomposing 2D contour data into !! polygonal regions between contour levels, supporting filled contour !! visualization across PNG, PDF, and ASCII backends. !! !! = Key Types = !! - contour_region_t: Represents a single contour region with boundary !! - contour_polygon_t: Polygon boundary with points !! !! = Core Functions = !! - extract_contour_regions: Main polygon decomposition function !! - get_region_boundaries: Extract boundary points for regions !! !! Author: fortplot contributors use iso_fortran_env, only: wp => real64 use fortplot_constants, only: EPSILON_GEOMETRY use fortplot_contour_algorithms, only: calculate_marching_squares_config, & get_contour_lines implicit none private ! Public interface public :: contour_region_t, contour_polygon_t public :: extract_contour_regions type :: contour_polygon_t !! Polygon boundary with ordered vertex points real(wp), allocatable :: x(:) ! X coordinates of boundary vertices real(wp), allocatable :: y(:) ! Y coordinates of boundary vertices logical :: is_closed = .false. ! Whether polygon is closed end type contour_polygon_t type :: contour_region_t !! Region between contour levels with boundary polygons real(wp) :: level_min = 0.0_wp ! Lower contour level real(wp) :: level_max = 0.0_wp ! Upper contour level type(contour_polygon_t), allocatable :: boundaries(:) ! Boundary polygons integer :: region_id = 0 ! Unique region identifier end type contour_region_t contains function extract_contour_regions(x_grid, y_grid, z_grid, levels) result(regions) !! Extract polygonal regions between contour levels using marching squares !! !! Given: 2D grid data and contour levels !! When: Processing grid to find regions between levels !! Then: Returns array of regions with boundary polygons !! !! This is the main function for contour polygon decomposition that !! enables filled contour rendering across all backends. real(wp), contiguous, intent(in) :: x_grid(:) ! X coordinates of grid points real(wp), contiguous, intent(in) :: y_grid(:) ! Y coordinates of grid points real(wp), contiguous, intent(in) :: z_grid(:, :) ! Z values at grid points real(wp), contiguous, intent(in) :: levels(:) ! Contour levels to extract type(contour_region_t), allocatable :: regions(:) integer :: n_levels, n_regions integer :: i n_levels = size(levels) n_regions = n_levels + 1 ! One more region than levels ! Allocate regions array allocate(regions(n_regions)) ! Create regions between and around contour levels do i = 1, n_regions regions(i)%region_id = i if (i == 1) then ! Region below first level regions(i)%level_min = minval(z_grid) - 1.0_wp regions(i)%level_max = levels(1) else if (i == n_regions) then ! Region above last level regions(i)%level_min = levels(n_levels) regions(i)%level_max = maxval(z_grid) + 1.0_wp else ! Region between two levels regions(i)%level_min = levels(max(1, i - 1)) regions(i)%level_max = levels(min(n_levels, i)) end if ! Extract boundary polygons for this region call extract_region_boundaries(x_grid, y_grid, z_grid, & regions(i)%level_min, & regions(i)%level_max, & regions(i)%boundaries) end do end function extract_contour_regions subroutine extract_region_boundaries(x_grid, y_grid, z_grid, level_min, & level_max, boundaries) !! Extract closed boundary rings for the band [level_min, level_max]. !! !! The band boundary is the union of the level_min and level_max !! isolines. Each isoline is traced with standard single-level marching !! squares (shared with the line renderer) and the resulting segments !! are chained into closed rings. Tracing each level independently avoids !! the stray grid-aligned fragments a direct isoband walk produced. real(wp), contiguous, intent(in) :: x_grid(:) real(wp), contiguous, intent(in) :: y_grid(:) real(wp), contiguous, intent(in) :: z_grid(:, :) real(wp), intent(in) :: level_min real(wp), intent(in) :: level_max type(contour_polygon_t), allocatable, intent(out) :: boundaries(:) real(wp), allocatable :: contour_x(:), contour_y(:) integer :: contour_count allocate(contour_x(0)) allocate(contour_y(0)) contour_count = 0 call collect_level_segments(level_min) call collect_level_segments(level_max) call finalize_boundaries(contour_x, contour_y, contour_count, boundaries) contains subroutine collect_level_segments(level) real(wp), intent(in) :: level integer :: nx, ny, i, j, config, num_lines real(wp) :: xc(4), yc(4), zc(4) real(wp) :: line_points(8) nx = size(x_grid) ny = size(y_grid) do j = 1, ny - 1 do i = 1, nx - 1 xc(1) = x_grid(i) ; yc(1) = y_grid(j) xc(2) = x_grid(i + 1); yc(2) = y_grid(j) xc(3) = x_grid(i + 1); yc(3) = y_grid(j + 1) xc(4) = x_grid(i) ; yc(4) = y_grid(j + 1) zc(1) = z_grid(j , i ) zc(2) = z_grid(j , i + 1) zc(3) = z_grid(j + 1, i + 1) zc(4) = z_grid(j + 1, i ) call calculate_marching_squares_config(zc(1), zc(2), zc(3), & zc(4), level, config) call get_contour_lines(config, xc(1), yc(1), xc(2), yc(2), & xc(3), yc(3), xc(4), yc(4), zc(1), & zc(2), zc(3), zc(4), level, & line_points, num_lines) if (num_lines >= 1) then call add_pair(line_points(1), line_points(2), & line_points(3), line_points(4)) end if if (num_lines >= 2) then call add_pair(line_points(5), line_points(6), & line_points(7), line_points(8)) end if end do end do end subroutine collect_level_segments subroutine add_pair(x1, y1, x2, y2) real(wp), intent(in) :: x1, y1, x2, y2 real(wp), allocatable :: tx(:), ty(:) integer :: n n = size(contour_x) allocate(tx(n+2)) allocate(ty(n+2)) if (n > 0) then tx(1:n) = contour_x ty(1:n) = contour_y end if tx(n+1) = x1; ty(n+1) = y1 tx(n+2) = x2; ty(n+2) = y2 call move_alloc(tx, contour_x) call move_alloc(ty, contour_y) contour_count = contour_count + 2 end subroutine add_pair end subroutine extract_region_boundaries subroutine finalize_boundaries(contour_x, contour_y, contour_count, boundaries) !! Create one or more boundary polygons by chaining contour segments real(wp), contiguous, intent(in) :: contour_x(:), contour_y(:) integer, intent(in) :: contour_count type(contour_polygon_t), allocatable, intent(out) :: boundaries(:) integer :: nseg, k, i_poly, max_pts real(wp), allocatable :: xs(:), ys(:) ! segment endpoints flattened [x1,y1,x2,y2] per segment logical, allocatable :: used(:) real(wp), allocatable :: px(:), py(:) integer :: npts, s real(wp) :: sx, sy, cx, cy, nxp, nyp ! Internal helpers are defined in the CONTAINS block below ! No segments collected if (contour_count < 2) then allocate(boundaries(0)) return end if ! Interpret consecutive point pairs as independent segments nseg = contour_count / 2 allocate(xs(2*nseg)) allocate(ys(2*nseg)) allocate(used(nseg)) used = .false. do k = 1, nseg xs(2*k-1) = contour_x(2*k-1) ys(2*k-1) = contour_y(2*k-1) xs(2*k ) = contour_x(2*k ) ys(2*k ) = contour_y(2*k ) end do ! We do not know polygon count upfront allocate(boundaries(0)) i_poly = 0 max_pts = contour_count + 1 allocate(px(max_pts)) allocate(py(max_pts)) do ! Find next unused segment to start a chain s = 0 do k = 1, nseg if (.not. used(k)) then s = k exit end if end do if (s == 0) exit ! no more segments ! Initialize chain with segment s sx = xs(2*s-1); sy = ys(2*s-1) cx = xs(2*s ); cy = ys(2*s ) used(s) = .true. npts = 2 px(1) = sx; py(1) = sy px(2) = cx; py(2) = cy ! Greedily connect subsequent segments by matching endpoints do ! Closed loop formed? if (points_close(cx, cy, sx, sy)) then npts = npts + 1 px(npts) = sx py(npts) = sy exit end if ! Search for a segment connecting to current endpoint s = 0 do k = 1, nseg if (used(k)) cycle if (points_close(xs(2*k-1), ys(2*k-1), cx, cy)) then nxp = xs(2*k ); nyp = ys(2*k ) s = k exit else if (points_close(xs(2*k ), ys(2*k ), cx, cy)) then nxp = xs(2*k-1); nyp = ys(2*k-1) s = k exit end if end do if (s == 0) then ! Could not close chain; discard open polyline npts = 0 exit end if ! Append next point and continue used(s) = .true. cx = nxp; cy = nyp npts = npts + 1 if (npts > max_pts) exit px(npts) = cx py(npts) = cy end do ! If a closed polygon was formed with >= 4 points (including closure), normalize orientation and store it if (npts >= 4) then call normalize_and_append(boundaries, px, py, npts) end if end do if (.not. allocated(boundaries)) then allocate(boundaries(0)) end if contains logical function points_close(ax, ay, bx, by) real(wp), intent(in) :: ax, ay, bx, by real(wp) :: dx, dy dx = ax - bx dy = ay - by points_close = (abs(dx) <= EPSILON_GEOMETRY .and. abs(dy) <= EPSILON_GEOMETRY) end function points_close subroutine normalize_and_append(arr, vx, vy, n) type(contour_polygon_t), allocatable, intent(inout) :: arr(:) real(wp), contiguous, intent(in) :: vx(:), vy(:) integer, intent(in) :: n type(contour_polygon_t), allocatable :: tmp(:) integer :: oldn real(wp) :: area integer :: i oldn = size(arr) allocate(tmp(oldn + 1)) if (oldn > 0) tmp(1:oldn) = arr call move_alloc(tmp, arr) allocate(arr(oldn + 1)%x(n)) allocate(arr(oldn + 1)%y(n)) ! Compute signed area to enforce CCW ordering for consistency area = 0.0_wp do i = 1, n-1 area = area + (vx(i) * vy(i+1) - vx(i+1) * vy(i)) end do area = area + (vx(n) * vy(1) - vx(1) * vy(n)) if (area < 0.0_wp) then ! Reverse to ensure CCW do i = 1, n arr(oldn + 1)%x(i) = vx(n - i + 1) arr(oldn + 1)%y(i) = vy(n - i + 1) end do else arr(oldn + 1)%x(1:n) = vx(1:n) arr(oldn + 1)%y(1:n) = vy(1:n) end if arr(oldn + 1)%is_closed = .true. end subroutine normalize_and_append end subroutine finalize_boundaries end module fortplot_contour_regions