Revert the optimization of the bilinear weights.
authorSteinar H. Gunderson <sgunderson@bigfoot.com>
Mon, 23 Feb 2015 00:19:03 +0000 (01:19 +0100)
committerSteinar H. Gunderson <sgunderson@bigfoot.com>
Mon, 23 Feb 2015 00:19:03 +0000 (01:19 +0100)
For the case where the resampling changed every frame (e.g. a zoom),
it just consumed too much CPU to be worth it, especially in memory
management; this is painful because it was an elegant solution to
a tricky problem, but it just has to go for now.

Also drop out to fp32 at the first sight of too-high error.

resample_effect.cpp

index f438a87..aa3a403 100644 (file)
@@ -224,112 +224,6 @@ double compute_sum_sq_error(const Tap<float>* weights, unsigned num_weights,
        return sum_sq_error;
 }
 
        return sum_sq_error;
 }
 
-// Given a predefined, fixed set of bilinear weight positions, try to optimize
-// their weights through some linear algebra. This can do a better job than
-// the weight calculation in combine_samples() because it can look at the entire
-// picture (an effective weight can sometimes be affected by multiple samples).
-// It will also optimize weights for non-combined samples, which is useful when
-// a sample happens in-between texels for numerical reasons.
-//
-// The math goes as follows: The desired result is a weighted sum, where the
-// weights are the coefficients in <weights>:
-//
-//   y = sum(c_j x_j, j)
-//
-// We try to approximate this by a different set of coefficients, which have
-// weights d_i and are placed at some fraction to the right of a source texel x_j.
-// This means it will influence two texels (x_j and x_{j+1}); generalizing this,
-// let us define that w_ij means the amount texel <j> influences bilinear weight
-// <i> (keeping in mind that w_ij = 0 for all but at most two different j).
-// This means the actually computed result is:
-//
-//   y' = sum(d_i w_ij x_j, j)
-//
-// We assume w_ij fixed and wish to find {d_i} so that y' gets as close to y
-// as possible. Specifically, let us consider the sum of squred errors of the
-// coefficients:
-//
-//   ε² = sum((sum( d_i w_ij, i ) - c_j)², j)
-//
-// The standard trick, which also applies just fine here, is to differentiate
-// the error with respect to each variable we wish to optimize, and set each
-// such expression to zero. Solving this equation set (which we can do efficiently
-// by letting Eigen invert a sparse matrix for us) yields the minimum possible
-// error. To see the form each such equation takes, pick any value k and
-// differentiate the expression by d_k:
-//
-//   ∂(ε²)/∂(d_k) = sum(2(sum( d_i w_ij, i ) - c_j) w_kj, j)
-//
-// Setting this expression equal to zero, dropping the irrelevant factor 2 and
-// rearranging yields:
-//
-//   sum(w_kj sum( d_i w_ij, i ), j) = sum(w_kj c_j, j)
-//
-// where again, we remember where the sums over j are over at most two elements,
-// since w_kj is nonzero for at most two values of j.
-template<class T>
-void optimize_sum_sq_error(const Tap<float>* weights, unsigned num_weights,
-                           Tap<T>* bilinear_weights, unsigned num_bilinear_weights,
-                           unsigned size)
-{
-       // Find the range of the desired weights.
-       int c_lower_pos = lrintf(weights[0].pos * size - 0.5);
-       int c_upper_pos = lrintf(weights[num_weights - 1].pos * size - 0.5) + 1;
-
-       SparseMatrix<float> A(num_bilinear_weights, num_bilinear_weights);
-       SparseVector<float> b(num_bilinear_weights);
-
-       // Convert each bilinear weight to the (x, frac) form for less junk in the code below.
-       int* pos = new int[num_bilinear_weights];
-       float* fracs = new float[num_bilinear_weights];
-       for (unsigned i = 0; i < num_bilinear_weights; ++i) {
-               const float pixel_pos = to_fp64(bilinear_weights[i].pos) * size - 0.5f;
-               const float f = pixel_pos - floor(pixel_pos);
-               pos[i] = int(floor(pixel_pos));
-               fracs[i] = lrintf(f / movit_texel_subpixel_precision) * movit_texel_subpixel_precision;
-       }
-
-       // The index ordering is a bit unusual to fit better with the
-       // notation in the derivation above.
-       for (unsigned k = 0; k < num_bilinear_weights; ++k) {
-               for (int j = pos[k]; j <= pos[k] + 1; ++j) {
-                       const float w_kj = (j == pos[k]) ? (1.0f - fracs[k]) : fracs[k];
-                       for (unsigned i = 0; i < num_bilinear_weights; ++i) {
-                               float w_ij;
-                               if (j == pos[i]) {
-                                       w_ij = 1.0f - fracs[i];
-                               } else if (j == pos[i] + 1) {
-                                       w_ij = fracs[i];
-                               } else {
-                                       // w_ij = 0
-                                       continue;
-                               }
-                               A.coeffRef(i, k) += w_kj * w_ij;
-                       }
-                       float c_j;
-                       if (j >= c_lower_pos && j < c_upper_pos) {
-                               c_j = weights[j - c_lower_pos].weight;
-                       } else {
-                               c_j = 0.0f;
-                       }
-                       b.coeffRef(k) += w_kj * c_j;
-               }
-       }
-       delete[] pos;
-       delete[] fracs;
-
-       A.makeCompressed();
-       SparseQR<SparseMatrix<float>, COLAMDOrdering<int> > qr(A);
-       assert(qr.info() == Success);
-       SparseMatrix<float> new_weights = qr.solve(b);
-       assert(qr.info() == Success);
-
-       for (unsigned i = 0; i < num_bilinear_weights; ++i) {
-               bilinear_weights[i].weight = from_fp64<T>(new_weights.coeff(i, 0));
-       }
-       normalize_sum(bilinear_weights, num_bilinear_weights);
-}
-
 }  // namespace
 
 ResampleEffect::ResampleEffect()
 }  // namespace
 
 ResampleEffect::ResampleEffect()
@@ -616,34 +510,28 @@ void SingleResamplePassEffect::update_texture(GLuint glsl_program_num, const str
 
        // Now make use of the bilinear filtering in the GPU to reduce the number of samples
        // we need to make. Try fp16 first; if it's not accurate enough, we go to fp32.
 
        // Now make use of the bilinear filtering in the GPU to reduce the number of samples
        // we need to make. Try fp16 first; if it's not accurate enough, we go to fp32.
+       // Our tolerance level for total error is a bit higher than the one for invididual
+       // samples, since one would assume overall errors in the shape don't matter as much.
+       const float max_error = 2.0f / (255.0f * 255.0f);
        Tap<fp16_int_t> *bilinear_weights_fp16;
        src_bilinear_samples = combine_many_samples(weights, src_size, src_samples, dst_samples, &bilinear_weights_fp16);
        Tap<float> *bilinear_weights_fp32 = NULL;
        bool fallback_to_fp32 = false;
        double max_sum_sq_error_fp16 = 0.0;
        for (unsigned y = 0; y < dst_samples; ++y) {
        Tap<fp16_int_t> *bilinear_weights_fp16;
        src_bilinear_samples = combine_many_samples(weights, src_size, src_samples, dst_samples, &bilinear_weights_fp16);
        Tap<float> *bilinear_weights_fp32 = NULL;
        bool fallback_to_fp32 = false;
        double max_sum_sq_error_fp16 = 0.0;
        for (unsigned y = 0; y < dst_samples; ++y) {
-               optimize_sum_sq_error(
-                       weights + y * src_samples, src_samples,
-                       bilinear_weights_fp16 + y * src_bilinear_samples, src_bilinear_samples,
-                       src_size);
                double sum_sq_error_fp16 = compute_sum_sq_error(
                        weights + y * src_samples, src_samples,
                        bilinear_weights_fp16 + y * src_bilinear_samples, src_bilinear_samples,
                        src_size);
                max_sum_sq_error_fp16 = std::max(max_sum_sq_error_fp16, sum_sq_error_fp16);
                double sum_sq_error_fp16 = compute_sum_sq_error(
                        weights + y * src_samples, src_samples,
                        bilinear_weights_fp16 + y * src_bilinear_samples, src_bilinear_samples,
                        src_size);
                max_sum_sq_error_fp16 = std::max(max_sum_sq_error_fp16, sum_sq_error_fp16);
+               if (max_sum_sq_error_fp16 > max_error) {
+                       break;
+               }
        }
 
        }
 
-       // Our tolerance level for total error is a bit higher than the one for invididual
-       // samples, since one would assume overall errors in the shape don't matter as much.
-       if (max_sum_sq_error_fp16 > 2.0f / (255.0f * 255.0f)) {
+       if (max_sum_sq_error_fp16 > max_error) {
                fallback_to_fp32 = true;
                src_bilinear_samples = combine_many_samples(weights, src_size, src_samples, dst_samples, &bilinear_weights_fp32);
                fallback_to_fp32 = true;
                src_bilinear_samples = combine_many_samples(weights, src_size, src_samples, dst_samples, &bilinear_weights_fp32);
-               for (unsigned y = 0; y < dst_samples; ++y) {
-                       optimize_sum_sq_error(
-                               weights + y * src_samples, src_samples,
-                               bilinear_weights_fp32 + y * src_bilinear_samples, src_bilinear_samples,
-                               src_size);
-               }
        }
 
        // Encode as a two-component texture. Note the GL_REPEAT.
        }
 
        // Encode as a two-component texture. Note the GL_REPEAT.