Merge pull request #494 from mronian/add-hist-matching

Adds Histogram Matching

docs/src/function_reference.md

@@ -285,6 +285,7 @@ findlocalminima
 imhist
 histeq
 adjust_gamma
+histmatch
 ```
 
 # Filtering kernels

src/Images.jl

@@ -227,6 +227,7 @@ export # types
     imhist,
     histeq,
     adjust_gamma,
+    histmatch,
     imlaplacian,
     imlineardiffusion,
     imlog,
@@ -292,7 +293,7 @@ Algorithms:
     - Resizing: `restrict`, `imresize` (not yet exported)
     - Filtering: `imfilter`, `imfilter_fft`, `imfilter_gaussian`, `imfilter_LoG`, `imROF`, `ncc`, `padarray`
     - Filtering kernels: `ando[345]`, `guassian2d`, `imaverage`, `imdog`, `imlaplacian`, `prewitt`, `sobel`
-    - Exposure : `imhist`, `histeq`, `adjust_gamma`
+    - Exposure : `imhist`, `histeq`, `adjust_gamma`, `histmatch`
     - Gradients: `backdiffx`, `backdiffy`, `forwarddiffx`, `forwarddiffy`, `imgradients`
     - Edge detection: `imedge`, `imgradients`, `thin_edges`, `magnitude`, `phase`, `magnitudephase`, `orientation`, `canny`
     - Corner detection: `imcorner`

src/algorithms.jl

@@ -1432,9 +1432,9 @@ 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)
+imhist{T<:Colorant}(img::AbstractArray{T}, nbins::Integer = 400) = imhist(convert(Array{Gray}, data(img)), nbins)
 
-function imhist{T<:Union{Gray,Number}}(img::AbstractArray{T}, nbins = 400)
+function imhist{T<:Union{Gray,Number}}(img::AbstractArray{T}, nbins::Integer = 400)
     minval = minfinite(img)
     maxval = maxfinite(img)
     imhist(img, nbins, minval, maxval)
@@ -1456,8 +1456,12 @@ maximum values present in the image are taken.
 `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})
+function imhist(img::AbstractArray, nbins::Integer, minval::Union{Gray,Real}, maxval::Union{Gray,Real})
     edges = StatsBase.histrange([Float64(minval), Float64(maxval)], nbins, :left)
+    imhist(img, edges)
+end
+
+function imhist(img::AbstractArray, edges::Range)
     histogram = zeros(Int, length(edges)+1)
     o = Base.Order.Forward
     G = graytype(eltype(img))
@@ -1473,7 +1477,6 @@ function imhist(img::AbstractArray, nbins, minval::Union{Gray,Real}, maxval::Uni
     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)
@@ -1576,6 +1579,47 @@ end
 
 adjust_gamma{T<:Number}(img::AbstractArray{T}, gamma::Number, minval::Number, maxval::Number) = map(i -> _gamma_pixel_rescale(i, gamma, minval, maxval), img)
 
+"""
+```
+hist_matched_img = histmatch(img, oimg, nbins)
+```
+
+Returns a grayscale histogram matched image with a granularity of `nbins` number of bins. `img` is the image to be 
+matched and `oimg` is the image having the desired histogram to be matched to. 
+
+"""
+histmatch(img::AbstractImage, oimg::AbstractArray, nbins::Integer = 400) = shareproperties(img, histmatch(data(img), oimg, nbins))
+
+_hist_match_pixel{T<:Union{Gray, Number}}(pixel::T, bins, lookup_table) = T(bins[lookup_table[searchsortedlast(bins, pixel)]])
+
+function _hist_match_pixel{T<:Color}(pixel::T, bins, lookup_table)
+    yiq = convert(YIQ, pixel)
+    y = _hist_match_pixel(yiq.y, bins, lookup_table)
+    convert(T, YIQ(y, yiq.i, yiq.q))
+end
+
+_hist_match_pixel{T<:TransparentColor}(pixel::T, bins, lookup_table) = base_colorant_type(T)(_hist_match_pixel(color(pixel), bins, lookup_table), alpha(pixel))
+
+function histmatch{T<:Colorant}(img::AbstractArray{T}, oimg::AbstractArray, nbins::Integer = 400)
+    el_gray = graytype(eltype(img))
+    oedges, ohist = imhist(oimg, nbins, zero(el_gray), one(el_gray))
+    _histmatch(img, oedges, ohist)
+end
+
+function _histmatch(img::AbstractArray, oedges::Range, ohist::AbstractArray{Int})
+    bins, histogram = imhist(img, oedges)
+    ohist[1] = 0
+    ohist[end] = 0
+    histogram[1] = 0 
+    histogram[end] = 0
+    cdf = cumsum(histogram)
+    cdf /= cdf[end]
+    ocdf = cumsum(ohist)
+    ocdf /= ocdf[end]
+    lookup_table = [indmin(abs(ocdf-val)) for val in cdf]
+    map(i -> _hist_match_pixel(i, bins, lookup_table), img)
+end
+
 # image gradients
 
 # forward and backward differences

test/algorithms.jl

@@ -229,6 +229,11 @@ facts("Algorithms") do
         @fact img == ret --> true
         @fact eltype(ret) == eltype(img) --> true
 
+        img = ones(Images.RGB{Images.U16}, 10, 10)
+        ret = histeq(img, 100)
+        @fact img == ret --> true
+        @fact eltype(ret) == eltype(img) --> true
+
         img = ones(Images.RGB{Float64}, 10, 10)
         ret = histeq(img, 100)
         @fact all(map((i, r) -> isapprox(i, r), img, ret)) --> true
@@ -314,7 +319,6 @@ facts("Algorithms") do
         img = ones(Images.ARGB{Images.U8}, 10, 10)
         ret = adjust_gamma(img, 1)
         @fact img == ret --> true
-        @fact eltype(ret) == eltype(img) --> true
 
         #Working
 
@@ -343,6 +347,50 @@ facts("Algorithms") do
         @fact alpha(r) --> alpha(b)
         @fact isapprox(r.val, b.val ^ 2) --> true
 
+        #Histogram Matching
+        #DataTypes
+        img = ones(Images.Gray{Float64}, 10, 10)
+        ret = histmatch(img, img)
+        @fact all(ret .== zero(eltype(img))) --> true
+        @fact eltype(ret) == eltype(img) --> true
+
+        img = ones(Images.Gray{Images.U8}, 10, 10)
+        ret = histmatch(img, img)
+        @fact all(ret .== zero(eltype(img))) --> true
+        @fact eltype(ret) == eltype(img) --> true
+
+        img = ones(Images.Gray{Images.U16}, 10, 10)
+        ret = histmatch(img, img)
+        @fact eltype(ret) == eltype(img) --> true
+
+        img = ones(Images.AGray{Images.U8}, 10, 10)
+        ret = histmatch(img, img)
+        @fact eltype(ret) == eltype(img) --> true
+
+        img = ones(Images.RGB{Images.U8}, 10, 10)
+        ret = histmatch(img, img)
+        @fact eltype(ret) == eltype(img) --> true
+
+        img = ones(Images.RGB{Images.U16}, 10, 10)
+        ret = histmatch(img, img)
+        @fact eltype(ret) == eltype(img) --> true
+
+        img = ones(Images.RGB{Float64}, 10, 10)
+        ret = histmatch(img, img)
+        @fact all(map((i, r) -> isapprox(zero(RGB), r, atol = 0.001), img, ret)) --> true
+        @fact eltype(ret) == eltype(img) --> true
+
+        img = ones(Images.ARGB{Images.U8}, 10, 10)
+        ret = histmatch(img, img)
+        @fact eltype(ret) == eltype(img) --> true
+
+        img = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
+        edges, hist = imhist(img, 2)
+        himg = Images._histmatch(img, edges, hist)
+        @fact himg == [0, 0, 0, 0, 0, 5, 5, 5, 5, 5] --> true
+        edges, hist = imhist(img, 5)
+        himg = Images._histmatch(img, edges, hist)
+        @fact himg == [0, 0, 2, 2, 4, 4, 6, 6, 8, 8] --> true
     end
 
     context("Array padding") do