#=----------------------------------------------------------------------------

Copyright (c) 2015 Peter Kovesi
pk@peterkovesi.com

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

----------------------------------------------------------------------------=#

import Images
export relief

#---------------------------------------------------------------------
"""
relief -  Generates relief shaded image

```
Usage:  shadeimg = relief(img, azimuth, elevation, gradscale, rgbimg)

Arguments: img - Image/heightmap to be relief shaded. ::Array{T<:Real,2}
       azimuth - Of light direction in degrees. Zero azimuth points
                 upwards and increases clockwise. Defaults to 45.
     elevation - Of light direction in degrees. Defaults to 45.
     gradscale - Scaling to apply to the surface gradients.  If the shading
                 is excessive decrease the scaling. Try successive doubling
                 or halving to find a good value.
        rgbimg - Optional RGB image to which the shading pattern derived
                 from 'img' is applied.   Alternatively, rgbimg can be a
                 colour map of type ::Array{ColorTypes.RGBA{Float64},1}
                 obtained from cmap().  This colour map is applied to the input
                 image/heightmap in order to obtain a RGB image to which
                 the shading pattern is applied.
```

Lambertian shading is used to form the relief image.  This obtained from the
cosine of the angle between the surface normal and light direction.  Note that
shadows are ignored.  Thus a small feature that might otherwise be in the
shadow of a nearby large structure is rendered as if the large feature was not
there.

See also: cmap, applycolourmap
"""
function relief(img::Array{T1,2}, az::Real=45, el::Real=45,
                gradscale::Real=1, rgbimg::Array{T2,3}=Array{Float64}(undef, 0,0,0)) where {T1<:Real,T2<:Real}
    # April    2014 Original MATLAB version
    # December 2015 Ported to Julia

    (rows, cols) = size(img)

    if !isempty(rgbimg)
        (r,c,chan) = size(rgbimg)
        if rows != r || cols != c || chan != 3
            error("rgbimg must be the same size as img")
        end
    end

    # Obtain surface normals of img
    loggrad = "lin"
    (n1, n2, n3) = surfacenormals(img, gradscale, loggrad)

    # Convert azimuth and elevation to a lighting direction vector.  Note that
    # the vector is constructed so that an azimuth of 0 points upwards and
    # increases clockwise.
    az = az/180*pi
    el = el/180*pi
    I = [cos(el)*sin(az), cos(el)*cos(az), sin(el)]
    I = I./norm(I)  # Ensure I is a unit vector

    # Generate Lambertian shading via the dot product between surface normals
    # and the light direction.  Note that the product with n2 is negated to
    # account for the image +ve y increasing downwards.
    shading = I[1]*n1 - I[2]*n2 + I[3]*n3

    # Remove -ve shading values which are generated by surface normals pointing
    # away from the light source.
    shading[shading .< 0] .= 0

    # If no RGB image has been supplied just return the raw shading image
    if isempty(rgbimg)
        shadeimg = shading

    else  #  Apply shading to the RGB image supplied
        shadeimg = zeros(size(rgbimg))

        for n = 1:3
            shadeimg[:,:,n]  = rgbimg[:,:,n].*shading
        end
    end

    # ** Resolve issue with RGB image being double or uint8

    return shadeimg
end


# Case where a colour map is supplied rather than a RGB img
function relief(img::Array{T1,2}, az::Real, el::Real,
                gradscale::Real,
                rgbmap::Array{ColorTypes.RGB{Float64},1}) where T1<:Real

    # Apply colour map to input image and use this as rgbimg
    rgbimg = applycolourmap(img, rgbmap)
    return relief(img, az, el, gradscale, rgbimg)
end


# Case where a colour map is supplied rather than a RGB img
function relief(img::Array{T1,2}, az::Real, el::Real,
                gradscale::Real,
                rgbmap::Array{ColorTypes.RGBA{Float64},1}) where T1<:Real

    # Apply colour map to input image and use this as rgbimg
    rgbimg = applycolourmap(img, rgbmap)
    return relief(img, az, el, gradscale, rgbimg)
end


#---------------------------------------------------------------------------
# Compute image/heightmap surface normals
#
# (n1, n2, n3) = surfacenormals(img, gradscale, loggrad)

function surfacenormals(img::Array{T,2}, gradscale::Real, loggrad::String) where T<:Real

    loggrad = lowercase(loggrad)

    # Compute partial derivatives of z.
    # p = dz/dx, q = dz/dy
    (q,p) = Images.imgradients(img, KernelFactors.ando3)
    p = p*gradscale
    q = q*gradscale

    # If specified take logs of gradient.
    # Note that taking the log of the surface gradients will produce a surface
    # that is not integrable (probably only of theoretical interest)
    if loggrad == "log"
        p = sign(p).*log1p(abs(p))
        q = sign(q).*log1p(abs(q))
    elseif loggrad == "loglog"
        p = sign(p).*log1p(log1p(abs(p)))
        q = sign(q).*log1p(log1p(abs(q)))
    elseif loggrad == "logloglog"
        p = sign(p).*log1p(log1p(log1p(abs(p))))
        q = sign(q).*log1p(log1p(log1p(abs(q))))
    end

    # Generate surface unit normal vectors. Components stored in n1, n2
    # and n3
    mag = sqrt.(1 .+ p.^2 .+ q.^2)
    n1 = -p./mag
    n2 = -q./mag
    n3 =  1 ./ mag

    return n1, n2, n3
end
