Piecewise bilinear reconstruction of JPEG-style DC
==================================================

Ok, I have no idea if I got this right...

Suppose we constrain the reconstruction impulses
described here:
  
   http://cbloomrants.blogspot.com/2009/06/06-16-09-inverse-box-sampling.html

to be piecewise linear with corners at low-res-pixel
centers. This will turn out to mean that you only
have to reconstruct the centers themselves (e.g.
separably), and then bilerp those centers, which
should make it more computationally feasible.

So you effectively want a piecewise linear approximation to
your equation (or to the windowed sinc), one that produces
the DC=1/0 result.

Interestingly, this is perfectly constrained -- I think it
has to produce the same result as my iterated thing, since
both produce an end result that's bilerpable. But this way
lets us just solve for it.

Suppose at the center of the centermost (low-res) pixel
we're scaling by (1+k), then at the edge of the pixel
we want to be scaling by (1-k) so the average over the
pixel is 1.

That means that if we follow that line to the center of the
next pixel, we find ourselves at (1-k*3), so that's our
weight for the next pixel.

The average of the next (and remaining) pixel needs to
be *0*. The left half average is (1-k + 1-k*3)/2, or (1-k*2).
So the right half needs to have an average of (-1+k*2). What
value at the next pixel crossing matches that?
   (1-k*3 + j)/2 = (-1+k*2)
   1-k*3+j = -2 + k*4
   1+j = -2 + k*7
     j = -3 + k*7

If we now trace the line through j into the next pixel,
that places the center point at:
   next = 1-k*3 + (j - (1-k*3))*2
        = 1-k*3 + (-3+k*7 -1 + k*3))*2
        = 1-k*3 + (-4 + k*10)*2
        = 1-k*3 - 8 + k*20
        = -7 + k*17

Find the new edge point by a similar process:
   pseudo-avg = (k*7-3 + k*17-7) = (k*24-10)
   edge pixel:  (k*17-7 + j) = -k*24 + 10
                 j = -k*41 + 17
   new center:  k*17-7 + (j - (k*17-7))*2
                 41 - 99 * k

Index of integer sequences is no help yet (and we kind of want fractional
sequences). Now, this sequence is going to "converge" in the sense that
we can ask what happens if we solve for the value of k that makes the
limit of these weights be 0. I assume that setting k to this value is
clearly the right thing to do, since we don't want the weights to diverge.

Also, we can stop early if we force a *crossing* point to go to 0 by
changing a weight. E.g, we constrain:

    edge pixel -k*41+17 = 0, k = 17/41ths

This gives us a filter with these weights:

   58/41, -10/41, 2/41, ** -2/41 (only go halfway to this, stop at edge)

The problem is that this requires us to stop the bilerp at
the edge (when it hits 0), which means we need to explicitly
generate the half-way points and bilerp from them. On the other
hand, our filter width is small enough that we don't need to
keep track of that many halfway points at once--in fact, just
one row's worth at a time (since we can generate them directly
from the other data). Mind, this is a continuous solution and
might not be accurate when solved discretely (although we could
pretty easily analyze and tweak this).

Of course, if you were going to interpolate from halfway-points,
you'd probably be better off tweaking the values at those halfway
points to better approximate the sinc-ish curve, instead of having
all but the last fall halfway between adjacent things, i.e. the
weights for the halfway points sugggest by the above are:

   58/41, 24/41, -10/41, -4/41, 2/41, 0, 0, 0, 0, 0,

but these should be tweaked since we have more actual degrees of
freedom.

If we stick to pure bilerp of the centers, then we have to carry
things out indefinitely, so let me symbolically solve how to generate
each successive center value instead of doing it by hand:

   compute center(i) from center(i-2) and center(i-1):
       crossing point into pixel i-1 from i-2:
           (center(i-2) + center(i-1))/2
       average of left half of region of i-1
           (center(i-1) + (center(i-2) + center(i-1))/2)/2
       average of right half of region i-1 to give average of 0:
           -(center(i-1) + (center(i-2) + center(i-1))/2)/2
       crossing point from pixel i-1 to pixel i:
           -(2*center(i-1) + (center(i-2) + center(i-1))/2)
       center of pixel i
           -6*center(i-1) - center(i-2)         

   let center(i) = ai + k*bi

   then
        ai = -6*a(i-1) - a(i-2)
        bi = -6*b(i-1) - b(i-2)

    i      a         b   |a/b|
    0      1         1
    1      1        -3   
    2     -7        17   0.411764
    3     41       -99   0.414141
    4   -239       577   0.414211438
    5   1393     -3363   0.41421349985
    6   8119     19601   0.41421356053

              sqrt(2) == 1.414213562373

(Online Encyclopeda of Integer Sequences indicates that the abs() of
the second sequence is the 'Chebyshev polynomials of the first kind evaluated at 3'.
The first sequence is the 'Newman-Shanks-Williams' sequence.)

Thus, the solution to the unbounded version of the problem, ignoring
edge fixup, is to use k = sqrt(2) - 1, and thus the weights are:

    i     weight = a+b*k
    0    1.414213562373095
    1   -0.242640687119285
    2    0.041630560342617
    3   -0.007142674936419
    4    0.001225489275899
    5   -0.000210260718974
    6    0.000036075037946
    7   -0.000006189508704
    8    0.000001062014280
    9   -0.000000182576977
   10    0.000000033447584

(interestingly, solving f(x)=-6*f(x-2)+f(x-1) with doubles diverges at this point already):
   11   -0.000000018108524
   12    0.000000075203562
   13   -0.000000433112847
   14    0.000002523473518
   15   -0.000014707728264
   16    0.000085722896063
   17   -0.000499629648113
   18    0.002912054992613
   19   -0.016972700307564     
   
Note that |weight(4)| < 1/256, so it should be possible to stop somewhere
around there naturally (maybe go another step or two). But it needs to be
tweaked to handle being discrete instead of continuous.

Edge fixup is probably best solved by computing the bilerp for all non-edge
pixels, then going back along the edge 8x8 blocks and reflecting the bilerped
data into the unfilled regions--this is an easy way to guarantee matched DC and
smoothness, even if it's pretty bogus (and totally wrong for wrapped textures).

Note this won't actually produce the same result as the iterated solution I
showed Charles previously, because the iterated solution actually clamped
to 0..255 while solving, so e.g. would actually reproduce a b&w checkerboard
that exactly fell on block boundaries, whereas the above approach will massively,
massively overshoot trying to smooth it out, since it's allowed to use arbitrary
numbers instead of clamping to 0..255.

The simplest approximation to this I can think of would be be to use half-pixel
bilerping with weights something like:

     x
    0.0    1.5
    0.5    0.5
    1.0   -0.25
    1.5    0
    2.0    0

The general equation here that gives you something tweakable:

   x      
   0.0    (1.0 + k) 
   0.5    (1.0 - k)
   1.0    -j
   1.5    0

Force DC of 0.5 .. 1.5 to be 0:
  (1.0 - k - j) = - (-j + 0)
  1 - k - j = j
  1-k = 2*j
  j = (1-k)/2

General weights to use with one free control factor:

    x       weight
   0.0      1+k
   0.5      1-k
   1.0      -(1-k)/2
   1.5      0
   2.0      0

Rewrite it in terms of weight_0:
   0.0      w
   0.5      2-w
   1.0      -(2-w)/2 = w/2-1
   1.5      0
   2.0      0

Sample values:

    i      w=1    w=1.2    w=1.4
   0.0      1      1.2      1.4 
   0.5      1      0.8      0.6 
   1.0     -0.5   -0.4     -0.3     
   1.5      0      0.0      0.0

Now, w isn't actually tunable; if the DCs are constant, we
want a constant output, we don't want something oscillating
up and down around it. So we want to have the constraint
that the sum of the weights at a given point is 1.
So weight(0)+2*weight(1) = 1, hence
     w + 2*(w/2-1) = 1, w+w-2 = 1, 2w=3, w=3/2
which is exactly the one I came up with above:

     x
    0.0    1.5
    0.5    0.5
    1.0   -0.25
    1.5    0
    2.0    0

Again, this needs to be tweaked to produce the correct DC
in the discrete case.

To be explicit about 'half-pixel bilerping', here's the 8x8
case with 1d pixels, using the above simplified approach:

   Given input values D[i] with 0 <= i < n:
   
   Construct a new array H[j] with 0 <= j < 2*n-1
      H(j) =
        j even: 1.5*D(j/2) + -0.25 * (D[j/2-1  ] + D[j/2+1  ])
        j odd:                0.50 * (D[j/2-1/2] + D[j/2+1/2])
    
   Now, construct the output samples by lerping H[j]:
   
   Let s(A,x) be a piecewise-linear sampling function:
      s(A,x) =  A[floor(x)] + frac(x) * (A[floor(x)+1] - A[floor(x)])
   
   Construct a new array O[k] with 4 <= k < 8*n-3
      O[k] = s(H, (k+0.5)/4)
   
   From a block-by-block standpoint, to solve the 1d block #b:
   
   Let p0 =    0.5  * D[b-1] +  0.5 * D[b]
   Let p1 =   -0.25 * D[b-1] +  1.0 * D[b] +  -0.25 * D[b+1]
   Let p2 =                     0.5 * D[b] +   0.50 * D[b+1]
   
   Let o0 = lerp(1/8, p0, p1)
   Let o1 = lerp(3/8, p0, p1)
   Let o2 = lerp(5/8, p0, p1)
   Let o3 = lerp(7/8, p0, p1)
   Let o4 = lerp(7/8, p1, p2)
   Let o5 = lerp(5/8, p1, p2)
   Let o6 = lerp(3/8, p1, p2)
   Let o7 = lerp(1/8, p1, p2)

Again, this hasn't been corrected for discreteness, but the above
formulation actually gives us the tools we need to solve it; we
just need to tweak the weights 0.5 and -0.25 to make sure
they produce DCs of 1 and 0 as appropriate -- i.e. tweaked to make
sure the average of all the o's is exactly D[b].