|
| 1 | +""" |
| 2 | + restrict(img[, region]) -> imgr |
| 3 | +
|
| 4 | +Reduce the size of `img` by approximately two-fold along the dimensions listed in |
| 5 | +`region`, or all spatial coordinates if `region` is not specified. |
| 6 | +The term `restrict` is taken from the coarsening operation of algebraic multigrid |
| 7 | +methods; it is the adjoint of "prolongation" (which is essentially interpolation). |
| 8 | +`restrict` anti-aliases the image as it goes, so is better than a naive summation |
| 9 | +over 2x2 blocks. |
| 10 | +The implementation of `restrict` has been tuned for performance, and should be |
| 11 | +a fast method for constructing [pyramids](https://en.wikipedia.org/wiki/Pyramid_(image_processing)). |
| 12 | +
|
| 13 | +If `l` is the size of `img` along a particular dimension, `restrict` produces an |
| 14 | +array of size `(l+1)÷2` for odd `l`, |
| 15 | +and `l÷2 + 1` for even `l`. See the example below for an explanation. |
| 16 | +
|
| 17 | +See also [`imresize`](@ref). |
| 18 | +
|
| 19 | +# Example |
| 20 | +
|
| 21 | +```julia |
| 22 | +a_course = [0, 1, 0.3] |
| 23 | +``` |
| 24 | +If we were to interpolate this at the halfway points, we'd get |
| 25 | +```julia |
| 26 | +a_fine = [0, 0.5, 1, 0.65, 0.3] |
| 27 | +``` |
| 28 | +Note that `a_fine` is obtained from `a_course` via the *prolongation* operator `P` as |
| 29 | +`P*a_course`, where |
| 30 | +```julia |
| 31 | +P = [1 0 0; # this line "copies over" the first point |
| 32 | + 0.5 0.5 0; # this line takes the mean of the first and second point |
| 33 | + 0 1 0; # copy the second point |
| 34 | + 0 0.5 0.5; # take the mean of the second and third |
| 35 | + 0 0 1] # copy the third |
| 36 | +``` |
| 37 | +`restrict` is the adjoint of prolongation. Consequently, |
| 38 | +```julia |
| 39 | +julia> restrict(a_fine) |
| 40 | +3-element Array{Float64,1}: |
| 41 | + 0.125 |
| 42 | + 0.7875 |
| 43 | + 0.3125 |
| 44 | +
|
| 45 | +julia> (P'*a_fine)/2 |
| 46 | +3-element Array{Float64,1}: |
| 47 | + 0.125 |
| 48 | + 0.7875 |
| 49 | + 0.3125 |
| 50 | +``` |
| 51 | +where the division by 2 approximately preserves the mean intensity of the input. |
| 52 | +
|
| 53 | +As we see here, for odd-length `a_fine`, restriction is the adjoint of interpolation at half-grid points. |
| 54 | +When `length(a_fine)` is even, restriction is the adjoint of interpolation at 1/4 and 3/4-grid points. |
| 55 | +This turns out to be the origin of the `l->l÷2 + 1` behavior. |
| 56 | +
|
| 57 | +One consequence of this definition is that the edges move towards zero: |
| 58 | +```julia |
| 59 | +julia> restrict(ones(11)) |
| 60 | +6-element Array{Float64,1}: |
| 61 | + 0.75 |
| 62 | + 1.0 |
| 63 | + 1.0 |
| 64 | + 1.0 |
| 65 | + 1.0 |
| 66 | + 0.75 |
| 67 | +``` |
| 68 | +In some applications (e.g., image registration), you may find it useful to trim the edges. |
| 69 | +""" |
| 70 | +restrict(img::AbstractArray, ::Tuple{}) = img |
| 71 | + |
| 72 | +restrict(A::AbstractArray, region::Vector{Int}) = restrict(A, (region...,)) |
| 73 | +restrict(A::AbstractArray) = restrict(A, coords_spatial(A)) |
| 74 | +function restrict(A::AbstractArray, region::Dims) |
| 75 | + restrict(restrict(A, region[1]), Base.tail(region)) |
| 76 | +end |
| 77 | + |
| 78 | +function restrict(A::AbstractArray{T,N}, dim::Integer) where {T,N} |
| 79 | + indsA = axes(A) |
| 80 | + newinds = map(UnitRange, ntuple(i->i==dim ? restrict_indices(indsA[dim]) : indsA[i], Val(N))) |
| 81 | + out = similar(Array{restrict_eltype(first(A)), N}, newinds) |
| 82 | + restrict!(out, A, dim) |
| 83 | + out |
| 84 | +end |
| 85 | + |
| 86 | +function restrict_eltype(A::AbstractArray) |
| 87 | + # infer the restrict_eltype on `eltype(A)` while preserving the container type |
| 88 | + # TODO: maybe there's more efficient way than the `similar` here.. |
| 89 | + typeof(similar(A, _restrict_eltype(eltype(A)), ntuple(_->1, ndims(A)))) |
| 90 | +end |
| 91 | +restrict_eltype(x) = _restrict_eltype(typeof(x)) |
| 92 | + |
| 93 | +for CT in (:AbstractGray, :AbstractRGB, :TransparentGray, :TransparentRGB) |
| 94 | + @eval _restrict_eltype(::Type{C}) where C<:$CT = __restrict_eltype(C) |
| 95 | +end |
| 96 | +_restrict_eltype(::Type{T}) where T = typeof(one(T)/4 + one(T)/2) |
| 97 | +_restrict_eltype(::Type{C}) where C<:Color = __restrict_eltype(RGB{eltype(C)}) |
| 98 | +_restrict_eltype(::Type{C}) where C<:Colorant = __restrict_eltype(ARGB{eltype(C)}) |
| 99 | +__restrict_eltype(::Type{C}) where C = base_colorant_type(C){promote_type(eltype(C), Float32)} |
| 100 | + |
| 101 | +function restrict!(out::AbstractArray{T,N}, A::AbstractArray, dim) where {T,N} |
| 102 | + if dim > N |
| 103 | + return copyto!(out, A) |
| 104 | + end |
| 105 | + indsout, indsA = axes(out), axes(A) |
| 106 | + ndims(out) == ndims(A) || throw(DimensionMismatch("input and output must have the same number of dimensions")) |
| 107 | + for d = 1:length(indsA) |
| 108 | + target = d==dim ? restrict_indices(indsA[d]) : indsA[d] |
| 109 | + indsout[d] == target || error("input and output must have corresponding indices; to be consistent with the input indices,\ndimension $d should be $target, got $(indsout[d])") |
| 110 | + end |
| 111 | + indspre, indspost = indsA[1:dim-1], indsA[dim+1:end] |
| 112 | + _restrict!(out, indsout[dim], A, indspre, indsA[dim], indspost) |
| 113 | +end |
| 114 | + |
| 115 | +@generated function _restrict!(out, indout, A, |
| 116 | + indspre::NTuple{Npre,AbstractUnitRange}, |
| 117 | + indA, |
| 118 | + indspost::NTuple{Npost,AbstractUnitRange}) where {Npre,Npost} |
| 119 | + Ipre = [Symbol(:ipre_, d) for d = 1:Npre] |
| 120 | + Ipost = [Symbol(:ipost_, d) for d = 1:Npost] |
| 121 | + quote |
| 122 | + $(Expr(:meta, :noinline)) |
| 123 | + T = eltype(out) |
| 124 | + if isodd(length(indA)) |
| 125 | + half = convert(eltype(T), 0.5) |
| 126 | + quarter = convert(eltype(T), 0.25) |
| 127 | + @nloops $Npost ipost d->indspost[d] begin |
| 128 | + iout = first(indout) |
| 129 | + @nloops $Npre ipre d->indspre[d] begin |
| 130 | + out[$(Ipre...), iout, $(Ipost...)] = zero(T) |
| 131 | + end |
| 132 | + ispeak = true |
| 133 | + for iA in indA |
| 134 | + if ispeak |
| 135 | + @inbounds @nloops $Npre ipre d->indspre[d] begin |
| 136 | + out[$(Ipre...), iout, $(Ipost...)] += |
| 137 | + half*convert(T, A[$(Ipre...), iA, $(Ipost...)]) |
| 138 | + end |
| 139 | + else |
| 140 | + @inbounds @nloops $Npre ipre d->indspre[d] begin |
| 141 | + tmp = quarter*convert(T, A[$(Ipre...), iA, $(Ipost...)]) |
| 142 | + out[$(Ipre...), iout, $(Ipost...)] += tmp |
| 143 | + out[$(Ipre...), iout+1, $(Ipost...)] = tmp |
| 144 | + end |
| 145 | + end |
| 146 | + ispeak = !ispeak |
| 147 | + iout += ispeak |
| 148 | + end |
| 149 | + end |
| 150 | + else |
| 151 | + threeeighths = convert(eltype(T), 0.375) |
| 152 | + oneeighth = convert(eltype(T), 0.125) |
| 153 | + z = zero(T) |
| 154 | + fill!(out, zero(T)) |
| 155 | + @nloops $Npost ipost d->indspost[d] begin |
| 156 | + peakfirst = true |
| 157 | + iout = first(indout) |
| 158 | + for iA in indA |
| 159 | + if peakfirst |
| 160 | + @inbounds @nloops $Npre ipre d->indspre[d] begin |
| 161 | + tmp = convert(T, A[$(Ipre...), iA, $(Ipost...)]) |
| 162 | + out[$(Ipre...), iout, $(Ipost...)] += threeeighths*tmp |
| 163 | + out[$(Ipre...), iout+1, $(Ipost...)] += oneeighth*tmp |
| 164 | + end |
| 165 | + else |
| 166 | + @inbounds @nloops $Npre ipre d->indspre[d] begin |
| 167 | + tmp = convert(T, A[$(Ipre...), iA, $(Ipost...)]) |
| 168 | + out[$(Ipre...), iout, $(Ipost...)] += oneeighth*tmp |
| 169 | + out[$(Ipre...), iout+1, $(Ipost...)] += threeeighths*tmp |
| 170 | + end |
| 171 | + end |
| 172 | + peakfirst = !peakfirst |
| 173 | + iout += peakfirst |
| 174 | + end |
| 175 | + end |
| 176 | + end |
| 177 | + out |
| 178 | + end |
| 179 | +end |
| 180 | + |
| 181 | +# If we're restricting along dimension 1, there are some additional efficiencies possible |
| 182 | +@generated function _restrict!(out, indout, A, ::NTuple{0,AbstractUnitRange}, |
| 183 | + indA, indspost::NTuple{Npost,AbstractUnitRange}) where Npost |
| 184 | + Ipost = [Symbol(:ipost_, d) for d = 1:Npost] |
| 185 | + quote |
| 186 | + $(Expr(:meta, :noinline)) |
| 187 | + T = eltype(out) |
| 188 | + if isodd(length(indA)) |
| 189 | + half = convert(eltype(T), 0.5) |
| 190 | + quarter = convert(eltype(T), 0.25) |
| 191 | + @inbounds @nloops $Npost ipost d->indspost[d] begin |
| 192 | + iout, iA = first(indout), first(indA) |
| 193 | + nxt = convert(T, A[iA+1, $(Ipost...)]) |
| 194 | + out[iout, $(Ipost...)] = half*convert(T, A[iA, $(Ipost...)]) + quarter*nxt |
| 195 | + for iA in first(indA)+2:2:last(indA)-2 |
| 196 | + prv = nxt |
| 197 | + nxt = convert(T, A[iA+1, $(Ipost...)]) |
| 198 | + out[iout+=1, $(Ipost...)] = quarter*(prv+nxt) + half*convert(T, A[iA, $(Ipost...)]) |
| 199 | + end |
| 200 | + out[iout+1, $(Ipost...)] = quarter*nxt + half*convert(T, A[last(indA), $(Ipost...)]) |
| 201 | + end |
| 202 | + else |
| 203 | + threeeighths = convert(eltype(T), 0.375) |
| 204 | + oneeighth = convert(eltype(T), 0.125) |
| 205 | + z = zero(T) |
| 206 | + @inbounds @nloops $Npost ipost d->indspost[d] begin |
| 207 | + c = d = z |
| 208 | + iA = first(indA) |
| 209 | + for iout = first(indout):last(indout)-1 |
| 210 | + a, b = c, d |
| 211 | + c, d = convert(T, A[iA, $(Ipost...)]), convert(T, A[iA+1, $(Ipost...)]) |
| 212 | + iA += 2 |
| 213 | + out[iout, $(Ipost...)] = oneeighth*(a+d) + threeeighths*(b+c) |
| 214 | + end |
| 215 | + out[last(indout), $(Ipost...)] = oneeighth*c + threeeighths*d |
| 216 | + end |
| 217 | + end |
| 218 | + out |
| 219 | + end |
| 220 | +end |
| 221 | + |
| 222 | +restrict_size(len::Integer) = isodd(len) ? (len+1)>>1 : (len>>1)+1 |
| 223 | +function restrict_indices(r::AbstractUnitRange) |
| 224 | + f, l = first(r), last(r) |
| 225 | + isodd(f) && return (f+1)>>1:restrict_size(l) |
| 226 | + f>>1 : (isodd(l) ? (l+1)>>1 : l>>1) |
| 227 | +end |
| 228 | +restrict_indices(r::Base.OneTo) = Base.OneTo(restrict_size(length(r))) |
| 229 | +function restrict_indices(r::UnitRange) |
| 230 | + f, l = first(r), last(r) |
| 231 | + isodd(f) && return (f+1)>>1:restrict_size(l) |
| 232 | + f>>1 : (isodd(l) ? (l+1)>>1 : l>>1) |
| 233 | +end |
0 commit comments