init_piecewise2D Function

    type(piecewise2D) function init_piecewise2D(cell_width, cell_height, image)
        !! Initalise the piecewise2D type with a given cell_width, cell_height and input image

        !> Input cell width
        real(kind=wp), intent(in) :: cell_width
        !> Input cell height
        real(kind=wp), intent(in) :: cell_height
        !> Input image
        real(kind=wp), intent(in) :: image(:,:)
        
        real(kind=wp), allocatable :: HC1D(:), imagenew(:,:)
        
        integer :: width, height, w2, h2
        integer(kind=int64) :: i
        integer(kind=int32) :: x, y
        
        width = size(image, 1)
        height = size(image, 2)
        
        ! need to pad image for z-order to work...
        w2 = nextpwr2(width)
        h2 = nextpwr2(height)
        
        allocate(imagenew(w2, h2))
        imagenew = 0.
        
        init_piecewise2D%xoffset = (h2 - height)/2
        init_piecewise2D%yoffset = (w2 - width)/2
        
        imagenew(init_piecewise2D%xoffset:init_piecewise2D%xoffset+width-1, &
        init_piecewise2D%yoffset:init_piecewise2D%yoffset+height-1) = image
        
        allocate(init_piecewise2D%cdf(w2 * h2))
        allocate(HC1D(w2 * h2))
        
        HC1D = 0.
        
        do i = 0, (h2*w2)-1
            call decode(i, x, y)
            HC1D(i+1) = imagenew(x+1, y+1)
        end do
        
        init_piecewise2D%cdf(1) = HC1D(1)
        do i = 2, size(HC1D)
            init_piecewise2D%cdf(i) = init_piecewise2D%cdf(i-1) + HC1D(i)
        end do
        
        init_piecewise2D%cell_height = cell_height
        init_piecewise2D%cell_width = cell_width
        init_piecewise2D%cdf = init_piecewise2D%cdf/ init_piecewise2D%cdf(size(init_piecewise2D%cdf))
        
    end function init_piecewise2D