Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

A few changes to histogram functionality #516

Merged
merged 4 commits into from
Jul 19, 2016
Merged

A few changes to histogram functionality #516

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
9549935
Move hist code to avoid interrupting `restrict`
timholy Jul 16, 2016
3d136db
Faster imhist
timholy Jul 16, 2016
97f0e0d
Simplify imhist and support on-the-fly grayscale conversion
timholy Jul 16, 2016
dfce601
Simplify `histeq`
timholy Jul 19, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
260 changes: 119 additions & 141 deletions src/algorithms.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -228,6 +228,13 @@ function minfinite{T}(A::AbstractArray{T})
end
ret
end
function minfinite(f, A::AbstractArray)
ret = sentinel_min(typeof(f(first(A))))
for a in A
ret = minfinite_scalar(f(a), ret)
end
ret
end

"""
`m = maxfinite(A)` calculates the maximum value in `A`, ignoring any values that are not finite (Inf or NaN).
Expand All @@ -239,6 +246,13 @@ function maxfinite{T}(A::AbstractArray{T})
end
ret
end
function maxfinite(f, A::AbstractArray)
ret = sentinel_max(typeof(f(first(A))))
for a in A
ret = maxfinite_scalar(f(a), ret)
end
ret
end

"""
`m = maxabsfinite(A)` calculates the maximum absolute value in `A`, ignoring any values that are not finite (Inf or NaN).
Expand Down Expand Up @@ -437,6 +451,9 @@ accum(::Type{Float32}) = Float32
accum{T<:Real}(::Type{T}) = Float64
accum{C<:Colorant}(::Type{C}) = base_colorant_type(C){accum(eltype(C))}

graytype{T<:Number}(::Type{T}) = T
graytype{C<:AbstractGray}(::Type{C}) = C
graytype{C<:Colorant}(::Type{C}) = Gray{eltype(C)}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was really needed :)


# normalized by Array size
"`s = ssdn(A, B)` computes the sum-of-squared differences over arrays/images A and B, normalized by array size"
Expand Down Expand Up @@ -1199,144 +1216,6 @@ if isdefined(:restrict)
import Grid.restrict
end

imhist{T<:Colorant}(img::AbstractArray{T}, nbins=400) = imhist(convert(Array{Gray}, data(img)), nbins)
imhist{T<:Colorant}(img::AbstractArray{T}, nbins, minval, maxval) = imhist(convert(Array{Gray}, data(img)), nbins, minval, maxval)

function imhist{T<:Union{Gray,Number}}(img::AbstractArray{T}, nbins = 400)
minval = minfinite(img)
maxval = maxfinite(img)
imhist(img, nbins, minval, maxval)
end

"""
```
edges, count = imhist(img, nbins)
edges, count = imhist(img, nbins, minval, maxval)
```

Generates a histogram for the image over nbins spread between `(minval, maxval]`.
If `minval` and `maxval` are not given, then the minimum and
maximum values present in the image are taken.

`edges` is a vector that specifies how the range is divided;
`count[i+1]` is the number of values `x` that satisfy `edges[i] <= x < edges[i+1]`.
`count[1]` is the number satisfying `x < edges[1]`, and
`count[end]` is the number satisfying `x >= edges[end]`. Consequently,
`length(count) == length(edges)+1`.
"""
function imhist{T<:Union{Gray,Number}}(img::AbstractArray{T}, nbins, minval::T, maxval::T)
edges = StatsBase.histrange([Float64(minval), Float64(maxval)], nbins, :left)
histogram = zeros(Int, length(edges)+1)
for val in img
if val>=edges[end]
histogram[end] += 1
continue
end
index = searchsortedlast(edges, val)
histogram[index+1] += 1
end
edges, histogram
end

imhist{T<:Union{Gray,Number}}(img::AbstractArray{T}, nbins, minval, maxval) = imhist(img, nbins, convert(T, minval), convert(T, maxval))

Y_MIN = 16
Y_MAX = 235

function _prep_image_for_histeq{D<:Union{U8, U16}}(img::AbstractArray, dtype::Type{D})
img_shape = size(img)
img = convert(Array{base_colorant_type(eltype(img)){dtype}}, img)
img
end

histeq{T<:Colorant, D<:Union{U8, U16}}(img::AbstractArray{T}, nbins, dtype::Type{D} = U8) = histeq(_prep_image_for_histeq(img, dtype), nbins, zero(YCbCr), zero(YCbCr))
histeq{T<:Colorant, D<:Union{U8, U16}}(img::AbstractArray{T}, nbins, minval, maxval, dtype::Type{D} = U8) = histeq(_prep_image_for_histeq(img, dtype), nbins, minval, maxval)

function recompose_y(c::Color, eq_Y::Float32)
c_ycbcr = convert(YCbCr, color(c))
ret_c = YCbCr(eq_Y, c_ycbcr.cb, c_ycbcr.cr)
ret_c
end

function recompose_y(c::TransparentColor, eq_Y::Float32)
c_ycbcr = convert(YCbCrA, color(c))
ret_c = YCbCrA(eq_Y, c_ycbcr.cb, c_ycbcr.cr, c_ycbcr.alpha)
ret_c
end

function histeq{T<:Colorant}(img::AbstractArray{T}, nbins, minval, maxval)
Y = map(c -> convert(YCbCr, color(c)).y, img)
if maxval == zero(YCbCr{Float32})
eq_Y = histeq(Y, nbins, Y_MIN, Y_MAX)
else
eq_Y = histeq(Y, nbins, minval, maxval)
end
hist_equalised_img = map((eq, c) -> convert(T, recompose_y(c, eq)), eq_Y, img)
hist_equalised_img
end

function histeq{D<:Union{U8, U16}}(img::AbstractImage, nbins, dtype::Type{D} = U8)
hist_equalised_img = histeq(data(img), nbins, dtype)
hist_equalised_img = shareproperties(img, hist_equalised_img)
hist_equalised_img
end

function histeq{D<:Union{U8, U16}}(img::AbstractImage, nbins, minval, maxval, dtype::Type{D} = U8)
hist_equalised_img = histeq(data(img), nbins, minval, maxval, dtype)
hist_equalised_img = shareproperties(img, hist_equalised_img)
hist_equalised_img
end

function histeq{T<:TransparentGray}(img::AbstractArray{T}, args...)
opaque_img = convert(Array{base_color_type(eltype(img))}, img)
hist_equalised_img = histeq(opaque_img, args...)
converted_img = map((eq, o) -> convert(T, AGray(eq, alpha(o))), hist_equalised_img, img)
converted_img
end

histeq{T<:Gray, D<:Union{U8, U16}}(img::AbstractArray{T}, nbins, dtype::Type{D} = U8) = histeq(_prep_image_for_histeq(img, dtype), nbins, ColorVectorSpace.typemin(Gray{dtype}), ColorVectorSpace.typemax(Gray{dtype}))

function _histeq_pixel_rescale{T<:Union{Gray,Number}}(pixel::T, cdf, minval::T, maxval::T, cdf_length::Integer)
bin_pixel = max(1,Int(ceil((pixel-minval)*cdf_length/(maxval-minval))))
rescaled_pixel = minval + ((cdf[bin_pixel]-cdf[1])*(maxval-minval)/(cdf[end]-cdf[1]))
converted_pixel = convert(T, rescaled_pixel)
end

"""
```
hist_equalised_img = histeq(img, nbins, dtype = U8)
hist_equalised_img = histeq(img, nbins, minval, maxval, dtype = U8)
```

Returns a histogram equalised image with a granularity of approximately `nbins`
number of bins. An optional `dtype` argument (defaulting to U8) can be specified
to choose the number of bits of the returned image.

The `histeq` function can handle a variety of input types. The returned image depends
on the input type. If the input is an `Image` then the resulting image is of the same type
and has the same properties.

For coloured images, the input is converted to YCbCr type and the Y channel is equalised. This
is the combined with the Cb and Cr channels and the resulting image converted to the same type
as the input.

If minval and maxval are specified then only the intensities in the range
(minval, maxval) are equalised.

"""
function histeq{T<:Union{Gray,Number}}(img::AbstractArray{T}, nbins::Int, minval::T, maxval::T)
bins, histogram = imhist(img, nbins, minval, maxval)
cdf = cumsum(histogram[2:end-1])
img_shape = size(img)
cdf_length = size(cdf)[1]
hist_equalised_img = map(p -> _histeq_pixel_rescale(p, cdf, minval, maxval, cdf_length), img)
hist_equalised_img
end

histeq{T<:Number}(img::AbstractArray{T}, nbins) = histeq(img, nbins, minfinite(img), maxfinite(img))

histeq{T<:Union{Gray,Number}}(img::AbstractArray{T}, nbins, minval, maxval) = histeq(img, Int(nbins), convert(T, minval), convert(T, maxval))

"""
`imgr = restrict(img[, region])` performs two-fold reduction in size
along the dimensions listed in `region`, or all spatial coordinates if
Expand Down Expand Up @@ -1552,6 +1431,105 @@ This assumes the input `img` has intensities between 0 and 1.
"""
imstretch(img::AbstractArray, m::Number, slope::Number) = _imstretch(float(img), m, slope)


imhist{T<:Colorant}(img::AbstractArray{T}, nbins=400) = imhist(convert(Array{Gray}, data(img)), nbins)

function imhist{T<:Union{Gray,Number}}(img::AbstractArray{T}, nbins = 400)
minval = minfinite(img)
maxval = maxfinite(img)
imhist(img, nbins, minval, maxval)
end

"""
```
edges, count = imhist(img, nbins)
edges, count = imhist(img, nbins, minval, maxval)
```

Generates a histogram for the image over nbins spread between `(minval, maxval]`.
If `minval` and `maxval` are not given, then the minimum and
maximum values present in the image are taken.

`edges` is a vector that specifies how the range is divided;
`count[i+1]` is the number of values `x` that satisfy `edges[i] <= x < edges[i+1]`.
`count[1]` is the number satisfying `x < edges[1]`, and
`count[end]` is the number satisfying `x >= edges[end]`. Consequently,
`length(count) == length(edges)+1`.
"""
function imhist(img::AbstractArray, nbins, minval::Union{Gray,Real}, maxval::Union{Gray,Real})
edges = StatsBase.histrange([Float64(minval), Float64(maxval)], nbins, :left)
histogram = zeros(Int, length(edges)+1)
o = Base.Order.Forward
G = graytype(eltype(img))
for v in img
val = convert(G, v)
if val>=edges[end]
histogram[end] += 1
continue
end
index = searchsortedlast(edges, val, o)
histogram[index+1] += 1
end
edges, histogram
end


function _histeq_pixel_rescale{T<:Union{Gray,Number}}(pixel::T, cdf, minval, maxval)
n = length(cdf)
bin_pixel = clamp(ceil(Int, (pixel-minval)*length(cdf)/(maxval-minval)), 1, n)
rescaled_pixel = minval + ((cdf[bin_pixel]-cdf[1])*(maxval-minval)/(cdf[end]-cdf[1]))
convert(T, rescaled_pixel)
end
function _histeq_pixel_rescale{C<:Color}(pixel::C, cdf, minval, maxval)
yiq = convert(YIQ, pixel)
y = _histeq_pixel_rescale(yiq.y, cdf, minval, maxval)
convert(C, YIQ(y, yiq.i, yiq.q))
end
function _histeq_pixel_rescale{C<:TransparentColor}(pixel::C, cdf, minval, maxval)
base_colorant_type(C)(_histeq_pixel_rescale(color(pixel), cdf, minval, maxval), alpha(pixel))
end

"""
```
hist_equalised_img = histeq(img, nbins)
hist_equalised_img = histeq(img, nbins, minval, maxval)
```

Returns a histogram equalised image with a granularity of approximately `nbins`
number of bins.

The `histeq` function can handle a variety of input types. The returned image depends
on the input type. If the input is an `Image` then the resulting image is of the same type
and has the same properties.

For coloured images, the input is converted to YCbCr type and the Y channel is equalised. This
is the combined with the Cb and Cr channels and the resulting image converted to the same type
as the input.

If minval and maxval are specified then intensities are equalized to the range
(minval, maxval). The default values are 0 and 1.

"""
function histeq(img::AbstractArray, nbins::Integer, minval::Union{Number,Gray}, maxval::Union{Number,Gray})
bins, histogram = imhist(img, nbins, minval, maxval)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@timholy I was doing the CLAHE implementation and realised that in case of color images, the histogram calculated will be of a converted Gray image but the values considered while scaling would be the Y channel of YIQ colorspace. We should either abandon the YIQ conversion and use grayscale itself, or we should use something like my previous implementation of histeq which would pass the entire Y channel for histogram calculation.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

julia> c = rand(RGB{U8})
RGB{U8}(0.533,0.839,0.298)

julia> Gray(c)
Gray{U8}(0.686)

julia> YIQ(c)
YIQ{Float32}(0.68606275f0,-0.0083590355f0,-0.2330596f0)

That's why I switched to YIQ, the Y component is identical to Gray up to the Float32/U8 distinction. (Do you think that matters?)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh. If that it is identical then its fine. I was worried it was taking two different spaces and using them to calculate the rescaled pixel. Thanks!

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm glad you're paying such close attention, that's great!

cdf = cumsum(histogram[2:end-1])
img_shape = size(img)
minval == maxval && return map(identity, img)
hist_equalised_img = map(p -> _histeq_pixel_rescale(p, cdf, minval, maxval), img)
hist_equalised_img
end

function histeq(img::AbstractArray, nbins::Integer)
T = graytype(eltype(img))
histeq(img, nbins, zero(T), one(T))
end

function histeq(img::AbstractImage, nbins::Integer, minval::Union{Number,Gray}, maxval::Union{Number,Gray})
newimg = histeq(data(img), nbins, minval, maxval)
shareproperties(img, newimg)
end
histeq(img::AbstractImage, nbins::Integer) = shareproperties(img, histeq(data(img), nbins))

# image gradients

# forward and backward differences
Expand Down Expand Up @@ -1965,8 +1943,8 @@ integral_img = integral_image(img)
```

Returns the integral image of an image. The integral image is calculated by assigning
to each pixel the sum of all pixels above it and to its left, i.e. the rectangle from
(1, 1) to the pixel.
to each pixel the sum of all pixels above it and to its left, i.e. the rectangle from
(1, 1) to the pixel.
"""
function integral_image(img::AbstractArray)
integral_img = Array{accum(eltype(img))}(size(img))
Expand All @@ -1976,4 +1954,4 @@ function integral_image(img::AbstractArray)
cumsum!(integral_img, integral_img, sd[i])
end
integral_img
end
end
Empty file modified test/algorithms.jl
100755 → 100644
View file
Open in desktop
Empty file.
Empty file modified test/map.jl
100755 → 100644
View file
Open in desktop
Empty file.