type :: contour_set_t
    ! Contour identification
    character(len=50) :: label = ''             ! Dataset label
    logical :: filled = .false.                 ! Line vs filled contours

    ! Grid data (regular grid support)
    real(wp), allocatable :: x(:), y(:)         ! 1D coordinate arrays
    real(wp), allocatable :: z(:,:)             ! 2D data grid (z(ny, nx))
    integer :: nx, ny                           ! Grid dimensions

    ! Contour levels and colors
    real(wp), allocatable :: levels(:)          ! Contour level values
    type(color_t), allocatable :: colors(:)     ! Level colors
    type(colormap_t) :: cmap                    ! Colormap for filled contours
    type(normalization_t) :: norm               ! Value normalization

    ! Generated contour paths
    type(contour_path_t), allocatable :: paths(:)    ! All contour paths
    integer :: n_paths = 0

    ! Labeling system
    type(contour_label_t), allocatable :: labels(:)  ! Text labels
    logical :: label_inline = .true.            ! Inline vs separate labels
    character(len=20) :: label_format = '%g'    ! Printf-style format

    ! Rendering properties
    real(wp) :: linewidth = 1.0
    character(len=20) :: linestyle = '-'
    real(wp) :: alpha = 1.0
contains
    procedure :: generate_contours
    procedure :: add_labels
    procedure :: optimize_paths
    procedure :: get_level_color
end type

type :: contour_path_t
    real(wp), allocatable :: vertices(:,:)      ! Path vertices (2, n_vertices)
    integer, allocatable :: codes(:)            ! Path codes (MOVETO, LINETO, etc.)
    real(wp) :: level                           ! Contour level value
    logical :: is_closed = .false.              ! Closed vs open path
    type(bbox_t) :: bounds                      ! Path bounding box
end type

type :: contour_label_t
    character(len=32) :: text = ''              ! Label text
    real(wp) :: x, y                           ! Label position
    real(wp) :: rotation = 0.0                 ! Text rotation angle
    type(font_properties_t) :: font             ! Font properties
    logical :: inline = .true.                  ! Inline vs external label
end type

type :: level_generator_t
    integer :: max_levels = 10                  ! Maximum number of levels
    logical :: symmetric = .false.              ! Symmetric around zero
    real(wp) :: extend_factor = 0.05           ! Extension beyond data range
contains
    procedure :: auto_levels                    ! Automatic level calculation
    procedure :: nice_levels                    ! "Nice" level spacing
    procedure :: log_levels                     ! Logarithmic level spacing
end type

function auto_levels(generator, z_data, n_levels) result(levels)
    ! Implements matplotlib's MaxNLocator-style algorithm
    real(wp), intent(in) :: z_data(:,:)
    integer, intent(in), optional :: n_levels
    real(wp), allocatable :: levels(:)

    ! Calculate data range
    z_min = minval(z_data, mask=is_finite(z_data))
    z_max = maxval(z_data, mask=is_finite(z_data))

    if (present(n_levels)) then
        target_levels = n_levels
    else
        target_levels = generator%max_levels
    end if

    ! Generate "nice" level spacing
    range_val = z_max - z_min
    if (range_val == 0.0) then
        levels = [z_min]
        return
    end if

    ! Find optimal step size (matplotlib MaxNLocator algorithm)
    log_range = log10(range_val)
    step_size = 10.0_wp**(floor(log_range))

    ! Refine step size for target number of levels
    nice_steps = [1.0_wp, 2.0_wp, 2.5_wp, 5.0_wp, 10.0_wp]
    do i = 1, size(nice_steps)
        test_step = step_size * nice_steps(i)
        n_steps = ceiling(range_val / test_step)
        if (n_steps <= target_levels) then
            final_step = test_step
            exit
        end if
    end do

    ! Generate level array
    level_min = floor(z_min / final_step) * final_step
    level_max = ceil(z_max / final_step) * final_step

    n_actual = nint((level_max - level_min) / final_step) + 1
    allocate(levels(n_actual))

    do i = 1, n_actual
        levels(i) = level_min + (i - 1) * final_step
    end do

    ! Filter levels within data range
    levels = pack(levels, levels >= z_min .and. levels <= z_max)
end function

type :: contour_labeler_t
    real(wp) :: label_spacing = 10.0            ! Minimum spacing between labels (points)
    logical :: smart_placement = .true.         ! Use intelligent placement
    real(wp) :: straightness_threshold = 0.2    ! Path straightness for labeling
    integer :: max_labels_per_contour = 5       ! Limit labels per contour
contains
    procedure :: locate_label_positions
    procedure :: calculate_label_rotation
    procedure :: check_label_collision
    procedure :: format_label_text
end type

subroutine locate_label_positions(labeler, contour_path, label_positions)
    ! Implements matplotlib's locate_label algorithm
    type(contour_path_t), intent(in) :: contour_path
    type(label_position_t), allocatable :: label_positions(:)

    ! Find straight segments suitable for labeling
    call find_straight_segments(contour_path, segments)

    ! Score segments by suitability
    do i = 1, size(segments)
        segments(i)%score = calculate_segment_score(segments(i))
    end do

    ! Sort by score and select best positions
    call sort_segments_by_score(segments)

    ! Place labels with spacing constraints
    allocate(label_positions(0))
    last_position = -huge(1.0_wp)

    do i = 1, size(segments)
        if (segments(i)%score > labeler%straightness_threshold) then
            position = segments(i)%center_position
            if (abs(position - last_position) > labeler%label_spacing) then
                label_positions = [label_positions, create_label_position(segments(i))]
                last_position = position
                if (size(label_positions) >= labeler%max_labels_per_contour) exit
            end if
        end if
    end do
end subroutine

function calculate_segment_score(segment) result(score)
    ! Score based on straightness and length
    type(path_segment_t), intent(in) :: segment
    real(wp) :: score

    ! Calculate path curvature
    curvature = calculate_path_curvature(segment%vertices)

    ! Calculate segment length
    segment_length = calculate_path_length(segment%vertices)

    ! Combine factors (straight, long segments score higher)
    score = (1.0 - curvature) * min(segment_length / 20.0, 1.0)
end function

subroutine calculate_label_rotation(path_segment, rotation_angle)
    ! Calculate text rotation to follow contour direction
    real(wp), intent(in) :: path_segment(:,:)   ! Segment vertices
    real(wp), intent(out) :: rotation_angle

    ! Calculate tangent direction at segment center
    n_vertices = size(path_segment, 2)
    center_idx = n_vertices / 2

    if (center_idx > 1 .and. center_idx < n_vertices) then
        dx = path_segment(1, center_idx + 1) - path_segment(1, center_idx - 1)
        dy = path_segment(2, center_idx + 1) - path_segment(2, center_idx - 1)

        rotation_angle = atan2(dy, dx) * 180.0 / pi

        ! Keep text readable (never upside down)
        if (abs(rotation_angle) > 90.0) then
            rotation_angle = rotation_angle + sign(180.0_wp, -rotation_angle)
        end if
    else
        rotation_angle = 0.0
    end if
end subroutine

type :: path_optimizer_t
    real(wp) :: smoothing_factor = 0.1          ! Path smoothing strength
    real(wp) :: simplification_tolerance = 0.5  ! Douglas-Peucker tolerance
    logical :: remove_artifacts = .true.        ! Remove marching squares artifacts
contains
    procedure :: smooth_path
    procedure :: simplify_path
    procedure :: remove_marching_artifacts
end type

subroutine smooth_path(optimizer, path, smoothed_path)
    ! Apply smoothing to reduce marching squares artifacts
    type(contour_path_t), intent(in) :: path
    type(contour_path_t), intent(out) :: smoothed_path

    if (size(path%vertices, 2) < 3) then
        smoothed_path = path
        return
    end if

    ! Apply moving average smoothing
    allocate(smoothed_path%vertices(2, size(path%vertices, 2)))
    smoothed_path%codes = path%codes
    smoothed_path%level = path%level

    ! Smooth interior vertices (preserve endpoints)
    smoothed_path%vertices(:, 1) = path%vertices(:, 1)

    do i = 2, size(path%vertices, 2) - 1
        alpha = optimizer%smoothing_factor
        smoothed_path%vertices(:, i) = &
            (1.0 - 2.0*alpha) * path%vertices(:, i) + &
            alpha * (path%vertices(:, i-1) + path%vertices(:, i+1))
    end do

    smoothed_path%vertices(:, size(path%vertices, 2)) = &
        path%vertices(:, size(path%vertices, 2))
end subroutine

subroutine simplify_path(optimizer, path, simplified_path)
    ! Douglas-Peucker line simplification
    type(contour_path_t), intent(in) :: path
    type(contour_path_t), intent(out) :: simplified_path

    call douglas_peucker_simplify(path%vertices, optimizer%simplification_tolerance, &
                                 simplified_vertices, keep_indices)

    simplified_path%vertices = simplified_vertices
    simplified_path%codes = path%codes(keep_indices)
    simplified_path%level = path%level
    simplified_path%is_closed = path%is_closed
end subroutine

subroutine generate_contour_paths(cs, grid_x, grid_y, grid_z, levels)
    ! Enhanced version of existing marching squares algorithm
    type(contour_set_t), intent(inout) :: cs
    real(wp), intent(in) :: grid_x(:), grid_y(:), grid_z(:,:)
    real(wp), intent(in) :: levels(:)

    allocate(cs%paths(0))
    cs%n_paths = 0

    ! Generate contours for each level
    do level_idx = 1, size(levels)
        current_level = levels(level_idx)

        ! Trace contour using existing marching squares
        call trace_contour_level(grid_x, grid_y, grid_z, current_level, level_paths)

        ! Add paths to contour set
        do path_idx = 1, size(level_paths)
            level_paths(path_idx)%level = current_level
            cs%paths = [cs%paths, level_paths(path_idx)]
            cs%n_paths = cs%n_paths + 1
        end do
    end do

    ! Apply path optimization
    if (cs%optimize_paths_flag) then
        call cs%optimize_paths()
    end if
end subroutine

subroutine trace_contour_level(x, y, z, level, paths)
    ! Wrapper around existing marching squares implementation
    ! Converts fortplot's current algorithm to path objects

    ! Use existing marching squares algorithm
    call generate_contour_lines(x, y, z, level, line_segments)

    ! Convert line segments to connected paths
    call connect_line_segments(line_segments, connected_paths)

    ! Convert to path objects
    allocate(paths(size(connected_paths)))
    do i = 1, size(connected_paths)
        call convert_to_path_object(connected_paths(i), paths(i))
    end do
end subroutine

type :: colorbar_t
    type(contour_set_t), pointer :: parent_contours => null()
    type(axes_t), pointer :: axes => null()       ! Colorbar axes

    ! Appearance
    character(len=20) :: orientation = 'vertical'  ! 'vertical', 'horizontal'
    real(wp) :: fraction = 0.15                   ! Fraction of parent axes
    real(wp) :: pad = 0.05                        ! Padding from parent axes
    character(len=20) :: aspect = 'auto'          ! Aspect ratio

    ! Ticks and labels
    type(tick_locator_t) :: tick_locator
    type(tick_formatter_t) :: tick_formatter
    character(len=100) :: label = ''              ! Colorbar label
contains
    procedure :: create_colorbar_axes
    procedure :: draw_color_gradient
    procedure :: add_tick_labels
end type

subroutine create_colorbar_for_contours(contours, colorbar)
    type(contour_set_t), intent(in) :: contours
    type(colorbar_t), intent(out) :: colorbar

    ! Link to parent contours
    colorbar%parent_contours => contours

    ! Create colorbar axes
    call colorbar%create_colorbar_axes(contours%parent_axes)

    ! Draw color gradient
    call colorbar%draw_color_gradient(contours%levels, contours%colors)

    ! Add tick labels
    call colorbar%add_tick_labels(contours%levels)
end subroutine

! Enhanced API maintaining backward compatibility
subroutine add_contour_enhanced(fig, x, y, z, levels, colors, labels)
    type(figure_t), intent(inout) :: fig
    real(wp), intent(in) :: x(:), y(:), z(:,:)
    real(wp), intent(in), optional :: levels(:)
    type(color_t), intent(in), optional :: colors(:)
    logical, intent(in), optional :: labels

    type(contour_set_t) :: contour_set

    ! Generate levels if not provided
    if (present(levels)) then
        contour_set%levels = levels
    else
        contour_set%levels = auto_generate_levels(z, 10)
    end if

    ! Generate contour paths using existing algorithm
    call contour_set%generate_contours(x, y, z)

    ! Add labels if requested
    if (present(labels) .and. labels) then
        call contour_set%add_labels()
    end if

    ! Store in figure (maintain existing storage system)
    call fig%add_contour_set(contour_set)
end subroutine

! Level filtering for performance
function filter_levels_in_range(levels, z_min, z_max) result(filtered_levels)
    ! Only compute contours that intersect data range
    real(wp), intent(in) :: levels(:), z_min, z_max
    real(wp), allocatable :: filtered_levels(:)

    filtered_levels = pack(levels, levels >= z_min .and. levels <= z_max)
end function

! Path caching for re-rendering
type :: contour_cache_t
    real(wp), allocatable :: cached_levels(:)
    type(contour_path_t), allocatable :: cached_paths(:)
    logical :: is_valid = .false.
contains
    procedure :: invalidate, lookup_paths
end type