|
| 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 | + |
0 commit comments