Skip to content

Instantly share code, notes, and snippets.

@andyferris
Last active June 25, 2016 16:40
Show Gist options
  • Save andyferris/9e35ebed4692f4a53a632a17f4608a92 to your computer and use it in GitHub Desktop.
Save andyferris/9e35ebed4692f4a53a632a17f4608a92 to your computer and use it in GitHub Desktop.
Possible image transformations design
module ImageTransformations
using CoordinateTransformations
import CoordinateTransformations.transform
export TransformedArray, ImageTransformation
abstract Filter
immutable NearestNeighbourFilter <: Filter; end
immutable LinearFilter <: Filter; end
immutable BoxFilter <: Filter; end
# Other filters might include Gaussian, sharpening or sampling filters, or even
# ones which account for non-linear colour spaces and so-forth.
"""
immutable TransformedArray{T, N, A <: AbstractArray, Trans <: AbstractTransformation, Filt <: Filter} <: AbstractArray{T, N}
TransformedArray(a::AbstractArray, trans, filter = LinearFilter())
A lazy view of an array which represents some geometric transformation of the
underlying data, such as a rotation of a 2D image. A transformation object
`trans` performs the corresponding point-by-point geometric transformation (e.g.
maps 2D points to 2D points). The `filter` specifies various sampling
techniques, such as the simple `NearestNeighbourFilter()`, the interpolating
`LinearFilter()` or the anti-aliasing `BoxFilter()` (not implemented yet).
"""
immutable TransformedArray{T, N, A <: AbstractArray, Trans <: AbstractTransformation, Filt <: Filter} <: AbstractArray{T, N}
data::A
tform::Trans
filter::Filter
# Maybe data about relative size and number of pixels in new image?
end
# Bare constructor, same T, N
TransformedArray(a::AbstractArray, tform, filter = LinearFilter()) = TransformedArray{eltype(a),ndims(a),typeof(a),typeof(tform),typeof(filter)}(a, tform, filter)
# Manually specify a different T or N
# Different T might arrise because of the filtering process. This might be quite
# useful.. the filter could convert between colour spaces or bit depths or whatever.
# Also, different N might occur for an slice (possibly diagonal) or an integrating filter (marginal projection).
# (It seems difficult to predict these changes automatically. )
(::Type{TransformedArray{T}}){T}(a::AbstractArray, tform, filter = LinearFilter()) = TransformedArray{T,ndims(a),typeof(a),typeof(tform),typeof(filter)}(a, tform, filter)
(::Type{TransformedArray{T,N}}){T,N}(a::AbstractArray, tform, filter = LinearFilter()) = TransformedArray{T,N,typeof(a),typeof(tform),typeof(filter)}(a, tform, filter)
Base.size(A::TransformedArray) = size(A.data)
function Base.getindex{T}(A::TransformedArray{T}, inds::Number...)
# Do we need to convert inds? E.g. Int to Float32?
x = transform(A.tform, Point(inds...))._ # Transform to-and-from Point (FixedSizeArrays)
@inbounds return convert(T, apply_filter(A.filter, A.data, x)) # Works for any inds, so inbounds by default...
end
Base.similar{T}(A::TransformedArray, ::Type{T}, dims::Dims) = Array{T}(dims)
function apply_filter(::NearestNeighbourFilter, A::AbstractArray, x::Tuple)
A[map(i -> (tmp = round(x[i]); tmp < 1 ? 1 : (tmp > size(A,i) ? size(A,i) : tmp)), 1:ndims(A))...]
end
# Linear interpolation
function apply_filter{IndT}(::LinearFilter, A::Vector, x::Tuple{IndT})
if x <= 1
return one(IndT) * A[1] # type-stability (maybe beter `convert(promote_type(IndT, eltype(A)), A[1])` ?)
elseif x >= length(A)
return one(IndT) * A[end] # type-stability
else
i1 = floor(x[1])
i2 = i1 + 1
c2 = x - i1
c1 = 1 - c2
return c1*A[i1] + c2*A[i2]
end
end
# Bi-linear interpolation
function apply_filter{IndT}(::LinearFilter, A::Matrix, xy::NTuple{2,IndT})
(x, y) = xy
if x <= 1
if y <= 1
return one(IndT) * A[1, 1]
elseif y >= size(A, 2)
return one(IndT) * A[1, end]
else
j1 = floor(y)
j2 = j1 + 1
c2y = y - j1
c1y = 1 - c2y
return c1y*A[1,j1] + c2y*A[1, j2]
end
elseif x >= size(A,1)
if y <= 1
return one(IndT) * A[end, 1]
elseif y >= size(A, 2)
return one(IndT) * A[end, end]
else
j1 = floor(y)
j2 = j1 + 1
c2y = y - j1
c1y = 1 - c2y
return c1y*A[end,j1] + c2y*A[end, j2]
end
else
i1 = floor(x)
i2 = i1 + 1
c2x = x - i1
c1x = 1 - c2x
if y <= 1
return c1x*A[i1, 1] + c2x*A[i2, 1]
elseif y >= size(A, 2)
return c1x*A[i1, end] + c2x*A[i2, end]
else
j1 = floor(y)
j2 = j1 + 1
c2y = y - j1
c1y = 1 - c2y
z1 = c1x*A[i1, j1] + c2x*A[i2, j1]
z2 = c1x*A[i1, j2] + c2x*A[i2, j2]
return c1y*z1 + c2y*z2
end
end
end
# Tri-linear interp, and maybe a generated function for generic multi-linear case?
# Or is there a fancy recusion algorithm?
# I haven't implemented a box filter, but it would use derivative information
# from the geometric transformation to create a rotated, rectangular filter of
# the original pixels
# Finally, an image transformation that computes the a whole image (non-lazy):
"""
An abstract transformation that applies a geometric transformation and a filter
"""
immutable ImageTransformation{NewSize <: Tuple, Trans <: AbstractTransformation, Filt <: Filter} <: AbstractTransformation{AbstractArray, AbstractArray}
tform::Trans
newsize::NewSize # A tuple of ranges, e.g. (1:pixels_x, 1:pixels_y)
filter::Filt
ImageTransformation(tform, sizes, filter = LinearFilter()) = new(tfrom,sizes,filter)
end
transform(trans::ImageTransformation, a::Array) = TransformedArray(a, trans.tform, trans.filter)[trans.newsizes...]
end # module
@andyferris
Copy link
Author

andyferris commented Jun 24, 2016

There seems to be a problem with the LinearFilter() when 1 < x < size(A,1), but maybe otherwise it is semi-working? It seems to be working for me now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment