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].