dis_otsu_thresholds Subroutine

  subroutine dis_otsu_thresholds(values, lower_bound, upper_bound, nbins, nclasses, thr, &
                                 nclasses_eff, degenerate)
    !================================================!
    !! Multi-class Otsu thresholding of a scalar distribution in
    !! [`lower_bound`, `upper_bound`].
    !! Histograms `values` into `nbins` equal-width bins over that FIXED range
    !! and finds the `nclasses`-1 ascending thresholds that maximise the between-class
    !! variance (exhaustive search, strict tie-break so the first-found maximum in
    !! ascending enumeration wins). Returns bin-centre thresholds in thr(1:nclasses_eff-1).
    !! The number of classes cannot exceed the number of populated bins (a cut
    !! must fall in a gap between clusters), so the effective class count is
    !! `nclasses_eff` = min(nclasses, populated bins); the caller is expected to
    !! note any reduction. Sets `degenerate` = .true. (thresholds undefined) only
    !! when fewer than three bins are populated.
    !================================================!
    implicit none

    real(kind=dp), intent(in) :: values(:)
    real(kind=dp), intent(in) :: lower_bound, upper_bound
    integer, intent(in) :: nbins, nclasses
    real(kind=dp), intent(out) :: thr(nclasses - 1)
    integer, intent(out) :: nclasses_eff
    logical, intent(out) :: degenerate

    integer :: n, i, k, npop, ncut, a, b
    integer, allocatable :: counts(:)
    real(kind=dp), allocatable :: pcum(:), scum(:), hmat(:, :)
    integer, allocatable :: cuts(:), best_cuts(:)
    real(kind=dp) :: d, denom, best_var

    degenerate = .false.
    nclasses_eff = 0
    thr = 0.0_dp

    n = size(values)
    d = (upper_bound - lower_bound)/real(nbins, dp)

    ! Fixed-range histogram; a value at upper_bound lands in the last bin.
    allocate (counts(0:nbins - 1))
    counts = 0
    do i = 1, n
      k = int((values(i) - lower_bound)/d)
      if (k < 0) k = 0
      if (k > nbins - 1) k = nbins - 1
      counts(k) = counts(k) + 1
    end do

    npop = count(counts > 0)
    if (npop < 3) then
      degenerate = .true.
      return
    end if

    ! Cannot resolve more classes than there are populated bins.
    nclasses_eff = min(nclasses, npop)
    ncut = nclasses_eff - 1

    ! Shortcut: as many classes as populated bins -> centres of the first
    ! nclasses_eff-1 populated bins.
    if (npop == nclasses_eff) then
      i = 0
      do k = 0, nbins - 1
        if (counts(k) > 0) then
          i = i + 1
          if (i > ncut) exit
          thr(i) = lower_bound + (real(k, dp) + 0.5_dp)*d
        end if
      end do
      return
    end if

    ! Cumulative count and index-moment (0-based bin index weighting), indexed
    ! from -1 so that segment [a,b] uses pcum(b)-pcum(a-1).
    allocate (pcum(-1:nbins - 1), scum(-1:nbins - 1))
    pcum(-1) = 0.0_dp
    scum(-1) = 0.0_dp
    do k = 0, nbins - 1
      pcum(k) = pcum(k - 1) + real(counts(k), dp)
      scum(k) = scum(k - 1) + real(k, dp)*real(counts(k), dp)
    end do

    ! Between-class variance table for every bin segment [a,b].
    allocate (hmat(0:nbins - 1, 0:nbins - 1))
    hmat = 0.0_dp
    do a = 0, nbins - 1
      do b = a, nbins - 1
        denom = pcum(b) - pcum(a - 1)
        if (denom > 0.0_dp) hmat(a, b) = (scum(b) - scum(a - 1))**2/denom
      end do
    end do

    allocate (cuts(ncut), best_cuts(ncut))
    best_var = -1.0_dp
    best_cuts = 0
    call search(1, 0)

    do i = 1, ncut
      thr(i) = lower_bound + (real(best_cuts(i), dp) + 0.5_dp)*d
    end do

  contains

    recursive subroutine search(level, lo)
      !! Enumerate ascending cut indices; `cuts(level)` is the last bin of the
      !! class below cut `level`. Objective = sum of between-class variance over
      !! the class segments; strict > keeps the first ascending maximum.
      integer, intent(in) :: level, lo
      integer :: idx, seg, aa, bb
      real(kind=dp) :: sigma

      if (level <= ncut) then
        do idx = lo, nbins - ncut + level - 2
          cuts(level) = idx
          call search(level + 1, idx + 1)
        end do
      else
        sigma = 0.0_dp
        aa = 0
        do seg = 1, ncut
          bb = cuts(seg)
          sigma = sigma + hmat(aa, bb)
          aa = cuts(seg) + 1
        end do
        sigma = sigma + hmat(aa, nbins - 1)
        if (sigma > best_var) then
          best_var = sigma
          best_cuts = cuts
        end if
      end if
    end subroutine search

  end subroutine dis_otsu_thresholds