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.

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

Histogram Equalisation #458

Merged
merged 3 commits into from
Jun 14, 2016
Merged

Histogram Equalisation #458

merged 3 commits into from
Jun 14, 2016

Conversation

mronian
Copy link
Contributor

@mronian mronian commented Mar 22, 2016

This commit adds a histogram function to help with histogram equalisation.

@mronian
Copy link
Contributor Author

mronian commented Mar 22, 2016

@bjarthur For Histogram Equalisation, I will need the maximum value for each datatype. For eg. in 8 bit Uint it would be 255 so I will have bins from 0-255 for the histogram. What is the best way to find this for all the types of Images?

raw(img) gives me the values in the original datatype. I'd like to get the values which Images.jl uses (between (0,1)) into a Float Array.

Update I did not know I could use the FixedPointNumbers datatype directly. It works now.

@coveralls
Copy link

Coverage Status

Coverage increased (+0.02%) to 75.342% when pulling 4b7c03c on mronian:add-hist-eq into cf1efb3 on timholy:master.

@timholy
Copy link
Member

timholy commented Mar 23, 2016

This was an excellent choice for a first pull request: it's a simple function, which gives you a chance to learn more about julia, Images, and the process of contributing to packages. 👍 for style.

I'll comment on the details in the code. (I was preparing a response to your previous version, and I noticed you just pushed a new version that addressed some of my concerns. Good timing!)

The histogram is calculated on the flattened image. For multi dimensional images, use the function
separately on each channel to obtain a histogram for each.
"""
function imhist{T}(img::AbstractArray{T,2}, nbins::Integer)
Copy link
Member

Choose a reason for hiding this comment

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

Is there a strong reason to limit this to 2d? Why couldn't I supply a 3d grayscale image?

Copy link
Member

Choose a reason for hiding this comment

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

What about an RGB image? What should this do in that case?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

No reason as such. I will extend it to N dims.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

For an RGB Image, there can be multiple ways. We could calculate the histogram for each channel or calculate the histogram of RGB values itself.

@timholy
Copy link
Member

timholy commented Mar 23, 2016

And now for a more general comment: do we need a separate function for this? Couldn't we just call hist?

@timholy
Copy link
Member

timholy commented Mar 23, 2016

(Note that if you try hist on an image with element type Gray, you'll get an error. I'll be happy to coach you through how to interpret the error, where to supply the fix, etc.)

@mronian
Copy link
Contributor Author

mronian commented Mar 23, 2016

Yes. hist was giving errors for the Datatypes involved so I thought of doing it this way. I would like to learn how to use it with Images.

I am getting this error for Gray datatype.

ERROR: MethodError: no method matching histrange(::Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}, ::Int64)
Closest candidates are:
  histrange{T<:AbstractFloat,N}(::AbstractArray{T<:AbstractFloat,N}, ::Integer)
  histrange{T<:Integer,N}(::AbstractArray{T<:Integer,N}, ::Integer)
 in hist(::Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}, ::Int64) at ./statistics.jl:732
 in eval(::Module, ::Any) at ./boot.jl:267

I'll work on the other comments and resubmit my PR.

Thanks a lot!

@timholy
Copy link
Member

timholy commented Mar 23, 2016

So, what that error means is that you might need to define a new histrange method. You do so like this:

Base.histrange{T}(A::AbstractArray{Gray{T}}, n) = ...

where ... is something I'm leaving for you to figure out. We should consider whether such a function goes here or in ColorVectorSpace.

One of the challenges we might face is whether to only compute histograms over [0,1], or whether it should be over [Gray(typemin(T)), Gray(typemax(T))] (what about T = Float64, though?), or possibly [minfinite(A), maxfinite(A)]. Thoughts?

@mronian
Copy link
Contributor Author

mronian commented Mar 23, 2016

I think we can add it to ColorVectorSpace as it has many other operations defined over the various Colors. I'll open an issue there and submit a PR.

I tried the hist(img[:], edg::AbstractVector) function. It works in this case. The other two ones where number of bins are given or just the data is given don't work due to the absence of a histrange function over ColorTypes. The hist(A, n) function anyways gives a range over [minfinite(A), maxfinite(A)] so that case is covered if we add the missing histrange function. For histeq we should use the hist(A, edg::AbstractVector) function with bins between [0,1] i.e. (1/n:1/n:1) where n is the number of bins desired. The fact that all Images are expressed using ColorTypes makes it even more easier as all Images are between [0,1].

@mronian
Copy link
Contributor Author

mronian commented Mar 23, 2016

I was going through Color.jl ColorTypes.jl and came to know that apart from RGB and XYZ, the others dont support the FixedPointNumbers package. This means we will have to take the histogram from ColorType(typemin(T)) to ColorType(typemax(T)) depending on the colorspace. Am I Right?

@mronian
Copy link
Contributor Author

mronian commented Mar 23, 2016

Submitted PR JuliaGraphics/ColorVectorSpace.jl#33.

@mronian
Copy link
Contributor Author

mronian commented Mar 24, 2016

PR JuliaGraphics/ColorVectorSpace.jl#35 adds zero() function for initialising cdf.

@mronian mronian force-pushed the add-hist-eq branch 2 times, most recently from 1eaf059 to 9b1f773 Compare March 24, 2016 13:40
@coveralls
Copy link

Coverage Status

Coverage increased (+0.2%) to 75.514% when pulling 9b1f773 on mronian:add-hist-eq into cf1efb3 on timholy:master.

@timholy
Copy link
Member

timholy commented Mar 24, 2016

I was going through Color.jl ColorTypes.jl and came to know that apart from RGB and XYZ, the others dont support the FixedPointNumbers package.

Most types require Fractional, which is defined here. So they take fractional.

The problem with typemin and typemax is they give Inf for things like Float32 and Float64. That probably won't be useful.

Here are some considerations: RGB{U8} doesn't allow ranges outside of 0, 1, but RGB{Float32} does. This is sometimes useful if, for example, you want to accumulate a value (e.g., in taking a mean). On the other hand, it's not meaningful to have values outside the range 0, 1, because they won't be displayed any differently than their clamped values. So on balance I favor restricting the range to 0, 1 for RGB.

For other color types, you're right that the range might need to be different. Overall it seems best to use the gamut, but I'm not always sure how you find out what that is in all cases. (It's pretty obvious in some, but in others less so.) I suppose one option is to define histrange only for those cases where the answer is clear, and leave it undefined for the rest.

```

"""
function histeq{T,N}(img::AbstractArray{T,N}, nbins::Int)
Copy link
Member

Choose a reason for hiding this comment

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

Might want histeq{T<:Union{Real,Gray},N} to restrict to grayscale?

Copy link
Member

Choose a reason for hiding this comment

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

Oh, I see, you're doing it by conversion to grayscale, so never mind. To avoid creating too many intermediaries, you might want to dispatch to two functions:

histeq{T<:Real,N}(...
    # this one calls maxfinite and minfinite, and does the scaling
    histeq_rescale(imgv)
end

histeq{T<:Colorant,N}(...
    imgv = [convert(Gray, p) for p in img]
    histeq_rescale(imgv)
end

or something.

@mronian
Copy link
Contributor Author

mronian commented Mar 25, 2016

I am getting this error sometimes if I use the julia notebook to display images returned by the histeq function.

WARNING: MethodError(Images.mapinfo_writemime_,(
Gray Images.Image with:
  data: 256x256 Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}
  properties:
    timedim: 0
    spatialorder:  y x
    pixelspacing:  1.0 1.0,))
 in #mapinfo_writemime#67(::Int64, ::Any, ::Images.Image{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2,Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}}) at /home/jarvis/.julia/v0.5/Images/src/writemime.jl:40
 in #save_#10(::Array{Any,1}, ::ImageMagick.#save_, ::FileIO.Stream{FileIO.DataFormat{:PNG},Base64EncodePipe}, ::Images.Image{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2,Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}}) at /home/jarvis/.julia/v0.5/ImageMagick/src/ImageMagick.jl:151
 [inlined code] from ./boot.jl:331
 in #save#6(::Array{Any,1}, ::Any, ::FileIO.Stream{FileIO.DataFormat{:PNG},Base64EncodePipe}, ::Images.Image{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2,Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}}, ::Vararg{Images.Image{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2,Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}}}) at /home/jarvis/.julia/v0.5/ImageMagick/src/ImageMagick.jl:67
 [inlined code] from ./boot.jl:331
 in #save#18(::Array{Any,1}, ::Any, ::FileIO.Stream{FileIO.DataFormat{:PNG},Base64EncodePipe}, ::Images.Image{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2,Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}}, ::Vararg{Images.Image{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2,Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}}}) at /home/jarvis/.julia/v0.5/FileIO/src/loadsave.jl:95
 [inlined code] from ./boot.jl:331
 in #writemime#64(::Images.MapNone{ColorTypes.RGB{FixedPointNumbers.UFixed{UInt8,8}}}, ::Int64, ::Int64, ::Any, ::Base64EncodePipe, ::MIME{symbol("image/png")}, ::Images.Image{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2,Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}}) at /home/jarvis/.julia/v0.5/Images/src/writemime.jl:30
 in base64encode(::Any, ::MIME{symbol("image/png")}, ::Vararg{Any}) at ./base64.jl:159
 in stringmime(::MIME{symbol("image/png")}, ::Images.Image{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2,Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}}) at ./multimedia.jl:84
 in display_dict(::Images.Image{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2,Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}}) at /home/jarvis/.julia/v0.5/IJulia/src/execute_request.jl:32
 in execute_request_0x535c5df2(::ZMQ.Socket, ::IJulia.Msg) at /home/jarvis/.julia/v0.5/IJulia/src/execute_request.jl:212
 [inlined code] from ./dict.jl:766
 in eventloop(::ZMQ.Socket) at /home/jarvis/.julia/v0.5/IJulia/src/IJulia.jl:141
 in (::IJulia.##24#30)() at ./task.jl:431
MethodError: no method matching mapinfo_writemime_(::Images.Image{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2,Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}})

@mronian
Copy link
Contributor Author

mronian commented Mar 25, 2016

There is a problem with scaling the image to values between [0-1] for the histogram. If there are intensities equal to 0, the histogram function will miss out on the 0 if we use the range 0:1/nbins:1. I feel it would be better to consider the bin centres instead i.e. take the range from 1/(2nbins):1/nbins:(1+1/(2nbins)). Is there a better way to do this?

@timholy
Copy link
Member

timholy commented Mar 28, 2016

Bin centers seem like a natural way to solve this. If you like that strategy, go for it.

If you prefer edge-based, then (given your observation) I have no objection to the "code duplication" of writing your own histogram algorithm. Your observation here is convincing incentive. So I'd say choose whichever path you like better.

colons = ntuple(d->Colon(), ndims(img))
Y = separate(YCbCr )[colons..., 1]
if maxval == zero(YCbCr)
eq_Y = histeq(Y, nbins, minfinite(Y), maxfinite(Y))
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Instead of minfinite and maxfinite, I'd like to take the entire range of the Y Channel possible. This depends on the type of conversion used though, so I wanted to check with you. In this it uses the range (16,235) for the Y Channel in 8 bits. Is this the correct one?

Copy link
Member

Choose a reason for hiding this comment

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

Best way to check would be to look at the convert methods for YCbCr in Colors.jl. This could surely use documentation, if it isn't documented already.

@mronian
Copy link
Contributor Author

mronian commented Jun 3, 2016

Finally Done! It works very well with multiple types of Images and ranges. The YCbCr conversion method used is the one from wikipedia. One slight change I have made is that while calculating the corrected intensity, it takes hist_length instead of nbin since the histogram calculated by StatsBase doesnt always return an histogram of length exactly nbins which was causing out-of-bounds exceptions.

@mronian
Copy link
Contributor Author

mronian commented Jun 3, 2016

Need to add some tests now. How are things like this which are very subjective and depend on the approach used usually tested? I mean, if someone were to modify the function a bit, the values returned would be different. Should I put tests to check if the return types etc are correct instead of testing the data returned?

hist_equalised_img = histeq(img, nbins)
hist_equalised_img = histeq(img, nbins, minval, maxval)
```
ColorVectorSpace.
Copy link
Member

Choose a reason for hiding this comment

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

What is this line for?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Thats crept in by mistake when I was switching between the terminal and my editor. Updating now.

@timholy
Copy link
Member

timholy commented Jun 4, 2016

Looking very good!

@mronian
Copy link
Contributor Author

mronian commented Jun 4, 2016

@timholy I forgot about Transparent Images. Is there a way to strip the entire alpha channel of a Colorant and then rejoin it after equalisation is done? I know about the alpha function but I think it can only be used on an element. Also, is there a way to know the non-transparent version of a colorspace given a TransparentColor element? This will help to handle all the different types with one function.

EDIT : Found base_color_type. It works beautifully !

I still need something for the alpha conversion though.

@coveralls
Copy link

Coverage Status

Coverage decreased (-2.4%) to 73.369% when pulling 6cd4793 on mronian:add-hist-eq into f50b32b on timholy:master.

@timholy
Copy link
Member

timholy commented Jun 4, 2016

Actually, in general it looks like a better pattern would be

Y = [convert(YCbCr, color(c)).y for c in img]   # this strips out the alpha
eq_Y = ...
histeq_img = reshape([recompose_y(c, eq_Y[i]) for (i,c) in enumerate(img)], size(img))

and define separate recompose_y methods for transparent vs opaque colors. This way you only end up creating three arrays: the output array, one temporary Y, and one temporary eq_Y. This might make up for the fact that recompose_y probably has to compute the YCbCr conversion all over again. (From a performance perspective, in principle the convert(YCbCr, color(c)).y might allow the compiler to skip calculation of the other channels of YCbCr, so this may not be as much redundant work as it seems. Having that work in practice may depend on inlining, however.)

@timholy
Copy link
Member

timholy commented Jun 4, 2016

In general, lots of "convert this array to that array" operations are characteristic of high-level languages that don't have good performance: it's how you "batch" work so that the (fast) underlying C code can do the real work. In julia, there's no need for that: focus on designing a good API for handling individual elements, chain all of your operations together on the elements, and then rebuild arrays as the "outermost" step.

@mronian
Copy link
Contributor Author

mronian commented Jun 6, 2016

@timholy I've made the changes. Please Review :)

@coveralls
Copy link

Coverage Status

Coverage decreased (-1.5%) to 74.259% when pulling 724735f on mronian:add-hist-eq into f50b32b on timholy:master.

@coveralls
Copy link

Coverage Status

Coverage decreased (-1.4%) to 74.292% when pulling 724735f on mronian:add-hist-eq into f50b32b on timholy:master.

@mronian mronian mentioned this pull request Jun 8, 2016
4 tasks
@timholy
Copy link
Member

timholy commented Jun 14, 2016

Sorry for the long delay, the gigantic JuliaLang/julia#16260 was all-consuming.

This looks good, I'm merging. At some point we should add some tests, though.

@timholy timholy merged commit 7581b75 into JuliaImages:master Jun 14, 2016
@mronian mronian mentioned this pull request Jun 14, 2016
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants