bezier Function

    pure function bezier(colors, levels) result(map)
        integer, intent(in) :: colors(:,:)
        integer, intent(in), optional :: levels
        integer, allocatable :: map(:,:)
        integer :: order, i, j, levels_, fact_order
        real(wp) :: r, g, b, coeff, t, omt

        ! Set default value for levels
        if (present(levels)) then
            levels_ = levels
        else
            levels_ = 256
        end if

        ! Order of the Bezier curve
        order = size(colors, 1) - 1
        if (order < 1) error stop "Error: At least two control colors are required for Bezier interpolation."

        allocate(map(levels_, 3))
        fact_order = factorial(order)
        do i = 1, levels_
            t = real(i-1, wp) / real(levels_-1, wp)
            omt = 1.0_wp - t
            r = 0.0_wp
            g = 0.0_wp
            b = 0.0_wp
            do j = 0, order
                coeff = real(fact_order, wp)/real(factorial(j)*factorial(order-j), wp)*t**j*omt**(order-j)
                r = r + real(colors(j+1,1), wp) * coeff
                g = g + real(colors(j+1,2), wp) * coeff
                b = b + real(colors(j+1,3), wp) * coeff
            end do
            map(i,1) = min(255, max(0, nint(r)))
            map(i,2) = min(255, max(0, nint(g)))
            map(i,3) = min(255, max(0, nint(b)))
        end do
    end function bezier