Commit 58f39e00 authored by Éric Thiébaut's avatar Éric Thiébaut
Browse files

Add an algorithm to find aligned boxes

parent e90617f1
#
# Geometry.jl --
#
# Methods for calibrating the geometry of the wavefront sensor.
#
#-------------------------------------------------------------------------------
#
# This file if part of the TAO software (https://github.com/emmt/TAO) licensed
# under the MIT license.
#
# Copyright (C) 2020, Éric Thiébaut.
#
module Geometry
using Printf
using TwoDimensional, ArrayTools, MappedArrays
"""
`Subtractable` is the union of real numbers that can be subtracted.
"""
const Subtractable = Union{Signed,AbstractFloat}
"""
```julia
findalignedboxes(img, msk, siz=26:40)
```
yields an array of aligned bounding boxes found in image `img` accroding to the
layout given by the true values in mask `msk`. Optional argument `siz` is to
specify a range of box sizes to try.
Set keyword `quiet=true` to suppress information messages.
Keyword `margin` can be used to specify an inner margin in pixels for the
boxes. A positive inner margin yields smaller boxes while a negative inner
margin yields large boxes. The default value is `margin=0`.
Keyword `plot` can be set with a reference to `PyPlot` module to plot results.
"""
function findalignedboxes(A::AbstractMatrix{<:Real},
mask::AbstractMatrix{Bool},
boxsizes::AbstractVector{<:Integer} = 26:40;
margin::Integer = 0,
quiet::Bool = false,
plot = nothing)
@assert length(boxsizes) > 0 && minimum(boxsizes) > 0
T = Float64
# Find number of sub-pupils in total and along each axis.
cnt = count(mask)
imin, imax, jmin, jmax = Tuple(BoundingBox(mask))
nx = imax - imin + 1
ny = jmax - jmin + 1
quiet || @printf("Mask has %d cells ∈ [%d:%d]×[%d:%d] = %d×%d cells\n",
cnt, imin, imax, jmin, jmax, nx, ny)
# Compute weights for the matching filters.
Wx = zeros(Float64, nx)
Wy = zeros(Float64, ny)
for j in jmin:jmax, i in imin:imax
if mask[i,j]
Wx[i - imin + 1] += 1
Wy[j - jmin + 1] += 1
end
end
# Integrate gradients along each axis. First integrate the image along
# a direction, second compute the gradient along the other direction
# --- Michel's idea ;-) ---. If carefully done, the last operation can
# be done in-place.
X, Y = axes(A)
xmin, xmax = axis_limits(X)
ymin, ymax = axis_limits(Y)
Gx = similar(Array{T}, X)
Gy = similar(Array{T}, Y)
fill!(Gx, 0)
@inbounds for y in ymin:ymax
s = zero(T)
@simd for x in xmin:xmax
a = T(A[x,y])
Gx[x] += a
s += a
end
Gy[y] = s
end
if plot !== nothing
plot.figure(2)
plot.clf()
plot.title("Integrated Profiles")
plot.plot(axes(Gx,1),Gx,label="1st dim.")
plot.plot(axes(Gy,1),Gy,label="2nd dim.")
plot.legend()
end
differentiate!(Gx)
differentiate!(Gy)
if plot !== nothing
plot.figure(3)
plot.clf()
plot.title("Derivatives of Integrated Profiles")
plot.plot(axes(Gx,1),Gx,label="1st dim.")
plot.plot(axes(Gy,1),Gy,label="2nd dim.")
plot.legend()
end
# Allocate workspaces for applying a matching filter on the gradients
# and store the list of the local maxima.
Qx = similar(Gx) # quality factors (as given by the matching filter)
Qy = similar(Gy) # quality factors (as given by the matching filter)
Sx = similar(Gx) # scratch array
Sy = similar(Gy) # scratch array
Ix = Int[]
Iy = Int[]
best_boxsize = 0
best_score = typemin(T)
best_xorigins = Array{Int}(undef, nx)
best_yorigins = Array{Int}(undef, ny)
for boxsize in of_eltype(Int, boxsizes)
matchopposite!(Qx, Gx, boxsize - 1)
matchopposite!(Qy, Gy, boxsize - 1)
findmaxima!(Ix, Qx, boxsize - 1, Sx)
findmaxima!(Iy, Qy, boxsize - 1, Sy)
if length(Ix) nx && length(Iy) ny
score1 = zero(T)
@inbounds for k in 1:nx
score1 += Wx[k]*Qx[Ix[k]]
end
score2 = zero(T)
@inbounds for k in 1:ny
score2 += Wy[k]*Qy[Iy[k]]
end
score = score1 + score2
if best_score < score
best_score = score
best_boxsize = boxsize
copyto!(best_xorigins, 1, Ix, 1, nx)
copyto!(best_yorigins, 1, Iy, 1, ny)
end
end
end
best_score > typemin(T) || error("algorithm failed: no solution found")
sort!(best_xorigins)
sort!(best_yorigins)
if !quiet
@printf("Best solution:\n - score: %.3e\n - boxes size: %d pixels\n",
best_score, best_boxsize)
print(" - X-origins: ")
print_values(best_xorigins, newline=true)
print(" - Y-origins: ")
print_values(best_yorigins, newline=true)
end
# Build list of bounding-boxes.
B = Array{BoundingBox{Int}}(undef, cnt)
ioff = imin - 1
joff = jmin - 1
k = 0
stp = best_boxsize - 1
mrg = Int(margin)
for j in 1:ny
y0 = best_yorigins[j]
y1 = y0 + stp
for i in 1:nx
if mask[i+ioff,j+joff]
k += 1
x0 = best_xorigins[i]
x1 = x0 + stp
B[k] = BoundingBox(x0,x1, y0,y1) - mrg
end
end
end
@assert k == cnt
if plot !== nothing
plot.figure(1)
plot.clf()
plot.title("Flat Image")
plot.imshow(permutedims(A), interpolation="nearest",
cmap="viridis", aspect="equal", origin="lower",
extent=(xmin-0.5,xmax+0.5,ymin-0.5,ymax+0.5))
plot.colorbar()
for i in eachindex(B)
x0, x1, y0, y1 = Tuple(B[i])
plot.plot([x0,x1,x1,x0,x0], [y0,y0,y1,y1,y0], "plum")
end
end
#return best_score, best_boxsize, best_xorigins, best_yorigins
return B
end
print_values(A::AbstractVector; kwds...) = print_values(stdout, A; kwds...)
function print_values(io::IO, A::AbstractVector; newline::Bool=false)
sep = ""
for x in A
print(io, sep, x)
sep = ", "
end
newline && println(io)
end
"""
```julia
differentiate(A) -> G
```
yields the discrete derivative of vector `A`. The result is a vector of
floating-point values of same size as `A`.
The output vector may be provided:
```julia
differentiate!(dst, src) -> dst
```
overwrites the contents of `dst` with the discrete derivative of vector `src`.
The operation can be done in-place, that is with `src` and `dst` the same
object. In that case, just call:
```julia
differentiate!(A) -> A
```
When the output array is provided and has integer elements, twice the discrete
derivatives are computed.
"""
differentiate(A::AbstractVector{T}) where {T<:Real} =
differentiate!(similar(Array{float(T)}, axes(A)), A)
function differentiate!(dst::AbstractVector{T},
src::AbstractVector{<:Real}) where {T<:Subtractable}
if dst === src
differentiate!(dst)
else
I = axes(dst,1)
axes(src,1) == I || throw(DimensionMismatch("arrays must have same axes"))
imin, imax = axis_limits(I)
if imax > imin
dst[imin] = diff1(T, src[imin+1], src[imin])
@simd for i in imin+1:imax-1
dst[i] = diff2(T, src[i+1], src[i-1])
end
dst[imax] = diff1(T, src[imax], src[imax-1])
else
fill!(dst, zero(T))
end
end
return dst
end
# In-place derivation along X.
function differentiate!(A::AbstractVector{T}) where {T<:Subtractable}
imin, imax = axis_limits(axes(A,1))
if imax > imin
a2, a3 = A[imin], A[imin]
A[imin] = diff1(T, a3, a2)
@simd for i in imin+1:imax-1
a1, a2, a3 = a2, a3, A[i+1]
A[i] = diff2(T, a3, a1)
end
A[imax] = diff1(T, a3, a2)
else
fill!(A, zero(T))
end
return A
end
@doc @doc(differentiate) differentiate!
@inline two(::Type{T}) where {T<:Real} = T(2)
@inline half(::Type{T}) where {T<:AbstractFloat} = one(T)/two(T)
# Finite difference for a step of 1 sample, T is the type of the result.
@inline diff1(::Type{T}, a, b) where {T<:AbstractFloat} = (T(a) - T(b))
@inline diff1(::Type{T}, a, b) where {T<:Signed} = (T(a) - T(b))*two(T)
# Finite difference for a step of 2 samples, T is the type of the result.
@inline diff2(::Type{T}, a, b) where {T<:AbstractFloat} = (T(a) - T(b))*half(T)
@inline diff2(::Type{T}, a, b) where {T<:Signed} = (T(a) - T(b))
"""
```julia
matchopposite(A, off) -> B
```
yields vector `B` whose elements are given by:
```
B[i] = A[i] - A[i+off] (∀i)
```
and assuming `A[i+off] = 0` if `i+off` is beyond limits.
This linear filter is useful to detect opposite features of a given size (here
`off+1`).
A destination array can be provided to avoid allocating memory:
```julia
matchopposite!(dst, src, off) -> dst
```
overwrites `dst` with the result. The operation can be applied *in-place*,
that is with `src` and `dst` being the same object.
"""
matchopposite(A::AbstractVector{T}, off::Int) where {T<:Subtractable} =
matchopposite!(similar(Array{T}, axes(A)), A, off)
function matchopposite!(dst::AbstractVector{<:Subtractable},
src::AbstractVector{<:Subtractable},
off::Int)
off 1 || throw(ArgumentError("offset must be at least 1"))
(I = axes(dst,1)) == axes(src,1) || throw_not_same_axes()
imin, imax = axis_limits(I)
@inbounds @simd for i in imin:imax-off
dst[i] = src[i] - src[i+off]
end
@inbounds @simd for i in imax-off+1:imax
dst[i] = src[i]
end
return dst
end
@doc @doc(matchopposite) matchopposite!
throw_not_same_axes() = throw(DimensionMismatch("arrays must have same axes"))
"""
```julia
findmaxima!(I, A, gap, ws=similar(A)) -> I
```
overwrites `I` with the list of the indices of the most significant local
maxima in vector `A`. Local maxima are detected in decreasing order of
significance and subsequent maxima must be at a distance greater than `gap`
from any previous ones.
Optional argument `ws` is a workspace array that may be provided to avoid
allocating memory. Argument `A` can be used workspace, but its contents
will be set to zero at all positions but those corresponding to the local
maxima found.
"""
function findmaxima!(I::AbstractVector{Int},
A::AbstractVector,
gap::Integer, args...)
return _find!(I, zero, identity, >, A, Int(gap), args...)
end
"""
```julia
findminima!(I, A, gap, ws=similar(A)) -> I
```
overwrites `I` with the list of the indices of the most significant local
minima in vector `A`. Local minima are detected in decreasing order of
significance and subsequent minima must be at a distance greater than `gap`
from any previous ones.
Optional argument `ws` is a workspace array that may be provided to avoid
allocating memory. Argument `A` can be used workspace, but its contents
will be set to zero at all positions but those corresponding to the local
minima found.
"""
function findminima!(I::AbstractVector{Int},
A::AbstractVector,
gap::Integer, args...)
return _find!(I, zero, identity, <, A, Int(gap), args...)
end
"""
```julia
findextrema!(I, A, gap, ws=similar(A)) -> I
```
overwrites `I` with the list of the indices of the most significant local
extrema in vector `A`. Local extrema are detected in decreasing order of
significance and subsequent extrema must be at a distance greater than
`gap` from any previous ones.
Optional argument `ws` is a workspace array that may be provided to avoid
allocating memory. Argument `A` can be used as workspace if its contents
is no longer needed.
"""
function findextrema!(I::AbstractVector{Int},
A::AbstractVector,
gap::Integer, args...)
return _find!(I, zero, abs, >, A, Int(gap), args...)
end
function _find!(I::AbstractVector{Int},
none::Function,
conv::Function,
comp::Function,
A::AbstractVector{T},
gap::Int,
ws::AbstractVector{T} = similar(A)) where {T}
gap 1 || throw(ArgumentError("offset must be at least 1"))
imin, imax = axis_limits(same_axes(A, ws)[1])
ws === A || copyto!(ws, A)
resize!(I, 0)
clr = none(T)::T
while true
# Find next most significant value.
flag = false
ibst = 0 # best index so far
vbst = clr # best value so far
@inbounds for i in imin:imax
val = conv(ws[i])
if comp(val, vbst)
ibst = i
vbst = val
flag = true
end
end
if ! flag
return I
end
# Register index.
push!(I, ibst)
# Clear neighborhood.
@inbounds @simd for i in max(imin,ibst-gap):min(imax,ibst+gap)
ws[i] = clr
end
end
end
end # module
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment