4dc3392fe24b186eb31468ddddcd148779ed1f1b
[libav.git] / libavcodec / snow.c
1 /*
2 * Copyright (C) 2004 Michael Niedermayer <michaelni@gmx.at>
3 *
4 * This library 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 of the License, or (at your option) any later version.
8 *
9 * This library 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 this library; if not, write to the Free Software
16 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 */
18
19 #include "avcodec.h"
20 #include "common.h"
21 #include "dsputil.h"
22 #include "cabac.h"
23
24 #include "mpegvideo.h"
25
26 #undef NDEBUG
27 #include <assert.h>
28
29 #define MAX_DECOMPOSITIONS 8
30 #define MAX_PLANES 4
31 #define DWTELEM int
32 #define QROOT 8
33
34 static const int8_t quant3[256]={
35 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
36 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
37 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
38 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
39 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
40 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
41 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
42 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
43 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
44 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
45 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
46 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
47 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
48 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
49 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
50 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 0,
51 };
52 static const int8_t quant3b[256]={
53 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
54 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
55 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
56 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
57 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
58 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
59 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
60 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
61 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
62 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
63 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
64 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
65 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
66 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
67 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
68 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
69 };
70 static const int8_t quant5[256]={
71 0, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
72 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
73 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
74 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
75 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
76 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
77 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
78 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
79 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
80 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
81 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
82 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
83 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
84 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
85 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
86 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-1,-1,-1,
87 };
88 static const int8_t quant7[256]={
89 0, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
90 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
91 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
92 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
93 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
94 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
95 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
96 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
97 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
98 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
99 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
100 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
101 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
102 -3,-3,-3,-3,-3,-3,-3,-3,-3,-2,-2,-2,-2,-2,-2,-2,
103 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
104 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-1,-1,
105 };
106 static const int8_t quant9[256]={
107 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3,
108 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
109 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
110 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
111 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
112 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
113 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
114 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
115 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
116 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
117 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
118 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
119 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
120 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
121 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-3,-3,-3,-3,
122 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-2,-2,-2,-2,-1,-1,
123 };
124 static const int8_t quant11[256]={
125 0, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4,
126 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
127 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
128 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
129 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
130 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
131 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
132 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
133 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
134 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
135 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
136 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
137 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
138 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-4,-4,
139 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
140 -4,-4,-4,-4,-4,-3,-3,-3,-3,-3,-3,-3,-2,-2,-2,-1,
141 };
142 static const int8_t quant13[256]={
143 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4,
144 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
145 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
146 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
147 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
148 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
149 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
150 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
151 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,
152 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,
153 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,
154 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,
155 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-5,
156 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
157 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
158 -4,-4,-4,-4,-4,-4,-4,-4,-4,-3,-3,-3,-3,-2,-2,-1,
159 };
160
161 #define OBMC_MAX 64
162 #if 0 //64*cubic
163 static const uint8_t obmc32[1024]={
164 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
165 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0,
166 0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 2, 2, 1, 1, 1, 1, 0, 0, 0,
167 0, 0, 1, 1, 2, 2, 3, 4, 4, 5, 6, 6, 7, 7, 8, 8, 8, 8, 7, 7, 6, 6, 5, 4, 4, 3, 2, 2, 1, 1, 0, 0,
168 0, 0, 1, 2, 2, 3, 4, 6, 7, 8, 9,10,11,12,12,12,12,12,12,11,10, 9, 8, 7, 6, 4, 3, 2, 2, 1, 0, 0,
169 0, 1, 1, 2, 3, 5, 6, 8,10,11,13,14,15,16,17,18,18,17,16,15,14,13,11,10, 8, 6, 5, 3, 2, 1, 1, 0,
170 0, 1, 1, 3, 4, 6, 8,10,13,15,17,19,20,22,22,23,23,22,22,20,19,17,15,13,10, 8, 6, 4, 3, 1, 1, 0,
171 0, 1, 2, 4, 6, 8,10,13,16,19,21,23,25,27,28,29,29,28,27,25,23,21,19,16,13,10, 8, 6, 4, 2, 1, 0,
172 0, 1, 2, 4, 7,10,13,16,19,22,25,28,31,33,34,35,35,34,33,31,28,25,22,19,16,13,10, 7, 4, 2, 1, 0,
173 0, 1, 3, 5, 8,11,15,19,22,26,30,33,36,38,40,41,41,40,38,36,33,30,26,22,19,15,11, 8, 5, 3, 1, 0,
174 0, 1, 3, 6, 9,12,17,21,25,30,34,38,41,44,45,46,46,45,44,41,38,34,30,25,21,17,12, 9, 6, 3, 1, 0,
175 0, 1, 3, 6,10,14,19,23,28,33,38,42,45,48,51,52,52,51,48,45,42,38,33,28,23,19,14,10, 6, 3, 1, 0,
176 0, 1, 4, 7,11,15,20,25,31,36,41,45,49,52,55,56,56,55,52,49,45,41,36,31,25,20,15,11, 7, 4, 1, 0,
177 0, 2, 4, 7,12,16,22,27,33,38,44,48,52,56,58,60,60,58,56,52,48,44,38,33,27,22,16,12, 7, 4, 2, 0,
178 0, 1, 4, 8,12,17,22,28,34,40,45,51,55,58,61,62,62,61,58,55,51,45,40,34,28,22,17,12, 8, 4, 1, 0,
179 0, 2, 4, 8,12,18,23,29,35,41,46,52,56,60,62,64,64,62,60,56,52,46,41,35,29,23,18,12, 8, 4, 2, 0,
180 0, 2, 4, 8,12,18,23,29,35,41,46,52,56,60,62,64,64,62,60,56,52,46,41,35,29,23,18,12, 8, 4, 2, 0,
181 0, 1, 4, 8,12,17,22,28,34,40,45,51,55,58,61,62,62,61,58,55,51,45,40,34,28,22,17,12, 8, 4, 1, 0,
182 0, 2, 4, 7,12,16,22,27,33,38,44,48,52,56,58,60,60,58,56,52,48,44,38,33,27,22,16,12, 7, 4, 2, 0,
183 0, 1, 4, 7,11,15,20,25,31,36,41,45,49,52,55,56,56,55,52,49,45,41,36,31,25,20,15,11, 7, 4, 1, 0,
184 0, 1, 3, 6,10,14,19,23,28,33,38,42,45,48,51,52,52,51,48,45,42,38,33,28,23,19,14,10, 6, 3, 1, 0,
185 0, 1, 3, 6, 9,12,17,21,25,30,34,38,41,44,45,46,46,45,44,41,38,34,30,25,21,17,12, 9, 6, 3, 1, 0,
186 0, 1, 3, 5, 8,11,15,19,22,26,30,33,36,38,40,41,41,40,38,36,33,30,26,22,19,15,11, 8, 5, 3, 1, 0,
187 0, 1, 2, 4, 7,10,13,16,19,22,25,28,31,33,34,35,35,34,33,31,28,25,22,19,16,13,10, 7, 4, 2, 1, 0,
188 0, 1, 2, 4, 6, 8,10,13,16,19,21,23,25,27,28,29,29,28,27,25,23,21,19,16,13,10, 8, 6, 4, 2, 1, 0,
189 0, 1, 1, 3, 4, 6, 8,10,13,15,17,19,20,22,22,23,23,22,22,20,19,17,15,13,10, 8, 6, 4, 3, 1, 1, 0,
190 0, 1, 1, 2, 3, 5, 6, 8,10,11,13,14,15,16,17,18,18,17,16,15,14,13,11,10, 8, 6, 5, 3, 2, 1, 1, 0,
191 0, 0, 1, 2, 2, 3, 4, 6, 7, 8, 9,10,11,12,12,12,12,12,12,11,10, 9, 8, 7, 6, 4, 3, 2, 2, 1, 0, 0,
192 0, 0, 1, 1, 2, 2, 3, 4, 4, 5, 6, 6, 7, 7, 8, 8, 8, 8, 7, 7, 6, 6, 5, 4, 4, 3, 2, 2, 1, 1, 0, 0,
193 0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 2, 2, 1, 1, 1, 1, 0, 0, 0,
194 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0,
195 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
196 //error:0.000022
197 };
198 static const uint8_t obmc16[256]={
199 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
200 0, 1, 1, 2, 4, 5, 5, 6, 6, 5, 5, 4, 2, 1, 1, 0,
201 0, 1, 4, 6, 9,11,13,15,15,13,11, 9, 6, 4, 1, 0,
202 0, 2, 6,11,15,20,24,26,26,24,20,15,11, 6, 2, 0,
203 0, 4, 9,15,23,29,34,38,38,34,29,23,15, 9, 4, 0,
204 0, 5,11,20,29,38,45,49,49,45,38,29,20,11, 5, 0,
205 1, 5,13,24,34,45,53,57,57,53,45,34,24,13, 5, 1,
206 1, 6,15,26,38,49,57,62,62,57,49,38,26,15, 6, 1,
207 1, 6,15,26,38,49,57,62,62,57,49,38,26,15, 6, 1,
208 1, 5,13,24,34,45,53,57,57,53,45,34,24,13, 5, 1,
209 0, 5,11,20,29,38,45,49,49,45,38,29,20,11, 5, 0,
210 0, 4, 9,15,23,29,34,38,38,34,29,23,15, 9, 4, 0,
211 0, 2, 6,11,15,20,24,26,26,24,20,15,11, 6, 2, 0,
212 0, 1, 4, 6, 9,11,13,15,15,13,11, 9, 6, 4, 1, 0,
213 0, 1, 1, 2, 4, 5, 5, 6, 6, 5, 5, 4, 2, 1, 1, 0,
214 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
215 //error:0.000033
216 };
217 #elif 1 // 64*linear
218 static const uint8_t obmc32[1024]={
219 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0,
220 0, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 2, 2, 2, 1, 1, 1, 0,
221 0, 1, 2, 2, 3, 3, 4, 5, 5, 6, 7, 7, 8, 8, 9,10,10, 9, 8, 8, 7, 7, 6, 5, 5, 4, 3, 3, 2, 2, 1, 0,
222 0, 1, 2, 3, 4, 5, 6, 7, 7, 8, 9,10,11,12,13,14,14,13,12,11,10, 9, 8, 7, 7, 6, 5, 4, 3, 2, 1, 0,
223 1, 2, 3, 4, 5, 6, 7, 8,10,11,12,13,14,15,16,17,17,16,15,14,13,12,11,10, 8, 7, 6, 5, 4, 3, 2, 1,
224 1, 2, 3, 5, 6, 8, 9,10,12,13,14,16,17,19,20,21,21,20,19,17,16,14,13,12,10, 9, 8, 6, 5, 3, 2, 1,
225 1, 2, 4, 6, 7, 9,11,12,14,15,17,19,20,22,24,25,25,24,22,20,19,17,15,14,12,11, 9, 7, 6, 4, 2, 1,
226 1, 3, 5, 7, 8,10,12,14,16,18,20,22,23,25,27,29,29,27,25,23,22,20,18,16,14,12,10, 8, 7, 5, 3, 1,
227 1, 3, 5, 7,10,12,14,16,18,20,22,24,27,29,31,33,33,31,29,27,24,22,20,18,16,14,12,10, 7, 5, 3, 1,
228 1, 4, 6, 8,11,13,15,18,20,23,25,27,30,32,34,37,37,34,32,30,27,25,23,20,18,15,13,11, 8, 6, 4, 1,
229 1, 4, 7, 9,12,14,17,20,22,25,28,30,33,35,38,41,41,38,35,33,30,28,25,22,20,17,14,12, 9, 7, 4, 1,
230 1, 4, 7,10,13,16,19,22,24,27,30,33,36,39,42,45,45,42,39,36,33,30,27,24,22,19,16,13,10, 7, 4, 1,
231 2, 5, 8,11,14,17,20,23,27,30,33,36,39,42,45,48,48,45,42,39,36,33,30,27,23,20,17,14,11, 8, 5, 2,
232 2, 5, 8,12,15,19,22,25,29,32,35,39,42,46,49,52,52,49,46,42,39,35,32,29,25,22,19,15,12, 8, 5, 2,
233 2, 5, 9,13,16,20,24,27,31,34,38,42,45,49,53,56,56,53,49,45,42,38,34,31,27,24,20,16,13, 9, 5, 2,
234 2, 6,10,14,17,21,25,29,33,37,41,45,48,52,56,60,60,56,52,48,45,41,37,33,29,25,21,17,14,10, 6, 2,
235 2, 6,10,14,17,21,25,29,33,37,41,45,48,52,56,60,60,56,52,48,45,41,37,33,29,25,21,17,14,10, 6, 2,
236 2, 5, 9,13,16,20,24,27,31,34,38,42,45,49,53,56,56,53,49,45,42,38,34,31,27,24,20,16,13, 9, 5, 2,
237 2, 5, 8,12,15,19,22,25,29,32,35,39,42,46,49,52,52,49,46,42,39,35,32,29,25,22,19,15,12, 8, 5, 2,
238 2, 5, 8,11,14,17,20,23,27,30,33,36,39,42,45,48,48,45,42,39,36,33,30,27,23,20,17,14,11, 8, 5, 2,
239 1, 4, 7,10,13,16,19,22,24,27,30,33,36,39,42,45,45,42,39,36,33,30,27,24,22,19,16,13,10, 7, 4, 1,
240 1, 4, 7, 9,12,14,17,20,22,25,28,30,33,35,38,41,41,38,35,33,30,28,25,22,20,17,14,12, 9, 7, 4, 1,
241 1, 4, 6, 8,11,13,15,18,20,23,25,27,30,32,34,37,37,34,32,30,27,25,23,20,18,15,13,11, 8, 6, 4, 1,
242 1, 3, 5, 7,10,12,14,16,18,20,22,24,27,29,31,33,33,31,29,27,24,22,20,18,16,14,12,10, 7, 5, 3, 1,
243 1, 3, 5, 7, 8,10,12,14,16,18,20,22,23,25,27,29,29,27,25,23,22,20,18,16,14,12,10, 8, 7, 5, 3, 1,
244 1, 2, 4, 6, 7, 9,11,12,14,15,17,19,20,22,24,25,25,24,22,20,19,17,15,14,12,11, 9, 7, 6, 4, 2, 1,
245 1, 2, 3, 5, 6, 8, 9,10,12,13,14,16,17,19,20,21,21,20,19,17,16,14,13,12,10, 9, 8, 6, 5, 3, 2, 1,
246 1, 2, 3, 4, 5, 6, 7, 8,10,11,12,13,14,15,16,17,17,16,15,14,13,12,11,10, 8, 7, 6, 5, 4, 3, 2, 1,
247 0, 1, 2, 3, 4, 5, 6, 7, 7, 8, 9,10,11,12,13,14,14,13,12,11,10, 9, 8, 7, 7, 6, 5, 4, 3, 2, 1, 0,
248 0, 1, 2, 2, 3, 3, 4, 5, 5, 6, 7, 7, 8, 8, 9,10,10, 9, 8, 8, 7, 7, 6, 5, 5, 4, 3, 3, 2, 2, 1, 0,
249 0, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 2, 2, 2, 1, 1, 1, 0,
250 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0,
251 //error:0.000020
252 };
253 static const uint8_t obmc16[256]={
254 0, 1, 1, 2, 2, 3, 3, 4, 4, 3, 3, 2, 2, 1, 1, 0,
255 1, 2, 4, 5, 7, 8,10,11,11,10, 8, 7, 5, 4, 2, 1,
256 1, 4, 6, 9,11,14,16,19,19,16,14,11, 9, 6, 4, 1,
257 2, 5, 9,12,16,19,23,26,26,23,19,16,12, 9, 5, 2,
258 2, 7,11,16,20,25,29,34,34,29,25,20,16,11, 7, 2,
259 3, 8,14,19,25,30,36,41,41,36,30,25,19,14, 8, 3,
260 3,10,16,23,29,36,42,49,49,42,36,29,23,16,10, 3,
261 4,11,19,26,34,41,49,56,56,49,41,34,26,19,11, 4,
262 4,11,19,26,34,41,49,56,56,49,41,34,26,19,11, 4,
263 3,10,16,23,29,36,42,49,49,42,36,29,23,16,10, 3,
264 3, 8,14,19,25,30,36,41,41,36,30,25,19,14, 8, 3,
265 2, 7,11,16,20,25,29,34,34,29,25,20,16,11, 7, 2,
266 2, 5, 9,12,16,19,23,26,26,23,19,16,12, 9, 5, 2,
267 1, 4, 6, 9,11,14,16,19,19,16,14,11, 9, 6, 4, 1,
268 1, 2, 4, 5, 7, 8,10,11,11,10, 8, 7, 5, 4, 2, 1,
269 0, 1, 1, 2, 2, 3, 3, 4, 4, 3, 3, 2, 2, 1, 1, 0,
270 //error:0.000015
271 };
272 #else //64*cos
273 static const uint8_t obmc32[1024]={
274 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
275 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
276 0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 1, 1, 1, 1, 0, 0, 0,
277 0, 0, 1, 1, 1, 2, 2, 3, 4, 5, 5, 6, 7, 7, 7, 7, 7, 7, 7, 7, 6, 5, 5, 4, 3, 2, 2, 1, 1, 1, 0, 0,
278 0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 9,10,11,11,12,12,12,12,11,11,10, 9, 7, 6, 5, 4, 3, 2, 1, 1, 0, 0,
279 0, 0, 1, 2, 3, 5, 6, 8, 9,11,12,14,15,16,17,17,17,17,16,15,14,12,11, 9, 8, 6, 5, 3, 2, 1, 0, 0,
280 0, 1, 1, 2, 4, 6, 8,10,12,15,17,19,20,21,22,23,23,22,21,20,19,17,15,12,10, 8, 6, 4, 2, 1, 1, 0,
281 0, 1, 2, 3, 5, 8,10,13,16,19,21,24,26,27,28,29,29,28,27,26,24,21,19,16,13,10, 8, 5, 3, 2, 1, 0,
282 0, 1, 2, 4, 6, 9,12,16,19,23,26,29,31,33,34,35,35,34,33,31,29,26,23,19,16,12, 9, 6, 4, 2, 1, 0,
283 0, 1, 3, 5, 7,11,15,19,23,26,30,34,37,39,40,41,41,40,39,37,34,30,26,23,19,15,11, 7, 5, 3, 1, 0,
284 0, 1, 3, 5, 9,12,17,21,26,30,35,38,42,44,46,47,47,46,44,42,38,35,30,26,21,17,12, 9, 5, 3, 1, 0,
285 0, 1, 3, 6, 9,14,19,24,29,34,38,43,46,49,51,52,52,51,49,46,43,38,34,29,24,19,14, 9, 6, 3, 1, 0,
286 0, 1, 3, 6,11,15,20,26,31,37,42,46,50,53,56,57,57,56,53,50,46,42,37,31,26,20,15,11, 6, 3, 1, 0,
287 0, 1, 3, 7,11,16,21,27,33,39,44,49,53,57,59,60,60,59,57,53,49,44,39,33,27,21,16,11, 7, 3, 1, 0,
288 0, 1, 4, 7,12,17,22,28,34,40,46,51,56,59,61,63,63,61,59,56,51,46,40,34,28,22,17,12, 7, 4, 1, 0,
289 0, 1, 4, 7,12,17,23,29,35,41,47,52,57,60,63,64,64,63,60,57,52,47,41,35,29,23,17,12, 7, 4, 1, 0,
290 0, 1, 4, 7,12,17,23,29,35,41,47,52,57,60,63,64,64,63,60,57,52,47,41,35,29,23,17,12, 7, 4, 1, 0,
291 0, 1, 4, 7,12,17,22,28,34,40,46,51,56,59,61,63,63,61,59,56,51,46,40,34,28,22,17,12, 7, 4, 1, 0,
292 0, 1, 3, 7,11,16,21,27,33,39,44,49,53,57,59,60,60,59,57,53,49,44,39,33,27,21,16,11, 7, 3, 1, 0,
293 0, 1, 3, 6,11,15,20,26,31,37,42,46,50,53,56,57,57,56,53,50,46,42,37,31,26,20,15,11, 6, 3, 1, 0,
294 0, 1, 3, 6, 9,14,19,24,29,34,38,43,46,49,51,52,52,51,49,46,43,38,34,29,24,19,14, 9, 6, 3, 1, 0,
295 0, 1, 3, 5, 9,12,17,21,26,30,35,38,42,44,46,47,47,46,44,42,38,35,30,26,21,17,12, 9, 5, 3, 1, 0,
296 0, 1, 3, 5, 7,11,15,19,23,26,30,34,37,39,40,41,41,40,39,37,34,30,26,23,19,15,11, 7, 5, 3, 1, 0,
297 0, 1, 2, 4, 6, 9,12,16,19,23,26,29,31,33,34,35,35,34,33,31,29,26,23,19,16,12, 9, 6, 4, 2, 1, 0,
298 0, 1, 2, 3, 5, 8,10,13,16,19,21,24,26,27,28,29,29,28,27,26,24,21,19,16,13,10, 8, 5, 3, 2, 1, 0,
299 0, 1, 1, 2, 4, 6, 8,10,12,15,17,19,20,21,22,23,23,22,21,20,19,17,15,12,10, 8, 6, 4, 2, 1, 1, 0,
300 0, 0, 1, 2, 3, 5, 6, 8, 9,11,12,14,15,16,17,17,17,17,16,15,14,12,11, 9, 8, 6, 5, 3, 2, 1, 0, 0,
301 0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 9,10,11,11,12,12,12,12,11,11,10, 9, 7, 6, 5, 4, 3, 2, 1, 1, 0, 0,
302 0, 0, 1, 1, 1, 2, 2, 3, 4, 5, 5, 6, 7, 7, 7, 7, 7, 7, 7, 7, 6, 5, 5, 4, 3, 2, 2, 1, 1, 1, 0, 0,
303 0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 1, 1, 1, 1, 0, 0, 0,
304 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
305 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
306 //error:0.000022
307 };
308 static const uint8_t obmc16[256]={
309 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0,
310 0, 0, 1, 2, 3, 4, 5, 5, 5, 5, 4, 3, 2, 1, 0, 0,
311 0, 1, 3, 6, 8,11,13,14,14,13,11, 8, 6, 3, 1, 0,
312 0, 2, 6,10,15,20,24,26,26,24,20,15,10, 6, 2, 0,
313 0, 3, 8,16,23,30,35,38,38,35,30,23,16, 8, 3, 0,
314 1, 4,11,20,30,39,46,49,49,46,39,30,20,11, 4, 1,
315 1, 5,13,24,35,46,54,58,58,54,46,35,24,13, 5, 1,
316 0, 5,14,26,38,49,58,63,63,58,49,38,26,14, 5, 0,
317 0, 5,14,26,38,49,58,63,63,58,49,38,26,14, 5, 0,
318 1, 5,13,24,35,46,54,58,58,54,46,35,24,13, 5, 1,
319 1, 4,11,20,30,39,46,49,49,46,39,30,20,11, 4, 1,
320 0, 3, 8,16,23,30,35,38,38,35,30,23,16, 8, 3, 0,
321 0, 2, 6,10,15,20,24,26,26,24,20,15,10, 6, 2, 0,
322 0, 1, 3, 6, 8,11,13,14,14,13,11, 8, 6, 3, 1, 0,
323 0, 0, 1, 2, 3, 4, 5, 5, 5, 5, 4, 3, 2, 1, 0, 0,
324 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0,
325 //error:0.000022
326 };
327 #endif
328
329
330 typedef struct SubBand{
331 int level;
332 int stride;
333 int width;
334 int height;
335 int qlog; ///< log(qscale)/log[2^(1/6)]
336 DWTELEM *buf;
337 struct SubBand *parent;
338 uint8_t state[/*7*2*/ 7 + 512][32];
339 }SubBand;
340
341 typedef struct Plane{
342 int width;
343 int height;
344 SubBand band[MAX_DECOMPOSITIONS][4];
345 }Plane;
346
347 typedef struct SnowContext{
348 // MpegEncContext m; // needed for motion estimation, should not be used for anything else, the idea is to make the motion estimation eventually independant of MpegEncContext, so this will be removed then (FIXME/XXX)
349
350 AVCodecContext *avctx;
351 CABACContext c;
352 DSPContext dsp;
353 AVFrame input_picture;
354 AVFrame current_picture;
355 AVFrame last_picture;
356 AVFrame mconly_picture;
357 // uint8_t q_context[16];
358 uint8_t header_state[32];
359 int keyframe;
360 int version;
361 int spatial_decomposition_type;
362 int temporal_decomposition_type;
363 int spatial_decomposition_count;
364 int temporal_decomposition_count;
365 DWTELEM *spatial_dwt_buffer;
366 DWTELEM *pred_buffer;
367 int colorspace_type;
368 int chroma_h_shift;
369 int chroma_v_shift;
370 int spatial_scalability;
371 int qlog;
372 int mv_scale;
373 int qbias;
374 #define QBIAS_SHIFT 3
375 int b_width; //FIXME remove?
376 int b_height; //FIXME remove?
377 Plane plane[MAX_PLANES];
378 SubBand mb_band;
379 SubBand mv_band[2];
380
381 uint16_t *mb_type;
382 uint8_t *mb_mean;
383 uint32_t *dummy;
384 int16_t (*motion_val8)[2];
385 int16_t (*motion_val16)[2];
386 MpegEncContext m; // needed for motion estimation, should not be used for anything else, the idea is to make the motion estimation eventually independant of MpegEncContext, so this will be removed then (FIXME/XXX)
387 }SnowContext;
388
389 #define QEXPSHIFT 7 //FIXME try to change this to 0
390 static const uint8_t qexp[8]={
391 128, 140, 152, 166, 181, 197, 215, 235
392 // 64, 70, 76, 83, 91, 99, 108, 117
393 // 32, 35, 38, 41, 45, 49, 54, 59
394 // 16, 17, 19, 21, 23, 25, 27, 29
395 // 8, 9, 10, 10, 11, 12, 13, 15
396 };
397
398 static inline int mirror(int v, int m){
399 if (v<0) return -v;
400 else if(v>m) return 2*m-v;
401 else return v;
402 }
403
404 static inline void put_symbol(CABACContext *c, uint8_t *state, int v, int is_signed){
405 int i;
406
407 if(v){
408 const int a= ABS(v);
409 const int e= av_log2(a);
410 #if 1
411 const int el= FFMIN(e, 10);
412 put_cabac(c, state+0, 0);
413
414 for(i=0; i<el; i++){
415 put_cabac(c, state+1+i, 1); //1..10
416 }
417 for(; i<e; i++){
418 put_cabac(c, state+1+9, 1); //1..10
419 }
420 put_cabac(c, state+1+FFMIN(i,9), 0);
421
422 for(i=e-1; i>=el; i--){
423 put_cabac(c, state+22+9, (a>>i)&1); //22..31
424 }
425 for(; i>=0; i--){
426 put_cabac(c, state+22+i, (a>>i)&1); //22..31
427 }
428
429 if(is_signed)
430 put_cabac(c, state+11 + el, v < 0); //11..21
431 #else
432
433 put_cabac(c, state+0, 0);
434 if(e<=9){
435 for(i=0; i<e; i++){
436 put_cabac(c, state+1+i, 1); //1..10
437 }
438 put_cabac(c, state+1+i, 0);
439
440 for(i=e-1; i>=0; i--){
441 put_cabac(c, state+22+i, (a>>i)&1); //22..31
442 }
443
444 if(is_signed)
445 put_cabac(c, state+11 + e, v < 0); //11..21
446 }else{
447 for(i=0; i<e; i++){
448 put_cabac(c, state+1+FFMIN(i,9), 1); //1..10
449 }
450 put_cabac(c, state+1+FFMIN(i,9), 0);
451
452 for(i=e-1; i>=0; i--){
453 put_cabac(c, state+22+FFMIN(i,9), (a>>i)&1); //22..31
454 }
455
456 if(is_signed)
457 put_cabac(c, state+11 + FFMIN(e,10), v < 0); //11..21
458 }
459 #endif
460 }else{
461 put_cabac(c, state+0, 1);
462 }
463 }
464
465 static inline int get_symbol(CABACContext *c, uint8_t *state, int is_signed){
466 if(get_cabac(c, state+0))
467 return 0;
468 else{
469 int i, e, a, el;
470 //FIXME try to merge loops with FFMIN() maybe they are equally fast and they are surly cuter
471 for(e=0; e<10; e++){
472 if(get_cabac(c, state + 1 + e)==0) // 1..10
473 break;
474 }
475 el= e;
476
477 if(e==10){
478 while(get_cabac(c, state + 1 + 9)) //10
479 e++;
480 }
481 a= 1;
482 for(i=e-1; i>=el; i--){
483 a += a + get_cabac(c, state+22+9); //31
484 }
485 for(; i>=0; i--){
486 a += a + get_cabac(c, state+22+i); //22..31
487 }
488
489 if(is_signed && get_cabac(c, state+11 + el)) //11..21
490 return -a;
491 else
492 return a;
493 }
494 }
495
496 static always_inline void lift(DWTELEM *dst, DWTELEM *src, DWTELEM *ref, int dst_step, int src_step, int ref_step, int width, int mul, int add, int shift, int highpass, int inverse){
497 const int mirror_left= !highpass;
498 const int mirror_right= (width&1) ^ highpass;
499 const int w= (width>>1) - 1 + (highpass & width);
500 int i;
501
502 #define LIFT(src, ref, inv) ((src) + ((inv) ? - (ref) : + (ref)))
503 if(mirror_left){
504 dst[0] = LIFT(src[0], ((mul*2*ref[0]+add)>>shift), inverse);
505 dst += dst_step;
506 src += src_step;
507 }
508
509 for(i=0; i<w; i++){
510 dst[i*dst_step] = LIFT(src[i*src_step], ((mul*(ref[i*ref_step] + ref[(i+1)*ref_step])+add)>>shift), inverse);
511 }
512
513 if(mirror_right){
514 dst[w*dst_step] = LIFT(src[w*src_step], ((mul*2*ref[w*ref_step]+add)>>shift), inverse);
515 }
516 }
517
518 static always_inline void lift5(DWTELEM *dst, DWTELEM *src, DWTELEM *ref, int dst_step, int src_step, int ref_step, int width, int mul, int add, int shift, int highpass, int inverse){
519 const int mirror_left= !highpass;
520 const int mirror_right= (width&1) ^ highpass;
521 const int w= (width>>1) - 1 + (highpass & width);
522 int i;
523
524 if(mirror_left){
525 int r= 3*2*ref[0];
526 r += r>>4;
527 r += r>>8;
528 dst[0] = LIFT(src[0], ((r+add)>>shift), inverse);
529 dst += dst_step;
530 src += src_step;
531 }
532
533 for(i=0; i<w; i++){
534 int r= 3*(ref[i*ref_step] + ref[(i+1)*ref_step]);
535 r += r>>4;
536 r += r>>8;
537 dst[i*dst_step] = LIFT(src[i*src_step], ((r+add)>>shift), inverse);
538 }
539
540 if(mirror_right){
541 int r= 3*2*ref[w*ref_step];
542 r += r>>4;
543 r += r>>8;
544 dst[w*dst_step] = LIFT(src[w*src_step], ((r+add)>>shift), inverse);
545 }
546 }
547
548
549 static void inplace_lift(int *dst, int width, int *coeffs, int n, int shift, int start, int inverse){
550 int x, i;
551
552 for(x=start; x<width; x+=2){
553 int64_t sum=0;
554
555 for(i=0; i<n; i++){
556 int x2= x + 2*i - n + 1;
557 if (x2< 0) x2= -x2;
558 else if(x2>=width) x2= 2*width-x2-2;
559 sum += coeffs[i]*(int64_t)dst[x2];
560 }
561 if(inverse) dst[x] -= (sum + (1<<shift)/2)>>shift;
562 else dst[x] += (sum + (1<<shift)/2)>>shift;
563 }
564 }
565
566 static void inplace_liftV(int *dst, int width, int height, int stride, int *coeffs, int n, int shift, int start, int inverse){
567 int x, y, i;
568 for(y=start; y<height; y+=2){
569 for(x=0; x<width; x++){
570 int64_t sum=0;
571
572 for(i=0; i<n; i++){
573 int y2= y + 2*i - n + 1;
574 if (y2< 0) y2= -y2;
575 else if(y2>=height) y2= 2*height-y2-2;
576 sum += coeffs[i]*(int64_t)dst[x + y2*stride];
577 }
578 if(inverse) dst[x + y*stride] -= (sum + (1<<shift)/2)>>shift;
579 else dst[x + y*stride] += (sum + (1<<shift)/2)>>shift;
580 }
581 }
582 }
583
584 #define SCALEX 1
585 #define LX0 0
586 #define LX1 1
587
588 #if 0 // more accurate 9/7
589 #define N1 2
590 #define SHIFT1 14
591 #define COEFFS1 (int[]){-25987,-25987}
592 #define N2 2
593 #define SHIFT2 19
594 #define COEFFS2 (int[]){-27777,-27777}
595 #define N3 2
596 #define SHIFT3 15
597 #define COEFFS3 (int[]){28931,28931}
598 #define N4 2
599 #define SHIFT4 15
600 #define COEFFS4 (int[]){14533,14533}
601 #elif 1 // 13/7 CRF
602 #define N1 4
603 #define SHIFT1 4
604 #define COEFFS1 (int[]){1,-9,-9,1}
605 #define N2 4
606 #define SHIFT2 4
607 #define COEFFS2 (int[]){-1,5,5,-1}
608 #define N3 0
609 #define SHIFT3 1
610 #define COEFFS3 NULL
611 #define N4 0
612 #define SHIFT4 1
613 #define COEFFS4 NULL
614 #elif 1 // 3/5
615 #define LX0 1
616 #define LX1 0
617 #define SCALEX 0.5
618 #define N1 2
619 #define SHIFT1 1
620 #define COEFFS1 (int[]){1,1}
621 #define N2 2
622 #define SHIFT2 2
623 #define COEFFS2 (int[]){-1,-1}
624 #define N3 0
625 #define SHIFT3 0
626 #define COEFFS3 NULL
627 #define N4 0
628 #define SHIFT4 0
629 #define COEFFS4 NULL
630 #elif 1 // 11/5
631 #define N1 0
632 #define SHIFT1 1
633 #define COEFFS1 NULL
634 #define N2 2
635 #define SHIFT2 2
636 #define COEFFS2 (int[]){-1,-1}
637 #define N3 2
638 #define SHIFT3 0
639 #define COEFFS3 (int[]){-1,-1}
640 #define N4 4
641 #define SHIFT4 7
642 #define COEFFS4 (int[]){-5,29,29,-5}
643 #define SCALEX 4
644 #elif 1 // 9/7 CDF
645 #define N1 2
646 #define SHIFT1 7
647 #define COEFFS1 (int[]){-203,-203}
648 #define N2 2
649 #define SHIFT2 12
650 #define COEFFS2 (int[]){-217,-217}
651 #define N3 2
652 #define SHIFT3 7
653 #define COEFFS3 (int[]){113,113}
654 #define N4 2
655 #define SHIFT4 9
656 #define COEFFS4 (int[]){227,227}
657 #define SCALEX 1
658 #elif 1 // 7/5 CDF
659 #define N1 0
660 #define SHIFT1 1
661 #define COEFFS1 NULL
662 #define N2 2
663 #define SHIFT2 2
664 #define COEFFS2 (int[]){-1,-1}
665 #define N3 2
666 #define SHIFT3 0
667 #define COEFFS3 (int[]){-1,-1}
668 #define N4 2
669 #define SHIFT4 4
670 #define COEFFS4 (int[]){3,3}
671 #elif 1 // 9/7 MN
672 #define N1 4
673 #define SHIFT1 4
674 #define COEFFS1 (int[]){1,-9,-9,1}
675 #define N2 2
676 #define SHIFT2 2
677 #define COEFFS2 (int[]){1,1}
678 #define N3 0
679 #define SHIFT3 1
680 #define COEFFS3 NULL
681 #define N4 0
682 #define SHIFT4 1
683 #define COEFFS4 NULL
684 #else // 13/7 CRF
685 #define N1 4
686 #define SHIFT1 4
687 #define COEFFS1 (int[]){1,-9,-9,1}
688 #define N2 4
689 #define SHIFT2 4
690 #define COEFFS2 (int[]){-1,5,5,-1}
691 #define N3 0
692 #define SHIFT3 1
693 #define COEFFS3 NULL
694 #define N4 0
695 #define SHIFT4 1
696 #define COEFFS4 NULL
697 #endif
698 static void horizontal_decomposeX(int *b, int width){
699 int temp[width];
700 const int width2= width>>1;
701 const int w2= (width+1)>>1;
702 int A1,A2,A3,A4, x;
703
704 inplace_lift(b, width, COEFFS1, N1, SHIFT1, LX1, 0);
705 inplace_lift(b, width, COEFFS2, N2, SHIFT2, LX0, 0);
706 inplace_lift(b, width, COEFFS3, N3, SHIFT3, LX1, 0);
707 inplace_lift(b, width, COEFFS4, N4, SHIFT4, LX0, 0);
708
709 for(x=0; x<width2; x++){
710 temp[x ]= b[2*x ];
711 temp[x+w2]= b[2*x + 1];
712 }
713 if(width&1)
714 temp[x ]= b[2*x ];
715 memcpy(b, temp, width*sizeof(int));
716 }
717
718 static void horizontal_composeX(int *b, int width){
719 int temp[width];
720 const int width2= width>>1;
721 int A1,A2,A3,A4, x;
722 const int w2= (width+1)>>1;
723
724 memcpy(temp, b, width*sizeof(int));
725 for(x=0; x<width2; x++){
726 b[2*x ]= temp[x ];
727 b[2*x + 1]= temp[x+w2];
728 }
729 if(width&1)
730 b[2*x ]= temp[x ];
731
732 inplace_lift(b, width, COEFFS4, N4, SHIFT4, LX0, 1);
733 inplace_lift(b, width, COEFFS3, N3, SHIFT3, LX1, 1);
734 inplace_lift(b, width, COEFFS2, N2, SHIFT2, LX0, 1);
735 inplace_lift(b, width, COEFFS1, N1, SHIFT1, LX1, 1);
736 }
737
738 static void spatial_decomposeX(int *buffer, int width, int height, int stride){
739 int x, y;
740
741 for(y=0; y<height; y++){
742 for(x=0; x<width; x++){
743 buffer[y*stride + x] *= SCALEX;
744 }
745 }
746
747 for(y=0; y<height; y++){
748 horizontal_decomposeX(buffer + y*stride, width);
749 }
750
751 inplace_liftV(buffer, width, height, stride, COEFFS1, N1, SHIFT1, LX1, 0);
752 inplace_liftV(buffer, width, height, stride, COEFFS2, N2, SHIFT2, LX0, 0);
753 inplace_liftV(buffer, width, height, stride, COEFFS3, N3, SHIFT3, LX1, 0);
754 inplace_liftV(buffer, width, height, stride, COEFFS4, N4, SHIFT4, LX0, 0);
755 }
756
757 static void spatial_composeX(int *buffer, int width, int height, int stride){
758 int x, y;
759
760 inplace_liftV(buffer, width, height, stride, COEFFS4, N4, SHIFT4, LX0, 1);
761 inplace_liftV(buffer, width, height, stride, COEFFS3, N3, SHIFT3, LX1, 1);
762 inplace_liftV(buffer, width, height, stride, COEFFS2, N2, SHIFT2, LX0, 1);
763 inplace_liftV(buffer, width, height, stride, COEFFS1, N1, SHIFT1, LX1, 1);
764
765 for(y=0; y<height; y++){
766 horizontal_composeX(buffer + y*stride, width);
767 }
768
769 for(y=0; y<height; y++){
770 for(x=0; x<width; x++){
771 buffer[y*stride + x] /= SCALEX;
772 }
773 }
774 }
775
776 static void horizontal_decompose53i(int *b, int width){
777 int temp[width];
778 const int width2= width>>1;
779 int A1,A2,A3,A4, x;
780 const int w2= (width+1)>>1;
781
782 for(x=0; x<width2; x++){
783 temp[x ]= b[2*x ];
784 temp[x+w2]= b[2*x + 1];
785 }
786 if(width&1)
787 temp[x ]= b[2*x ];
788 #if 0
789 A2= temp[1 ];
790 A4= temp[0 ];
791 A1= temp[0+width2];
792 A1 -= (A2 + A4)>>1;
793 A4 += (A1 + 1)>>1;
794 b[0+width2] = A1;
795 b[0 ] = A4;
796 for(x=1; x+1<width2; x+=2){
797 A3= temp[x+width2];
798 A4= temp[x+1 ];
799 A3 -= (A2 + A4)>>1;
800 A2 += (A1 + A3 + 2)>>2;
801 b[x+width2] = A3;
802 b[x ] = A2;
803
804 A1= temp[x+1+width2];
805 A2= temp[x+2 ];
806 A1 -= (A2 + A4)>>1;
807 A4 += (A1 + A3 + 2)>>2;
808 b[x+1+width2] = A1;
809 b[x+1 ] = A4;
810 }
811 A3= temp[width-1];
812 A3 -= A2;
813 A2 += (A1 + A3 + 2)>>2;
814 b[width -1] = A3;
815 b[width2-1] = A2;
816 #else
817 lift(b+w2, temp+w2, temp, 1, 1, 1, width, -1, 0, 1, 1, 0);
818 lift(b , temp , b+w2, 1, 1, 1, width, 1, 2, 2, 0, 0);
819 #endif
820 }
821
822 static void vertical_decompose53iH0(int *b0, int *b1, int *b2, int width){
823 int i;
824
825 for(i=0; i<width; i++){
826 b1[i] -= (b0[i] + b2[i])>>1;
827 }
828 }
829
830 static void vertical_decompose53iL0(int *b0, int *b1, int *b2, int width){
831 int i;
832
833 for(i=0; i<width; i++){
834 b1[i] += (b0[i] + b2[i] + 2)>>2;
835 }
836 }
837
838 static void spatial_decompose53i(int *buffer, int width, int height, int stride){
839 int x, y;
840 DWTELEM *b0= buffer + mirror(-2-1, height-1)*stride;
841 DWTELEM *b1= buffer + mirror(-2 , height-1)*stride;
842
843 for(y=-2; y<height; y+=2){
844 DWTELEM *b2= buffer + mirror(y+1, height-1)*stride;
845 DWTELEM *b3= buffer + mirror(y+2, height-1)*stride;
846
847 {START_TIMER
848 if(b1 <= b3) horizontal_decompose53i(b2, width);
849 if(y+2 < height) horizontal_decompose53i(b3, width);
850 STOP_TIMER("horizontal_decompose53i")}
851
852 {START_TIMER
853 if(b1 <= b3) vertical_decompose53iH0(b1, b2, b3, width);
854 if(b0 <= b2) vertical_decompose53iL0(b0, b1, b2, width);
855 STOP_TIMER("vertical_decompose53i*")}
856
857 b0=b2;
858 b1=b3;
859 }
860 }
861
862 #define lift5 lift
863 #if 1
864 #define W_AM 3
865 #define W_AO 0
866 #define W_AS 1
867
868 #define W_BM 1
869 #define W_BO 8
870 #define W_BS 4
871
872 #undef lift5
873 #define W_CM 9999
874 #define W_CO 2
875 #define W_CS 2
876
877 #define W_DM 15
878 #define W_DO 16
879 #define W_DS 5
880 #elif 0
881 #define W_AM 55
882 #define W_AO 16
883 #define W_AS 5
884
885 #define W_BM 3
886 #define W_BO 32
887 #define W_BS 6
888
889 #define W_CM 127
890 #define W_CO 64
891 #define W_CS 7
892
893 #define W_DM 7
894 #define W_DO 8
895 #define W_DS 4
896 #elif 0
897 #define W_AM 97
898 #define W_AO 32
899 #define W_AS 6
900
901 #define W_BM 63
902 #define W_BO 512
903 #define W_BS 10
904
905 #define W_CM 13
906 #define W_CO 8
907 #define W_CS 4
908
909 #define W_DM 15
910 #define W_DO 16
911 #define W_DS 5
912
913 #else
914
915 #define W_AM 203
916 #define W_AO 64
917 #define W_AS 7
918
919 #define W_BM 217
920 #define W_BO 2048
921 #define W_BS 12
922
923 #define W_CM 113
924 #define W_CO 64
925 #define W_CS 7
926
927 #define W_DM 227
928 #define W_DO 128
929 #define W_DS 9
930 #endif
931 static void horizontal_decompose97i(int *b, int width){
932 int temp[width];
933 const int w2= (width+1)>>1;
934
935 lift (temp+w2, b +1, b , 1, 2, 2, width, -W_AM, W_AO, W_AS, 1, 0);
936 lift (temp , b , temp+w2, 1, 2, 1, width, -W_BM, W_BO, W_BS, 0, 0);
937 lift5(b +w2, temp+w2, temp , 1, 1, 1, width, W_CM, W_CO, W_CS, 1, 0);
938 lift (b , temp , b +w2, 1, 1, 1, width, W_DM, W_DO, W_DS, 0, 0);
939 }
940
941
942 static void vertical_decompose97iH0(int *b0, int *b1, int *b2, int width){
943 int i;
944
945 for(i=0; i<width; i++){
946 b1[i] -= (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
947 }
948 }
949
950 static void vertical_decompose97iH1(int *b0, int *b1, int *b2, int width){
951 int i;
952
953 for(i=0; i<width; i++){
954 #ifdef lift5
955 b1[i] += (W_CM*(b0[i] + b2[i])+W_CO)>>W_CS;
956 #else
957 int r= 3*(b0[i] + b2[i]);
958 r+= r>>4;
959 r+= r>>8;
960 b1[i] += (r+W_CO)>>W_CS;
961 #endif
962 }
963 }
964
965 static void vertical_decompose97iL0(int *b0, int *b1, int *b2, int width){
966 int i;
967
968 for(i=0; i<width; i++){
969 b1[i] -= (W_BM*(b0[i] + b2[i])+W_BO)>>W_BS;
970 }
971 }
972
973 static void vertical_decompose97iL1(int *b0, int *b1, int *b2, int width){
974 int i;
975
976 for(i=0; i<width; i++){
977 b1[i] += (W_DM*(b0[i] + b2[i])+W_DO)>>W_DS;
978 }
979 }
980
981 static void spatial_decompose97i(int *buffer, int width, int height, int stride){
982 int x, y;
983 DWTELEM *b0= buffer + mirror(-4-1, height-1)*stride;
984 DWTELEM *b1= buffer + mirror(-4 , height-1)*stride;
985 DWTELEM *b2= buffer + mirror(-4+1, height-1)*stride;
986 DWTELEM *b3= buffer + mirror(-4+2, height-1)*stride;
987
988 for(y=-4; y<height; y+=2){
989 DWTELEM *b4= buffer + mirror(y+3, height-1)*stride;
990 DWTELEM *b5= buffer + mirror(y+4, height-1)*stride;
991
992 {START_TIMER
993 if(b3 <= b5) horizontal_decompose97i(b4, width);
994 if(y+4 < height) horizontal_decompose97i(b5, width);
995 if(width>400){
996 STOP_TIMER("horizontal_decompose97i")
997 }}
998
999 {START_TIMER
1000 if(b3 <= b5) vertical_decompose97iH0(b3, b4, b5, width);
1001 if(b2 <= b4) vertical_decompose97iL0(b2, b3, b4, width);
1002 if(b1 <= b3) vertical_decompose97iH1(b1, b2, b3, width);
1003 if(b0 <= b2) vertical_decompose97iL1(b0, b1, b2, width);
1004
1005 if(width>400){
1006 STOP_TIMER("vertical_decompose97i")
1007 }}
1008
1009 b0=b2;
1010 b1=b3;
1011 b2=b4;
1012 b3=b5;
1013 }
1014 }
1015
1016 static void spatial_dwt(SnowContext *s, int *buffer, int width, int height, int stride){
1017 int level;
1018
1019 for(level=0; level<s->spatial_decomposition_count; level++){
1020 switch(s->spatial_decomposition_type){
1021 case 0: spatial_decompose97i(buffer, width>>level, height>>level, stride<<level); break;
1022 case 1: spatial_decompose53i(buffer, width>>level, height>>level, stride<<level); break;
1023 case 2: spatial_decomposeX (buffer, width>>level, height>>level, stride<<level); break;
1024 }
1025 }
1026 }
1027
1028 static void horizontal_compose53i(int *b, int width){
1029 int temp[width];
1030 const int width2= width>>1;
1031 const int w2= (width+1)>>1;
1032 int A1,A2,A3,A4, x;
1033
1034 #if 0
1035 A2= temp[1 ];
1036 A4= temp[0 ];
1037 A1= temp[0+width2];
1038 A1 -= (A2 + A4)>>1;
1039 A4 += (A1 + 1)>>1;
1040 b[0+width2] = A1;
1041 b[0 ] = A4;
1042 for(x=1; x+1<width2; x+=2){
1043 A3= temp[x+width2];
1044 A4= temp[x+1 ];
1045 A3 -= (A2 + A4)>>1;
1046 A2 += (A1 + A3 + 2)>>2;
1047 b[x+width2] = A3;
1048 b[x ] = A2;
1049
1050 A1= temp[x+1+width2];
1051 A2= temp[x+2 ];
1052 A1 -= (A2 + A4)>>1;
1053 A4 += (A1 + A3 + 2)>>2;
1054 b[x+1+width2] = A1;
1055 b[x+1 ] = A4;
1056 }
1057 A3= temp[width-1];
1058 A3 -= A2;
1059 A2 += (A1 + A3 + 2)>>2;
1060 b[width -1] = A3;
1061 b[width2-1] = A2;
1062 #else
1063 lift(temp , b , b+w2, 1, 1, 1, width, 1, 2, 2, 0, 1);
1064 lift(temp+w2, b+w2, temp, 1, 1, 1, width, -1, 0, 1, 1, 1);
1065 #endif
1066 for(x=0; x<width2; x++){
1067 b[2*x ]= temp[x ];
1068 b[2*x + 1]= temp[x+w2];
1069 }
1070 if(width&1)
1071 b[2*x ]= temp[x ];
1072 }
1073
1074 static void vertical_compose53iH0(int *b0, int *b1, int *b2, int width){
1075 int i;
1076
1077 for(i=0; i<width; i++){
1078 b1[i] += (b0[i] + b2[i])>>1;
1079 }
1080 }
1081
1082 static void vertical_compose53iL0(int *b0, int *b1, int *b2, int width){
1083 int i;
1084
1085 for(i=0; i<width; i++){
1086 b1[i] -= (b0[i] + b2[i] + 2)>>2;
1087 }
1088 }
1089
1090 static void spatial_compose53i(int *buffer, int width, int height, int stride){
1091 int x, y;
1092 DWTELEM *b0= buffer + mirror(-1-1, height-1)*stride;
1093 DWTELEM *b1= buffer + mirror(-1 , height-1)*stride;
1094
1095 for(y=-1; y<=height; y+=2){
1096 DWTELEM *b2= buffer + mirror(y+1, height-1)*stride;
1097 DWTELEM *b3= buffer + mirror(y+2, height-1)*stride;
1098
1099 {START_TIMER
1100 if(b1 <= b3) vertical_compose53iL0(b1, b2, b3, width);
1101 if(b0 <= b2) vertical_compose53iH0(b0, b1, b2, width);
1102 STOP_TIMER("vertical_compose53i*")}
1103
1104 {START_TIMER
1105 if(y-1 >= 0) horizontal_compose53i(b0, width);
1106 if(b0 <= b2) horizontal_compose53i(b1, width);
1107 STOP_TIMER("horizontal_compose53i")}
1108
1109 b0=b2;
1110 b1=b3;
1111 }
1112 }
1113
1114
1115 static void horizontal_compose97i(int *b, int width){
1116 int temp[width];
1117 const int w2= (width+1)>>1;
1118
1119 lift (temp , b , b +w2, 1, 1, 1, width, W_DM, W_DO, W_DS, 0, 1);
1120 lift5(temp+w2, b +w2, temp , 1, 1, 1, width, W_CM, W_CO, W_CS, 1, 1);
1121 lift (b , temp , temp+w2, 2, 1, 1, width, -W_BM, W_BO, W_BS, 0, 1);
1122 lift (b+1 , temp+w2, b , 2, 1, 2, width, -W_AM, W_AO, W_AS, 1, 1);
1123 }
1124
1125 static void vertical_compose97iH0(int *b0, int *b1, int *b2, int width){
1126 int i;
1127
1128 for(i=0; i<width; i++){
1129 b1[i] += (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
1130 }
1131 }
1132
1133 static void vertical_compose97iH1(int *b0, int *b1, int *b2, int width){
1134 int i;
1135
1136 for(i=0; i<width; i++){
1137 #ifdef lift5
1138 b1[i] -= (W_CM*(b0[i] + b2[i])+W_CO)>>W_CS;
1139 #else
1140 int r= 3*(b0[i] + b2[i]);
1141 r+= r>>4;
1142 r+= r>>8;
1143 b1[i] -= (r+W_CO)>>W_CS;
1144 #endif
1145 }
1146 }
1147
1148 static void vertical_compose97iL0(int *b0, int *b1, int *b2, int width){
1149 int i;
1150
1151 for(i=0; i<width; i++){
1152 b1[i] += (W_BM*(b0[i] + b2[i])+W_BO)>>W_BS;
1153 }
1154 }
1155
1156 static void vertical_compose97iL1(int *b0, int *b1, int *b2, int width){
1157 int i;
1158
1159 for(i=0; i<width; i++){
1160 b1[i] -= (W_DM*(b0[i] + b2[i])+W_DO)>>W_DS;
1161 }
1162 }
1163
1164 static void spatial_compose97i(int *buffer, int width, int height, int stride){
1165 int x, y;
1166 DWTELEM *b0= buffer + mirror(-3-1, height-1)*stride;
1167 DWTELEM *b1= buffer + mirror(-3 , height-1)*stride;
1168 DWTELEM *b2= buffer + mirror(-3+1, height-1)*stride;
1169 DWTELEM *b3= buffer + mirror(-3+2, height-1)*stride;
1170
1171 for(y=-3; y<=height; y+=2){
1172 DWTELEM *b4= buffer + mirror(y+3, height-1)*stride;
1173 DWTELEM *b5= buffer + mirror(y+4, height-1)*stride;
1174
1175 if(stride == width && y+4 < height && 0){
1176 int x;
1177 for(x=0; x<width/2; x++)
1178 b5[x] += 64*2;
1179 for(; x<width; x++)
1180 b5[x] += 169*2;
1181 }
1182
1183 {START_TIMER
1184 if(b3 <= b5) vertical_compose97iL1(b3, b4, b5, width);
1185 if(b2 <= b4) vertical_compose97iH1(b2, b3, b4, width);
1186 if(b1 <= b3) vertical_compose97iL0(b1, b2, b3, width);
1187 if(b0 <= b2) vertical_compose97iH0(b0, b1, b2, width);
1188 if(width>400){
1189 STOP_TIMER("vertical_compose97i")}}
1190
1191 {START_TIMER
1192 if(y-1>= 0) horizontal_compose97i(b0, width);
1193 if(b0 <= b2) horizontal_compose97i(b1, width);
1194 if(width>400 && b0 <= b2){
1195 STOP_TIMER("horizontal_compose97i")}}
1196
1197 b0=b2;
1198 b1=b3;
1199 b2=b4;
1200 b3=b5;
1201 }
1202 }
1203
1204 static void spatial_idwt(SnowContext *s, int *buffer, int width, int height, int stride){
1205 int level;
1206
1207 for(level=s->spatial_decomposition_count-1; level>=0; level--){
1208 switch(s->spatial_decomposition_type){
1209 case 0: spatial_compose97i(buffer, width>>level, height>>level, stride<<level); break;
1210 case 1: spatial_compose53i(buffer, width>>level, height>>level, stride<<level); break;
1211 case 2: spatial_composeX (buffer, width>>level, height>>level, stride<<level); break;
1212 }
1213 }
1214 }
1215
1216 static const int hilbert[16][2]={
1217 {0,0}, {1,0}, {1,1}, {0,1},
1218 {0,2}, {0,3}, {1,3}, {1,2},
1219 {2,2}, {2,3}, {3,3}, {3,2},
1220 {3,1}, {2,1}, {2,0}, {3,0},
1221 };
1222 #if 0
1223 -o o-
1224 | |
1225 o-o
1226
1227 -o-o o-o-
1228 | |
1229 o-o o-o
1230 | |
1231 o o-o o
1232 | | | |
1233 o-o o-o
1234
1235 0112122312232334122323342334
1236 0123456789ABCDEF0123456789AB
1237 RLLRMRRLLRRMRLLMLRRLMLLRRLLM
1238
1239 4 B F 14 1B
1240 4 11 15 20 27
1241
1242 -o o-o-o o-o-o o-
1243 | | | | | |
1244 o-o o-o o-o o-o
1245 | |
1246 o-o o-o o-o o-o
1247 | | | | | |
1248 o o-o-o o-o-o o
1249 | |
1250 o-o o-o-o-o o-o
1251 | | | |
1252 o-o o-o o-o o-o
1253 | | | |
1254 o o-o o o o-o o
1255 | | | | | | | |
1256 o-o o-o o-o o-o
1257
1258 #endif
1259
1260 #define SVI(a, i, x, y) \
1261 {\
1262 a[i][0]= x;\
1263 a[i][1]= y;\
1264 i++;\
1265 }
1266
1267 static int sig_cmp(const void *a, const void *b){
1268 const int16_t* da = (const int16_t *) a;
1269 const int16_t* db = (const int16_t *) b;
1270
1271 if(da[1] != db[1]) return da[1] - db[1];
1272 else return da[0] - db[0];
1273 }
1274
1275
1276 static void encode_subband(SnowContext *s, SubBand *b, DWTELEM *src, DWTELEM *parent, int stride, int orientation){
1277 const int level= b->level;
1278 const int w= b->width;
1279 const int h= b->height;
1280 int x, y;
1281
1282 #if 1
1283 if(orientation==3 && parent && 0){
1284 int16_t candidate[w*h][2];
1285 uint8_t state[w*h];
1286 int16_t boarder[3][w*h*4][2];
1287 int16_t significant[w*h][2];
1288 int candidate_count=0;
1289 int boarder_count[3]={0,0,0};
1290 int significant_count=0;
1291 int rle_pos=0;
1292 int v, last_v;
1293 int primary= orientation==1;
1294
1295 memset(candidate, 0, sizeof(candidate));
1296 memset(state, 0, sizeof(state));
1297 memset(boarder, 0, sizeof(boarder));
1298
1299 for(y=0; y<h; y++){
1300 for(x=0; x<w; x++){
1301 if(parent[(x>>1) + (y>>1)*2*stride])
1302 SVI(candidate, candidate_count, x, y)
1303 }
1304 }
1305
1306 for(;;){
1307 while(candidate_count && !boarder_count[0] && !boarder_count[1] && !boarder_count[2]){
1308 candidate_count--;
1309 x= candidate[ candidate_count][0];
1310 y= candidate[ candidate_count][1];
1311 if(state[x + y*w])
1312 continue;
1313 state[x + y*w]= 1;
1314 v= !!src[x + y*stride];
1315 put_cabac(&s->c, &b->state[0][0], v);
1316 if(v){
1317 SVI(significant, significant_count, x,y)
1318 if(x && !state[x - 1 + y *w]) SVI(boarder[0],boarder_count[0],x-1,y )
1319 if(y && !state[x + (y-1)*w]) SVI(boarder[1],boarder_count[1],x ,y-1)
1320 if(x+1<w && !state[x + 1 + y *w]) SVI(boarder[0],boarder_count[0],x+1,y )
1321 if(y+1<h && !state[x + (y+1)*w]) SVI(boarder[1],boarder_count[1],x ,y+1)
1322 if(x && y && !state[x - 1 + (y-1)*w]) SVI(boarder[2],boarder_count[2],x-1,y-1)
1323 if(x && y+1<h && !state[x - 1 + (y+1)*w]) SVI(boarder[2],boarder_count[2],x-1,y+1)
1324 if(x+1<w && y+1<h && !state[x + 1 + (y+1)*w]) SVI(boarder[2],boarder_count[2],x+1,y+1)
1325 if(x+1<w && y && !state[x + 1 + (y-1)*w]) SVI(boarder[2],boarder_count[2],x+1,y-1)
1326 }
1327 }
1328 while(!boarder_count[0] && !boarder_count[1] && !boarder_count[2] && rle_pos < w*h){
1329 int run=0;
1330 for(; rle_pos < w*h;){
1331 x= rle_pos % w; //FIXME speed
1332 y= rle_pos / w;
1333 rle_pos++;
1334 if(state[x + y*w])
1335 continue;
1336 state[x + y*w]= 1;
1337 v= !!src[x + y*stride];
1338 if(v){
1339 put_symbol(&s->c, b->state[1], run, 0);
1340 SVI(significant, significant_count, x,y)
1341 if(x && !state[x - 1 + y *w]) SVI(boarder[0],boarder_count[0],x-1,y )
1342 if(y && !state[x + (y-1)*w]) SVI(boarder[1],boarder_count[1],x ,y-1)
1343 if(x+1<w && !state[x + 1 + y *w]) SVI(boarder[0],boarder_count[0],x+1,y )
1344 if(y+1<h && !state[x + (y+1)*w]) SVI(boarder[1],boarder_count[1],x ,y+1)
1345 if(x && y && !state[x - 1 + (y-1)*w]) SVI(boarder[2],boarder_count[2],x-1,y-1)
1346 if(x && y+1<h && !state[x - 1 + (y+1)*w]) SVI(boarder[2],boarder_count[2],x-1,y+1)
1347 if(x+1<w && y+1<h && !state[x + 1 + (y+1)*w]) SVI(boarder[2],boarder_count[2],x+1,y+1)
1348 if(x+1<w && y && !state[x + 1 + (y-1)*w]) SVI(boarder[2],boarder_count[2],x+1,y-1)
1349 break;
1350 //FIXME note only right & down can be boarders
1351 }
1352 run++;
1353 }
1354 }
1355 if(!boarder_count[0] && !boarder_count[1] && !boarder_count[2])
1356 break;
1357
1358 while(boarder_count[0] || boarder_count[1] || boarder_count[2]){
1359 int index;
1360
1361 if (boarder_count[ primary]) index= primary;
1362 else if(boarder_count[1-primary]) index=1-primary;
1363 else index=2;
1364
1365 boarder_count[index]--;
1366 x= boarder[index][ boarder_count[index] ][0];
1367 y= boarder[index][ boarder_count[index] ][1];
1368 if(state[x + y*w]) //FIXME maybe check earlier
1369 continue;
1370 state[x + y*w]= 1;
1371 v= !!src[x + y*stride];
1372 put_cabac(&s->c, &b->state[0][index+1], v);
1373 if(v){
1374 SVI(significant, significant_count, x,y)
1375 if(x && !state[x - 1 + y *w]) SVI(boarder[0],boarder_count[0],x-1,y )
1376 if(y && !state[x + (y-1)*w]) SVI(boarder[1],boarder_count[1],x ,y-1)
1377 if(x+1<w && !state[x + 1 + y *w]) SVI(boarder[0],boarder_count[0],x+1,y )
1378 if(y+1<h && !state[x + (y+1)*w]) SVI(boarder[1],boarder_count[1],x ,y+1)
1379 if(x && y && !state[x - 1 + (y-1)*w]) SVI(boarder[2],boarder_count[2],x-1,y-1)
1380 if(x && y+1<h && !state[x - 1 + (y+1)*w]) SVI(boarder[2],boarder_count[2],x-1,y+1)
1381 if(x+1<w && y+1<h && !state[x + 1 + (y+1)*w]) SVI(boarder[2],boarder_count[2],x+1,y+1)
1382 if(x+1<w && y && !state[x + 1 + (y-1)*w]) SVI(boarder[2],boarder_count[2],x+1,y-1)
1383 }
1384 }
1385 }
1386 //FIXME sort significant coeffs maybe
1387 if(1){
1388 qsort(significant, significant_count, sizeof(int16_t[2]), sig_cmp);
1389 }
1390
1391 last_v=1;
1392 while(significant_count){
1393 int context= 3 + quant7[last_v&0xFF]; //use significance of suroundings
1394 significant_count--;
1395 x= significant[significant_count][0];//FIXME try opposit direction
1396 y= significant[significant_count][1];
1397 v= src[x + y*stride];
1398 put_symbol(&s->c, b->state[context + 2], v, 1); //FIXME try to avoid first bit, try this with the old code too!!
1399 last_v= v;
1400 }
1401 }
1402 #endif
1403 if(1){
1404 int w, h;
1405 int run=0;
1406 int last_run=0;
1407 int last_v=1;
1408 int last_x=0;
1409 int runs[b->width*b->height];
1410 int run_index=0;
1411
1412 if(orientation==1){
1413 w= b->height;
1414 h= b->width;
1415 }else{
1416 w= b->width;
1417 h= b->height;
1418 }
1419
1420 for(y=0; y<h; y++){
1421 for(x=0; x<w; x++){
1422 int v, p=0;
1423 int /*ll=0, */l=0, lt=0, t=0, rt=0;
1424 if(orientation==1) v= src[y + x*stride];
1425 else v= src[x + y*stride];
1426
1427 if(y){
1428 if(orientation==1) t= src[y - 1 + x*stride];
1429 else t= src[x + (y-1)*stride];
1430 if(x){
1431 if(orientation==1) lt= src[y - 1 + (x-1)*stride];
1432 else lt= src[x - 1 + (y-1)*stride];
1433 }
1434 if(x + 1 < w){
1435 if(orientation==1) rt= src[y - 1 + (x+1)*stride];
1436 else rt= src[x + 1 + (y-1)*stride];
1437 }
1438 }
1439 if(x){
1440 if(orientation==1) l= src[y + (x-1)*stride];
1441 else l= src[x - 1 + y*stride];
1442 /*if(x > 1){
1443 if(orientation==1) ll= src[y + (x-2)*stride];
1444 else ll= src[x - 2 + y*stride];
1445 }*/
1446 }
1447 if(parent){
1448 int px= (orientation==1 ? y : x)>>1;
1449 int py= (orientation==1 ? x : y)>>1;
1450 if(px<b->parent->width && py<b->parent->height)
1451 p= parent[px + py*2*stride];
1452 }
1453 if(!(/*ll|*/l|lt|t|rt|p)){
1454 if(v){
1455 runs[run_index++]= run;
1456 run=0;
1457 }else{
1458 run++;
1459 }
1460 }
1461 }
1462 }
1463 runs[run_index++]= run;
1464 run_index=0;
1465 run= runs[run_index++];
1466
1467 put_symbol(&s->c, b->state[1], run, 0);
1468
1469 for(y=0; y<h; y++){
1470 for(x=0; x<w; x++){
1471 int v, p=0;
1472 int /*ll=0, */l=0, lt=0, t=0, rt=0;
1473 if(orientation==1) v= src[y + x*stride];
1474 else v= src[x + y*stride];
1475
1476 if(y){
1477 if(orientation==1) t= src[y - 1 + x*stride];
1478 else t= src[x + (y-1)*stride];
1479 if(x){
1480 if(orientation==1) lt= src[y - 1 + (x-1)*stride];
1481 else lt= src[x - 1 + (y-1)*stride];
1482 }
1483 if(x + 1 < w){
1484 if(orientation==1) rt= src[y - 1 + (x+1)*stride];
1485 else rt= src[x + 1 + (y-1)*stride];
1486 }
1487 }
1488 if(x){
1489 if(orientation==1) l= src[y + (x-1)*stride];
1490 else l= src[x - 1 + y*stride];
1491 /*if(x > 1){
1492 if(orientation==1) ll= src[y + (x-2)*stride];
1493 else ll= src[x - 2 + y*stride];
1494 }*/
1495 }
1496 if(parent){
1497 int px= (orientation==1 ? y : x)>>1;
1498 int py= (orientation==1 ? x : y)>>1;
1499 if(px<b->parent->width && py<b->parent->height)
1500 p= parent[px + py*2*stride];
1501 }
1502 if(/*ll|*/l|lt|t|rt|p){
1503 int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));
1504
1505 put_cabac(&s->c, &b->state[0][context], !!v);
1506 }else{
1507 if(!run){
1508 run= runs[run_index++];
1509 put_symbol(&s->c, b->state[1], run, 0);
1510 assert(v);
1511 }else{
1512 run--;
1513 assert(!v);
1514 }
1515 }
1516 if(v){
1517 int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));
1518
1519 put_symbol(&s->c, b->state[context + 2], ABS(v)-1, 0);
1520 put_cabac(&s->c, &b->state[0][16 + 1 + 3 + quant3b[l&0xFF] + 3*quant3b[t&0xFF]], v<0);
1521 }
1522 }
1523 }
1524 return;
1525 }
1526 if(1){
1527 int w, h;
1528 int run=0;
1529 int last_run=0;
1530 int last_v=1;
1531 int last_x=0;
1532
1533 if(orientation==1){
1534 w= b->height;
1535 h= b->width;
1536 }else{
1537 w= b->width;
1538 h= b->height;
1539 }
1540
1541 for(y=0; y<h; y++){
1542 for(x=0; x<w; x++){
1543 int v, p;
1544 if(orientation==1) v= src[y + x*stride];
1545 else v= src[x + y*stride];
1546
1547 if(v){
1548 int sec=0, pri=0/*, lt=0, rt=0, pri2=0*/;
1549
1550 if(y){
1551 if(orientation==1) sec= src[y - 1 + x*stride];
1552 else sec= src[x + (y-1)*stride];
1553 /* if(x){
1554 if(orientation==1) lt= src[y - 1 + (x-1)*stride];
1555 else lt= src[x - 1 + (y-1)*stride];
1556 }*/
1557 }
1558 /* if (orientation==1 && y + 1 < h && x) rt= src[y + 1 + (x-1)*stride];
1559 else if (orientation!=1 && x + 1 < w && y) rt= src[x + 1 + (y-1)*stride];*/
1560 if(x){
1561 if(orientation==1) pri= src[y + (x-1)*stride];
1562 else pri= src[x - 1 + y*stride];
1563 /* if(x > 1){
1564 if(orientation==1) pri2= src[y + (x-2)*stride];
1565 else pri2= src[x - 2 + y*stride];
1566 }*/
1567 }
1568
1569 put_symbol(&s->c, b->state[ABS(quant7[last_v&0xFF])], run, 0);
1570 // context= 3*7 + 3 + quant11[last_v&0xFF] /*+ 7*quant3[(run-1)&0xFF]*/;
1571 // if(parent && orientation!=1) context += 3*7+7*quant3[parent[(x>>1) + (y>>1)*2*stride]&0xFF];
1572
1573 // put_symbol(&s->c, b->state[ABS(quant7[pri&0xFF]) + 4*ABS(quant7[sec&0xFF]) + 8], ABS(v)-1, 0);
1574 put_symbol(&s->c, b->state[quant13[(ABS(pri) + ABS(sec))&0x7F] + 8], ABS(v)-1, 0);
1575
1576 // context= quant3b[sec&0xFF] + 3*quant3b[pri&0xFF];
1577 // put_cabac(&s->c, &b->state[7][ABS(context)], (v^context)<0);
1578 put_cabac(&s->c, &b->state[7][1 + 3 + quant3b[sec&0xFF] + 3*quant3b[pri&0xFF]], v<0);
1579 last_run=run;
1580 last_v=v;
1581
1582 run=0;
1583 }else{
1584 run++;
1585 }
1586 }
1587 }
1588 if(run){
1589 put_symbol(&s->c, b->state[ABS(quant7[last_v&0xFF])], run, 0);
1590 }
1591 return;
1592 }
1593
1594 {
1595 int run=0;
1596 int mode=1;
1597 for(y=0; y<h; y++){
1598 for(x=0; x<w; x++){
1599 // int context=5*11*11 + 5*11+5;
1600 int context=3*7 + 3 + 1;
1601 int v= src[x + y*stride];
1602
1603 if(mode){
1604 if(v){
1605 put_symbol(&s->c, b->state[0], run, 0);
1606 mode=0;
1607 run=0;
1608 }else{
1609 run++;
1610 continue;
1611 }
1612 }
1613 if(x /*&& orientation!=1*/) context += quant7[src[x-1 + y*stride ]&0xFF];
1614 if(y /*&& orientation!=2*/) context += 7*quant7[src[x + y*stride-stride]&0xFF];
1615
1616 // if(x>1) context += 11*11*quant3 [src[x-2 + y*stride ]&0xFF];
1617 // if(parent) context += 11*11*quant3[parent[(x>>1) + (y>>1)*2*stride]&0xFF];
1618
1619 put_symbol(&s->c, b->state[context], v, 1);
1620 if(v==0){
1621 mode=1;
1622 }
1623 }
1624 }
1625 if(mode && run)
1626 put_symbol(&s->c, b->state[0], run, 0);
1627 }
1628 }
1629
1630 static inline void decode_subbandX(SnowContext *s, SubBand *b, DWTELEM *src, DWTELEM *parent, int stride, int orientation){
1631 const int level= b->level;
1632 const int w= b->width;
1633 const int h= b->height;
1634 int x,y;
1635
1636 START_TIMER
1637
1638 if(1){
1639 int w,h;
1640 int run;
1641 int last_run;
1642
1643 for(y=0; y<b->height; y++)
1644 memset(&src[y*stride], 0, b->width*sizeof(DWTELEM));
1645
1646 if(orientation==1){
1647 w= b->height;
1648 h= b->width;
1649 }else{
1650 w= b->width;
1651 h= b->height;
1652 }
1653 x=y=0;
1654 run= get_symbol(&s->c, b->state[1], 0);
1655 for(y=0; y<h; y++){
1656 for(x=0; x<w; x++){
1657 int v, p=0;
1658 int /*ll=0, */l=0, lt=0, t=0, rt=0;
1659
1660 if(y){
1661 if(orientation==1) t= src[y - 1 + x*stride];
1662 else t= src[x + (y-1)*stride];
1663 if(x){
1664 if(orientation==1) lt= src[y - 1 + (x-1)*stride];
1665 else lt= src[x - 1 + (y-1)*stride];
1666 }
1667 if(x + 1 < w){
1668 if(orientation==1) rt= src[y - 1 + (x+1)*stride];
1669 else rt= src[x + 1 + (y-1)*stride];
1670 }
1671 }
1672 if(x){
1673 if(orientation==1) l= src[y + (x-1)*stride];
1674 else l= src[x - 1 + y*stride];
1675 /*if(x > 1){
1676 if(orientation==1) ll= src[y + (x-2)*stride];
1677 else ll= src[x - 2 + y*stride];
1678 }*/
1679 }
1680 if(parent){
1681 int px= (orientation==1 ? y : x)>>1;
1682 int py= (orientation==1 ? x : y)>>1;
1683 if(px<b->parent->width && py<b->parent->height)
1684 p= parent[px + py*2*stride];
1685 }
1686 if(/*ll|*/l|lt|t|rt|p){
1687 int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));
1688
1689 v=get_cabac(&s->c, &b->state[0][context]);
1690 }else{
1691 if(!run){
1692 run= get_symbol(&s->c, b->state[1], 0);
1693 //FIXME optimize this here
1694 //FIXME try to store a more naive run
1695 v=1;
1696 }else{
1697 run--;
1698 v=0;
1699 }
1700 }
1701 if(v){
1702 int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));
1703 v= get_symbol(&s->c, b->state[context + 2], 0) + 1;
1704 if(get_cabac(&s->c, &b->state[0][16 + 1 + 3 + quant3b[l&0xFF] + 3*quant3b[t&0xFF]]))
1705 v= -v;
1706 if(orientation==1) src[y + x*stride]= v;
1707 else src[x + y*stride]= v;
1708 }
1709 }
1710 }
1711 if(level+1 == s->spatial_decomposition_count){
1712 STOP_TIMER("decode_subband")
1713 }
1714
1715 return;
1716 }
1717 if(1){
1718 int w,h;
1719 int last_run=0;
1720 int last_v=1;
1721
1722 for(y=0; y<b->height; y++)
1723 memset(&src[y*stride], 0, b->width*sizeof(DWTELEM));
1724
1725 if(orientation==1){
1726 w= b->height;
1727 h= b->width;
1728 }else{
1729 w= b->width;
1730 h= b->height;
1731 }
1732 x=y=0;
1733 for(;;){
1734 int sec=0, pri=0/*, lt=0, rt=0, pri2=0*/;
1735
1736 last_run= get_symbol(&s->c, b->state[ABS(quant7[last_v&0xFF])], 0);
1737
1738 // context= 3*7 + 3 + quant7[last_v&0xFF]/* + 7*quant7[(last_run-1)&0xFF]*/;
1739 x += last_run;
1740 while(x>=w){
1741 x-= w;
1742 if(++y>=h) break;
1743 /* if(orientation!=1)
1744 memset(&src[y*stride], 0, w*sizeof(DWTELEM));*/
1745 }
1746 if(y>=h) break;
1747 if(y){
1748 if(orientation==1) sec= src[y - 1 + x*stride];
1749 else sec= src[x + (y-1)*stride];
1750 /* if(x){
1751 if(orientation==1) lt= src[y - 1 + (x-1)*stride];
1752 else lt= src[x - 1 + (y-1)*stride];
1753 }*/
1754 }
1755 /* if (orientation==1 && y + 1 < h && x) rt= src[y + 1 + (x-1)*stride];
1756 else if (orientation!=1 && x + 1 < w && y) rt= src[x + 1 + (y-1)*stride];*/
1757 if(x){
1758 if(orientation==1) pri= src[y + (x-1)*stride];
1759 else pri= src[x - 1 + y*stride];
1760 /* if(x > 1){
1761 if(orientation==1) pri2= src[y + (x-2)*stride];
1762 else pri2= src[x - 2 + y*stride];
1763 }*/
1764 }
1765
1766 last_v= get_symbol(&s->c, b->state[quant13[(ABS(pri) + ABS(sec))&0x7F] + 8], 0) + 1;
1767 if(get_cabac(&s->c, &b->state[7][1 + 3 + quant3b[sec&0xFF] + 3*quant3b[pri&0xFF]]))
1768 last_v= -last_v;
1769
1770 if(orientation==1) src[y + x*stride]= last_v;
1771 else src[x + y*stride]= last_v;
1772
1773 if(++x==w){
1774 x=0;
1775 if(++y>=h) break;
1776 /* if(orientation!=1)
1777 memset(&src[y*stride], 0, w*sizeof(DWTELEM));*/
1778 }
1779 }
1780
1781 if(level+1 == s->spatial_decomposition_count){
1782 STOP_TIMER("decode_subband")
1783 }
1784
1785 return;
1786 }
1787 }
1788
1789 static void decode_subband(SnowContext *s, SubBand *b, DWTELEM *src, DWTELEM *parent, int stride, int orientation){
1790 if(orientation==1) decode_subbandX(s, b, src, parent, stride, 1);
1791 else decode_subbandX(s, b, src, parent, stride, 0);
1792 }
1793
1794 static void reset_contexts(SnowContext *s){
1795 int plane_index, level, orientation;
1796
1797 for(plane_index=0; plane_index<2; plane_index++){
1798 for(level=0; level<s->spatial_decomposition_count; level++){
1799 for(orientation=level ? 1:0; orientation<4; orientation++){
1800 memset(s->plane[plane_index].band[level][orientation].state, 0, sizeof(s->plane[plane_index].band[level][orientation].state));
1801 }
1802 }
1803 }
1804 memset(s->mb_band.state, 0, sizeof(s->mb_band.state));
1805 memset(s->mv_band[0].state, 0, sizeof(s->mv_band[0].state));
1806 memset(s->mv_band[1].state, 0, sizeof(s->mv_band[1].state));
1807 memset(s->header_state, 0, sizeof(s->header_state));
1808 }
1809
1810 static void mc_block(uint8_t *dst, uint8_t *src, uint8_t *tmp, int stride, int b_w, int b_h, int dx, int dy){
1811 int x, y;
1812
1813 for(y=0; y < b_h+5; y++){
1814 for(x=0; x < b_w; x++){
1815 int a0= src[x + y*stride];
1816 int a1= src[x + 1 + y*stride];
1817 int a2= src[x + 2 + y*stride];
1818 int a3= src[x + 3 + y*stride];
1819 int a4= src[x + 4 + y*stride];
1820 int a5= src[x + 5 + y*stride];
1821 // int am= 9*(a1+a2) - (a0+a3);
1822 int am= 20*(a2+a3) - 5*(a1+a4) + (a0+a5);
1823 // int am= 18*(a2+a3) - 2*(a1+a4);
1824 // int aL= (-7*a0 + 105*a1 + 35*a2 - 5*a3)>>3;
1825 // int aR= (-7*a3 + 105*a2 + 35*a1 - 5*a0)>>3;
1826
1827 // if(b_w==16) am= 8*(a1+a2);
1828
1829 if(dx<8) tmp[x + y*stride]= (32*a2*( 8-dx) + am* dx + 128)>>8;
1830 else tmp[x + y*stride]= ( am*(16-dx) + 32*a3*(dx-8) + 128)>>8;
1831
1832 /* if (dx< 4) tmp[x + y*stride]= (16*a1*( 4-dx) + aL* dx + 32)>>6;
1833 else if(dx< 8) tmp[x + y*stride]= ( aL*( 8-dx) + am*(dx- 4) + 32)>>6;
1834 else if(dx<12) tmp[x + y*stride]= ( am*(12-dx) + aR*(dx- 8) + 32)>>6;
1835 else tmp[x + y*stride]= ( aR*(16-dx) + 16*a2*(dx-12) + 32)>>6;*/
1836 }
1837 }
1838 for(y=0; y < b_h; y++){
1839 for(x=0; x < b_w; x++){
1840 int a0= tmp[x + y *stride];
1841 int a1= tmp[x + (y + 1)*stride];
1842 int a2= tmp[x + (y + 2)*stride];
1843 int a3= tmp[x + (y + 3)*stride];
1844 int a4= tmp[x + (y + 4)*stride];
1845 int a5= tmp[x + (y + 5)*stride];
1846 int am= 20*(a2+a3) - 5*(a1+a4) + (a0+a5);
1847 // int am= 18*(a2+a3) - 2*(a1+a4);
1848 /* int aL= (-7*a0 + 105*a1 + 35*a2 - 5*a3)>>3;
1849 int aR= (-7*a3 + 105*a2 + 35*a1 - 5*a0)>>3;*/
1850
1851 // if(b_w==16) am= 8*(a1+a2);
1852
1853 if(dy<8) dst[x + y*stride]= (32*a2*( 8-dy) + am* dy + 128)>>8;
1854 else dst[x + y*stride]= ( am*(16-dy) + 32*a3*(dy-8) + 128)>>8;
1855
1856 /* if (dy< 4) tmp[x + y*stride]= (16*a1*( 4-dy) + aL* dy + 32)>>6;
1857 else if(dy< 8) tmp[x + y*stride]= ( aL*( 8-dy) + am*(dy- 4) + 32)>>6;
1858 else if(dy<12) tmp[x + y*stride]= ( am*(12-dy) + aR*(dy- 8) + 32)>>6;
1859 else tmp[x + y*stride]= ( aR*(16-dy) + 16*a2*(dy-12) + 32)>>6;*/
1860 }
1861 }
1862 }
1863
1864 #define mcb(dx,dy,b_w)\
1865 static void mc_block ## dx ## dy(uint8_t *dst, uint8_t *src, int stride){\
1866 uint8_t tmp[stride*(b_w+5)];\
1867 mc_block(dst, src-2-2*stride, tmp, stride, b_w, b_w, dx, dy);\
1868 }
1869
1870 mcb( 0, 0,16)
1871 mcb( 4, 0,16)
1872 mcb( 8, 0,16)
1873 mcb(12, 0,16)
1874 mcb( 0, 4,16)
1875 mcb( 4, 4,16)
1876 mcb( 8, 4,16)
1877 mcb(12, 4,16)
1878 mcb( 0, 8,16)
1879 mcb( 4, 8,16)
1880 mcb( 8, 8,16)
1881 mcb(12, 8,16)
1882 mcb( 0,12,16)
1883 mcb( 4,12,16)
1884 mcb( 8,12,16)
1885 mcb(12,12,16)
1886
1887 #define mca(dx,dy,b_w)\
1888 static void mc_block_hpel ## dx ## dy(uint8_t *dst, uint8_t *src, int stride, int h){\
1889 uint8_t tmp[stride*(b_w+5)];\
1890 assert(h==b_w);\
1891 mc_block(dst, src-2-2*stride, tmp, stride, b_w, b_w, dx, dy);\
1892 }
1893
1894 mca( 0, 0,16)
1895 mca( 8, 0,16)
1896 mca( 0, 8,16)
1897 mca( 8, 8,16)
1898
1899 static void add_xblock(DWTELEM *dst, uint8_t *src, uint8_t *obmc, int s_x, int s_y, int b_w, int b_h, int mv_x, int mv_y, int w, int h, int dst_stride, int src_stride, int obmc_stride, int mb_type, int add){
1900 uint8_t tmp[src_stride*(b_h+5)]; //FIXME move to context to gurantee alignment
1901 int x,y;
1902
1903 if(s_x<0){
1904 obmc -= s_x;
1905 b_w += s_x;
1906 s_x=0;
1907 }else if(s_x + b_w > w){
1908 b_w = w - s_x;
1909 }
1910 if(s_y<0){
1911 obmc -= s_y*obmc_stride;
1912 b_h += s_y;
1913 s_y=0;
1914 }else if(s_y + b_h> h){
1915 b_h = h - s_y;
1916 }
1917
1918 dst += s_x + s_y*dst_stride;
1919
1920 if(mb_type==1){
1921 src += s_x + s_y*src_stride;
1922 for(y=0; y < b_h; y++){
1923 for(x=0; x < b_w; x++){
1924 if(add) dst[x + y*dst_stride] += obmc[x + y*obmc_stride] * 128 * (256/OBMC_MAX);
1925 else dst[x + y*dst_stride] -= obmc[x + y*obmc_stride] * 128 * (256/OBMC_MAX);
1926 }
1927 }
1928 }else{
1929 int dx= mv_x&15;
1930 int dy= mv_y&15;
1931 // int dxy= (mv_x&1) + 2*(mv_y&1);
1932
1933 s_x += (mv_x>>4) - 2;
1934 s_y += (mv_y>>4) - 2;
1935 src += s_x + s_y*src_stride;
1936 //use dsputil
1937
1938 if( (unsigned)s_x >= w - b_w - 4
1939 || (unsigned)s_y >= h - b_h - 4){
1940 ff_emulated_edge_mc(tmp + 32, src, src_stride, b_w+5, b_h+5, s_x, s_y, w, h);
1941 src= tmp + 32;
1942 }
1943
1944 if(mb_type==0){
1945 mc_block(tmp, src, tmp + 64+8, src_stride, b_w, b_h, dx, dy);
1946 }else{
1947 int sum=0;
1948 for(y=0; y < b_h; y++){
1949 for(x=0; x < b_w; x++){
1950 sum += src[x+ y*src_stride];
1951 }
1952 }
1953 sum= (sum + b_h*b_w/2) / (b_h*b_w);
1954 for(y=0; y < b_h; y++){
1955 for(x=0; x < b_w; x++){
1956 tmp[x + y*src_stride]= sum;
1957 }
1958 }
1959 }
1960
1961 for(y=0; y < b_h; y++){
1962 for(x=0; x < b_w; x++){
1963 if(add) dst[x + y*dst_stride] += obmc[x + y*obmc_stride] * tmp[x + y*src_stride] * (256/OBMC_MAX);
1964 else dst[x + y*dst_stride] -= obmc[x + y*obmc_stride] * tmp[x + y*src_stride] * (256/OBMC_MAX);
1965 }
1966 }
1967 }
1968 }
1969
1970 static void predict_plane(SnowContext *s, DWTELEM *buf, int plane_index, int add){
1971 Plane *p= &s->plane[plane_index];
1972 const int mb_w= s->mb_band.width;
1973 const int mb_h= s->mb_band.height;
1974 const int mb_stride= s->mb_band.stride;
1975 int x, y, mb_x, mb_y;
1976 int scale = plane_index ? s->mv_scale : 2*s->mv_scale;
1977 int block_w = plane_index ? 8 : 16;
1978 uint8_t *obmc = plane_index ? obmc16 : obmc32;
1979 int obmc_stride= plane_index ? 16 : 32;
1980 int ref_stride= s->last_picture.linesize[plane_index];
1981 uint8_t *ref = s->last_picture.data[plane_index];
1982 int w= p->width;
1983 int h= p->height;
1984
1985 if(s->avctx->debug&512){
1986 for(y=0; y<h; y++){
1987 for(x=0; x<w; x++){
1988 if(add) buf[x + y*w]+= 128*256;
1989 else buf[x + y*w]-= 128*256;
1990 }
1991 }
1992
1993 return;
1994 }
1995 for(mb_y=-1; mb_y<=mb_h; mb_y++){
1996 for(mb_x=-1; mb_x<=mb_w; mb_x++){
1997 int index= clip(mb_x, 0, mb_w-1) + clip(mb_y, 0, mb_h-1)*mb_stride;
1998
1999 add_xblock(buf, ref, obmc,
2000 block_w*mb_x - block_w/2,
2001 block_w*mb_y - block_w/2,
2002 2*block_w, 2*block_w,
2003 s->mv_band[0].buf[index]*scale, s->mv_band[1].buf[index]*scale,
2004 w, h,
2005 w, ref_stride, obmc_stride,
2006 s->mb_band.buf[index], add);
2007
2008 }
2009 }
2010 }
2011
2012 static void quantize(SnowContext *s, SubBand *b, DWTELEM *src, int stride, int bias){
2013 const int level= b->level;
2014 const int w= b->width;
2015 const int h= b->height;
2016 const int qlog= clip(s->qlog + b->qlog, 0, 128);
2017 const int qmul= qexp[qlog&7]<<(qlog>>3);
2018 int x,y;
2019
2020 assert(QROOT==8);
2021
2022 bias= bias ? 0 : (3*qmul)>>3;
2023
2024 if(!bias){
2025 for(y=0; y<h; y++){
2026 for(x=0; x<w; x++){
2027 int i= src[x + y*stride];
2028 //FIXME use threshold
2029 //FIXME optimize
2030 //FIXME bias
2031 if(i>=0){
2032 i<<= QEXPSHIFT;
2033 i/= qmul;
2034 src[x + y*stride]= i;
2035 }else{
2036 i= -i;
2037 i<<= QEXPSHIFT;
2038 i/= qmul;
2039 src[x + y*stride]= -i;
2040 }
2041 }
2042 }
2043 }else{
2044 for(y=0; y<h; y++){
2045 for(x=0; x<w; x++){
2046 int i= src[x + y*stride];
2047
2048 //FIXME use threshold
2049 //FIXME optimize
2050 //FIXME bias
2051 if(i>=0){
2052 i<<= QEXPSHIFT;
2053 i= (i + bias) / qmul;
2054 src[x + y*stride]= i;
2055 }else{
2056 i= -i;
2057 i<<= QEXPSHIFT;
2058 i= (i + bias) / qmul;
2059 src[x + y*stride]= -i;
2060 }
2061 }
2062 }
2063 }
2064 }
2065
2066 static void dequantize(SnowContext *s, SubBand *b, DWTELEM *src, int stride){
2067 const int level= b->level;
2068 const int w= b->width;
2069 const int h= b->height;
2070 const int qlog= clip(s->qlog + b->qlog, 0, 128);
2071 const int qmul= qexp[qlog&7]<<(qlog>>3);
2072 const int qadd= (s->qbias*qmul)>>QBIAS_SHIFT;
2073 int x,y;
2074
2075 assert(QROOT==8);
2076
2077 for(y=0; y<h; y++){
2078 for(x=0; x<w; x++){
2079 int i= src[x + y*stride];
2080 if(i<0){
2081 src[x + y*stride]= -((-i*qmul + qadd)>>(QEXPSHIFT)); //FIXME try different bias
2082 }else if(i>0){
2083 src[x + y*stride]= (( i*qmul + qadd)>>(QEXPSHIFT));
2084 }
2085 }
2086 }
2087 }
2088
2089 static void decorrelate(SnowContext *s, SubBand *b, DWTELEM *src, int stride, int inverse, int use_median){
2090 const int w= b->width;
2091 const int h= b->height;
2092 int x,y;
2093
2094 for(y=h-1; y>=0; y--){
2095 for(x=w-1; x>=0; x--){
2096 int i= x + y*stride;
2097
2098 if(x){
2099 if(use_median){
2100 if(y && x+1<w) src[i] -= mid_pred(src[i - 1], src[i - stride], src[i - stride + 1]);
2101 else src[i] -= src[i - 1];
2102 }else{
2103 if(y) src[i] -= mid_pred(src[i - 1], src[i - stride], src[i - 1] + src[i - stride] - src[i - 1 - stride]);
2104 else src[i] -= src[i - 1];
2105 }
2106 }else{
2107 if(y) src[i] -= src[i - stride];
2108 }
2109 }
2110 }
2111 }
2112
2113 static void correlate(SnowContext *s, SubBand *b, DWTELEM *src, int stride, int inverse, int use_median){
2114 const int w= b->width;
2115 const int h= b->height;
2116 int x,y;
2117
2118 for(y=0; y<h; y++){
2119 for(x=0; x<w; x++){
2120 int i= x + y*stride;
2121
2122 if(x){
2123 if(use_median){
2124 if(y && x+1<w) src[i] += mid_pred(src[i - 1], src[i - stride], src[i - stride + 1]);
2125 else src[i] += src[i - 1];
2126 }else{
2127 if(y) src[i] += mid_pred(src[i - 1], src[i - stride], src[i - 1] + src[i - stride] - src[i - 1 - stride]);
2128 else src[i] += src[i - 1];
2129 }
2130 }else{
2131 if(y) src[i] += src[i - stride];
2132 }
2133 }
2134 }
2135 }
2136
2137 static void encode_header(SnowContext *s){
2138 int plane_index, level, orientation;
2139
2140 put_cabac(&s->c, s->header_state, s->keyframe); // state clearing stuff?
2141 if(s->keyframe){
2142 put_symbol(&s->c, s->header_state, s->version, 0);
2143 put_symbol(&s->c, s->header_state, s->temporal_decomposition_type, 0);
2144 put_symbol(&s->c, s->header_state, s->temporal_decomposition_count, 0);
2145 put_symbol(&s->c, s->header_state, s->spatial_decomposition_count, 0);
2146 put_symbol(&s->c, s->header_state, s->colorspace_type, 0);
2147 put_symbol(&s->c, s->header_state, s->b_width, 0);
2148 put_symbol(&s->c, s->header_state, s->b_height, 0);
2149 put_symbol(&s->c, s->header_state, s->chroma_h_shift, 0);
2150 put_symbol(&s->c, s->header_state, s->chroma_v_shift, 0);
2151 put_cabac(&s->c, s->header_state, s->spatial_scalability);
2152 // put_cabac(&s->c, s->header_state, s->rate_scalability);
2153
2154 for(plane_index=0; plane_index<2; plane_index++){
2155 for(level=0; level<s->spatial_decomposition_count; level++){
2156 for(orientation=level ? 1:0; orientation<4; orientation++){
2157 if(orientation==2) continue;
2158 put_symbol(&s->c, s->header_state, s->plane[plane_index].band[level][orientation].qlog, 1);
2159 }
2160 }
2161 }
2162 }
2163 put_symbol(&s->c, s->header_state, s->spatial_decomposition_type, 0);
2164 put_symbol(&s->c, s->header_state, s->qlog, 1);
2165 put_symbol(&s->c, s->header_state, s->mv_scale, 0);
2166 put_symbol(&s->c, s->header_state, s->qbias, 1);
2167 }
2168
2169 static int decode_header(SnowContext *s){
2170 int plane_index, level, orientation;
2171
2172 s->keyframe= get_cabac(&s->c, s->header_state);
2173 if(s->keyframe){
2174 s->version= get_symbol(&s->c, s->header_state, 0);
2175 if(s->version>0){
2176 av_log(s->avctx, AV_LOG_ERROR, "version %d not supported", s->version);
2177 return -1;
2178 }
2179 s->temporal_decomposition_type= get_symbol(&s->c, s->header_state, 0);
2180 s->temporal_decomposition_count= get_symbol(&s->c, s->header_state, 0);
2181 s->spatial_decomposition_count= get_symbol(&s->c, s->header_state, 0);
2182 s->colorspace_type= get_symbol(&s->c, s->header_state, 0);
2183 s->b_width= get_symbol(&s->c, s->header_state, 0);
2184 s->b_height= get_symbol(&s->c, s->header_state, 0);
2185 s->chroma_h_shift= get_symbol(&s->c, s->header_state, 0);
2186 s->chroma_v_shift= get_symbol(&s->c, s->header_state, 0);
2187 s->spatial_scalability= get_cabac(&s->c, s->header_state);
2188 // s->rate_scalability= get_cabac(&s->c, s->header_state);
2189
2190 for(plane_index=0; plane_index<3; plane_index++){
2191 for(level=0; level<s->spatial_decomposition_count; level++){
2192 for(orientation=level ? 1:0; orientation<4; orientation++){
2193 int q;
2194 if (plane_index==2) q= s->plane[1].band[level][orientation].qlog;
2195 else if(orientation==2) q= s->plane[plane_index].band[level][1].qlog;
2196 else q= get_symbol(&s->c, s->header_state, 1);
2197 s->plane[plane_index].band[level][orientation].qlog= q;
2198 }
2199 }
2200 }
2201 }
2202
2203 s->spatial_decomposition_type= get_symbol(&s->c, s->header_state, 0);
2204 if(s->spatial_decomposition_type > 2){
2205 av_log(s->avctx, AV_LOG_ERROR, "spatial_decomposition_type %d not supported", s->spatial_decomposition_type);
2206 return -1;
2207 }
2208
2209 s->qlog= get_symbol(&s->c, s->header_state, 1);
2210 s->mv_scale= get_symbol(&s->c, s->header_state, 0);
2211 s->qbias= get_symbol(&s->c, s->header_state, 1);
2212
2213 return 0;
2214 }
2215
2216 static int common_init(AVCodecContext *avctx){
2217 SnowContext *s = avctx->priv_data;
2218 int width, height;
2219 int level, orientation, plane_index, dec;
2220
2221 s->avctx= avctx;
2222
2223 dsputil_init(&s->dsp, avctx);
2224
2225 #define mcf(dx,dy)\
2226 s->dsp.put_qpel_pixels_tab [0][dy+dx/4]=\
2227 s->dsp.put_no_rnd_qpel_pixels_tab[0][dy+dx/4]=\
2228 mc_block ## dx ## dy;
2229
2230 mcf( 0, 0)
2231 mcf( 4, 0)
2232 mcf( 8, 0)
2233 mcf(12, 0)
2234 mcf( 0, 4)
2235 mcf( 4, 4)
2236 mcf( 8, 4)
2237 mcf(12, 4)
2238 mcf( 0, 8)
2239 mcf( 4, 8)
2240 mcf( 8, 8)
2241 mcf(12, 8)
2242 mcf( 0,12)
2243 mcf( 4,12)
2244 mcf( 8,12)
2245 mcf(12,12)
2246
2247 #define mcfh(dx,dy)\
2248 s->dsp.put_pixels_tab [0][dy/4+dx/8]=\
2249 s->dsp.put_no_rnd_pixels_tab[0][dy/4+dx/8]=\
2250 mc_block_hpel ## dx ## dy;
2251
2252 mcfh(0, 0)
2253 mcfh(8, 0)
2254 mcfh(0, 8)
2255 mcfh(8, 8)
2256
2257 dec= s->spatial_decomposition_count= 5;
2258 s->spatial_decomposition_type= avctx->prediction_method; //FIXME add decorrelator type r transform_type
2259
2260 s->chroma_h_shift= 1; //FIXME XXX
2261 s->chroma_v_shift= 1;
2262
2263 // dec += FFMAX(s->chroma_h_shift, s->chroma_v_shift);
2264
2265 s->b_width = (s->avctx->width +(1<<dec)-1)>>dec;
2266 s->b_height= (s->avctx->height+(1<<dec)-1)>>dec;
2267
2268 s->spatial_dwt_buffer= av_mallocz(s->b_width*s->b_height*sizeof(DWTELEM)<<(2*dec));
2269 s->pred_buffer= av_mallocz(s->b_width*s->b_height*sizeof(DWTELEM)<<(2*dec));
2270
2271 s->mv_scale= (s->avctx->flags & CODEC_FLAG_QPEL) ? 2 : 4;
2272
2273 for(plane_index=0; plane_index<3; plane_index++){
2274 int w= s->avctx->width;
2275 int h= s->avctx->height;
2276
2277 if(plane_index){
2278 w>>= s->chroma_h_shift;
2279 h>>= s->chroma_v_shift;
2280 }
2281 s->plane[plane_index].width = w;
2282 s->plane[plane_index].height= h;
2283 av_log(NULL, AV_LOG_DEBUG, "%d %d\n", w, h);
2284 for(level=s->spatial_decomposition_count-1; level>=0; level--){
2285 for(orientation=level ? 1 : 0; orientation<4; orientation++){
2286 SubBand *b= &s->plane[plane_index].band[level][orientation];
2287
2288 b->buf= s->spatial_dwt_buffer;
2289 b->level= level;
2290 b->stride= s->plane[plane_index].width << (s->spatial_decomposition_count - level);
2291 b->width = (w + !(orientation&1))>>1;
2292 b->height= (h + !(orientation>1))>>1;
2293
2294 if(orientation&1) b->buf += (w+1)>>1;
2295 if(orientation>1) b->buf += b->stride>>1;
2296
2297 if(level)
2298 b->parent= &s->plane[plane_index].band[level-1][orientation];
2299 }
2300 w= (w+1)>>1;
2301 h= (h+1)>>1;
2302 }
2303 }
2304
2305 //FIXME init_subband() ?
2306 s->mb_band.stride= s->mv_band[0].stride= s->mv_band[1].stride=
2307 s->mb_band.width = s->mv_band[0].width = s->mv_band[1].width = (s->avctx->width + 15)>>4;
2308 s->mb_band.height= s->mv_band[0].height= s->mv_band[1].height= (s->avctx->height+ 15)>>4;
2309 s->mb_band .buf= av_mallocz(s->mb_band .stride * s->mb_band .height*sizeof(DWTELEM));
2310 s->mv_band[0].buf= av_mallocz(s->mv_band[0].stride * s->mv_band[0].height*sizeof(DWTELEM));
2311 s->mv_band[1].buf= av_mallocz(s->mv_band[1].stride * s->mv_band[1].height*sizeof(DWTELEM));
2312
2313 reset_contexts(s);
2314 /*
2315 width= s->width= avctx->width;
2316 height= s->height= avctx->height;
2317
2318 assert(width && height);
2319 */
2320 s->avctx->get_buffer(s->avctx, &s->mconly_picture);
2321
2322 return 0;
2323 }
2324
2325
2326 static void calculate_vissual_weight(SnowContext *s, Plane *p){
2327 int width = p->width;
2328 int height= p->height;
2329 int i, level, orientation, x, y;
2330
2331 for(level=0; level<s->spatial_decomposition_count; level++){
2332 for(orientation=level ? 1 : 0; orientation<4; orientation++){
2333 SubBand *b= &p->band[level][orientation];
2334 DWTELEM *buf= b->buf;
2335 int64_t error=0;
2336
2337 memset(s->spatial_dwt_buffer, 0, sizeof(int)*width*height);
2338 buf[b->width/2 + b->height/2*b->stride]= 256*256;
2339 spatial_idwt(s, s->spatial_dwt_buffer, width, height, width);
2340 for(y=