]> git.sesse.net Git - movit/blob - resample_effect.cpp
Some small cleanups after we got rid of GLSL 1.10; we can now unify 1.30 and ES 3...
[movit] / resample_effect.cpp
1 // Three-lobed Lanczos, the most common choice.
2 // Note that if you change this, the accuracy for LANCZOS_TABLE_SIZE
3 // needs to be recomputed.
4 #define LANCZOS_RADIUS 3.0
5
6 #include <epoxy/gl.h>
7 #include <assert.h>
8 #include <limits.h>
9 #include <math.h>
10 #include <stdio.h>
11 #include <algorithm>
12 #include <Eigen/Sparse>
13 #include <Eigen/SparseQR>
14 #include <Eigen/OrderingMethods>
15
16 #include "effect_chain.h"
17 #include "effect_util.h"
18 #include "fp16.h"
19 #include "init.h"
20 #include "resample_effect.h"
21 #include "util.h"
22
23 using namespace Eigen;
24 using namespace std;
25
26 namespace movit {
27
28 namespace {
29
30 template<class T>
31 struct Tap {
32         T weight;
33         T pos;
34 };
35
36 float sinc(float x)
37 {
38         if (fabs(x) < 1e-6) {
39                 return 1.0f - fabs(x);
40         } else {
41                 return sin(x) / x;
42         }
43 }
44
45 float lanczos_weight(float x)
46 {
47         if (fabs(x) > LANCZOS_RADIUS) {
48                 return 0.0f;
49         } else {
50                 return sinc(M_PI * x) * sinc((M_PI / LANCZOS_RADIUS) * x);
51         }
52 }
53
54 // The weight function can be expensive to compute over and over again
55 // (which will happen during e.g. a zoom), but it is also easy to interpolate
56 // linearly. We compute the right half of the function (in the range of
57 // 0..LANCZOS_RADIUS), with two guard elements for easier interpolation, and
58 // linearly interpolate to get our function.
59 //
60 // We want to scale the table so that the maximum error is always smaller
61 // than 1e-6. As per http://www-solar.mcs.st-andrews.ac.uk/~clare/Lectures/num-analysis/Numan_chap3.pdf,
62 // the error for interpolating a function linearly between points [a,b] is
63 //
64 //   e = 1/2 (x-a)(x-b) f''(u_x)
65 //
66 // for some point u_x in [a,b] (where f(x) is our Lanczos function; we're
67 // assuming LANCZOS_RADIUS=3 from here on). Obviously this is bounded by
68 // f''(x) over the entire range. Numeric optimization shows the maximum of
69 // |f''(x)| to be in x=1.09369819474562880, with the value 2.40067758733152381.
70 // So if the steps between consecutive values are called d, we get
71 //
72 //   |e| <= 1/2 (d/2)^2 2.4007
73 //   |e| <= 0.1367 d^2
74 //
75 // Solve for e = 1e-6 yields a step size of 0.0027, which to cover the range
76 // 0..3 needs 1109 steps. We round up to the next power of two, just to be sure.
77 #define LANCZOS_TABLE_SIZE 2048
78 bool lanczos_table_init_done = false;
79 float lanczos_table[LANCZOS_TABLE_SIZE + 2];
80
81 void init_lanczos_table()
82 {
83         for (unsigned i = 0; i < LANCZOS_TABLE_SIZE + 2; ++i) {
84                 lanczos_table[i] = lanczos_weight(float(i) * (LANCZOS_RADIUS / LANCZOS_TABLE_SIZE));
85         }
86         lanczos_table_init_done = true;
87 }
88
89 float lanczos_weight_cached(float x)
90 {
91         if (!lanczos_table_init_done) {
92                 // Could in theory race between two threads if we are unlucky,
93                 // but that is harmless, since they'll write the same data.
94                 init_lanczos_table();
95         }
96         x = fabs(x);
97         if (x > LANCZOS_RADIUS) {
98                 return 0.0f;
99         }
100         float table_pos = x * (LANCZOS_TABLE_SIZE / LANCZOS_RADIUS);
101         int table_pos_int = int(table_pos);  // Truncate towards zero.
102         float table_pos_frac = table_pos - table_pos_int;
103         assert(table_pos < LANCZOS_TABLE_SIZE + 2);
104         return lanczos_table[table_pos_int] +
105                 table_pos_frac * (lanczos_table[table_pos_int + 1] - lanczos_table[table_pos_int]);
106 }
107
108 // Euclid's algorithm, from Wikipedia.
109 unsigned gcd(unsigned a, unsigned b)
110 {
111         while (b != 0) {
112                 unsigned t = b;
113                 b = a % b;
114                 a = t;
115         }
116         return a;
117 }
118
119 template<class DestFloat>
120 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)
121 {
122         // Cut off near-zero values at both sides.
123         unsigned num_samples_saved = 0;
124         while (num_samples_saved < max_samples_saved &&
125                num_src_samples > 0 &&
126                fabs(src[0].weight) < 1e-6) {
127                 ++src;
128                 --num_src_samples;
129                 ++num_samples_saved;
130         }
131         while (num_samples_saved < max_samples_saved &&
132                num_src_samples > 0 &&
133                fabs(src[num_src_samples - 1].weight) < 1e-6) {
134                 --num_src_samples;
135                 ++num_samples_saved;
136         }
137
138         for (unsigned i = 0, j = 0; i < num_src_samples; ++i, ++j) {
139                 // Copy the sample directly; it will be overwritten later if we can combine.
140                 if (dst != NULL) {
141                         dst[j].weight = convert_float<float, DestFloat>(src[i].weight);
142                         dst[j].pos = convert_float<float, DestFloat>(src[i].pos);
143                 }
144
145                 if (i == num_src_samples - 1) {
146                         // Last sample; cannot combine.
147                         continue;
148                 }
149                 assert(num_samples_saved <= max_samples_saved);
150                 if (num_samples_saved == max_samples_saved) {
151                         // We could maybe save more here, but other rows can't, so don't bother.
152                         continue;
153                 }
154
155                 float w1 = src[i].weight;
156                 float w2 = src[i + 1].weight;
157                 if (w1 * w2 < 0.0f) {
158                         // Differing signs; cannot combine.
159                         continue;
160                 }
161
162                 float pos1 = src[i].pos;
163                 float pos2 = src[i + 1].pos;
164                 assert(pos2 > pos1);
165
166                 DestFloat pos, total_weight;
167                 float sum_sq_error;
168                 combine_two_samples(w1, w2, pos1, pos2, num_subtexels, inv_num_subtexels, &pos, &total_weight, &sum_sq_error);
169
170                 // If the interpolation error is larger than that of about sqrt(2) of
171                 // a level at 8-bit precision, don't combine. (You'd think 1.0 was enough,
172                 // but since the artifacts are not really random, they can get quite
173                 // visible. On the other hand, going to 0.25f, I can see no change at
174                 // all with 8-bit output, so it would not seem to be worth it.)
175                 if (sum_sq_error > 0.5f / (255.0f * 255.0f)) {
176                         continue;
177                 }
178
179                 // OK, we can combine this and the next sample.
180                 if (dst != NULL) {
181                         dst[j].weight = total_weight;
182                         dst[j].pos = pos;
183                 }
184
185                 ++i;  // Skip the next sample.
186                 ++num_samples_saved;
187         }
188         return num_samples_saved;
189 }
190
191 // Normalize so that the sum becomes one. Note that we do it twice;
192 // this sometimes helps a tiny little bit when we have many samples.
193 template<class T>
194 void normalize_sum(Tap<T>* vals, unsigned num)
195 {
196         for (int normalize_pass = 0; normalize_pass < 2; ++normalize_pass) {
197                 double sum = 0.0;
198                 for (unsigned i = 0; i < num; ++i) {
199                         sum += to_fp64(vals[i].weight);
200                 }
201                 double inv_sum = 1.0 / sum;
202                 for (unsigned i = 0; i < num; ++i) {
203                         vals[i].weight = from_fp64<T>(to_fp64(vals[i].weight) * inv_sum);
204                 }
205         }
206 }
207
208 // Make use of the bilinear filtering in the GPU to reduce the number of samples
209 // we need to make. This is a bit more complex than BlurEffect since we cannot combine
210 // two neighboring samples if their weights have differing signs, so we first need to
211 // figure out the maximum number of samples. Then, we downconvert all the weights to
212 // that number -- we could have gone for a variable-length system, but this is simpler,
213 // and the gains would probably be offset by the extra cost of checking when to stop.
214 //
215 // The greedy strategy for combining samples is optimal.
216 template<class DestFloat>
217 unsigned combine_many_samples(const Tap<float> *weights, unsigned src_size, unsigned src_samples, unsigned dst_samples, Tap<DestFloat> **bilinear_weights)
218 {
219         float num_subtexels = src_size / movit_texel_subpixel_precision;
220         float inv_num_subtexels = movit_texel_subpixel_precision / src_size;
221
222         unsigned max_samples_saved = UINT_MAX;
223         for (unsigned y = 0; y < dst_samples && max_samples_saved > 0; ++y) {
224                 unsigned num_samples_saved = combine_samples<DestFloat>(weights + y * src_samples, NULL, num_subtexels, inv_num_subtexels, src_samples, max_samples_saved);
225                 max_samples_saved = min(max_samples_saved, num_samples_saved);
226         }
227
228         // Now that we know the right width, actually combine the samples.
229         unsigned src_bilinear_samples = src_samples - max_samples_saved;
230         *bilinear_weights = new Tap<DestFloat>[dst_samples * src_bilinear_samples];
231         for (unsigned y = 0; y < dst_samples; ++y) {
232                 Tap<DestFloat> *bilinear_weights_ptr = *bilinear_weights + y * src_bilinear_samples;
233                 unsigned num_samples_saved = combine_samples(
234                         weights + y * src_samples,
235                         bilinear_weights_ptr,
236                         num_subtexels,
237                         inv_num_subtexels,
238                         src_samples,
239                         max_samples_saved);
240                 assert(num_samples_saved == max_samples_saved);
241                 normalize_sum(bilinear_weights_ptr, src_bilinear_samples);
242         }
243         return src_bilinear_samples;
244 }
245
246 // Compute the sum of squared errors between the ideal weights (which are
247 // assumed to fall exactly on pixel centers) and the weights that result
248 // from sampling at <bilinear_weights>. The primary reason for the difference
249 // is inaccuracy in the sampling positions, both due to limited precision
250 // in storing them (already inherent in sending them in as fp16_int_t)
251 // and in subtexel sampling precision (which we calculate in this function).
252 template<class T>
253 double compute_sum_sq_error(const Tap<float>* weights, unsigned num_weights,
254                             const Tap<T>* bilinear_weights, unsigned num_bilinear_weights,
255                             unsigned size)
256 {
257         // Find the effective range of the bilinear-optimized kernel.
258         // Due to rounding of the positions, this is not necessarily the same
259         // as the intended range (ie., the range of the original weights).
260         int lower_pos = int(floor(to_fp64(bilinear_weights[0].pos) * size - 0.5));
261         int upper_pos = int(ceil(to_fp64(bilinear_weights[num_bilinear_weights - 1].pos) * size - 0.5)) + 2;
262         lower_pos = min<int>(lower_pos, lrintf(weights[0].pos * size - 0.5));
263         upper_pos = max<int>(upper_pos, lrintf(weights[num_weights - 1].pos * size - 0.5) + 1);
264
265         float* effective_weights = new float[upper_pos - lower_pos];
266         for (int i = 0; i < upper_pos - lower_pos; ++i) {
267                 effective_weights[i] = 0.0f;
268         }
269
270         // Now find the effective weights that result from this sampling.
271         for (unsigned i = 0; i < num_bilinear_weights; ++i) {
272                 const float pixel_pos = to_fp64(bilinear_weights[i].pos) * size - 0.5f;
273                 const int x0 = int(floor(pixel_pos)) - lower_pos;
274                 const int x1 = x0 + 1;
275                 const float f = lrintf((pixel_pos - (x0 + lower_pos)) / movit_texel_subpixel_precision) * movit_texel_subpixel_precision;
276
277                 assert(x0 >= 0);
278                 assert(x1 >= 0);
279                 assert(x0 < upper_pos - lower_pos);
280                 assert(x1 < upper_pos - lower_pos);
281
282                 effective_weights[x0] += to_fp64(bilinear_weights[i].weight) * (1.0 - f);
283                 effective_weights[x1] += to_fp64(bilinear_weights[i].weight) * f;
284         }
285
286         // Subtract the desired weights to get the error.
287         for (unsigned i = 0; i < num_weights; ++i) {
288                 const int x = lrintf(weights[i].pos * size - 0.5f) - lower_pos;
289                 assert(x >= 0);
290                 assert(x < upper_pos - lower_pos);
291
292                 effective_weights[x] -= weights[i].weight;
293         }
294
295         double sum_sq_error = 0.0;
296         for (unsigned i = 0; i < num_weights; ++i) {
297                 sum_sq_error += effective_weights[i] * effective_weights[i];
298         }
299
300         delete[] effective_weights;
301         return sum_sq_error;
302 }
303
304 }  // namespace
305
306 ResampleEffect::ResampleEffect()
307         : input_width(1280),
308           input_height(720),
309           offset_x(0.0f), offset_y(0.0f),
310           zoom_x(1.0f), zoom_y(1.0f),
311           zoom_center_x(0.5f), zoom_center_y(0.5f)
312 {
313         register_int("width", &output_width);
314         register_int("height", &output_height);
315
316         // The first blur pass will forward resolution information to us.
317         hpass = new SingleResamplePassEffect(this);
318         CHECK(hpass->set_int("direction", SingleResamplePassEffect::HORIZONTAL));
319         vpass = new SingleResamplePassEffect(NULL);
320         CHECK(vpass->set_int("direction", SingleResamplePassEffect::VERTICAL));
321
322         update_size();
323 }
324
325 void ResampleEffect::rewrite_graph(EffectChain *graph, Node *self)
326 {
327         Node *hpass_node = graph->add_node(hpass);
328         Node *vpass_node = graph->add_node(vpass);
329         graph->connect_nodes(hpass_node, vpass_node);
330         graph->replace_receiver(self, hpass_node);
331         graph->replace_sender(self, vpass_node);
332         self->disabled = true;
333
334
335 // We get this information forwarded from the first blur pass,
336 // since we are not part of the chain ourselves.
337 void ResampleEffect::inform_input_size(unsigned input_num, unsigned width, unsigned height)
338 {
339         assert(input_num == 0);
340         assert(width != 0);
341         assert(height != 0);
342         input_width = width;
343         input_height = height;
344         update_size();
345 }
346
347 void ResampleEffect::update_size()
348 {
349         bool ok = true;
350         ok |= hpass->set_int("input_width", input_width);
351         ok |= hpass->set_int("input_height", input_height);
352         ok |= hpass->set_int("output_width", output_width);
353         ok |= hpass->set_int("output_height", input_height);
354
355         ok |= vpass->set_int("input_width", output_width);
356         ok |= vpass->set_int("input_height", input_height);
357         ok |= vpass->set_int("output_width", output_width);
358         ok |= vpass->set_int("output_height", output_height);
359
360         assert(ok);
361
362         // The offset added due to zoom may have changed with the size.
363         update_offset_and_zoom();
364 }
365
366 void ResampleEffect::update_offset_and_zoom()
367 {
368         bool ok = true;
369
370         // Zoom from the right origin. (zoom_center is given in normalized coordinates,
371         // i.e. 0..1.)
372         float extra_offset_x = zoom_center_x * (1.0f - 1.0f / zoom_x) * input_width;
373         float extra_offset_y = (1.0f - zoom_center_y) * (1.0f - 1.0f / zoom_y) * input_height;
374
375         ok |= hpass->set_float("offset", extra_offset_x + offset_x);
376         ok |= vpass->set_float("offset", extra_offset_y - offset_y);  // Compensate for the bottom-left origin.
377         ok |= hpass->set_float("zoom", zoom_x);
378         ok |= vpass->set_float("zoom", zoom_y);
379
380         assert(ok);
381 }
382
383 bool ResampleEffect::set_float(const string &key, float value) {
384         if (key == "width") {
385                 output_width = value;
386                 update_size();
387                 return true;
388         }
389         if (key == "height") {
390                 output_height = value;
391                 update_size();
392                 return true;
393         }
394         if (key == "top") {
395                 offset_y = value;
396                 update_offset_and_zoom();
397                 return true;
398         }
399         if (key == "left") {
400                 offset_x = value;
401                 update_offset_and_zoom();
402                 return true;
403         }
404         if (key == "zoom_x") {
405                 if (value <= 0.0f) {
406                         return false;
407                 }
408                 zoom_x = value;
409                 update_offset_and_zoom();
410                 return true;
411         }
412         if (key == "zoom_y") {
413                 if (value <= 0.0f) {
414                         return false;
415                 }
416                 zoom_y = value;
417                 update_offset_and_zoom();
418                 return true;
419         }
420         if (key == "zoom_center_x") {
421                 zoom_center_x = value;
422                 update_offset_and_zoom();
423                 return true;
424         }
425         if (key == "zoom_center_y") {
426                 zoom_center_y = value;
427                 update_offset_and_zoom();
428                 return true;
429         }
430         return false;
431 }
432
433 SingleResamplePassEffect::SingleResamplePassEffect(ResampleEffect *parent)
434         : parent(parent),
435           direction(HORIZONTAL),
436           input_width(1280),
437           input_height(720),
438           offset(0.0),
439           zoom(1.0),
440           last_input_width(-1),
441           last_input_height(-1),
442           last_output_width(-1),
443           last_output_height(-1),
444           last_offset(0.0 / 0.0),  // NaN.
445           last_zoom(0.0 / 0.0),  // NaN.
446           last_texture_width(-1), last_texture_height(-1)
447 {
448         register_int("direction", (int *)&direction);
449         register_int("input_width", &input_width);
450         register_int("input_height", &input_height);
451         register_int("output_width", &output_width);
452         register_int("output_height", &output_height);
453         register_float("offset", &offset);
454         register_float("zoom", &zoom);
455         register_uniform_sampler2d("sample_tex", &uniform_sample_tex);
456         register_uniform_int("num_samples", &uniform_num_samples);
457         register_uniform_float("num_loops", &uniform_num_loops);
458         register_uniform_float("slice_height", &uniform_slice_height);
459         register_uniform_float("sample_x_scale", &uniform_sample_x_scale);
460         register_uniform_float("sample_x_offset", &uniform_sample_x_offset);
461         register_uniform_float("whole_pixel_offset", &uniform_whole_pixel_offset);
462
463         glGenTextures(1, &texnum);
464 }
465
466 SingleResamplePassEffect::~SingleResamplePassEffect()
467 {
468         glDeleteTextures(1, &texnum);
469 }
470
471 string SingleResamplePassEffect::output_fragment_shader()
472 {
473         char buf[256];
474         sprintf(buf, "#define DIRECTION_VERTICAL %d\n", (direction == VERTICAL));
475         return buf + read_file("resample_effect.frag");
476 }
477
478 // Using vertical scaling as an example:
479 //
480 // Generally out[y] = w0 * in[yi] + w1 * in[yi + 1] + w2 * in[yi + 2] + ...
481 //
482 // Obviously, yi will depend on y (in a not-quite-linear way), but so will
483 // the weights w0, w1, w2, etc.. The easiest way of doing this is to encode,
484 // for each sample, the weight and the yi value, e.g. <yi, w0>, <yi + 1, w1>,
485 // and so on. For each y, we encode these along the x-axis (since that is spare),
486 // so out[0] will read from parameters <x,y> = <0,0>, <1,0>, <2,0> and so on.
487 //
488 // For horizontal scaling, we fill in the exact same texture;
489 // the shader just interprets it differently.
490 void SingleResamplePassEffect::update_texture(GLuint glsl_program_num, const string &prefix, unsigned *sampler_num)
491 {
492         unsigned src_size, dst_size;
493         if (direction == SingleResamplePassEffect::HORIZONTAL) {
494                 assert(input_height == output_height);
495                 src_size = input_width;
496                 dst_size = output_width;
497         } else if (direction == SingleResamplePassEffect::VERTICAL) {
498                 assert(input_width == output_width);
499                 src_size = input_height;
500                 dst_size = output_height;
501         } else {
502                 assert(false);
503         }
504
505         // For many resamplings (e.g. 640 -> 1280), we will end up with the same
506         // set of samples over and over again in a loop. Thus, we can compute only
507         // the first such loop, and then ask the card to repeat the texture for us.
508         // This is both easier on the texture cache and lowers our CPU cost for
509         // generating the kernel somewhat.
510         float scaling_factor;
511         if (fabs(zoom - 1.0f) < 1e-6) {
512                 num_loops = gcd(src_size, dst_size);
513                 scaling_factor = float(dst_size) / float(src_size);
514         } else {
515                 // If zooming is enabled (ie., zoom != 1), we turn off the looping.
516                 // We _could_ perhaps do it for rational zoom levels (especially
517                 // things like 2:1), but it doesn't seem to be worth it, given that
518                 // the most common use case would seem to be varying the zoom
519                 // from frame to frame.
520                 num_loops = 1;
521                 scaling_factor = zoom * float(dst_size) / float(src_size);
522         }
523         slice_height = 1.0f / num_loops;
524         unsigned dst_samples = dst_size / num_loops;
525
526         // Sample the kernel in the right place. A diagram with a triangular kernel
527         // (corresponding to linear filtering, and obviously with radius 1)
528         // for easier ASCII art drawing:
529         //
530         //                *
531         //               / \                      |
532         //              /   \                     |
533         //             /     \                    |
534         //    x---x---x   x   x---x---x---x
535         //
536         // Scaling up (in this case, 2x) means sampling more densely:
537         //
538         //                *
539         //               / \                      |
540         //              /   \                     |
541         //             /     \                    |
542         //   x-x-x-x-x-x x x x-x-x-x-x-x-x-x
543         //
544         // When scaling up, any destination pixel will only be influenced by a few
545         // (in this case, two) neighboring pixels, and more importantly, the number
546         // will not be influenced by the scaling factor. (Note, however, that the
547         // pixel centers have moved, due to OpenGL's center-pixel convention.)
548         // The only thing that changes is the weights themselves, as the sampling
549         // points are at different distances from the original pixels.
550         //
551         // Scaling down is a different story:
552         //
553         //                *
554         //               / \                      |
555         //              /   \                     |
556         //             /     \                    |
557         //    --x------ x     --x-------x--
558         //
559         // Again, the pixel centers have moved in a maybe unintuitive fashion,
560         // although when you consider that there are multiple source pixels around,
561         // it's not so bad as at first look:
562         //
563         //            *   *   *   *
564         //           / \ / \ / \ / \              |
565         //          /   X   X   X   \             |
566         //         /   / \ / \ / \   \            |
567         //    --x-------x-------x-------x--
568         //
569         // As you can see, the new pixels become averages of the two neighboring old
570         // ones (the situation for Lanczos is of course more complex).
571         //
572         // Anyhow, in this case we clearly need to look at more source pixels
573         // to compute the destination pixel, and how many depend on the scaling factor.
574         // Thus, the kernel width will vary with how much we scale.
575         float radius_scaling_factor = min(scaling_factor, 1.0f);
576         int int_radius = lrintf(LANCZOS_RADIUS / radius_scaling_factor);
577         int src_samples = int_radius * 2 + 1;
578         Tap<float> *weights = new Tap<float>[dst_samples * src_samples];
579         float subpixel_offset = offset - lrintf(offset);  // The part not covered by whole_pixel_offset.
580         assert(subpixel_offset >= -0.5f && subpixel_offset <= 0.5f);
581         for (unsigned y = 0; y < dst_samples; ++y) {
582                 // Find the point around which we want to sample the source image,
583                 // compensating for differing pixel centers as the scale changes.
584                 float center_src_y = (y + 0.5f) / scaling_factor - 0.5f;
585                 int base_src_y = lrintf(center_src_y);
586
587                 // Now sample <int_radius> pixels on each side around that point.
588                 for (int i = 0; i < src_samples; ++i) {
589                         int src_y = base_src_y + i - int_radius;
590                         float weight = lanczos_weight_cached(radius_scaling_factor * (src_y - center_src_y - subpixel_offset));
591                         weights[y * src_samples + i].weight = weight * radius_scaling_factor;
592                         weights[y * src_samples + i].pos = (src_y + 0.5) / float(src_size);
593                 }
594         }
595
596         // Now make use of the bilinear filtering in the GPU to reduce the number of samples
597         // we need to make. Try fp16 first; if it's not accurate enough, we go to fp32.
598         // Our tolerance level for total error is a bit higher than the one for invididual
599         // samples, since one would assume overall errors in the shape don't matter as much.
600         const float max_error = 2.0f / (255.0f * 255.0f);
601         Tap<fp16_int_t> *bilinear_weights_fp16;
602         src_bilinear_samples = combine_many_samples(weights, src_size, src_samples, dst_samples, &bilinear_weights_fp16);
603         Tap<float> *bilinear_weights_fp32 = NULL;
604         bool fallback_to_fp32 = false;
605         double max_sum_sq_error_fp16 = 0.0;
606         for (unsigned y = 0; y < dst_samples; ++y) {
607                 double sum_sq_error_fp16 = compute_sum_sq_error(
608                         weights + y * src_samples, src_samples,
609                         bilinear_weights_fp16 + y * src_bilinear_samples, src_bilinear_samples,
610                         src_size);
611                 max_sum_sq_error_fp16 = std::max(max_sum_sq_error_fp16, sum_sq_error_fp16);
612                 if (max_sum_sq_error_fp16 > max_error) {
613                         break;
614                 }
615         }
616
617         if (max_sum_sq_error_fp16 > max_error) {
618                 fallback_to_fp32 = true;
619                 src_bilinear_samples = combine_many_samples(weights, src_size, src_samples, dst_samples, &bilinear_weights_fp32);
620         }
621
622         // Encode as a two-component texture. Note the GL_REPEAT.
623         glActiveTexture(GL_TEXTURE0 + *sampler_num);
624         check_error();
625         glBindTexture(GL_TEXTURE_2D, texnum);
626         check_error();
627         if (last_texture_width == -1) {
628                 // Need to set this state the first time.
629                 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
630                 check_error();
631                 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
632                 check_error();
633                 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
634                 check_error();
635         }
636
637         GLenum type, internal_format;
638         void *pixels;
639         if (fallback_to_fp32) {
640                 type = GL_FLOAT;
641                 internal_format = GL_RG32F;
642                 pixels = bilinear_weights_fp32;
643         } else {
644                 type = GL_HALF_FLOAT;
645                 internal_format = GL_RG16F;
646                 pixels = bilinear_weights_fp16;
647         }
648
649         if (int(src_bilinear_samples) == last_texture_width &&
650             int(dst_samples) == last_texture_height &&
651             internal_format == last_texture_internal_format) {
652                 // Texture dimensions and type are unchanged; it is more efficient
653                 // to just update it rather than making an entirely new texture.
654                 glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, src_bilinear_samples, dst_samples, GL_RG, type, pixels);
655         } else {
656                 glTexImage2D(GL_TEXTURE_2D, 0, internal_format, src_bilinear_samples, dst_samples, 0, GL_RG, type, pixels);
657                 last_texture_width = src_bilinear_samples;
658                 last_texture_height = dst_samples;
659                 last_texture_internal_format = internal_format;
660         }
661         check_error();
662
663         delete[] weights;
664         delete[] bilinear_weights_fp16;
665         delete[] bilinear_weights_fp32;
666 }
667
668 void SingleResamplePassEffect::set_gl_state(GLuint glsl_program_num, const string &prefix, unsigned *sampler_num)
669 {
670         Effect::set_gl_state(glsl_program_num, prefix, sampler_num);
671
672         assert(input_width > 0);
673         assert(input_height > 0);
674         assert(output_width > 0);
675         assert(output_height > 0);
676
677         if (input_width != last_input_width ||
678             input_height != last_input_height ||
679             output_width != last_output_width ||
680             output_height != last_output_height ||
681             offset != last_offset ||
682             zoom != last_zoom) {
683                 update_texture(glsl_program_num, prefix, sampler_num);
684                 last_input_width = input_width;
685                 last_input_height = input_height;
686                 last_output_width = output_width;
687                 last_output_height = output_height;
688                 last_offset = offset;
689                 last_zoom = zoom;
690         }
691
692         glActiveTexture(GL_TEXTURE0 + *sampler_num);
693         check_error();
694         glBindTexture(GL_TEXTURE_2D, texnum);
695         check_error();
696
697         uniform_sample_tex = *sampler_num;
698         ++*sampler_num;
699         uniform_num_samples = src_bilinear_samples;
700         uniform_num_loops = num_loops;
701         uniform_slice_height = slice_height;
702
703         // Instructions for how to convert integer sample numbers to positions in the weight texture.
704         uniform_sample_x_scale = 1.0f / src_bilinear_samples;
705         uniform_sample_x_offset = 0.5f / src_bilinear_samples;
706
707         if (direction == SingleResamplePassEffect::VERTICAL) {
708                 uniform_whole_pixel_offset = lrintf(offset) / float(input_height);
709         } else {
710                 uniform_whole_pixel_offset = lrintf(offset) / float(input_width);
711         }
712
713         // We specifically do not want mipmaps on the input texture;
714         // they break minification.
715         Node *self = chain->find_node_for_effect(this);
716         glActiveTexture(chain->get_input_sampler(self, 0));
717         check_error();
718         glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
719         check_error();
720 }
721
722 }  // namespace movit