1/*
2 * This file is part of FFmpeg.
3 *
4 * FFmpeg is free software; you can redistribute it and/or
5 * modify it under the terms of the GNU Lesser General Public
6 * License as published by the Free Software Foundation; either
7 * version 2.1 of the License, or (at your option) any later version.
8 *
9 * FFmpeg is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12 * Lesser General Public License for more details.
13 *
14 * You should have received a copy of the GNU Lesser General Public
15 * License along with FFmpeg; if not, write to the Free Software
16 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17 */
18
19/**
20 * @file
21 *@brief IntraX8 frame subdecoder image manipulation routines
22 */
23
24#include "intrax8dsp.h"
25#include "libavutil/common.h"
26
27/*
28 * area positions, #3 is 1 pixel only, other are 8 pixels
29 *    |66666666|
30 *   3|44444444|55555555|
31 * - -+--------+--------+
32 * 1 2|XXXXXXXX|
33 * 1 2|XXXXXXXX|
34 * 1 2|XXXXXXXX|
35 * 1 2|XXXXXXXX|
36 * 1 2|XXXXXXXX|
37 * 1 2|XXXXXXXX|
38 * 1 2|XXXXXXXX|
39 * 1 2|XXXXXXXX|
40 * ^-start
41 */
42
43#define area1 (0)
44#define area2 (8)
45#define area3 (8 + 8)
46#define area4 (8 + 8 + 1)
47#define area5 (8 + 8 + 1 + 8)
48#define area6 (8 + 8 + 1 + 16)
49
50/**
51 Collect statistics and prepare the edge pixels required by the other spatial compensation functions.
52
53 * @param src pointer to the beginning of the processed block
54 * @param dst pointer to emu_edge, edge pixels are stored the way other compensation routines do.
55 * @param linesize byte offset between 2 vertical pixels in the source image
56 * @param range pointer to the variable where the edge pixel range is to be stored (max-min values)
57 * @param psum  pointer to the variable where the edge pixel sum is to be stored
58 * @param edges Informs this routine that the block is on an image border, so it has to interpolate the missing edge pixels.
59                and some of the edge pixels should be interpolated, the flag has the following meaning:
60                1   - mb_x==0 - first block in the row, interpolate area #1,#2,#3;
61                2   - mb_y==0 - first row, interpolate area #3,#4,#5,#6;
62        note:   1|2 - mb_x==mb_y==0 - first block, use 0x80 value for all areas;
63                4   - mb_x>= (mb_width-1) last block in the row, interpolate area #5;
64-*/
65static void x8_setup_spatial_compensation(uint8_t *src, uint8_t *dst,
66                                          ptrdiff_t stride, int *range,
67                                          int *psum, int edges)
68{
69    uint8_t *ptr;
70    int sum;
71    int i;
72    int min_pix, max_pix;
73    uint8_t c;
74
75    if ((edges & 3) == 3) {
76        *psum  = 0x80 * (8 + 1 + 8 + 2);
77        *range = 0;
78        memset(dst, 0x80, 16 + 1 + 16 + 8);
79        /* this triggers flat_dc for sure. flat_dc avoids all (other)
80         * prediction modes, but requires dc_level decoding. */
81        return;
82    }
83
84    min_pix = 256;
85    max_pix = -1;
86
87    sum = 0;
88
89    if (!(edges & 1)) { // (mb_x != 0) // there is previous block on this row
90        ptr = src - 1; // left column, area 2
91        for (i = 7; i >= 0; i--) {
92            c              = *(ptr - 1); // area1, same mb as area2, no need to check
93            dst[area1 + i] = c;
94            c              = *ptr;
95
96            sum           += c;
97            min_pix        = FFMIN(min_pix, c);
98            max_pix        = FFMAX(max_pix, c);
99            dst[area2 + i] = c;
100
101            ptr += stride;
102        }
103    }
104
105    if (!(edges & 2)) { // (mb_y != 0) // there is row above
106        ptr = src - stride; // top line
107        for (i = 0; i < 8; i++) {
108            c       = *(ptr + i);
109            sum    += c;
110            min_pix = FFMIN(min_pix, c);
111            max_pix = FFMAX(max_pix, c);
112        }
113        if (edges & 4) { // last block on the row?
114            memset(dst + area5, c, 8); // set with last pixel fr
115            memcpy(dst + area4, ptr, 8);
116        } else {
117            memcpy(dst + area4, ptr, 16); // both area4 and 5
118        }
119        // area6 always present in the above block
120        memcpy(dst + area6, ptr - stride, 8);
121    }
122    // now calculate the stuff we need
123    if (edges & 3) { // mb_x ==0 || mb_y == 0) {
124        int avg = (sum + 4) >> 3;
125
126        if (edges & 1) // (mb_x == 0) { // implies mb_y !=0
127            memset(dst + area1, avg, 8 + 8 + 1); // areas 1, 2, 3 are averaged
128        else // implies y == 0 x != 0
129            memset(dst + area3, avg, 1 + 16 + 8); // areas 3, 4, 5, 6
130
131        sum += avg * 9;
132    } else {
133        // the edge pixel, in the top line and left column
134        uint8_t c = *(src - 1 - stride);
135        dst[area3] = c;
136        sum       += c;
137        // edge pixel is not part of min/max
138    }
139    *range = max_pix - min_pix;
140    sum   += *(dst + area5) + *(dst + area5 + 1);
141    *psum  = sum;
142}
143
144static const uint16_t zero_prediction_weights[64 * 2] = {
145    640,  640, 669,  480, 708,  354, 748, 257,
146    792,  198, 760,  143, 808,  101, 772,  72,
147    480,  669, 537,  537, 598,  416, 661, 316,
148    719,  250, 707,  185, 768,  134, 745,  97,
149    354,  708, 416,  598, 488,  488, 564, 388,
150    634,  317, 642,  241, 716,  179, 706, 132,
151    257,  748, 316,  661, 388,  564, 469, 469,
152    543,  395, 571,  311, 655,  238, 660, 180,
153    198,  792, 250,  719, 317,  634, 395, 543,
154    469,  469, 507,  380, 597,  299, 616, 231,
155    161,  855, 206,  788, 266,  710, 340, 623,
156    411,  548, 455,  455, 548,  366, 576, 288,
157    122,  972, 159,  914, 211,  842, 276, 758,
158    341,  682, 389,  584, 483,  483, 520, 390,
159    110, 1172, 144, 1107, 193, 1028, 254, 932,
160    317,  846, 366,  731, 458,  611, 499, 499,
161};
162
163static void spatial_compensation_0(uint8_t *src, uint8_t *dst, ptrdiff_t stride)
164{
165    int i, j;
166    int x, y;
167    unsigned int p; // power divided by 2
168    int a;
169    uint16_t left_sum[2][8] = { { 0 } };
170    uint16_t  top_sum[2][8] = { { 0 } };
171
172    for (i = 0; i < 8; i++) {
173        a = src[area2 + 7 - i] << 4;
174        for (j = 0; j < 8; j++) {
175            p                   = abs(i - j);
176            left_sum[p & 1][j] += a >> (p >> 1);
177        }
178    }
179
180    for (i = 0; i < 8; i++) {
181        a = src[area4 + i] << 4;
182        for (j = 0; j < 8; j++) {
183            p                  = abs(i - j);
184            top_sum[p & 1][j] += a >> (p >> 1);
185        }
186    }
187    for (; i < 10; i++) {
188        a = src[area4 + i] << 4;
189        for (j = 5; j < 8; j++) {
190            p                  = abs(i - j);
191            top_sum[p & 1][j] += a >> (p >> 1);
192        }
193    }
194    for (; i < 12; i++) {
195        a = src[area4 + i] << 4;
196        for (j = 7; j < 8; j++) {
197            p                  = abs(i - j);
198            top_sum[p & 1][j] += a >> (p >> 1);
199        }
200    }
201
202    for (i = 0; i < 8; i++) {
203        top_sum[0][i]  +=  (top_sum[1][i] * 181 + 128) >> 8; // 181 is sqrt(2)/2
204        left_sum[0][i] += (left_sum[1][i] * 181 + 128) >> 8;
205    }
206    for (y = 0; y < 8; y++) {
207        for (x = 0; x < 8; x++)
208            dst[x] = ((uint32_t)  top_sum[0][x] * zero_prediction_weights[y * 16 + x * 2 + 0] +
209                      (uint32_t) left_sum[0][y] * zero_prediction_weights[y * 16 + x * 2 + 1] +
210                      0x8000) >> 16;
211        dst += stride;
212    }
213}
214
215static void spatial_compensation_1(uint8_t *src, uint8_t *dst, ptrdiff_t stride)
216{
217    int x, y;
218
219    for (y = 0; y < 8; y++) {
220        for (x = 0; x < 8; x++)
221            dst[x] = src[area4 + FFMIN(2 * y + x + 2, 15)];
222        dst += stride;
223    }
224}
225
226static void spatial_compensation_2(uint8_t *src, uint8_t *dst, ptrdiff_t stride)
227{
228    int x, y;
229
230    for (y = 0; y < 8; y++) {
231        for (x = 0; x < 8; x++)
232            dst[x] = src[area4 + 1 + y + x];
233        dst += stride;
234    }
235}
236
237static void spatial_compensation_3(uint8_t *src, uint8_t *dst, ptrdiff_t stride)
238{
239    int x, y;
240
241    for (y = 0; y < 8; y++) {
242        for (x = 0; x < 8; x++)
243            dst[x] = src[area4 + ((y + 1) >> 1) + x];
244        dst += stride;
245    }
246}
247
248static void spatial_compensation_4(uint8_t *src, uint8_t *dst, ptrdiff_t stride)
249{
250    int x, y;
251
252    for (y = 0; y < 8; y++) {
253        for (x = 0; x < 8; x++)
254            dst[x] = (src[area4 + x] + src[area6 + x] + 1) >> 1;
255        dst += stride;
256    }
257}
258
259static void spatial_compensation_5(uint8_t *src, uint8_t *dst, ptrdiff_t stride)
260{
261    int x, y;
262
263    for (y = 0; y < 8; y++) {
264        for (x = 0; x < 8; x++) {
265            if (2 * x - y < 0)
266                dst[x] = src[area2 + 9 + 2 * x - y];
267            else
268                dst[x] = src[area4 + x - ((y + 1) >> 1)];
269        }
270        dst += stride;
271    }
272}
273
274static void spatial_compensation_6(uint8_t *src, uint8_t *dst, ptrdiff_t stride)
275{
276    int x, y;
277
278    for (y = 0; y < 8; y++) {
279        for (x = 0; x < 8; x++)
280            dst[x] = src[area3 + x - y];
281        dst += stride;
282    }
283}
284
285static void spatial_compensation_7(uint8_t *src, uint8_t *dst, ptrdiff_t stride)
286{
287    int x, y;
288
289    for (y = 0; y < 8; y++) {
290        for (x = 0; x < 8; x++) {
291            if (x - 2 * y > 0)
292                dst[x] = (src[area3 - 1 + x - 2 * y] + src[area3 + x - 2 * y] + 1) >> 1;
293            else
294                dst[x] = src[area2 + 8 - y + (x >> 1)];
295        }
296        dst += stride;
297    }
298}
299
300static void spatial_compensation_8(uint8_t *src, uint8_t *dst, ptrdiff_t stride)
301{
302    int x, y;
303
304    for (y = 0; y < 8; y++) {
305        for (x = 0; x < 8; x++)
306            dst[x] = (src[area1 + 7 - y] + src[area2 + 7 - y] + 1) >> 1;
307        dst += stride;
308    }
309}
310
311static void spatial_compensation_9(uint8_t *src, uint8_t *dst, ptrdiff_t stride)
312{
313    int x, y;
314
315    for (y = 0; y < 8; y++) {
316        for (x = 0; x < 8; x++)
317            dst[x] = src[area2 + 6 - FFMIN(x + y, 6)];
318        dst += stride;
319    }
320}
321
322static void spatial_compensation_10(uint8_t *src, uint8_t *dst, ptrdiff_t stride)
323{
324    int x, y;
325
326    for (y = 0; y < 8; y++) {
327        for (x = 0; x < 8; x++)
328            dst[x] = (src[area2 + 7 - y] * (8 - x) + src[area4 + x] * x + 4) >> 3;
329        dst += stride;
330    }
331}
332
333static void spatial_compensation_11(uint8_t *src, uint8_t *dst, ptrdiff_t stride)
334{
335    int x, y;
336
337    for (y = 0; y < 8; y++) {
338        for (x = 0; x < 8; x++)
339            dst[x] = (src[area2 + 7 - y] * y + src[area4 + x] * (8 - y) + 4) >> 3;
340        dst += stride;
341    }
342}
343
344static void x8_loop_filter(uint8_t *ptr, const ptrdiff_t a_stride,
345                           const ptrdiff_t b_stride, int quant)
346{
347    int i, t;
348    int p0, p1, p2, p3, p4, p5, p6, p7, p8, p9;
349    int ql = (quant + 10) >> 3;
350
351    for (i = 0; i < 8; i++, ptr += b_stride) {
352        p0 = ptr[-5 * a_stride];
353        p1 = ptr[-4 * a_stride];
354        p2 = ptr[-3 * a_stride];
355        p3 = ptr[-2 * a_stride];
356        p4 = ptr[-1 * a_stride];
357        p5 = ptr[0];
358        p6 = ptr[1 * a_stride];
359        p7 = ptr[2 * a_stride];
360        p8 = ptr[3 * a_stride];
361        p9 = ptr[4 * a_stride];
362
363        t = (FFABS(p1 - p2) <= ql) +
364            (FFABS(p2 - p3) <= ql) +
365            (FFABS(p3 - p4) <= ql) +
366            (FFABS(p4 - p5) <= ql);
367
368        // You need at least 1 to be able to reach a total score of 6.
369        if (t > 0) {
370            t += (FFABS(p5 - p6) <= ql) +
371                 (FFABS(p6 - p7) <= ql) +
372                 (FFABS(p7 - p8) <= ql) +
373                 (FFABS(p8 - p9) <= ql) +
374                 (FFABS(p0 - p1) <= ql);
375            if (t >= 6) {
376                int min, max;
377
378                min = max = p1;
379                min = FFMIN(min, p3);
380                max = FFMAX(max, p3);
381                min = FFMIN(min, p5);
382                max = FFMAX(max, p5);
383                min = FFMIN(min, p8);
384                max = FFMAX(max, p8);
385                if (max - min < 2 * quant) { // early stop
386                    min = FFMIN(min, p2);
387                    max = FFMAX(max, p2);
388                    min = FFMIN(min, p4);
389                    max = FFMAX(max, p4);
390                    min = FFMIN(min, p6);
391                    max = FFMAX(max, p6);
392                    min = FFMIN(min, p7);
393                    max = FFMAX(max, p7);
394                    if (max - min < 2 * quant) {
395                        ptr[-2 * a_stride] = (4 * p2 + 3 * p3 + 1 * p7 + 4) >> 3;
396                        ptr[-1 * a_stride] = (3 * p2 + 3 * p4 + 2 * p7 + 4) >> 3;
397                        ptr[0]             = (2 * p2 + 3 * p5 + 3 * p7 + 4) >> 3;
398                        ptr[1 * a_stride]  = (1 * p2 + 3 * p6 + 4 * p7 + 4) >> 3;
399                        continue;
400                    }
401                }
402            }
403        }
404        {
405            int x, x0, x1, x2;
406            int m;
407
408            x0 = (2 * p3 - 5 * p4 + 5 * p5 - 2 * p6 + 4) >> 3;
409            if (FFABS(x0) < quant) {
410                x1 = (2 * p1 - 5 * p2 + 5 * p3 - 2 * p4 + 4) >> 3;
411                x2 = (2 * p5 - 5 * p6 + 5 * p7 - 2 * p8 + 4) >> 3;
412
413                x = FFABS(x0) - FFMIN(FFABS(x1), FFABS(x2));
414                m = p4 - p5;
415
416                if (x > 0 && (m ^ x0) < 0) {
417                    int32_t sign;
418
419                    sign = m >> 31;
420                    m    = (m ^ sign) - sign; // abs(m)
421                    m  >>= 1;
422
423                    x = 5 * x >> 3;
424
425                    if (x > m)
426                        x = m;
427
428                    x = (x ^ sign) - sign;
429
430                    ptr[-1 * a_stride] -= x;
431                    ptr[0]             += x;
432                }
433            }
434        }
435    }
436}
437
438static void x8_h_loop_filter(uint8_t *src, ptrdiff_t stride, int qscale)
439{
440    x8_loop_filter(src, stride, 1, qscale);
441}
442
443static void x8_v_loop_filter(uint8_t *src, ptrdiff_t stride, int qscale)
444{
445    x8_loop_filter(src, 1, stride, qscale);
446}
447
448av_cold void ff_intrax8dsp_init(IntraX8DSPContext *dsp)
449{
450    dsp->h_loop_filter              = x8_h_loop_filter;
451    dsp->v_loop_filter              = x8_v_loop_filter;
452    dsp->setup_spatial_compensation = x8_setup_spatial_compensation;
453    dsp->spatial_compensation[0]    = spatial_compensation_0;
454    dsp->spatial_compensation[1]    = spatial_compensation_1;
455    dsp->spatial_compensation[2]    = spatial_compensation_2;
456    dsp->spatial_compensation[3]    = spatial_compensation_3;
457    dsp->spatial_compensation[4]    = spatial_compensation_4;
458    dsp->spatial_compensation[5]    = spatial_compensation_5;
459    dsp->spatial_compensation[6]    = spatial_compensation_6;
460    dsp->spatial_compensation[7]    = spatial_compensation_7;
461    dsp->spatial_compensation[8]    = spatial_compensation_8;
462    dsp->spatial_compensation[9]    = spatial_compensation_9;
463    dsp->spatial_compensation[10]   = spatial_compensation_10;
464    dsp->spatial_compensation[11]   = spatial_compensation_11;
465}
466