# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1574093616 18000
# Node ID eef71dd7202ba39830ea2beecde068407aa7e445

Initialise

diff -r 000000000000 -r eef71dd7202b Project.toml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Project.toml	Mon Nov 18 11:13:36 2019 -0500
@@ -0,0 +1,4 @@
+name = "ImageTools"
+uuid = "b548cc0d-4ade-417e-bf62-0e39f9d2eee9"
+authors = ["tuomov "]
+version = "0.1.0"
diff -r 000000000000 -r eef71dd7202b src/Gradient.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Gradient.jl	Mon Nov 18 11:13:36 2019 -0500
@@ -0,0 +1,185 @@
+########################
+# Discretised gradients
+########################
+
+module Gradient
+
+##############
+# Our exports
+##############
+
+export ∇₂!, ∇₂ᵀ!, ∇₂fold!,
+       ∇₂_norm₂₂_est, ∇₂_norm₂₂_est²,
+       ∇₂_norm₂∞_est, ∇₂_norm₂∞_est²,
+       ∇₂c!,
+       ∇₃!, ∇₃ᵀ!,
+       vec∇₃!, vec∇₃ᵀ!
+
+##################
+# Helper routines
+##################
+
+@inline function imfold₂′!(f_aa!, f_a0!, f_ab!,
+                           f_0a!, f_00!, f_0b!,
+                           f_ba!, f_b0!, f_bb!,
+                           n, m, state)
+    # First row
+    state = f_aa!(state, (1, 1))
+    for j = 2:m-1
+        state = f_a0!(state, (1, j))
+    end
+    state = f_ab!(state, (1, m))
+
+    # Middle rows
+    for i=2:n-1
+        state = f_0a!(state, (i, 1))
+        for j = 2:m-1
+            state = f_00!(state, (i, j))
+        end
+        state = f_0b!(state, (i, m))
+    end
+
+    # Last row
+    state = f_ba!(state, (n, 1))
+    for  j =2:m-1
+        state = f_b0!(state, (n, j))
+    end
+    return f_bb!(state, (n, m))
+end
+
+#########################
+# 2D forward differences
+#########################
+
+∇₂_norm₂₂_est² = 8
+∇₂_norm₂₂_est = √∇₂_norm₂₂_est²
+∇₂_norm₂∞_est² = 2
+∇₂_norm₂∞_est = √∇₂_norm₂∞_est²
+
+function ∇₂!(u₁, u₂, u)
+    @. @views begin
+        u₁[1:(end-1), :] = u[2:end, :] - u[1:(end-1), :]
+        u₁[end, :, :] = 0
+
+        u₂[:, 1:(end-1)] = u[:, 2:end] - u[:, 1:(end-1)]
+        u₂[:, end] = 0
+    end
+    return u₁, u₂
+end
+
+function ∇₂!(v, u)
+    ∇₂!(@view(v[1, :, :]), @view(v[2, :, :]), u)
+end
+
+@inline function ∇₂fold!(f!::Function, u, state)
+    g! = (state, pt) -> begin
+        (i, j) = pt
+        g = @inbounds [u[i+1, j]-u[i, j], u[i, j+1]-u[i, j]]
+        return f!(g, state, pt)
+    end
+    gr! = (state, pt) -> begin
+        (i, j) = pt
+        g = @inbounds [u[i+1, j]-u[i, j], 0.0]
+        return f!(g, state, pt)
+    end
+    gb! = (state, pt) -> begin
+        (i, j) = pt
+        g = @inbounds [0.0, u[i, j+1]-u[i, j]]
+        return f!(g, state, pt)
+    end
+    g0! = (state, pt) -> begin
+        return f!([0.0, 0.0], state, pt)
+    end
+    return imfold₂′!(g!, g!, gr!,
+                     g!, g!, gr!,
+                     gb!, gb!, g0!,
+                     size(u, 1), size(u, 2), state)
+end
+
+function ∇₂ᵀ!(v, v₁, v₂)
+    @. @views begin
+        v[2:(end-1), :] = v₁[1:(end-2), :] - v₁[2:(end-1), :]
+        v[1, :] = -v₁[1, :]
+        v[end, :] = v₁[end-1, :]
+
+        v[:, 2:(end-1)] += v₂[:, 1:(end-2)] - v₂[:, 2:(end-1)]
+        v[:, 1] += -v₂[:, 1]
+        v[:, end] += v₂[:, end-1]
+    end
+    return v
+end
+
+function ∇₂ᵀ!(u, v)
+    ∇₂ᵀ!(u, @view(v[1, :, :]), @view(v[2, :, :]))
+end
+
+##################################################
+# 2D central differences (partial implementation)
+##################################################
+
+function ∇₂c!(v, u)
+    @. @views begin
+        v[1, 2:(end-1), :] = (u[3:end, :] - u[1:(end-2), :])/2
+        v[1, end, :] = (u[end, :] - u[end-1, :])/2
+        v[1, 1, :] = (u[2, :] - u[1, :])/2
+
+        v[2, :, 2:(end-1)] = (u[:, 3:end] - u[:, 1:(end-2)])/2
+        v[2, :, end] = (u[:, end] - u[:, end-1])/2
+        v[2, :, 1] = (u[:, 2] - u[:, 1])/2
+    end
+end
+
+#########################
+# 3D forward differences
+#########################
+
+function ∇₃!(u₁,u₂,u₃,u)
+    @. @views begin
+        u₁[1:(end-1), :, :] = u[2:end, :, :] - u[1:(end-1), :, :]
+        u₁[end, :, :] = 0
+
+        u₂[:, 1:(end-1), :] = u[:, 2:end, :] - u[:, 1:(end-1), :]
+        u₂[:, end, :] = 0
+
+        u₃[:, :, 1:(end-1)] = u[:, :, 2:end] - u[:, :, 1:(end-1)]
+        u₃[:, :, end] = 0
+    end
+    return u₁, u₂, u₃
+end
+
+function ∇₃ᵀ!(v,v₁,v₂,v₃)
+    @. @views begin
+        v[2:(end-1), :, :] = v₁[1:(end-2), :, :] - v₁[2:(end-1), :, :]
+        v[1, :, :] = -v₁[1, :, :]
+        v[end, :, :] = v₁[end-1, :, :]
+
+        v[:, 2:(end-1), :] += v₂[:, 1:(end-2), :] - v₂[:, 2:(end-1), :]
+        v[:, 1, :] += -v₂[:, 1, :]
+        v[:, end, :] += v₂[:, end-1, :]
+
+        v[:, :, 2:(end-1)] += v₃[:, :, 1:(end-2)] - v₃[:, :, 2:(end-1)]
+        v[:, :, 1] += -v₃[:, :, 1]
+        v[:, :, end] += v₃[:, :, end-1]
+    end
+    return v
+end
+
+###########################################
+# 3D forward differences for vector fields
+###########################################
+
+function vec∇₃!(u₁,u₂,u₃,u)
+    @. @views for j=1:size(u, 1)
+        ∇₃!(u₁[j, :, :, :],u₂[j, :, :, :],u₃[j, :, :, :],u[j, :, :, :])
+    end
+    return u₁, u₂, u₃
+end
+
+function vec∇₃ᵀ!(u,v₁,v₂,v₃)
+    @. @views for j=1:size(u, 1)
+        ∇₃ᵀ!(u[j, :, :, :],v₁[j, :, :, :],v₂[j, :, :, :],v₃[j, :, :, :])
+    end
+    return u
+end
+
+end # Module
diff -r 000000000000 -r eef71dd7202b src/Translate.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Translate.jl	Mon Nov 18 11:13:36 2019 -0500
@@ -0,0 +1,148 @@
+######################################
+# Image subpixel accuracy translation
+######################################
+
+module Translate
+
+##########
+# Exports
+##########
+
+export interpolate2d,
+       interpolate2d_quadrants,
+       extract_subimage!,
+       translate_image!,
+       DisplacementFull,
+       DisplacementConstant,
+       Displacement,
+       Image
+
+##################
+# Types
+##################
+
+# Two different types of displacement data supported:
+#  a) given in each pixel
+#  b) constant in space
+Image = Array{Float64,2}
+DisplacementFull = Array{Float64,3}
+DisplacementConstant = Array{Float64,1}
+Displacement = Union{DisplacementFull,DisplacementConstant}
+
+#############################
+# Base interpolation routine
+#############################
+
+@inline function interpolate2d_quadrants(v, (x, y))
+    (m, n) = size(v)
+    clipx = xʹ -> max(1, min(xʹ, m))
+    clipy = yʹ -> max(1, min(yʹ, n))
+
+    xfℤ = clipx(floor(Int, x))
+    xcℤ = clipx(ceil(Int, x))
+    yfℤ = clipy(floor(Int, y))
+    ycℤ = clipy(ceil(Int, y))
+       
+    xf = convert(Float64, xfℤ)
+    xc = convert(Float64, xcℤ)
+    yf = convert(Float64, yfℤ)
+    yc = convert(Float64, ycℤ)
+    xm = (xf+xc)/2
+    ym = (yf+yc)/2
+
+    vff = @inbounds v[xfℤ, yfℤ]
+    vfc = @inbounds v[xfℤ, ycℤ]
+    vcf = @inbounds v[xcℤ, yfℤ]
+    vcc = @inbounds v[xcℤ, ycℤ]
+    vmm = (vff+vfc+vcf+vcc)/4
+
+    if xfℤ==xcℤ
+        if yfℤ==ycℤ
+            # Completely degenerate case
+            v = vmm
+        else
+            # Degenerate x
+            v = vff+(y-yf)/(yc-yf)*(vfc-vff)
+        end
+    elseif yfℤ==ycℤ
+        # Degenerate y
+        v = vff + (x-xf)/(xc-xf)*(vcf-vff)
+    elseif y-ym ≥ x-xm
+        # top-left half
+        if (y-ym) + (x-xm) ≥ 0
+            # top quadrant
+            v = vfc + (x-xf)/(xc-xf)*(vcc-vfc) + (y-yc)/(ym-yc)*(vmm-(vcc+vfc)/2)
+        else
+            # left quadrant
+            v = vff + (y-yf)/(yc-yf)*(vfc-vff) + (x-xf)/(xm-xf)*(vmm-(vfc+vff)/2)
+        end
+    else
+        # bottom-left half
+        if (y-ym) + (x-xm) ≥ 0
+            # right quadrant
+            v = vcf + (y-yf)/(yc-yf)*(vcc-vcf) + (x-xc)/(xm-xc)*(vmm-(vcc+vcf)/2)
+        else
+            # bottom quadrant
+            v = vff + (x-xf)/(xc-xf)*(vcf-vff) + (y-yf)/(ym-yf)*(vmm-(vcf+vff)/2)
+        end
+    end
+
+    return v
+end
+
+interpolate2d = interpolate2d_quadrants
+
+##############
+# Translation
+##############
+
+@polly function translate_image!(x, z, u::DisplacementFull)
+    @assert(size(u, 1)==2 && size(x)==size(u)[2:end] && size(x)==size(z))
+
+    @inbounds @simd for i=1:size(x, 1)
+        @simd for j=1:size(x, 2)
+            pt = (i - u[1, i, j], j - u[2, i, j])
+            x[i, j] = interpolate2d(z, pt)
+        end
+    end
+end
+
+@polly function translate_image!(x, z, u::DisplacementConstant)
+    @assert(size(u)==(2,) && size(x)==size(z))
+
+    @inbounds @simd for i=1:size(x, 1)
+        @simd for j=1:size(x, 2)
+            pt = (i - u[1], j - u[2])
+            x[i, j] = interpolate2d(z, pt)
+        end
+    end
+end
+
+######################
+# Subimage extraction
+######################
+
+@polly function extract_subimage!(b, im, v::DisplacementConstant)
+    (imx, imy) = size(im)
+    (bx, by) = size(b)
+
+    # Translation from target to source coordinates
+    vxʹ = v[1] + (imx-bx)/2
+    vyʹ = v[2] + (imy-by)/2
+
+    # Target image indices within source image
+    px = ceil(Int, max(1, vxʹ + 1) - vxʹ)
+    py = ceil(Int, max(1, vyʹ + 1) - vyʹ)
+    qx = floor(Int, min(imx, vxʹ + bx) - vxʹ)
+    qy = floor(Int, min(imy, vyʹ + by) - vyʹ)
+    
+    b .= 0
+
+    @inbounds @simd for i=px:qx
+        for j=py:qy
+            b[i, j] = interpolate2d(im, (i+vxʹ, j+vyʹ))
+        end
+    end
+end
+
+end