[request] Array Interpolation via Stencils

[fosskers]

I assume that it's currently possible "expand" an image, say from 256x256 -> 512x512 via some manual traverse calls and simple interpolation.

Both for mapalgebra and for implementing a Haskell version of the hqx algorithm, it would be very handy if images could be upsampled/downsampled (I can never remember the difference) using stencils.

For instance, hxq wants to know the neighbourhood around each pixel. Based on the "shape" of the surrounding colours, it replaces the current pixel with 4 (a 2x2 grid) new pixels, effectively ballooning the image up by double. Currently, stencil operations assume "read some neighbourhood, write 1 pixel". Is this an unchangeable assumption?

Cheers, Colin

[lehins]

@fosskers That's a great question! Your assumption is correct, it would be fairly easy to implement upsampling (the ballooning up or adding new elements) right now using traverse or even backpermute. Same with downsampling, i.e. dropping current elements. Both of which I actually have an implementation for laying around in another repo, although only for 2D case of course. A naive version of Bilinear interpolation for 2D and linear interpolation for 1D could be added pretty quickly as well.

So the naive and slower way would be to first upsample an array, compute it, and only after that apply interpolation to it. This works only for scaling factors such as 2x, 3x, 4x, etc., but does not for fractional scaling, such as x1.7, for instance. Also it requires an intermediate in-memory representation, which, I am pretty sure, could be avoided.

The combination that you bring up of Interpolation via Stencils, I think, is an outstanding idea. I will ponder around on it and try to see what I can come up with.

[lehins]

Hey @fosskers, hope you are doing well. Just wanted to let you know that in the newest release of massiv-0.2.1 I added Stride functionality, which is the opposite of "balooning" up the image, but is still very useful, eg. when scaling down. Which means if you map a 3x3 stencil over an array and do computeWithStride on it with stride (3 :. 3), you will get O(n) computation.

I am still pondering on how to do the upsampling in a similarly efficient manner and don't have a solution yet, but do have some good ideas about it.

[fosskers]

Cool! So striding can help us shrink an image by tossing out intermediate pixels. Do we have any access to the ignored pixels during the shrinking process that would allow us to implement some sort of interpolation (i.e. the average colour of all the pixels in the Stride area)?

[lehins]

That's precisely where Stride can be useful. So here is an example you are talking about:

downSizeX3 :: (Load r' Ix2 Double, Manifest r' Ix2 Double, Mutable r Ix2 Double) =>
  Array r' Ix2 Double -> Array r Ix2 Double
downSizeX3 = computeWithStride (Stride 3) . mapStencil Edge average3x3Filter

average3x3Filter :: Stencil Ix2 Double Double
average3x3Filter = makeStencil (3 :. 3) (1 :. 1) $ \ get ->
  (  get (-1 :. -1) + get (-1 :. 0) + get (-1 :. 1) +
     get ( 0 :. -1) + get ( 0 :. 0) + get ( 0 :. 1) +
     get ( 1 :. -1) + get ( 1 :. 0) + get ( 1 :. 1)   ) / 9

[fosskers]

Oh, duh, functional programming. Do one thing and then do the other thing, hahaha.

In that case, I can add downsampling to mapalgebra right away.

[lehins]

In case if you want to scale by two:

average2x2Filter :: Stencil Ix2 Double Double
average2x2Filter = makeStencil (2 :. 2) (0 :. 0) $ \ get ->
  (  get ( 0 :. 0) + get ( 0 :. 1) +
     get ( 1 :. 0) + get ( 1 :. 1)   ) / 4

After that just use Stride 2