Skip to content

Commit f12ee3a

Browse files
zygmuntszpakzygmuntszpak
committed
Adds PiecewiseLinearStretching
1 parent d545b4e commit f12ee3a

File tree

6 files changed

+358
-0
lines changed

6 files changed

+358
-0
lines changed

‎docs/Project.toml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
[deps]
22
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
33
FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
4+
ImageContrastAdjustment = "f332f351-ec65-5f6a-b3d1-319c6670881a"
45
ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
56
ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
67
ImageShow = "4e3cecfd-b093-5904-9786-8bbb286a6a31"

‎docs/src/reference.md‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,3 +48,8 @@ Matching
4848
```@docs
4949
MidwayEqualization
5050
```
51+
52+
### PiecewiseLinearStretching
53+
```@docs
54+
PiecewiseLinearStretching
55+
```

‎src/ImageContrastAdjustment.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ include("algorithms/contrast_stretching.jl")
2424
include("algorithms/gamma_correction.jl")
2525
include("algorithms/matching.jl")
2626
include("algorithms/midway_equalization.jl")
27+
include("algorithms/piecewise_linear_stretching.jl")
2728
include("compat.jl")
2829

2930
export
@@ -35,6 +36,7 @@ export
3536
GammaCorrection,
3637
LinearStretching,
3738
ContrastStretching,
39+
PiecewiseLinearStretching,
3840
build_histogram,
3941
adjust_histogram,
4042
adjust_histogram!

‎src/algorithms/piecewise_linear_stretching.jl‎

Lines changed: 179 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,179 @@
1+
"""
2+
```
3+
PiecewiseLinearStretching <: AbstractHistogramAdjustmentAlgorithm
4+
PiecewiseLinearStretching(; src_intervals = (0.0 => 0.1, 0.1 => 0.9, 0.9 => 1.0),
5+
dst_intervals = (0.0 => 0.2, 0.2 => 1.0, 1.0 => 1.0))
6+
7+
adjust_histogram([T,] img, f::PiecewiseLinearStretching)
8+
adjust_histogram!([out,] img, f::PiecewiseLinearStretching)
9+
```
10+
11+
Returns an image for which the intensity range was partitioned into subintervals (specified
12+
by `src_intervals`) and mapped linearly to a new set of subintervals
13+
(specified by `dst_intervals`).
14+
15+
# Details
16+
17+
Piecewise linear stretching is a contrast enhancing transformation that is used to modify
18+
the dynamic range of an image. The concept involves partitioning the intensity
19+
range of an image, ``I = [A,B]`` into a sequence of ``m`` disjoint intervals
20+
21+
```math
22+
I_1 = [A_1,A_2], I_2 = (A_2, A_3], \\ldots, I_{m} = (A_m,A_{m+1}]
23+
```
24+
which encompass the entire range of intensity values such that ``\\sum_j |I_j | = | I |``.
25+
One subsequently constructs a concomitant new set of disjoint intervals
26+
```math
27+
I_{1}^{′} = [a_1,a_2], I_{2}^{′} = (a_2, a_3], \\ldots, I_{m}^{′} = (a_m,a_{m+1}]
28+
```
29+
such that ``\\sum_j |I_{j}^{′} | = | I^{′} |``, where ``I^{′} = [a,b]``. Finally, one introduces a
30+
sequence of mappings ``f_j : I_j \\to I_{j}^{′}`` ``(j = 1, \\ldots, m)`` given by
31+
```math
32+
f_j(x) = (x-A_{j}) \\frac{a_{j+1}-a_j}{A_{j+1}-A_{j}} + a_j.
33+
```
34+
which linearly map intensities from the interval ``I_j`` to the interval ``I_{j}^{′}`` .
35+
36+
# Options
37+
38+
Various options for the parameters of the `adjust_histogram` and
39+
`PiecewiseLinearStretching` type are described in more detail below.
40+
41+
## Choices for `img`
42+
43+
The function can handle a variety of input types. The returned image depends on the input
44+
type.
45+
46+
For colored images, the input is converted to the [YIQ](https://en.wikipedia.org/wiki/YIQ)
47+
type and the intensities of the Y channel are stretched to the specified range. The modified
48+
Y channel is then combined with the I and Q channels and the resulting image converted to
49+
the same type as the input.
50+
51+
## Choices for `src_intervals`
52+
53+
The `src_intervals` must be specified as an `ntuple` of pairs. For example,
54+
`src_intervals = (0.0 => 0.2, 0.2 => 0.8, 0.8 => 1.0)` partitions the unit
55+
interval into three subintervals, [0.0, 0.2], (0.2, 0.8] and (0.8, 1.0].
56+
57+
To specify a single interval you must include a trailing comma, e.g.
58+
`src_intervals = (0.0 => 1.0,)`.
59+
60+
The start and end of a subinterval in `src_intervals` cannot be equal. For example,
61+
`src_intervals = (0.0 => 0.5, 0.5 => 0.5, 0.5 => 1.0)` will result in an `ArgumentError`
62+
because of the pair `0.5 => 0.5`.
63+
64+
## Choices for `dst_intervals`
65+
66+
The `dst_intervals` must be specified as an `ntuple` of pairs and there needs to be a
67+
corresponding pair for each pair in `src_intervals`.
68+
69+
Unlike `src_intervals`, the start and end of an interval in `src_intervals` can be equal.
70+
For example, the combination `src_intervals = (0.0 => 0.5, 0.5 => 1.0)` and `dst_intervals =
71+
(0.0 => 0.0 , 1.0 => 1.0)` will produce a binary image.
72+
73+
Care must be taken to ensure that the subintervals specified in `dst_intervals` are
74+
representable within the type of image on which the piecewise linear stretching will be
75+
applied. For example, if you specify a negative output interval like `dst_intervals = (1.0
76+
=> -1.0),` and try to apply the transformation on an image with type `Gray{N0f8}`, you will
77+
receive a DomainError.
78+
79+
## Saturating pixels
80+
81+
If the specified `src_intervals` don't span the entire range of intensities in an image,
82+
then intensities that fall outside `src_intervals` are automatically set to the boundary
83+
values of the `dst_intervals`. For example, given the combination `src_intervals = (0.1 =>
84+
0.9, )` and `dst_intervals = (0.0 => 1.0)`, intensities below `0.1` are set to zero, and
85+
intensities above `0.9` are set to one.
86+
87+
# Example
88+
89+
```julia
90+
using ImageContrastAdjustment, TestImages
91+
92+
img = testimage("mandril_gray")
93+
# Stretches the contrast in `img` so that it spans the unit interval.
94+
f = PiecewiseLinearStretching(src_intervals = (0.0 => 0.2, 0.2 => 0.8, 0.8 => 1.0),
95+
dst_intervals = (0.0 => 0.4, 0.4 => 0.9, 0.9 => 1.0))
96+
imgo = adjust_histogram(img, f)
97+
```
98+
99+
# References
100+
1. W. Burger and M. J. Burge. *Digital Image Processing*. Texts in Computer Science, 2016. [doi:10.1007/978-1-4471-6684-9](https://doi.org/10.1007/978-1-4471-6684-9)
101+
"""
102+
@with_kw struct PiecewiseLinearStretching{T} <: AbstractHistogramAdjustmentAlgorithm where T <: NTuple{N, Pair{<: Union{Real, AbstractGray}, <: Union{Real, AbstractGray}}} where {N}
103+
src_intervals::T = (0.0 => 0.1, 0.1 => 0.9, 0.9 => 1.0)
104+
dst_intervals::T = (0.0 => 0.2, 0.2 => 1.0, 1.0 => 1.0)
105+
function PiecewiseLinearStretching(src_intervals::T₁, dst_intervals::T₂) where {T₁ <: NTuple{N₁, Pair{<: Union{Real, AbstractGray}, <: Union{Real, AbstractGray}}} where {N₁}, T₂ <: NTuple{N₂, Pair{<: Union{Real, AbstractGray}, <: Union{Real, AbstractGray}}} where {N₂}}
106+
length(src_intervals) == length(dst_intervals) || throw(ArgumentError("The dst_intervals must define an interval pair for each consecutive src_interval pair, i.e. length(src_intervals) == length(dst_intervals)."))
107+
for src_pair in src_intervals
108+
if first(src_pair) == last(src_pair)
109+
throw(ArgumentError("The start and end of an interval in src_intervals cannot be equal."))
110+
end
111+
end
112+
T = promote_type(T₁, T₂)
113+
new{T}(src_intervals, dst_intervals)
114+
end
115+
end
116+
117+
function (f::PiecewiseLinearStretching)(out::GenericGrayImage, img::GenericGrayImage)
118+
# Verify that the specified dst_intervals fall inside the permissible range
119+
# of the output type.
120+
T = eltype(out)
121+
for (a,b) in f.dst_intervals
122+
p = typemin(T)
123+
r = typemax(T)
124+
if !(p <= a <= r) || !(p <= b <= r)
125+
throw(DomainError("The specified interval [$a, $b] falls outside the permissible range [$p, $r] associated with $T"))
126+
end
127+
end
128+
# A linear mapping of a value x in the interval [A, B] to the interval [a,b] is given by
129+
# the formula f(x) = (x-A) * ((b-a)/(B-A)) + a. The operation r * x - o is equivalent to
130+
# (x-A) * ((b-a)/(B-A)) + a. Precalculate the r and o so that later on an inner loop
131+
# contains only multiplication and addition to get better performance.
132+
N = length(f.src_intervals)
133+
rs = zeros(Float64, N)
134+
os = zeros(Float64, N)
135+
for (src, dst) in zip(pairs(f.src_intervals), pairs(f.dst_intervals))
136+
n, (A,B) = src
137+
_, (a,b) = dst
138+
rs[n] = (b - a) / (B - A)
139+
os[n] = (A*b - B*a) / (B - A)
140+
end
141+
142+
# Determine which interval the source pixel fall into, and apply the accompanying linear
143+
# transformation taking into account NaN values.
144+
@inbounds for p in eachindex(img)
145+
val = img[p]
146+
if isnan(val)
147+
out[p] = val
148+
else
149+
is_inside_src_intervals = false
150+
for (n, (A,B)) in pairs(f.src_intervals)
151+
if (A <= val < B) || (A <= val <= B && n == N)
152+
newval = rs[n] * val - os[n]
153+
out[p] = T <: Integer ? round(Int, newval) : newval
154+
is_inside_src_intervals = true
155+
break
156+
end
157+
end
158+
if !is_inside_src_intervals
159+
# Saturate pixels that fall outside the specified src_intervals by setting
160+
# them equal to the start/end of the edge dst_intervals.
161+
a = first(first(f.dst_intervals))
162+
b = last(last(f.dst_intervals))
163+
A = first(first(f.src_intervals))
164+
newval = val < A ? a : b
165+
out[p] = T <: Integer ? round(Int, newval) : newval
166+
end
167+
end
168+
end
169+
return out
170+
end
171+
172+
function (f::PiecewiseLinearStretching)(out::AbstractArray{<:Color3}, img::AbstractArray{<:Color3})
173+
T = eltype(out)
174+
yiq = convert.(YIQ, img)
175+
yiq_view = channelview(yiq)
176+
adjust_histogram!(view(yiq_view,1,:,:), f)
177+
out .= convert.(T, yiq)
178+
end
179+

‎test/piecewise_linear_stretching.jl‎

Lines changed: 170 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,170 @@
1+
@testset "Piecewise Linear Stretching" begin
2+
3+
@testset "constructors" begin
4+
@test_throws ArgumentError PiecewiseLinearStretching(
5+
src_intervals = ((0.0 => 0.2),(0.2 => 0.2), (0.8 => 1.0)),
6+
dst_intervals = ((0.0 => 0.2),(-1.0 => 0.0),(1 => -2)))
7+
8+
@test_throws ArgumentError PiecewiseLinearStretching(
9+
src_intervals = ((0.0 => 0.2),(0.2 => 0.8), (0.8 => 1.0)),
10+
dst_intervals = ((0.0 => 0.2), (1 => -2)))
11+
end
12+
@testset "miscellaneous" begin
13+
# We should be able to invert the image.
14+
f = PiecewiseLinearStretching(src_intervals = (0.0 => 1.0,), dst_intervals = (1.0 => 0.0,))
15+
for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
16+
img = T.(testimage("mandril_gray"))
17+
ret = adjust_histogram(img, f)
18+
@test ret ≈ (1 .- img)
19+
end
20+
21+
# We should be able to binarize the image.
22+
f = PiecewiseLinearStretching(
23+
src_intervals = ((0.0 => 0.5), (0.5 => 1.0)),
24+
dst_intervals = ((0.0 => 0.0),(1.0 => 1.0)))
25+
for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
26+
img = T.(testimage("mandril_gray"))
27+
28+
ret = adjust_histogram(img, f)
29+
vals = sort(unique(ret))
30+
@test first(vals) == 0.0
31+
@test last(vals) == 1.0
32+
end
33+
34+
# Setting the src_intervals equal to the dst_intervals should result in an identity
35+
# mapping (no change to the image) provided that the src_intervals span the entire
36+
# permissible contrast range, i.e. the unit interval.
37+
f = PiecewiseLinearStretching(
38+
src_intervals = ((0.0 => 0.3), (0.3 => 0.8), (0.8 => 1.0)),
39+
dst_intervals = ((0.0 => 0.3), (0.3 => 0.8), (0.8 => 1.0)))
40+
for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
41+
img = T.(testimage("mandril_gray"))
42+
ret = adjust_histogram(img, f)
43+
@test ret ≈ img
44+
end
45+
46+
# If the src_intervals equal the dst_intervals, but the src_intervals don't span the
47+
# entire permissible contrast range, i.e. the unit interval, then the resulting
48+
# image will not match the original image because pixels that fall outside the
49+
# specified src_intervals are satuared to the edge value of src_intervals.
50+
f = PiecewiseLinearStretching(
51+
src_intervals = ((0.3 => 0.8),),
52+
dst_intervals = ((0.3 => 0.8),))
53+
for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
54+
img = T.(testimage("mandril_gray"))
55+
ret = adjust_histogram(img, f)
56+
@test !isapprox(ret, 1 .- img, atol=1e-1)
57+
end
58+
59+
for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
60+
#=
61+
Stretching an image consisting of a linear ramp should not change the image
62+
if the specified minimum and maximum values match the minimum and maximum of
63+
the image.
64+
=#
65+
img = T.(collect(reshape(1/100:1/100:1, 10, 10)))
66+
minval = minimum(img)
67+
maxval = maximum(img)
68+
f = PiecewiseLinearStretching(src_intervals = (minval => maxval,), dst_intervals = (minval => maxval,))
69+
ret = adjust_histogram(img, f)
70+
if T <: Gray{Float32} || T <: Gray{Float64}
71+
@test all(map((i, r) -> isapprox(i, r), img, ret))
72+
else
73+
@test ret == img
74+
end
75+
76+
# Verify that NaN is also handled correctly.
77+
if T <: Gray{Float32} || T <: Gray{Float64}
78+
img[10] = NaN
79+
ret = adjust_histogram(img, f)
80+
@test isapprox(first(img), first(ret))
81+
@test isapprox(last(img), last(ret))
82+
@test isnan(ret[10])
83+
end
84+
85+
# Verify that the smallest and largest values match the specified minval and maxval.
86+
img = T.(collect(reshape(1/100:1/100:1, 10, 10)))
87+
minval = minimum(img)
88+
maxval = maximum(img)
89+
f = PiecewiseLinearStretching(src_intervals = (minval => maxval,), dst_intervals = (0.0 => 1.0,))
90+
ret = adjust_histogram(img, f)
91+
@test isapprox(0, first(ret))
92+
@test isapprox(1, last(ret))
93+
@test isapprox(0, minimum(ret[.!isnan.(ret)]))
94+
@test isapprox(1, maximum(ret[.!isnan.(ret)]))
95+
96+
# Verify that the return type matches the input type.
97+
img = T.(testimage("mandril_gray"))
98+
minval = minimum(img)
99+
maxval = maximum(img)
100+
f = PiecewiseLinearStretching(src_intervals = (minval => maxval,), dst_intervals = (0.0 => 1.0,))
101+
ret = adjust_histogram(img, f)
102+
@test eltype(ret) == eltype(img)
103+
@test isapprox(0, minimum(ret))
104+
@test isapprox(1, maximum(ret))
105+
106+
f = PiecewiseLinearStretching(src_intervals = (0.0 => 1.0,), dst_intervals = (0.2 => 0.8,))
107+
ret = adjust_histogram(img, f)
108+
@test eltype(ret) == eltype(img)
109+
# We are mapping the interval [0, 1] to [0.2, 0.8] but since the
110+
# smallest intensity in the mandril image need not be zero, and the largest
111+
# need to be one, smallest and largest observed values in the
112+
# image need not be 0.2 and 0.8, but should still be within those bounds.
113+
@test minimum(ret) >= 0.2
114+
@test maximum(ret) <= 0.8
115+
116+
# Verify that results are correctly clamped to [0.2, 0.9] if it exceeds the range
117+
f = PiecewiseLinearStretching(src_intervals = (0.3 => 0.8,),
118+
dst_intervals = (0.2 => 0.9,))
119+
ret = adjust_histogram(img, f)
120+
@test eltype(ret) == eltype(img)
121+
@test minimum(ret) == T(0.2)
122+
@test maximum(ret) == T(0.9)
123+
end
124+
125+
img = Float32.(collect(reshape(1/100:1/100:1, 10, 10)))
126+
127+
f = PiecewiseLinearStretching(src_intervals = ((0.0 => 1.0),), dst_intervals = ((-1.0 => 0.0),))
128+
ret = adjust_histogram(img, f)
129+
@test eltype(ret) == eltype(img)
130+
# @test minimum(ret) ≈ -1.0
131+
# @test maximum(ret) ≈ 0.0
132+
@test isapprox(minimum(ret), -1.0, atol=1e-1)
133+
@test isapprox(maximum(ret), 0.0, atol=1e-1)
134+
@test_throws DomainError adjust_histogram(Gray{N0f8}.(img), f)
135+
136+
img = Gray{Float32}.(collect(reshape(1/100:1/100:1, 10, 10)))
137+
f = PiecewiseLinearStretching(src_intervals = ((0.0 => 0.2), (0.2 => 0.8), (0.8 => 1.0)), dst_intervals = ((0.0 => 0.2), (-1.0 => 0.0), (1 => -2)))
138+
ret = adjust_histogram(img, f)
139+
@test eltype(ret) == eltype(img)
140+
@test minimum(ret) ≈ -2.0
141+
@test maximum(ret) ≈ 1.0
142+
143+
for T in (RGB{N0f8}, RGB{N0f16}, RGB{Float32}, RGB{Float64})
144+
#=
145+
Create a color image that spans a narrow graylevel range. Then
146+
quantize the 256 bins down to 32 and determine how many bins have
147+
non-zero counts.
148+
=#
149+
imgg = Gray{Float32}.([i/255.0 for i = 64:128, j = 1:10])
150+
img = colorview(RGB,imgg,imgg,imgg)
151+
img = T.(img)
152+
_, counts_before = build_histogram(img,32, minval = 0, maxval = 1)
153+
nonzero_before = sum(counts_before .!= 0)
154+
155+
#=
156+
Stretch the histogram. Then quantize the 256 bins down to 32 and
157+
verify that all 32 bins have non-zero counts. This will confirm that
158+
the dynamic range of the original image has been increased.
159+
=#
160+
f = PiecewiseLinearStretching(src_intervals = (0.3 => 0.5,),
161+
dst_intervals = (0.0 => 1.0,))
162+
ret = adjust_histogram(img, f)
163+
edges, counts_after = build_histogram(ret, 32, minval = 0, maxval = 1)
164+
nonzero_after = sum(counts_after .!= 0)
165+
@test nonzero_before < nonzero_after
166+
@test nonzero_after == 32
167+
@test eltype(img) == eltype(ret)
168+
end
169+
end
170+
end

‎test/runtests.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,4 +11,5 @@ using Test, ImageCore, ImageFiltering, TestImages, LinearAlgebra
1111
include("gamma_adjustment.jl")
1212
include("linear_stretching.jl")
1313
include("contrast_stretching.jl")
14+
include("piecewise_linear_stretching.jl")
1415
end

0 commit comments

Comments
 (0)