// Three-lobed Lanczos, the most common choice.
+// Note that if you change this, the accuracy for LANCZOS_TABLE_SIZE
+// needs to be recomputed.
#define LANCZOS_RADIUS 3.0
#include <epoxy/gl.h>
}
}
-float lanczos_weight(float x, float a)
+float lanczos_weight(float x)
{
- if (fabs(x) > a) {
+ if (fabs(x) > LANCZOS_RADIUS) {
return 0.0f;
} else {
- return sinc(M_PI * x) * sinc(M_PI * x / a);
+ return sinc(M_PI * x) * sinc((M_PI / LANCZOS_RADIUS) * x);
}
}
+// The weight function can be expensive to compute over and over again
+// (which will happen during e.g. a zoom), but it is also easy to interpolate
+// linearly. We compute the right half of the function (in the range of
+// 0..LANCZOS_RADIUS), with two guard elements for easier interpolation, and
+// linearly interpolate to get our function.
+//
+// We want to scale the table so that the maximum error is always smaller
+// than 1e-6. As per http://www-solar.mcs.st-andrews.ac.uk/~clare/Lectures/num-analysis/Numan_chap3.pdf,
+// the error for interpolating a function linearly between points [a,b] is
+//
+// e = 1/2 (x-a)(x-b) f''(u_x)
+//
+// for some point u_x in [a,b] (where f(x) is our Lanczos function; we're
+// assuming LANCZOS_RADIUS=3 from here on). Obviously this is bounded by
+// f''(x) over the entire range. Numeric optimization shows the maximum of
+// |f''(x)| to be in x=1.09369819474562880, with the value 2.40067758733152381.
+// So if the steps between consecutive values are called d, we get
+//
+// |e| <= 1/2 (d/2)^2 2.4007
+// |e| <= 0.1367 d^2
+//
+// Solve for e = 1e-6 yields a step size of 0.0027, which to cover the range
+// 0..3 needs 1109 steps. We round up to the next power of two, just to be sure.
+//
+// You need to call lanczos_table_init_done before the first call to
+// lanczos_weight_cached.
+#define LANCZOS_TABLE_SIZE 2048
+bool lanczos_table_init_done = false;
+float lanczos_table[LANCZOS_TABLE_SIZE + 2];
+
+void init_lanczos_table()
+{
+ for (unsigned i = 0; i < LANCZOS_TABLE_SIZE + 2; ++i) {
+ lanczos_table[i] = lanczos_weight(float(i) * (LANCZOS_RADIUS / LANCZOS_TABLE_SIZE));
+ }
+ lanczos_table_init_done = true;
+}
+
+float lanczos_weight_cached(float x)
+{
+ x = fabs(x);
+ if (x > LANCZOS_RADIUS) {
+ return 0.0f;
+ }
+ float table_pos = x * (LANCZOS_TABLE_SIZE / LANCZOS_RADIUS);
+ int table_pos_int = int(table_pos); // Truncate towards zero.
+ float table_pos_frac = table_pos - table_pos_int;
+ assert(table_pos < LANCZOS_TABLE_SIZE + 2);
+ return lanczos_table[table_pos_int] +
+ table_pos_frac * (lanczos_table[table_pos_int + 1] - lanczos_table[table_pos_int]);
+}
+
// Euclid's algorithm, from Wikipedia.
unsigned gcd(unsigned a, unsigned b)
{
}
template<class DestFloat>
-unsigned combine_samples(const Tap<float> *src, Tap<DestFloat> *dst, unsigned src_size, unsigned num_src_samples, unsigned max_samples_saved)
+unsigned combine_samples(const Tap<float> *src, Tap<DestFloat> *dst, float num_subtexels, float inv_num_subtexels, unsigned num_src_samples, unsigned max_samples_saved)
{
// Cut off near-zero values at both sides.
unsigned num_samples_saved = 0;
float pos2 = src[i + 1].pos;
assert(pos2 > pos1);
- fp16_int_t pos, total_weight;
+ DestFloat pos, total_weight;
float sum_sq_error;
- combine_two_samples(w1, w2, pos1, pos2, src_size, &pos, &total_weight, &sum_sq_error);
+ combine_two_samples(w1, w2, pos1, pos2, num_subtexels, inv_num_subtexels, &pos, &total_weight, &sum_sq_error);
// If the interpolation error is larger than that of about sqrt(2) of
// a level at 8-bit precision, don't combine. (You'd think 1.0 was enough,
void normalize_sum(Tap<T>* vals, unsigned num)
{
for (int normalize_pass = 0; normalize_pass < 2; ++normalize_pass) {
- double sum = 0.0;
+ float sum = 0.0;
for (unsigned i = 0; i < num; ++i) {
- sum += to_fp64(vals[i].weight);
+ sum += to_fp32(vals[i].weight);
}
+ float inv_sum = 1.0 / sum;
for (unsigned i = 0; i < num; ++i) {
- vals[i].weight = from_fp64<T>(to_fp64(vals[i].weight) / sum);
+ vals[i].weight = from_fp32<T>(to_fp32(vals[i].weight) * inv_sum);
}
}
}
template<class DestFloat>
unsigned combine_many_samples(const Tap<float> *weights, unsigned src_size, unsigned src_samples, unsigned dst_samples, Tap<DestFloat> **bilinear_weights)
{
- int src_bilinear_samples = 0;
- for (unsigned y = 0; y < dst_samples; ++y) {
- unsigned num_samples_saved = combine_samples<DestFloat>(weights + y * src_samples, NULL, src_size, src_samples, UINT_MAX);
- src_bilinear_samples = max<int>(src_bilinear_samples, src_samples - num_samples_saved);
+ float num_subtexels = src_size / movit_texel_subpixel_precision;
+ float inv_num_subtexels = movit_texel_subpixel_precision / src_size;
+
+ unsigned max_samples_saved = UINT_MAX;
+ for (unsigned y = 0; y < dst_samples && max_samples_saved > 0; ++y) {
+ unsigned num_samples_saved = combine_samples<DestFloat>(weights + y * src_samples, NULL, num_subtexels, inv_num_subtexels, src_samples, max_samples_saved);
+ max_samples_saved = min(max_samples_saved, num_samples_saved);
}
// Now that we know the right width, actually combine the samples.
+ unsigned src_bilinear_samples = src_samples - max_samples_saved;
*bilinear_weights = new Tap<DestFloat>[dst_samples * src_bilinear_samples];
for (unsigned y = 0; y < dst_samples; ++y) {
Tap<DestFloat> *bilinear_weights_ptr = *bilinear_weights + y * src_bilinear_samples;
unsigned num_samples_saved = combine_samples(
weights + y * src_samples,
bilinear_weights_ptr,
- src_size,
+ num_subtexels,
+ inv_num_subtexels,
src_samples,
- src_samples - src_bilinear_samples);
- assert(int(src_samples) - int(num_samples_saved) == src_bilinear_samples);
+ max_samples_saved);
+ assert(num_samples_saved == max_samples_saved);
normalize_sum(bilinear_weights_ptr, src_bilinear_samples);
}
return src_bilinear_samples;
// Find the effective range of the bilinear-optimized kernel.
// Due to rounding of the positions, this is not necessarily the same
// as the intended range (ie., the range of the original weights).
- int lower_pos = int(floor(to_fp64(bilinear_weights[0].pos) * size - 0.5));
- int upper_pos = int(ceil(to_fp64(bilinear_weights[num_bilinear_weights - 1].pos) * size - 0.5)) + 2;
+ int lower_pos = int(floor(to_fp32(bilinear_weights[0].pos) * size - 0.5));
+ int upper_pos = int(ceil(to_fp32(bilinear_weights[num_bilinear_weights - 1].pos) * size - 0.5)) + 2;
lower_pos = min<int>(lower_pos, lrintf(weights[0].pos * size - 0.5));
upper_pos = max<int>(upper_pos, lrintf(weights[num_weights - 1].pos * size - 0.5) + 1);
// Now find the effective weights that result from this sampling.
for (unsigned i = 0; i < num_bilinear_weights; ++i) {
- const float pixel_pos = to_fp64(bilinear_weights[i].pos) * size - 0.5f;
+ const float pixel_pos = to_fp32(bilinear_weights[i].pos) * size - 0.5f;
const int x0 = int(floor(pixel_pos)) - lower_pos;
const int x1 = x0 + 1;
const float f = lrintf((pixel_pos - (x0 + lower_pos)) / movit_texel_subpixel_precision) * movit_texel_subpixel_precision;
assert(x0 < upper_pos - lower_pos);
assert(x1 < upper_pos - lower_pos);
- effective_weights[x0] += to_fp64(bilinear_weights[i].weight) * (1.0 - f);
- effective_weights[x1] += to_fp64(bilinear_weights[i].weight) * f;
+ effective_weights[x0] += to_fp32(bilinear_weights[i].weight) * (1.0 - f);
+ effective_weights[x1] += to_fp32(bilinear_weights[i].weight) * f;
}
// Subtract the desired weights to get the error.
register_float("offset", &offset);
register_float("zoom", &zoom);
register_uniform_sampler2d("sample_tex", &uniform_sample_tex);
- register_uniform_int("num_samples", &uniform_num_samples); // FIXME: What about GLSL pre-1.30?
+ register_uniform_int("num_samples", &uniform_num_samples);
register_uniform_float("num_loops", &uniform_num_loops);
register_uniform_float("slice_height", &uniform_slice_height);
register_uniform_float("sample_x_scale", &uniform_sample_x_scale);
register_uniform_float("whole_pixel_offset", &uniform_whole_pixel_offset);
glGenTextures(1, &texnum);
+
+ if (!lanczos_table_init_done) {
+ // Could in theory race between two threads if we are unlucky,
+ // but that is harmless, since they'll write the same data.
+ init_lanczos_table();
+ }
}
SingleResamplePassEffect::~SingleResamplePassEffect()
// Now sample <int_radius> pixels on each side around that point.
for (int i = 0; i < src_samples; ++i) {
int src_y = base_src_y + i - int_radius;
- float weight = lanczos_weight(radius_scaling_factor * (src_y - center_src_y - subpixel_offset), LANCZOS_RADIUS);
+ float weight = lanczos_weight_cached(radius_scaling_factor * (src_y - center_src_y - subpixel_offset));
weights[y * src_samples + i].weight = weight * radius_scaling_factor;
weights[y * src_samples + i].pos = (src_y + 0.5) / float(src_size);
}
// We specifically do not want mipmaps on the input texture;
// they break minification.
Node *self = chain->find_node_for_effect(this);
- glActiveTexture(chain->get_input_sampler(self, 0));
- check_error();
- glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
- check_error();
+ if (chain->has_input_sampler(self, 0)) {
+ glActiveTexture(chain->get_input_sampler(self, 0));
+ check_error();
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
+ check_error();
+ }
}
} // namespace movit
projects / movit / blobdiff