Possible image transformations design

ImageTransformations.jl

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

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.