houghtransform.jl

using Images

"""
```
lines = hough_transform_standard(
    img_edges::AbstractMatrix;
    stepsize=1,
    angles=range(0,stop=pi,length=minimum(size(img))),
    vote_threshold=minimum(size(img)) / stepsize -1,
    max_linecount=typemax(Int))
```

Returns a vector of tuples corresponding to the tuples of (r,t)
where r and t are parameters for normal form of line:
    `x * cos(t) + y * sin(t) = r`

-   `r` = length of perpendicular from (1,1) to the line
-   `t` = angle between perpendicular from (1,1) to the line and x-axis

The lines are generated by applying hough transform on the image.

Parameters:
-    `img_edges`      = Image to be transformed (eltype should be `Bool`)
-    `stepsize`       = Discrete step size for perpendicular length of line
-    `angles`         = List of angles for which the transform is computed
-    `vote_threshold` = Accumulator threshold for line detection
-    `max_linecount`  = Maximum no of lines to return

# Example
```julia
julia> using ImageFeatures

julia> img = fill(false,5,5); img[3,:] .= true; img
5×5 Array{Bool,2}:
 false  false  false  false  false
 false  false  false  false  false
  true   true   true   true   true
 false  false  false  false  false
 false  false  false  false  false

julia> hough_transform_standard(img)
1-element Array{Tuple{Float64,Float64},1}:
 (3.0, 1.5707963267948966)
```
"""
function hough_transform_standard(
    img_edges::AbstractMatrix{Bool};
    stepsize=1,
    angles=range(0,stop=pi,length=minimum(size(img_edges))),
    vote_threshold=minimum(size(img_edges)) / stepsize -1,
    max_linecount=typemax(Int))

    stepsize > 0 || error("Discrete step size must be positive")
    ρ = stepsize
    θ = angles

    #function to compute local maximum lines with values > threshold and return a vector containing them
    function findlocalmaxima!(validLines::AbstractVector{CartesianIndex{2}}, accumulator_matrix::Array{Int,2}, threshold)
        for val in CartesianIndices(size(accumulator_matrix))
            if  accumulator_matrix[val] >  threshold                             &&
                accumulator_matrix[val] >  accumulator_matrix[val[1],val[2] - 1] &&
                accumulator_matrix[val] >= accumulator_matrix[val[1],val[2] + 1] &&
                accumulator_matrix[val] >  accumulator_matrix[val[1] - 1,val[2]] &&
                accumulator_matrix[val] >= accumulator_matrix[val[1] + 1,val[2]]
                push!(validLines,val)
            end
        end
    end

    indsy, indsx = axes(img_edges)
    ρinv = 1 / ρ
    numangle = length(θ)
    numrho = round(Int,(2(length(indsx) + length(indsy)) + 1)*ρinv)

    accumulator_matrix = zeros(Int, numangle + 2, numrho + 2)

    #Pre-Computed sines and cosines in tables
    sinθ, cosθ = sin.(θ).*ρinv, cos.(θ).*ρinv

    #Hough Transform implementation
    constadd = round(Int,(numrho -1)/2)
    for pix in CartesianIndices(size(img_edges))
        if img_edges[pix]
            for i in 1:numangle
                dist = round(Int, pix[1] * sinθ[i] + pix[2] * cosθ[i])
                dist += constadd
                accumulator_matrix[i + 1, dist + 1] += 1
            end
        end
    end

    #Finding local maximum lines
    validLines = Vector{CartesianIndex{2}}(undef, 0)
    findlocalmaxima!(validLines, accumulator_matrix, vote_threshold)

    #Sorting by value in accumulator_matrix
    @noinline sort_by_votes(validLines, accumulator_matrix) = sort!(validLines, lt = (a,b)-> accumulator_matrix[a]>accumulator_matrix[b])
    sort_by_votes(validLines, accumulator_matrix)

    max_linecount = min(max_linecount, length(validLines))

    lines = Vector{Tuple{Float64,Float64}}(undef, 0)

    #Getting lines with Maximum value in accumulator_matrix && size(lines) < max_linecount
    for l in 1:max_linecount
        lrho = ((validLines[l][2]-1) - (numrho-1)*0.5)*ρ
        langle = θ[validLines[l][1]-1]
        push!(lines,(lrho,langle))
    end

    lines
end

"""
```
circle_centers, circle_radius = hough_circle_gradient(img_edges, img_phase, radii; scale=1, min_dist=minimum(radii), vote_threshold)
```
Returns two vectors, corresponding to circle centers and radius.

The circles are generated using a hough transform variant in which a non-zero point only votes for circle
centers perpendicular to the local gradient. In case of concentric circles, only the largest circle is detected.

Parameters:
-   `img_edges`    = edges of the image
-   `img_phase`    = phase of the gradient image
-   `radii`        = circle radius range
-   `scale`        = relative accumulator resolution factor
-   `min_dist`     = minimum distance between detected circle centers
-   `vote_threshold`   = accumulator threshold for circle detection

[`canny`](@ref) and [`phase`](@ref) can be used for obtaining img_edges and img_phase respectively.

# Example

```julia
julia> using Images, ImageFeatures, FileIO, ImageView

julia> img = load(download("http://docs.opencv.org/3.1.0/water_coins.jpg"));

julia> img = Gray.(img);

julia> img_edges = canny(img, (Percentile(99), Percentile(80)));

julia> dx, dy=imgradients(img, KernelFactors.ando5);

julia> img_phase = phase(dx, dy);

julia> centers, radii = hough_circle_gradient(img_edges, img_phase, 20:30);

julia> img_demo = Float64.(img_edges); for c in centers img_demo[c] = 2; end

julia> imshow(img_demo)
```
"""
function hough_circle_gradient(
        img_edges::AbstractArray{Bool,2},
        img_phase::AbstractArray{<:Number,2},
        radii::AbstractRange{<:Integer};
        scale::Number=1,
        min_dist::Number=minimum(radii),
        vote_threshold::Number=minimum(radii)*min(scale, length(radii)))

    rows,cols=size(img_edges)

    non_zeros=CartesianIndex{2}[]
    centers=CartesianIndex{2}[]
    circle_centers=CartesianIndex{2}[]
    circle_radius=Int[]
    accumulator_matrix=zeros(Int, Int(floor(rows/scale))+1, Int(floor(cols/scale))+1)

    function vote!(accumulator_matrix, x, y)
        fx = Int(floor(x))
        fy = Int(floor(y))

        for i in fx:fx+1
            for j in fy:fy+1
                if checkbounds(Bool, accumulator_matrix, i, j)
                    @inbounds accumulator_matrix[i, j] += 1
                end
            end
        end
    end

    for j in axes(img_edges, 2)
        for i in axes(img_edges, 1)
            if img_edges[i,j]
                sinθ = -cos(img_phase[i,j]);
                cosθ = sin(img_phase[i,j]);

                for r in radii
                    x=(i+r*sinθ)/scale
                    y=(j+r*cosθ)/scale
                    vote!(accumulator_matrix, x, y)

                    x=(i-r*sinθ)/scale
                    y=(j-r*cosθ)/scale
                    vote!(accumulator_matrix, x, y)
                end
                push!(non_zeros, CartesianIndex{2}(i,j));
            end
        end
    end

    for i in findlocalmaxima(accumulator_matrix)
        if accumulator_matrix[i]>vote_threshold
            push!(centers, i);
        end
    end

    @noinline sort_by_votes(centers, accumulator_matrix) = sort!(centers, lt=(a, b) -> accumulator_matrix[a]>accumulator_matrix[b])

    sort_by_votes(centers, accumulator_matrix)

    dist(a, b) = sqrt(sum(abs2, (a-b).I))

    f = CartesianIndex(map(r->first(r), axes(accumulator_matrix)))
    l = CartesianIndex(map(r->last(r), axes(accumulator_matrix)))
    radius_accumulator=Vector{Int}(undef, Int(floor(dist(f,l)/scale)+1))

    for center in centers
        center=(center-1*_oneunit(center))*scale
        fill!(radius_accumulator, 0)

        too_close=false
        for circle_center in circle_centers
            if dist(center, circle_center)< min_dist
                too_close=true
                break
            end
        end
        if too_close
            continue;
        end

        for point in non_zeros
            r=Int(floor(dist(center, point)/scale))
            if radii.start/scale<=r<=radii.stop/scale
                radius_accumulator[r+1]+=1
            end
        end

        voters, radius = findmax(radius_accumulator)
        radius=(radius-1)*scale;

        if voters>vote_threshold
            push!(circle_centers, center)
            push!(circle_radius, radius)
        end
    end
    return circle_centers, circle_radius
end

@deprecate hough_circle_gradient(
        img_edges, img_phase,
        scale, min_dist,
        vote_threshold, radii) hough_circle_gradient(img_edges, img_phase, radii; scale=scale, min_dist=min_dist, vote_threshold=vote_threshold)

@deprecate hough_transform_standard(image, ρ, θ, threshold, linesMax) hough_transform_standard(
    image, stepsize=ρ, angles=θ, vote_threshold=threshold, max_linecount=linesMax)