pure function bezier(colors, levels) result(map)
integer, dimension(:,:), intent(in) :: colors
integer, intent(in), optional :: levels
integer, dimension(:,:), allocatable :: map
real(wp), dimension(:,:), allocatable :: map_r
integer :: order, i, j, levels_
real(wp) :: t
! 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_r(levels_,3), map(levels_,3)) ! 3 for RGB
do i = 1,levels_
t = real(i-1, wp) / real(levels_-1, wp)
map_r(i,:) = 0.0_wp
do j = 0, order
map_r(i,:) = map_r(i,:) + real(colors(j+1,:), wp)*&
real(factorial(order), wp)/(real(factorial(j), wp)*real(factorial(order-j), wp)) * t**j * (1.0_wp-t)**(order-j)
end do
map(i,1) = min(255, max(0, nint(map_r(i,1))))
map(i,2) = min(255, max(0, nint(map_r(i,2))))
map(i,3) = min(255, max(0, nint(map_r(i,3))))
end do
end function bezier