add basic image warping

src/SemiSmoothNewton.jl

+module SemiSmoothNewton
+
+include("image.jl")
+include("mesh.jl")
+
+function run()
+
+end
+
+
+function clip_line(r0, r1, l0, l1)
+    d_l = -(l1[1] - l0[1])
+    d_r = +(l1[1] - l0[1])
+    d_b = -(l1[2] - l0[2])
+    d_t = +(l1[2] - l0[2])
+
+    t_l = -(r0[1] - l0[1]) / d_l
+    t_r = +(r1[1] - l0[1]) / d_r
+    t_b = -(r0[2] - l0[2]) / d_b
+    t_t = +(r1[2] - l0[2]) / d_t
+
+    t0 = max(
+	d_l <= 0 ? t_l : 0.,
+	d_r <= 0 ? t_r : 0.,
+	d_b <= 0 ? t_b : 0.,
+	d_t <= 0 ? t_t : 0.)
+
+    t1 = min(
+	d_l >= 0 ? t_l : 1.,
+	d_r >= 0 ? t_r : 1.,
+	d_b >= 0 ? t_b : 1.,
+	d_t >= 0 ? t_t : 1.)
+
+    # tolerance on non-intersection
+    t1 < t0 + eps() && return false
+
+    c0 = l0 + t0 * (l1 - l0)
+    c1 = l0 + t1 * (l1 - l0)
+
+    return c0, c1
+end
+
+
+
+# meshfunction
+#   mesh
+
+# mesh
+#   vertices: idx -> vertex
+#   cells: idx -> cell
+
+"project array data onto a grid function"
+function project_array!(u, img::Array)
+    vertices = u.mesh.vertices
+    mesh_box = (
+	reduce((a, b) -> min.(a, b), vertices),
+	reduce((a, b) -> max.(a, b), vertices))
+
+
+    for cell in u.mesh.cells
+	cell_box = (
+	    reduce((a, b) -> min.(a, b), cell.vertices),
+	    reduce((a, b) -> max.(a, b), cell.vertices))
+
+	u[cell]
+
+    end
+end
+
+
+
+end # module

src/image.jl

+eval_neumann(img, x) = img[clamp.(Tuple(x), axes(img))...]
+
+function interpolate_bilinear(img, x)
+    x0 = floor.(Int, x)
+    x1 = x0 .+ 1
+
+    val = zero(eltype(img))
+    for idx in CartesianIndices(ntuple(_->0:1, ndims(img)))
+	cornerbool = Bool.(Tuple(idx))
+	λ = ifelse.(cornerbool, x .- x0, x1 .- x)
+	corner = ifelse.(cornerbool, x1, x0)
+	val += prod(λ) * eval_neumann(img, CartesianIndex(corner))
+    end
+    return val
+end
+
+function warp_backwards(img, u)
+    axes(img) == axes(u)[2:end] ||
+	throw(ArgumentError("invalid dimensions"))
+
+    res = similar(img)
+    for I in CartesianIndices(img)
+	x = Tuple(I) .+ u[:, I]
+	res[I] = interpolate_bilinear(img, (x...,))
+    end
+    return res
+end

src/mesh.jl

+using WriteVTK
+
+# 2d, simplex grid
+struct Mesh
+    vertices::Array{Float64, 2}
+    cells::Array{NTuple{3, Int64}}
+end
+
+
+
+
+#function append_data!(vtkfile, f::GlobalMeshFunction{<:TriangleGrid, n}) where n
+#    fgrid = mesh(f)
+#    fdata = Array{Float64, 2}(undef, n, size(fgrid.vertices, 2))
+#    for i in axes(fdata, 2)
+#        fdata[:, i] .= f.f(fgrid.vertices[:, i])
+#    end
+#    vtk_point_data(vtkfile, fdata, "f")
+#end
+#
+#cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, grid.cells[:, i]) for i in axes(grid.cells, 2)]
+#vtkfile = vtk_grid(filename, grid.vertices, cells)
+#for (i, f) in enumerate(fs)
+#    append_data!(vtkfile, f)
+#end
+#vtk_save(vtkfile)