|
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) |
| 1 | +# restrict |
| 2 | +function restrict(img::Union{WarpedView, InvWarpedView}, dims::Integer) |
| 3 | + # preserve axes information but eagerly evaluate the transformation |
| 4 | + # TODO: this eager evaluation destroy the laziness of `WarpedView` |
| 5 | + # and `InvWarpedView`. |
| 6 | + restrict(OffsetArray(collect(img), axes(img)), dims) |
233 | 7 | end |
234 | 8 |
|
235 | 9 | # imresize |
|
0 commit comments