So after a good sleep the matrix math stuff I was reading seemed less scary.

There are still some things happening that I don’t understand but for all intents and purposes, I got it working!

Here’s the result for anyone interested. If I can get this fully working for my use-case, I’ll write up a post about it.

So this function takes an image and a list of destination points to skew it to. The destination points are just given as a flat list of repeated x, ys in the order: [bottom_left_x, bottom_right_y, bottom_right_x, bottom_right_y, top_right_x, top_right_y, top_left_x, top_left_y]. This implementation always uses the four corners of the original image as the source points.

This is obviously super spikey code that shouldn’t have that very specific point-order requirement:

defmodule Foo do
  def skew(%Vix.Vips.Image{} = image, points) do
    {width, height} = {Image.width(image), Image.height(image)}

    [sx1, sy1, sx2, sy2, sx3, sy3, sx4, sy4] = [0, 0, width, 0, width, height, 0, height]
    [dx1, dy1, dx2, dy2, dx3, dy3, dx4, dy4] = points

    src =
      [
        [sx1, sy1, 1, 0, 0, 0, -sx1 * dx1, -sy1 * dx1],
        [sx2, sy2, 1, 0, 0, 0, -sx2 * dx2, -sy2 * dx2],
        [sx3, sy3, 1, 0, 0, 0, -sx3 * dx3, -sy3 * dx3],
        [sx4, sy4, 1, 0, 0, 0, -sx4 * dx4, -sy4 * dx4],
        [0, 0, 0, sx1, sy1, 1, -sx1 * dy1, -sy1 * dy1],
        [0, 0, 0, sx2, sy2, 1, -sx2 * dy2, -sy2 * dy2],
        [0, 0, 0, sx3, sy3, 1, -sx3 * dy3, -sy3 * dy3],
        [0, 0, 0, sx4, sy4, 1, -sx4 * dy4, -sy4 * dy4],
      ]

    dest = [dx1, dx2, dx3, dx4, dy1, dy2, dy3, dy4]

    result = Nx.LinAlg.solve(Nx.tensor(src), Nx.tensor(dest))

    [t0, t1, t2, t3, t4, t5, t6, t7] = Nx.to_list(result)

    alias Vix.Vips.Operation, as: O

    with {:ok, index_image} <- O.xyz(width, height),
         {:ok, copy_image} <- O.copy(index_image),
         {:ok, band_0} <- O.extract_band(index_image, 0),
         {:ok, band_1} <- O.extract_band(copy_image, 1),
         {:ok, x_1} <- O.linear(band_0, [t0], [0]),
         {:ok, x_2} <- O.linear(band_1, [t1], [0]),
         {:ok, x_3} <- O.add(x_1, x_2),
         {:ok, x_a} <- O.linear(x_3, [1], [t2]),
         {:ok, x_4} <- O.linear(band_0, [t6], [0]),
         {:ok, x_5} <- O.linear(band_1, [t7], [0]),
         {:ok, x_6} <- O.add(x_4, x_5),
         {:ok, x_b} <- O.linear(x_6, [1], [1]),
         {:ok, x} <- O.divide(x_a, x_b),
         {:ok, y_1} <- O.linear(band_0, [t3], [0]),
         {:ok, y_2} <- O.linear(band_1, [t4], [0]),
         {:ok, y_3} <- O.add(y_1, y_2),
         {:ok, y_a} <- O.linear(y_3, [1], [t5]),
         {:ok, y_4} <- O.linear(band_0, [t6], [0]),
         {:ok, y_5} <- O.linear(band_1, [t7], [0]),
         {:ok, y_6} <- O.add(y_4, y_5),
         {:ok, y_b} <- O.linear(y_6, [1], [1]),
         {:ok, y} <- O.divide(y_a, y_b),
         {:ok, mapimage} <- O.bandjoin([x, y]) do
      O.mapim(image, mapimage)
    end
  end
end

Usage:

iex> {:ok, image} = Image.open("path/to/image.jpg")
iex> {:ok, result} = Foo.skew(image, [0, 0, Image.width(image), 0, Image.width(image), Image.height(image), 0, Image.height(image)]
iex> Image.write(result, "foo.jpg")

A bit of a trolling example but if it works you should get a completely unaltered image Otherwise, play around with the point values.

         {:ok, x_1} <- O.linear(band_0, [t0], [0]),
         {:ok, x_2} <- O.linear(band_1, [t1], [t2]),
         {:ok, x_a} <- O.add(x_1, x_2),
         {:ok, x_4} <- O.linear(band_0, [t6], [0]),
         {:ok, x_5} <- O.linear(band_1, [t7], [1]),
         {:ok, x_b} <- O.add(x_4, x_5),
         {:ok, x} <- O.divide(x_a, x_b),
         {:ok, y_1} <- O.linear(band_0, [t3], [0]),
         {:ok, y_2} <- O.linear(band_1, [t4], [t5]),
         {:ok, y_a} <- O.add(y_1, y_2),
         {:ok, y_4} <- O.linear(band_0, [t6], [0]),
         {:ok, y_5} <- O.linear(band_1, [t7], [1]),
         {:ok, y_b} <- O.add(y_4, y_5),
         {:ok, y} <- O.divide(y_a, y_b),

Slight simplification, rather than multiplying by 1 and adding a new value, you can add the constant in the previous call to linear (instead of adding 0).

You can import image operators for solving the equation, that should make it read better.

  defp generate_map(width, height, tensor) do
    alias Vix.Vips.Operation, as: O
    import Kernel, except: [+: 2, *: 2, /: 2]
    import Image.Math, only: [+: 2, *: 2, /: 2]

    [t0, t1, t2, t3, t4, t5, t6, t7] = Nx.to_list(tensor)
    index = O.xyz!(width, height)

    x =
      (
        t = index * [t0, t1]
        x_a = t[0] + t[1] + t2

        t = index * [t6, t7]
        x_b = t[0] + t[1] + 1

        x_a / x_b
      )

    y =
      (
        t = index * [t3, t4]
        y_a = t[0] + t[1] + t5

        t = index * [t6, t7]
        y_b = t[0] + t[1] + 1

        y_a / y_b
      )

    O.bandjoin!([x, y])
  end

But, maybe Nx is better suited for these types of array operations. You can switch between Vix and Nx, back and forth, easily. Using Vix for reading, and writing image in different formats and using Nx for something like this.

And yes, making a copy is not necessary.