Merge vertical_compose53iL0*()
[ffmpeg.git] / libavcodec / snow.c
1 /*
2  * Copyright (C) 2004 Michael Niedermayer <michaelni@gmx.at>
3  *
4  * This file is part of FFmpeg.
5  *
6  * FFmpeg is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU Lesser General Public
8  * License as published by the Free Software Foundation; either
9  * version 2.1 of the License, or (at your option) any later version.
10  *
11  * FFmpeg is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with FFmpeg; if not, write to the Free Software
18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19  */
20
21 #include "avcodec.h"
22 #include "dsputil.h"
23 #include "snow.h"
24
25 #include "rangecoder.h"
26 #include "mathops.h"
27
28 #include "mpegvideo.h"
29
30 #undef NDEBUG
31 #include <assert.h>
32
33 static const int8_t quant3[256]={
34  0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
35  1, 1, 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, 0,
50 };
51 static const int8_t quant3b[256]={
52  0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
53  1, 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 };
69 static const int8_t quant3bA[256]={
70  0, 0, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
71  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
72  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
73  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
74  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
75  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
76  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
77  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
78  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
79  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
80  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
81  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
82  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
83  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
84  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
85  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
86 };
87 static const int8_t quant5[256]={
88  0, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
89  2, 2, 2, 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, 2, 2, 2, 2, 2, 2, 2, 2,
92  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
93  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
94  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
95  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
96 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
97 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
98 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
99 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
100 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
101 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
102 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
103 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-1,-1,-1,
104 };
105 static const int8_t quant7[256]={
106  0, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
107  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
108  2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
109  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
110  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
111  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
112  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
113  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
114 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
115 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
116 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
117 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
118 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
119 -3,-3,-3,-3,-3,-3,-3,-3,-3,-2,-2,-2,-2,-2,-2,-2,
120 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
121 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-1,-1,
122 };
123 static const int8_t quant9[256]={
124  0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3,
125  3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 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, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
128  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
129  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
130  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
131  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
132 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
133 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
134 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
135 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
136 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
137 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
138 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-3,-3,-3,-3,
139 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-2,-2,-2,-2,-1,-1,
140 };
141 static const int8_t quant11[256]={
142  0, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4,
143  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
144  4, 4, 4, 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, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
147  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
148  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
149  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
150 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
151 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
152 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
153 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
154 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
155 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-4,-4,
156 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
157 -4,-4,-4,-4,-4,-3,-3,-3,-3,-3,-3,-3,-2,-2,-2,-1,
158 };
159 static const int8_t quant13[256]={
160  0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4,
161  4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
162  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
163  5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
164  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
165  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
166  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
167  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
168 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,
169 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,
170 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,
171 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,
172 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-5,
173 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
174 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
175 -4,-4,-4,-4,-4,-4,-4,-4,-4,-3,-3,-3,-3,-2,-2,-1,
176 };
177
178 #if 0 //64*cubic
179 static const uint8_t obmc32[1024]={
180   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,
181   0,  0,  0,  0,  0,  4,  4,  4,  4,  4,  4,  4,  4,  8,  8,  8,  8,  8,  8,  4,  4,  4,  4,  4,  4,  4,  4,  0,  0,  0,  0,  0,
182   0,  0,  0,  4,  4,  4,  4,  8,  8, 12, 12, 12, 16, 16, 16, 16, 16, 16, 16, 16, 12, 12, 12,  8,  8,  4,  4,  4,  4,  0,  0,  0,
183   0,  0,  4,  4,  8,  8, 12, 16, 16, 20, 24, 24, 28, 28, 32, 32, 32, 32, 28, 28, 24, 24, 20, 16, 16, 12,  8,  8,  4,  4,  0,  0,
184   0,  0,  4,  8,  8, 12, 16, 24, 28, 32, 36, 40, 44, 48, 48, 48, 48, 48, 48, 44, 40, 36, 32, 28, 24, 16, 12,  8,  8,  4,  0,  0,
185   0,  4,  4,  8, 12, 20, 24, 32, 40, 44, 52, 56, 60, 64, 68, 72, 72, 68, 64, 60, 56, 52, 44, 40, 32, 24, 20, 12,  8,  4,  4,  0,
186   0,  4,  4, 12, 16, 24, 32, 40, 52, 60, 68, 76, 80, 88, 88, 92, 92, 88, 88, 80, 76, 68, 60, 52, 40, 32, 24, 16, 12,  4,  4,  0,
187   0,  4,  8, 16, 24, 32, 40, 52, 64, 76, 84, 92,100,108,112,116,116,112,108,100, 92, 84, 76, 64, 52, 40, 32, 24, 16,  8,  4,  0,
188   0,  4,  8, 16, 28, 40, 52, 64, 76, 88,100,112,124,132,136,140,140,136,132,124,112,100, 88, 76, 64, 52, 40, 28, 16,  8,  4,  0,
189   0,  4, 12, 20, 32, 44, 60, 76, 88,104,120,132,144,152,160,164,164,160,152,144,132,120,104, 88, 76, 60, 44, 32, 20, 12,  4,  0,
190   0,  4, 12, 24, 36, 48, 68, 84,100,120,136,152,164,176,180,184,184,180,176,164,152,136,120,100, 84, 68, 48, 36, 24, 12,  4,  0,
191   0,  4, 12, 24, 40, 56, 76, 92,112,132,152,168,180,192,204,208,208,204,192,180,168,152,132,112, 92, 76, 56, 40, 24, 12,  4,  0,
192   0,  4, 16, 28, 44, 60, 80,100,124,144,164,180,196,208,220,224,224,220,208,196,180,164,144,124,100, 80, 60, 44, 28, 16,  4,  0,
193   0,  8, 16, 28, 48, 64, 88,108,132,152,176,192,208,224,232,240,240,232,224,208,192,176,152,132,108, 88, 64, 48, 28, 16,  8,  0,
194   0,  4, 16, 32, 48, 68, 88,112,136,160,180,204,220,232,244,248,248,244,232,220,204,180,160,136,112, 88, 68, 48, 32, 16,  4,  0,
195   1,  8, 16, 32, 48, 72, 92,116,140,164,184,208,224,240,248,255,255,248,240,224,208,184,164,140,116, 92, 72, 48, 32, 16,  8,  1,
196   1,  8, 16, 32, 48, 72, 92,116,140,164,184,208,224,240,248,255,255,248,240,224,208,184,164,140,116, 92, 72, 48, 32, 16,  8,  1,
197   0,  4, 16, 32, 48, 68, 88,112,136,160,180,204,220,232,244,248,248,244,232,220,204,180,160,136,112, 88, 68, 48, 32, 16,  4,  0,
198   0,  8, 16, 28, 48, 64, 88,108,132,152,176,192,208,224,232,240,240,232,224,208,192,176,152,132,108, 88, 64, 48, 28, 16,  8,  0,
199   0,  4, 16, 28, 44, 60, 80,100,124,144,164,180,196,208,220,224,224,220,208,196,180,164,144,124,100, 80, 60, 44, 28, 16,  4,  0,
200   0,  4, 12, 24, 40, 56, 76, 92,112,132,152,168,180,192,204,208,208,204,192,180,168,152,132,112, 92, 76, 56, 40, 24, 12,  4,  0,
201   0,  4, 12, 24, 36, 48, 68, 84,100,120,136,152,164,176,180,184,184,180,176,164,152,136,120,100, 84, 68, 48, 36, 24, 12,  4,  0,
202   0,  4, 12, 20, 32, 44, 60, 76, 88,104,120,132,144,152,160,164,164,160,152,144,132,120,104, 88, 76, 60, 44, 32, 20, 12,  4,  0,
203   0,  4,  8, 16, 28, 40, 52, 64, 76, 88,100,112,124,132,136,140,140,136,132,124,112,100, 88, 76, 64, 52, 40, 28, 16,  8,  4,  0,
204   0,  4,  8, 16, 24, 32, 40, 52, 64, 76, 84, 92,100,108,112,116,116,112,108,100, 92, 84, 76, 64, 52, 40, 32, 24, 16,  8,  4,  0,
205   0,  4,  4, 12, 16, 24, 32, 40, 52, 60, 68, 76, 80, 88, 88, 92, 92, 88, 88, 80, 76, 68, 60, 52, 40, 32, 24, 16, 12,  4,  4,  0,
206   0,  4,  4,  8, 12, 20, 24, 32, 40, 44, 52, 56, 60, 64, 68, 72, 72, 68, 64, 60, 56, 52, 44, 40, 32, 24, 20, 12,  8,  4,  4,  0,
207   0,  0,  4,  8,  8, 12, 16, 24, 28, 32, 36, 40, 44, 48, 48, 48, 48, 48, 48, 44, 40, 36, 32, 28, 24, 16, 12,  8,  8,  4,  0,  0,
208   0,  0,  4,  4,  8,  8, 12, 16, 16, 20, 24, 24, 28, 28, 32, 32, 32, 32, 28, 28, 24, 24, 20, 16, 16, 12,  8,  8,  4,  4,  0,  0,
209   0,  0,  0,  4,  4,  4,  4,  8,  8, 12, 12, 12, 16, 16, 16, 16, 16, 16, 16, 16, 12, 12, 12,  8,  8,  4,  4,  4,  4,  0,  0,  0,
210   0,  0,  0,  0,  0,  4,  4,  4,  4,  4,  4,  4,  4,  8,  8,  8,  8,  8,  8,  4,  4,  4,  4,  4,  4,  4,  4,  0,  0,  0,  0,  0,
211   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,
212 //error:0.000022
213 };
214 static const uint8_t obmc16[256]={
215   0,  0,  0,  0,  0,  0,  4,  4,  4,  4,  0,  0,  0,  0,  0,  0,
216   0,  4,  4,  8, 16, 20, 20, 24, 24, 20, 20, 16,  8,  4,  4,  0,
217   0,  4, 16, 24, 36, 44, 52, 60, 60, 52, 44, 36, 24, 16,  4,  0,
218   0,  8, 24, 44, 60, 80, 96,104,104, 96, 80, 60, 44, 24,  8,  0,
219   0, 16, 36, 60, 92,116,136,152,152,136,116, 92, 60, 36, 16,  0,
220   0, 20, 44, 80,116,152,180,196,196,180,152,116, 80, 44, 20,  0,
221   4, 20, 52, 96,136,180,212,228,228,212,180,136, 96, 52, 20,  4,
222   4, 24, 60,104,152,196,228,248,248,228,196,152,104, 60, 24,  4,
223   4, 24, 60,104,152,196,228,248,248,228,196,152,104, 60, 24,  4,
224   4, 20, 52, 96,136,180,212,228,228,212,180,136, 96, 52, 20,  4,
225   0, 20, 44, 80,116,152,180,196,196,180,152,116, 80, 44, 20,  0,
226   0, 16, 36, 60, 92,116,136,152,152,136,116, 92, 60, 36, 16,  0,
227   0,  8, 24, 44, 60, 80, 96,104,104, 96, 80, 60, 44, 24,  8,  0,
228   0,  4, 16, 24, 36, 44, 52, 60, 60, 52, 44, 36, 24, 16,  4,  0,
229   0,  4,  4,  8, 16, 20, 20, 24, 24, 20, 20, 16,  8,  4,  4,  0,
230   0,  0,  0,  0,  0,  0,  4,  4,  4,  4,  0,  0,  0,  0,  0,  0,
231 //error:0.000033
232 };
233 #elif 1 // 64*linear
234 static const uint8_t obmc32[1024]={
235   0,  0,  0,  0,  4,  4,  4,  4,  4,  4,  4,  4,  8,  8,  8,  8,  8,  8,  8,  8,  4,  4,  4,  4,  4,  4,  4,  4,  0,  0,  0,  0,
236   0,  4,  4,  4,  8,  8,  8, 12, 12, 16, 16, 16, 20, 20, 20, 24, 24, 20, 20, 20, 16, 16, 16, 12, 12,  8,  8,  8,  4,  4,  4,  0,
237   0,  4,  8,  8, 12, 12, 16, 20, 20, 24, 28, 28, 32, 32, 36, 40, 40, 36, 32, 32, 28, 28, 24, 20, 20, 16, 12, 12,  8,  8,  4,  0,
238   0,  4,  8, 12, 16, 20, 24, 28, 28, 32, 36, 40, 44, 48, 52, 56, 56, 52, 48, 44, 40, 36, 32, 28, 28, 24, 20, 16, 12,  8,  4,  0,
239   4,  8, 12, 16, 20, 24, 28, 32, 40, 44, 48, 52, 56, 60, 64, 68, 68, 64, 60, 56, 52, 48, 44, 40, 32, 28, 24, 20, 16, 12,  8,  4,
240   4,  8, 12, 20, 24, 32, 36, 40, 48, 52, 56, 64, 68, 76, 80, 84, 84, 80, 76, 68, 64, 56, 52, 48, 40, 36, 32, 24, 20, 12,  8,  4,
241   4,  8, 16, 24, 28, 36, 44, 48, 56, 60, 68, 76, 80, 88, 96,100,100, 96, 88, 80, 76, 68, 60, 56, 48, 44, 36, 28, 24, 16,  8,  4,
242   4, 12, 20, 28, 32, 40, 48, 56, 64, 72, 80, 88, 92,100,108,116,116,108,100, 92, 88, 80, 72, 64, 56, 48, 40, 32, 28, 20, 12,  4,
243   4, 12, 20, 28, 40, 48, 56, 64, 72, 80, 88, 96,108,116,124,132,132,124,116,108, 96, 88, 80, 72, 64, 56, 48, 40, 28, 20, 12,  4,
244   4, 16, 24, 32, 44, 52, 60, 72, 80, 92,100,108,120,128,136,148,148,136,128,120,108,100, 92, 80, 72, 60, 52, 44, 32, 24, 16,  4,
245   4, 16, 28, 36, 48, 56, 68, 80, 88,100,112,120,132,140,152,164,164,152,140,132,120,112,100, 88, 80, 68, 56, 48, 36, 28, 16,  4,
246   4, 16, 28, 40, 52, 64, 76, 88, 96,108,120,132,144,156,168,180,180,168,156,144,132,120,108, 96, 88, 76, 64, 52, 40, 28, 16,  4,
247   8, 20, 32, 44, 56, 68, 80, 92,108,120,132,144,156,168,180,192,192,180,168,156,144,132,120,108, 92, 80, 68, 56, 44, 32, 20,  8,
248   8, 20, 32, 48, 60, 76, 88,100,116,128,140,156,168,184,196,208,208,196,184,168,156,140,128,116,100, 88, 76, 60, 48, 32, 20,  8,
249   8, 20, 36, 52, 64, 80, 96,108,124,136,152,168,180,196,212,224,224,212,196,180,168,152,136,124,108, 96, 80, 64, 52, 36, 20,  8,
250   8, 24, 40, 56, 68, 84,100,116,132,148,164,180,192,208,224,240,240,224,208,192,180,164,148,132,116,100, 84, 68, 56, 40, 24,  8,
251   8, 24, 40, 56, 68, 84,100,116,132,148,164,180,192,208,224,240,240,224,208,192,180,164,148,132,116,100, 84, 68, 56, 40, 24,  8,
252   8, 20, 36, 52, 64, 80, 96,108,124,136,152,168,180,196,212,224,224,212,196,180,168,152,136,124,108, 96, 80, 64, 52, 36, 20,  8,
253   8, 20, 32, 48, 60, 76, 88,100,116,128,140,156,168,184,196,208,208,196,184,168,156,140,128,116,100, 88, 76, 60, 48, 32, 20,  8,
254   8, 20, 32, 44, 56, 68, 80, 92,108,120,132,144,156,168,180,192,192,180,168,156,144,132,120,108, 92, 80, 68, 56, 44, 32, 20,  8,
255   4, 16, 28, 40, 52, 64, 76, 88, 96,108,120,132,144,156,168,180,180,168,156,144,132,120,108, 96, 88, 76, 64, 52, 40, 28, 16,  4,
256   4, 16, 28, 36, 48, 56, 68, 80, 88,100,112,120,132,140,152,164,164,152,140,132,120,112,100, 88, 80, 68, 56, 48, 36, 28, 16,  4,
257   4, 16, 24, 32, 44, 52, 60, 72, 80, 92,100,108,120,128,136,148,148,136,128,120,108,100, 92, 80, 72, 60, 52, 44, 32, 24, 16,  4,
258   4, 12, 20, 28, 40, 48, 56, 64, 72, 80, 88, 96,108,116,124,132,132,124,116,108, 96, 88, 80, 72, 64, 56, 48, 40, 28, 20, 12,  4,
259   4, 12, 20, 28, 32, 40, 48, 56, 64, 72, 80, 88, 92,100,108,116,116,108,100, 92, 88, 80, 72, 64, 56, 48, 40, 32, 28, 20, 12,  4,
260   4,  8, 16, 24, 28, 36, 44, 48, 56, 60, 68, 76, 80, 88, 96,100,100, 96, 88, 80, 76, 68, 60, 56, 48, 44, 36, 28, 24, 16,  8,  4,
261   4,  8, 12, 20, 24, 32, 36, 40, 48, 52, 56, 64, 68, 76, 80, 84, 84, 80, 76, 68, 64, 56, 52, 48, 40, 36, 32, 24, 20, 12,  8,  4,
262   4,  8, 12, 16, 20, 24, 28, 32, 40, 44, 48, 52, 56, 60, 64, 68, 68, 64, 60, 56, 52, 48, 44, 40, 32, 28, 24, 20, 16, 12,  8,  4,
263   0,  4,  8, 12, 16, 20, 24, 28, 28, 32, 36, 40, 44, 48, 52, 56, 56, 52, 48, 44, 40, 36, 32, 28, 28, 24, 20, 16, 12,  8,  4,  0,
264   0,  4,  8,  8, 12, 12, 16, 20, 20, 24, 28, 28, 32, 32, 36, 40, 40, 36, 32, 32, 28, 28, 24, 20, 20, 16, 12, 12,  8,  8,  4,  0,
265   0,  4,  4,  4,  8,  8,  8, 12, 12, 16, 16, 16, 20, 20, 20, 24, 24, 20, 20, 20, 16, 16, 16, 12, 12,  8,  8,  8,  4,  4,  4,  0,
266   0,  0,  0,  0,  4,  4,  4,  4,  4,  4,  4,  4,  8,  8,  8,  8,  8,  8,  8,  8,  4,  4,  4,  4,  4,  4,  4,  4,  0,  0,  0,  0,
267  //error:0.000020
268 };
269 static const uint8_t obmc16[256]={
270   0,  4,  4,  8,  8, 12, 12, 16, 16, 12, 12,  8,  8,  4,  4,  0,
271   4,  8, 16, 20, 28, 32, 40, 44, 44, 40, 32, 28, 20, 16,  8,  4,
272   4, 16, 24, 36, 44, 56, 64, 76, 76, 64, 56, 44, 36, 24, 16,  4,
273   8, 20, 36, 48, 64, 76, 92,104,104, 92, 76, 64, 48, 36, 20,  8,
274   8, 28, 44, 64, 80,100,116,136,136,116,100, 80, 64, 44, 28,  8,
275  12, 32, 56, 76,100,120,144,164,164,144,120,100, 76, 56, 32, 12,
276  12, 40, 64, 92,116,144,168,196,196,168,144,116, 92, 64, 40, 12,
277  16, 44, 76,104,136,164,196,224,224,196,164,136,104, 76, 44, 16,
278  16, 44, 76,104,136,164,196,224,224,196,164,136,104, 76, 44, 16,
279  12, 40, 64, 92,116,144,168,196,196,168,144,116, 92, 64, 40, 12,
280  12, 32, 56, 76,100,120,144,164,164,144,120,100, 76, 56, 32, 12,
281   8, 28, 44, 64, 80,100,116,136,136,116,100, 80, 64, 44, 28,  8,
282   8, 20, 36, 48, 64, 76, 92,104,104, 92, 76, 64, 48, 36, 20,  8,
283   4, 16, 24, 36, 44, 56, 64, 76, 76, 64, 56, 44, 36, 24, 16,  4,
284   4,  8, 16, 20, 28, 32, 40, 44, 44, 40, 32, 28, 20, 16,  8,  4,
285   0,  4,  4,  8,  8, 12, 12, 16, 16, 12, 12,  8,  8,  4,  4,  0,
286 //error:0.000015
287 };
288 #else //64*cos
289 static const uint8_t obmc32[1024]={
290   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,
291   0,  0,  0,  0,  0,  0,  4,  4,  4,  4,  4,  4,  4,  4,  8,  4,  4,  8,  4,  4,  4,  4,  4,  4,  4,  4,  0,  0,  0,  0,  0,  0,
292   0,  0,  0,  4,  4,  4,  4,  8,  8, 12, 12, 12, 12, 16, 16, 16, 16, 16, 16, 12, 12, 12, 12,  8,  8,  4,  4,  4,  4,  0,  0,  0,
293   0,  0,  4,  4,  4,  8,  8, 12, 16, 20, 20, 24, 28, 28, 28, 28, 28, 28, 28, 28, 24, 20, 20, 16, 12,  8,  8,  4,  4,  4,  0,  0,
294   0,  0,  4,  4,  8, 12, 16, 20, 24, 28, 36, 40, 44, 44, 48, 48, 48, 48, 44, 44, 40, 36, 28, 24, 20, 16, 12,  8,  4,  4,  0,  0,
295   0,  0,  4,  8, 12, 20, 24, 32, 36, 44, 48, 56, 60, 64, 68, 68, 68, 68, 64, 60, 56, 48, 44, 36, 32, 24, 20, 12,  8,  4,  0,  0,
296   0,  4,  4,  8, 16, 24, 32, 40, 48, 60, 68, 76, 80, 84, 88, 92, 92, 88, 84, 80, 76, 68, 60, 48, 40, 32, 24, 16,  8,  4,  4,  0,
297   0,  4,  8, 12, 20, 32, 40, 52, 64, 76, 84, 96,104,108,112,116,116,112,108,104, 96, 84, 76, 64, 52, 40, 32, 20, 12,  8,  4,  0,
298   0,  4,  8, 16, 24, 36, 48, 64, 76, 92,104,116,124,132,136,140,140,136,132,124,116,104, 92, 76, 64, 48, 36, 24, 16,  8,  4,  0,
299   0,  4, 12, 20, 28, 44, 60, 76, 92,104,120,136,148,156,160,164,164,160,156,148,136,120,104, 92, 76, 60, 44, 28, 20, 12,  4,  0,
300   0,  4, 12, 20, 36, 48, 68, 84,104,120,140,152,168,176,184,188,188,184,176,168,152,140,120,104, 84, 68, 48, 36, 20, 12,  4,  0,
301   0,  4, 12, 24, 36, 56, 76, 96,116,136,152,172,184,196,204,208,208,204,196,184,172,152,136,116, 96, 76, 56, 36, 24, 12,  4,  0,
302   0,  4, 12, 24, 44, 60, 80,104,124,148,168,184,200,212,224,228,228,224,212,200,184,168,148,124,104, 80, 60, 44, 24, 12,  4,  0,
303   0,  4, 12, 28, 44, 64, 84,108,132,156,176,196,212,228,236,240,240,236,228,212,196,176,156,132,108, 84, 64, 44, 28, 12,  4,  0,
304   0,  4, 16, 28, 48, 68, 88,112,136,160,184,204,224,236,244,252,252,244,236,224,204,184,160,136,112, 88, 68, 48, 28, 16,  4,  0,
305   1,  4, 16, 28, 48, 68, 92,116,140,164,188,208,228,240,252,255,255,252,240,228,208,188,164,140,116, 92, 68, 48, 28, 16,  4,  1,
306   1,  4, 16, 28, 48, 68, 92,116,140,164,188,208,228,240,252,255,255,252,240,228,208,188,164,140,116, 92, 68, 48, 28, 16,  4,  1,
307   0,  4, 16, 28, 48, 68, 88,112,136,160,184,204,224,236,244,252,252,244,236,224,204,184,160,136,112, 88, 68, 48, 28, 16,  4,  0,
308   0,  4, 12, 28, 44, 64, 84,108,132,156,176,196,212,228,236,240,240,236,228,212,196,176,156,132,108, 84, 64, 44, 28, 12,  4,  0,
309   0,  4, 12, 24, 44, 60, 80,104,124,148,168,184,200,212,224,228,228,224,212,200,184,168,148,124,104, 80, 60, 44, 24, 12,  4,  0,
310   0,  4, 12, 24, 36, 56, 76, 96,116,136,152,172,184,196,204,208,208,204,196,184,172,152,136,116, 96, 76, 56, 36, 24, 12,  4,  0,
311   0,  4, 12, 20, 36, 48, 68, 84,104,120,140,152,168,176,184,188,188,184,176,168,152,140,120,104, 84, 68, 48, 36, 20, 12,  4,  0,
312   0,  4, 12, 20, 28, 44, 60, 76, 92,104,120,136,148,156,160,164,164,160,156,148,136,120,104, 92, 76, 60, 44, 28, 20, 12,  4,  0,
313   0,  4,  8, 16, 24, 36, 48, 64, 76, 92,104,116,124,132,136,140,140,136,132,124,116,104, 92, 76, 64, 48, 36, 24, 16,  8,  4,  0,
314   0,  4,  8, 12, 20, 32, 40, 52, 64, 76, 84, 96,104,108,112,116,116,112,108,104, 96, 84, 76, 64, 52, 40, 32, 20, 12,  8,  4,  0,
315   0,  4,  4,  8, 16, 24, 32, 40, 48, 60, 68, 76, 80, 84, 88, 92, 92, 88, 84, 80, 76, 68, 60, 48, 40, 32, 24, 16,  8,  4,  4,  0,
316   0,  0,  4,  8, 12, 20, 24, 32, 36, 44, 48, 56, 60, 64, 68, 68, 68, 68, 64, 60, 56, 48, 44, 36, 32, 24, 20, 12,  8,  4,  0,  0,
317   0,  0,  4,  4,  8, 12, 16, 20, 24, 28, 36, 40, 44, 44, 48, 48, 48, 48, 44, 44, 40, 36, 28, 24, 20, 16, 12,  8,  4,  4,  0,  0,
318   0,  0,  4,  4,  4,  8,  8, 12, 16, 20, 20, 24, 28, 28, 28, 28, 28, 28, 28, 28, 24, 20, 20, 16, 12,  8,  8,  4,  4,  4,  0,  0,
319   0,  0,  0,  4,  4,  4,  4,  8,  8, 12, 12, 12, 12, 16, 16, 16, 16, 16, 16, 12, 12, 12, 12,  8,  8,  4,  4,  4,  4,  0,  0,  0,
320   0,  0,  0,  0,  0,  0,  4,  4,  4,  4,  4,  4,  4,  4,  8,  4,  4,  8,  4,  4,  4,  4,  4,  4,  4,  4,  0,  0,  0,  0,  0,  0,
321   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,
322 //error:0.000022
323 };
324 static const uint8_t obmc16[256]={
325   0,  0,  0,  0,  0,  4,  4,  4,  4,  4,  4,  0,  0,  0,  0,  0,
326   0,  0,  4,  8, 12, 16, 20, 20, 20, 20, 16, 12,  8,  4,  0,  0,
327   0,  4, 12, 24, 32, 44, 52, 56, 56, 52, 44, 32, 24, 12,  4,  0,
328   0,  8, 24, 40, 60, 80, 96,104,104, 96, 80, 60, 40, 24,  8,  0,
329   0, 12, 32, 64, 92,120,140,152,152,140,120, 92, 64, 32, 12,  0,
330   4, 16, 44, 80,120,156,184,196,196,184,156,120, 80, 44, 16,  4,
331   4, 20, 52, 96,140,184,216,232,232,216,184,140, 96, 52, 20,  4,
332   0, 20, 56,104,152,196,232,252,252,232,196,152,104, 56, 20,  0,
333   0, 20, 56,104,152,196,232,252,252,232,196,152,104, 56, 20,  0,
334   4, 20, 52, 96,140,184,216,232,232,216,184,140, 96, 52, 20,  4,
335   4, 16, 44, 80,120,156,184,196,196,184,156,120, 80, 44, 16,  4,
336   0, 12, 32, 64, 92,120,140,152,152,140,120, 92, 64, 32, 12,  0,
337   0,  8, 24, 40, 60, 80, 96,104,104, 96, 80, 60, 40, 24,  8,  0,
338   0,  4, 12, 24, 32, 44, 52, 56, 56, 52, 44, 32, 24, 12,  4,  0,
339   0,  0,  4,  8, 12, 16, 20, 20, 20, 20, 16, 12,  8,  4,  0,  0,
340   0,  0,  0,  0,  0,  4,  4,  4,  4,  4,  4,  0,  0,  0,  0,  0,
341 //error:0.000022
342 };
343 #endif /* 0 */
344
345 //linear *64
346 static const uint8_t obmc8[64]={
347   4, 12, 20, 28, 28, 20, 12,  4,
348  12, 36, 60, 84, 84, 60, 36, 12,
349  20, 60,100,140,140,100, 60, 20,
350  28, 84,140,196,196,140, 84, 28,
351  28, 84,140,196,196,140, 84, 28,
352  20, 60,100,140,140,100, 60, 20,
353  12, 36, 60, 84, 84, 60, 36, 12,
354   4, 12, 20, 28, 28, 20, 12,  4,
355 //error:0.000000
356 };
357
358 //linear *64
359 static const uint8_t obmc4[16]={
360  16, 48, 48, 16,
361  48,144,144, 48,
362  48,144,144, 48,
363  16, 48, 48, 16,
364 //error:0.000000
365 };
366
367 static const uint8_t * const obmc_tab[4]={
368     obmc32, obmc16, obmc8, obmc4
369 };
370
371 static int scale_mv_ref[MAX_REF_FRAMES][MAX_REF_FRAMES];
372
373 typedef struct BlockNode{
374     int16_t mx;
375     int16_t my;
376     uint8_t ref;
377     uint8_t color[3];
378     uint8_t type;
379 //#define TYPE_SPLIT    1
380 #define BLOCK_INTRA   1
381 #define BLOCK_OPT     2
382 //#define TYPE_NOCOLOR  4
383     uint8_t level; //FIXME merge into type?
384 }BlockNode;
385
386 static const BlockNode null_block= { //FIXME add border maybe
387     .color= {128,128,128},
388     .mx= 0,
389     .my= 0,
390     .ref= 0,
391     .type= 0,
392     .level= 0,
393 };
394
395 #define LOG2_MB_SIZE 4
396 #define MB_SIZE (1<<LOG2_MB_SIZE)
397 #define ENCODER_EXTRA_BITS 4
398 #define HTAPS_MAX 8
399
400 typedef struct x_and_coeff{
401     int16_t x;
402     uint16_t coeff;
403 } x_and_coeff;
404
405 typedef struct SubBand{
406     int level;
407     int stride;
408     int width;
409     int height;
410     int qlog;        ///< log(qscale)/log[2^(1/6)]
411     DWTELEM *buf;
412     IDWTELEM *ibuf;
413     int buf_x_offset;
414     int buf_y_offset;
415     int stride_line; ///< Stride measured in lines, not pixels.
416     x_and_coeff * x_coeff;
417     struct SubBand *parent;
418     uint8_t state[/*7*2*/ 7 + 512][32];
419 }SubBand;
420
421 typedef struct Plane{
422     int width;
423     int height;
424     SubBand band[MAX_DECOMPOSITIONS][4];
425
426     int htaps;
427     int8_t hcoeff[HTAPS_MAX/2];
428     int diag_mc;
429     int fast_mc;
430
431     int last_htaps;
432     int8_t last_hcoeff[HTAPS_MAX/2];
433     int last_diag_mc;
434 }Plane;
435
436 typedef struct SnowContext{
437
438     AVCodecContext *avctx;
439     RangeCoder c;
440     DSPContext dsp;
441     AVFrame new_picture;
442     AVFrame input_picture;              ///< new_picture with the internal linesizes
443     AVFrame current_picture;
444     AVFrame last_picture[MAX_REF_FRAMES];
445     uint8_t *halfpel_plane[MAX_REF_FRAMES][4][4];
446     AVFrame mconly_picture;
447 //     uint8_t q_context[16];
448     uint8_t header_state[32];
449     uint8_t block_state[128 + 32*128];
450     int keyframe;
451     int always_reset;
452     int version;
453     int spatial_decomposition_type;
454     int last_spatial_decomposition_type;
455     int temporal_decomposition_type;
456     int spatial_decomposition_count;
457     int last_spatial_decomposition_count;
458     int temporal_decomposition_count;
459     int max_ref_frames;
460     int ref_frames;
461     int16_t (*ref_mvs[MAX_REF_FRAMES])[2];
462     uint32_t *ref_scores[MAX_REF_FRAMES];
463     DWTELEM *spatial_dwt_buffer;
464     IDWTELEM *spatial_idwt_buffer;
465     int colorspace_type;
466     int chroma_h_shift;
467     int chroma_v_shift;
468     int spatial_scalability;
469     int qlog;
470     int last_qlog;
471     int lambda;
472     int lambda2;
473     int pass1_rc;
474     int mv_scale;
475     int last_mv_scale;
476     int qbias;
477     int last_qbias;
478 #define QBIAS_SHIFT 3
479     int b_width;
480     int b_height;
481     int block_max_depth;
482     int last_block_max_depth;
483     Plane plane[MAX_PLANES];
484     BlockNode *block;
485 #define ME_CACHE_SIZE 1024
486     int me_cache[ME_CACHE_SIZE];
487     int me_cache_generation;
488     slice_buffer sb;
489
490     MpegEncContext m; // needed for motion estimation, should not be used for anything else, the idea is to eventually make the motion estimation independent of MpegEncContext, so this will be removed then (FIXME/XXX)
491
492     uint8_t *scratchbuf;
493 }SnowContext;
494
495 typedef struct {
496     IDWTELEM *b0;
497     IDWTELEM *b1;
498     IDWTELEM *b2;
499     IDWTELEM *b3;
500     int y;
501 } DWTCompose;
502
503 #define slice_buffer_get_line(slice_buf, line_num) ((slice_buf)->line[line_num] ? (slice_buf)->line[line_num] : slice_buffer_load_line((slice_buf), (line_num)))
504 //#define slice_buffer_get_line(slice_buf, line_num) (slice_buffer_load_line((slice_buf), (line_num)))
505
506 static void slice_buffer_init(slice_buffer * buf, int line_count, int max_allocated_lines, int line_width, IDWTELEM * base_buffer)
507 {
508     int i;
509
510     buf->base_buffer = base_buffer;
511     buf->line_count = line_count;
512     buf->line_width = line_width;
513     buf->data_count = max_allocated_lines;
514     buf->line = av_mallocz (sizeof(IDWTELEM *) * line_count);
515     buf->data_stack = av_malloc (sizeof(IDWTELEM *) * max_allocated_lines);
516
517     for(i = 0; i < max_allocated_lines; i++){
518         buf->data_stack[i] = av_malloc (sizeof(IDWTELEM) * line_width);
519     }
520
521     buf->data_stack_top = max_allocated_lines - 1;
522 }
523
524 static IDWTELEM * slice_buffer_load_line(slice_buffer * buf, int line)
525 {
526     IDWTELEM * buffer;
527
528     assert(buf->data_stack_top >= 0);
529 //  assert(!buf->line[line]);
530     if (buf->line[line])
531         return buf->line[line];
532
533     buffer = buf->data_stack[buf->data_stack_top];
534     buf->data_stack_top--;
535     buf->line[line] = buffer;
536
537     return buffer;
538 }
539
540 static void slice_buffer_release(slice_buffer * buf, int line)
541 {
542     IDWTELEM * buffer;
543
544     assert(line >= 0 && line < buf->line_count);
545     assert(buf->line[line]);
546
547     buffer = buf->line[line];
548     buf->data_stack_top++;
549     buf->data_stack[buf->data_stack_top] = buffer;
550     buf->line[line] = NULL;
551 }
552
553 static void slice_buffer_flush(slice_buffer * buf)
554 {
555     int i;
556     for(i = 0; i < buf->line_count; i++){
557         if (buf->line[i])
558             slice_buffer_release(buf, i);
559     }
560 }
561
562 static void slice_buffer_destroy(slice_buffer * buf)
563 {
564     int i;
565     slice_buffer_flush(buf);
566
567     for(i = buf->data_count - 1; i >= 0; i--){
568         av_freep(&buf->data_stack[i]);
569     }
570     av_freep(&buf->data_stack);
571     av_freep(&buf->line);
572 }
573
574 #ifdef __sgi
575 // Avoid a name clash on SGI IRIX
576 #undef qexp
577 #endif
578 #define QEXPSHIFT (7-FRAC_BITS+8) //FIXME try to change this to 0
579 static uint8_t qexp[QROOT];
580
581 static inline int mirror(int v, int m){
582     while((unsigned)v > (unsigned)m){
583         v=-v;
584         if(v<0) v+= 2*m;
585     }
586     return v;
587 }
588
589 static inline void put_symbol(RangeCoder *c, uint8_t *state, int v, int is_signed){
590     int i;
591
592     if(v){
593         const int a= FFABS(v);
594         const int e= av_log2(a);
595 #if 1
596         const int el= FFMIN(e, 10);
597         put_rac(c, state+0, 0);
598
599         for(i=0; i<el; i++){
600             put_rac(c, state+1+i, 1);  //1..10
601         }
602         for(; i<e; i++){
603             put_rac(c, state+1+9, 1);  //1..10
604         }
605         put_rac(c, state+1+FFMIN(i,9), 0);
606
607         for(i=e-1; i>=el; i--){
608             put_rac(c, state+22+9, (a>>i)&1); //22..31
609         }
610         for(; i>=0; i--){
611             put_rac(c, state+22+i, (a>>i)&1); //22..31
612         }
613
614         if(is_signed)
615             put_rac(c, state+11 + el, v < 0); //11..21
616 #else
617
618         put_rac(c, state+0, 0);
619         if(e<=9){
620             for(i=0; i<e; i++){
621                 put_rac(c, state+1+i, 1);  //1..10
622             }
623             put_rac(c, state+1+i, 0);
624
625             for(i=e-1; i>=0; i--){
626                 put_rac(c, state+22+i, (a>>i)&1); //22..31
627             }
628
629             if(is_signed)
630                 put_rac(c, state+11 + e, v < 0); //11..21
631         }else{
632             for(i=0; i<e; i++){
633                 put_rac(c, state+1+FFMIN(i,9), 1);  //1..10
634             }
635             put_rac(c, state+1+9, 0);
636
637             for(i=e-1; i>=0; i--){
638                 put_rac(c, state+22+FFMIN(i,9), (a>>i)&1); //22..31
639             }
640
641             if(is_signed)
642                 put_rac(c, state+11 + 10, v < 0); //11..21
643         }
644 #endif /* 1 */
645     }else{
646         put_rac(c, state+0, 1);
647     }
648 }
649
650 static inline int get_symbol(RangeCoder *c, uint8_t *state, int is_signed){
651     if(get_rac(c, state+0))
652         return 0;
653     else{
654         int i, e, a;
655         e= 0;
656         while(get_rac(c, state+1 + FFMIN(e,9))){ //1..10
657             e++;
658         }
659
660         a= 1;
661         for(i=e-1; i>=0; i--){
662             a += a + get_rac(c, state+22 + FFMIN(i,9)); //22..31
663         }
664
665         e= -(is_signed && get_rac(c, state+11 + FFMIN(e,10))); //11..21
666         return (a^e)-e;
667     }
668 }
669
670 static inline void put_symbol2(RangeCoder *c, uint8_t *state, int v, int log2){
671     int i;
672     int r= log2>=0 ? 1<<log2 : 1;
673
674     assert(v>=0);
675     assert(log2>=-4);
676
677     while(v >= r){
678         put_rac(c, state+4+log2, 1);
679         v -= r;
680         log2++;
681         if(log2>0) r+=r;
682     }
683     put_rac(c, state+4+log2, 0);
684
685     for(i=log2-1; i>=0; i--){
686         put_rac(c, state+31-i, (v>>i)&1);
687     }
688 }
689
690 static inline int get_symbol2(RangeCoder *c, uint8_t *state, int log2){
691     int i;
692     int r= log2>=0 ? 1<<log2 : 1;
693     int v=0;
694
695     assert(log2>=-4);
696
697     while(get_rac(c, state+4+log2)){
698         v+= r;
699         log2++;
700         if(log2>0) r+=r;
701     }
702
703     for(i=log2-1; i>=0; i--){
704         v+= get_rac(c, state+31-i)<<i;
705     }
706
707     return v;
708 }
709
710 static av_always_inline void
711 lift(DWTELEM *dst, DWTELEM *src, DWTELEM *ref,
712      int dst_step, int src_step, int ref_step,
713      int width, int mul, int add, int shift,
714      int highpass, int inverse){
715     const int mirror_left= !highpass;
716     const int mirror_right= (width&1) ^ highpass;
717     const int w= (width>>1) - 1 + (highpass & width);
718     int i;
719
720 #define LIFT(src, ref, inv) ((src) + ((inv) ? - (ref) : + (ref)))
721     if(mirror_left){
722         dst[0] = LIFT(src[0], ((mul*2*ref[0]+add)>>shift), inverse);
723         dst += dst_step;
724         src += src_step;
725     }
726
727     for(i=0; i<w; i++){
728         dst[i*dst_step] =
729             LIFT(src[i*src_step],
730                  ((mul*(ref[i*ref_step] + ref[(i+1)*ref_step])+add)>>shift),
731                  inverse);
732     }
733
734     if(mirror_right){
735         dst[w*dst_step] =
736             LIFT(src[w*src_step],
737                  ((mul*2*ref[w*ref_step]+add)>>shift),
738                  inverse);
739     }
740 }
741
742 static av_always_inline void
743 inv_lift(IDWTELEM *dst, IDWTELEM *src, IDWTELEM *ref,
744          int dst_step, int src_step, int ref_step,
745          int width, int mul, int add, int shift,
746          int highpass, int inverse){
747     const int mirror_left= !highpass;
748     const int mirror_right= (width&1) ^ highpass;
749     const int w= (width>>1) - 1 + (highpass & width);
750     int i;
751
752 #define LIFT(src, ref, inv) ((src) + ((inv) ? - (ref) : + (ref)))
753     if(mirror_left){
754         dst[0] = LIFT(src[0], ((mul*2*ref[0]+add)>>shift), inverse);
755         dst += dst_step;
756         src += src_step;
757     }
758
759     for(i=0; i<w; i++){
760         dst[i*dst_step] =
761             LIFT(src[i*src_step],
762                  ((mul*(ref[i*ref_step] + ref[(i+1)*ref_step])+add)>>shift),
763                  inverse);
764     }
765
766     if(mirror_right){
767         dst[w*dst_step] =
768             LIFT(src[w*src_step],
769                  ((mul*2*ref[w*ref_step]+add)>>shift),
770                  inverse);
771     }
772 }
773
774 #ifndef liftS
775 static av_always_inline void
776 liftS(DWTELEM *dst, DWTELEM *src, DWTELEM *ref,
777       int dst_step, int src_step, int ref_step,
778       int width, int mul, int add, int shift,
779       int highpass, int inverse){
780     const int mirror_left= !highpass;
781     const int mirror_right= (width&1) ^ highpass;
782     const int w= (width>>1) - 1 + (highpass & width);
783     int i;
784
785     assert(shift == 4);
786 #define LIFTS(src, ref, inv) \
787         ((inv) ? \
788             (src) + (((ref) + 4*(src))>>shift): \
789             -((-16*(src) + (ref) + add/4 + 1 + (5<<25))/(5*4) - (1<<23)))
790     if(mirror_left){
791         dst[0] = LIFTS(src[0], mul*2*ref[0]+add, inverse);
792         dst += dst_step;
793         src += src_step;
794     }
795
796     for(i=0; i<w; i++){
797         dst[i*dst_step] =
798             LIFTS(src[i*src_step],
799                   mul*(ref[i*ref_step] + ref[(i+1)*ref_step])+add,
800                   inverse);
801     }
802
803     if(mirror_right){
804         dst[w*dst_step] =
805             LIFTS(src[w*src_step], mul*2*ref[w*ref_step]+add, inverse);
806     }
807 }
808 static av_always_inline void
809 inv_liftS(IDWTELEM *dst, IDWTELEM *src, IDWTELEM *ref,
810           int dst_step, int src_step, int ref_step,
811           int width, int mul, int add, int shift,
812           int highpass, int inverse){
813     const int mirror_left= !highpass;
814     const int mirror_right= (width&1) ^ highpass;
815     const int w= (width>>1) - 1 + (highpass & width);
816     int i;
817
818     assert(shift == 4);
819 #define LIFTS(src, ref, inv) \
820     ((inv) ? \
821         (src) + (((ref) + 4*(src))>>shift): \
822         -((-16*(src) + (ref) + add/4 + 1 + (5<<25))/(5*4) - (1<<23)))
823     if(mirror_left){
824         dst[0] = LIFTS(src[0], mul*2*ref[0]+add, inverse);
825         dst += dst_step;
826         src += src_step;
827     }
828
829     for(i=0; i<w; i++){
830         dst[i*dst_step] =
831             LIFTS(src[i*src_step],
832                   mul*(ref[i*ref_step] + ref[(i+1)*ref_step])+add,
833                   inverse);
834     }
835
836     if(mirror_right){
837         dst[w*dst_step] =
838             LIFTS(src[w*src_step], mul*2*ref[w*ref_step]+add, inverse);
839     }
840 }
841 #endif /* ! liftS */
842
843 static void horizontal_decompose53i(DWTELEM *b, int width){
844     DWTELEM temp[width];
845     const int width2= width>>1;
846     int x;
847     const int w2= (width+1)>>1;
848
849     for(x=0; x<width2; x++){
850         temp[x   ]= b[2*x    ];
851         temp[x+w2]= b[2*x + 1];
852     }
853     if(width&1)
854         temp[x   ]= b[2*x    ];
855 #if 0
856     {
857     int A1,A2,A3,A4;
858     A2= temp[1       ];
859     A4= temp[0       ];
860     A1= temp[0+width2];
861     A1 -= (A2 + A4)>>1;
862     A4 += (A1 + 1)>>1;
863     b[0+width2] = A1;
864     b[0       ] = A4;
865     for(x=1; x+1<width2; x+=2){
866         A3= temp[x+width2];
867         A4= temp[x+1     ];
868         A3 -= (A2 + A4)>>1;
869         A2 += (A1 + A3 + 2)>>2;
870         b[x+width2] = A3;
871         b[x       ] = A2;
872
873         A1= temp[x+1+width2];
874         A2= temp[x+2       ];
875         A1 -= (A2 + A4)>>1;
876         A4 += (A1 + A3 + 2)>>2;
877         b[x+1+width2] = A1;
878         b[x+1       ] = A4;
879     }
880     A3= temp[width-1];
881     A3 -= A2;
882     A2 += (A1 + A3 + 2)>>2;
883     b[width -1] = A3;
884     b[width2-1] = A2;
885     }
886 #else
887     lift(b+w2, temp+w2, temp, 1, 1, 1, width, -1, 0, 1, 1, 0);
888     lift(b   , temp   , b+w2, 1, 1, 1, width,  1, 2, 2, 0, 0);
889 #endif /* 0 */
890 }
891
892 static void vertical_decompose53iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
893     int i;
894
895     for(i=0; i<width; i++){
896         b1[i] -= (b0[i] + b2[i])>>1;
897     }
898 }
899
900 static void vertical_decompose53iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
901     int i;
902
903     for(i=0; i<width; i++){
904         b1[i] += (b0[i] + b2[i] + 2)>>2;
905     }
906 }
907
908 static void spatial_decompose53i(DWTELEM *buffer, int width, int height, int stride){
909     int y;
910     DWTELEM *b0= buffer + mirror(-2-1, height-1)*stride;
911     DWTELEM *b1= buffer + mirror(-2  , height-1)*stride;
912
913     for(y=-2; y<height; y+=2){
914         DWTELEM *b2= buffer + mirror(y+1, height-1)*stride;
915         DWTELEM *b3= buffer + mirror(y+2, height-1)*stride;
916
917         if(y+1<(unsigned)height) horizontal_decompose53i(b2, width);
918         if(y+2<(unsigned)height) horizontal_decompose53i(b3, width);
919
920         if(y+1<(unsigned)height) vertical_decompose53iH0(b1, b2, b3, width);
921         if(y+0<(unsigned)height) vertical_decompose53iL0(b0, b1, b2, width);
922
923         b0=b2;
924         b1=b3;
925     }
926 }
927
928 static void horizontal_decompose97i(DWTELEM *b, int width){
929     DWTELEM temp[width];
930     const int w2= (width+1)>>1;
931
932     lift (temp+w2, b    +1, b      , 1, 2, 2, width,  W_AM, W_AO, W_AS, 1, 1);
933     liftS(temp   , b      , temp+w2, 1, 2, 1, width,  W_BM, W_BO, W_BS, 0, 0);
934     lift (b   +w2, temp+w2, temp   , 1, 1, 1, width,  W_CM, W_CO, W_CS, 1, 0);
935     lift (b      , temp   , b   +w2, 1, 1, 1, width,  W_DM, W_DO, W_DS, 0, 0);
936 }
937
938
939 static void vertical_decompose97iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
940     int i;
941
942     for(i=0; i<width; i++){
943         b1[i] -= (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
944     }
945 }
946
947 static void vertical_decompose97iH1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
948     int i;
949
950     for(i=0; i<width; i++){
951         b1[i] += (W_CM*(b0[i] + b2[i])+W_CO)>>W_CS;
952     }
953 }
954
955 static void vertical_decompose97iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
956     int i;
957
958     for(i=0; i<width; i++){
959 #ifdef liftS
960         b1[i] -= (W_BM*(b0[i] + b2[i])+W_BO)>>W_BS;
961 #else
962         b1[i] = (16*4*b1[i] - 4*(b0[i] + b2[i]) + W_BO*5 + (5<<27)) / (5*16) - (1<<23);
963 #endif
964     }
965 }
966
967 static void vertical_decompose97iL1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
968     int i;
969
970     for(i=0; i<width; i++){
971         b1[i] += (W_DM*(b0[i] + b2[i])+W_DO)>>W_DS;
972     }
973 }
974
975 static void spatial_decompose97i(DWTELEM *buffer, int width, int height, int stride){
976     int y;
977     DWTELEM *b0= buffer + mirror(-4-1, height-1)*stride;
978     DWTELEM *b1= buffer + mirror(-4  , height-1)*stride;
979     DWTELEM *b2= buffer + mirror(-4+1, height-1)*stride;
980     DWTELEM *b3= buffer + mirror(-4+2, height-1)*stride;
981
982     for(y=-4; y<height; y+=2){
983         DWTELEM *b4= buffer + mirror(y+3, height-1)*stride;
984         DWTELEM *b5= buffer + mirror(y+4, height-1)*stride;
985
986         if(y+3<(unsigned)height) horizontal_decompose97i(b4, width);
987         if(y+4<(unsigned)height) horizontal_decompose97i(b5, width);
988
989         if(y+3<(unsigned)height) vertical_decompose97iH0(b3, b4, b5, width);
990         if(y+2<(unsigned)height) vertical_decompose97iL0(b2, b3, b4, width);
991         if(y+1<(unsigned)height) vertical_decompose97iH1(b1, b2, b3, width);
992         if(y+0<(unsigned)height) vertical_decompose97iL1(b0, b1, b2, width);
993
994         b0=b2;
995         b1=b3;
996         b2=b4;
997         b3=b5;
998     }
999 }
1000
1001 void ff_spatial_dwt(DWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count){
1002     int level;
1003
1004     for(level=0; level<decomposition_count; level++){
1005         switch(type){
1006         case DWT_97: spatial_decompose97i(buffer, width>>level, height>>level, stride<<level); break;
1007         case DWT_53: spatial_decompose53i(buffer, width>>level, height>>level, stride<<level); break;
1008         }
1009     }
1010 }
1011
1012 static void horizontal_compose53i(IDWTELEM *b, int width){
1013     IDWTELEM temp[width];
1014     const int width2= width>>1;
1015     const int w2= (width+1)>>1;
1016     int x;
1017
1018     for(x=0; x<width2; x++){
1019         temp[2*x    ]= b[x   ];
1020         temp[2*x + 1]= b[x+w2];
1021     }
1022     if(width&1)
1023         temp[2*x    ]= b[x   ];
1024
1025     b[0] = temp[0] - ((temp[1]+1)>>1);
1026     for(x=2; x<width-1; x+=2){
1027         b[x  ] = temp[x  ] - ((temp[x-1] + temp[x+1]+2)>>2);
1028         b[x-1] = temp[x-1] + ((b   [x-2] + b   [x  ]+1)>>1);
1029     }
1030     if(width&1){
1031         b[x  ] = temp[x  ] - ((temp[x-1]+1)>>1);
1032         b[x-1] = temp[x-1] + ((b   [x-2] + b  [x  ]+1)>>1);
1033     }else
1034         b[x-1] = temp[x-1] + b[x-2];
1035 }
1036
1037 static void vertical_compose53iH0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
1038     int i;
1039
1040     for(i=0; i<width; i++){
1041         b1[i] += (b0[i] + b2[i])>>1;
1042     }
1043 }
1044
1045 static void vertical_compose53iL0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
1046     int i;
1047
1048     for(i=0; i<width; i++){
1049         b1[i] -= (b0[i] + b2[i] + 2)>>2;
1050     }
1051 }
1052
1053 static void spatial_compose53i_buffered_init(DWTCompose *cs, slice_buffer * sb, int height, int stride_line){
1054     cs->b0 = slice_buffer_get_line(sb, mirror(-1-1, height-1) * stride_line);
1055     cs->b1 = slice_buffer_get_line(sb, mirror(-1  , height-1) * stride_line);
1056     cs->y = -1;
1057 }
1058
1059 static void spatial_compose53i_init(DWTCompose *cs, IDWTELEM *buffer, int height, int stride){
1060     cs->b0 = buffer + mirror(-1-1, height-1)*stride;
1061     cs->b1 = buffer + mirror(-1  , height-1)*stride;
1062     cs->y = -1;
1063 }
1064
1065 static void spatial_compose53i_dy_buffered(DWTCompose *cs, slice_buffer * sb, int width, int height, int stride_line){
1066     int y= cs->y;
1067
1068     IDWTELEM *b0= cs->b0;
1069     IDWTELEM *b1= cs->b1;
1070     IDWTELEM *b2= slice_buffer_get_line(sb, mirror(y+1, height-1) * stride_line);
1071     IDWTELEM *b3= slice_buffer_get_line(sb, mirror(y+2, height-1) * stride_line);
1072
1073     if(y+1<(unsigned)height && y<(unsigned)height){
1074         int x;
1075
1076         for(x=0; x<width; x++){
1077             b2[x] -= (b1[x] + b3[x] + 2)>>2;
1078             b1[x] += (b0[x] + b2[x])>>1;
1079         }
1080     }else{
1081         if(y+1<(unsigned)height) vertical_compose53iL0(b1, b2, b3, width);
1082         if(y+0<(unsigned)height) vertical_compose53iH0(b0, b1, b2, width);
1083     }
1084
1085         if(y-1<(unsigned)height) horizontal_compose53i(b0, width);
1086         if(y+0<(unsigned)height) horizontal_compose53i(b1, width);
1087
1088     cs->b0 = b2;
1089     cs->b1 = b3;
1090     cs->y += 2;
1091 }
1092
1093 static void spatial_compose53i_dy(DWTCompose *cs, IDWTELEM *buffer, int width, int height, int stride){
1094     int y= cs->y;
1095     IDWTELEM *b0= cs->b0;
1096     IDWTELEM *b1= cs->b1;
1097     IDWTELEM *b2= buffer + mirror(y+1, height-1)*stride;
1098     IDWTELEM *b3= buffer + mirror(y+2, height-1)*stride;
1099
1100         if(y+1<(unsigned)height) vertical_compose53iL0(b1, b2, b3, width);
1101         if(y+0<(unsigned)height) vertical_compose53iH0(b0, b1, b2, width);
1102
1103         if(y-1<(unsigned)height) horizontal_compose53i(b0, width);
1104         if(y+0<(unsigned)height) horizontal_compose53i(b1, width);
1105
1106     cs->b0 = b2;
1107     cs->b1 = b3;
1108     cs->y += 2;
1109 }
1110
1111 static void av_unused spatial_compose53i(IDWTELEM *buffer, int width, int height, int stride){
1112     DWTCompose cs;
1113     spatial_compose53i_init(&cs, buffer, height, stride);
1114     while(cs.y <= height)
1115         spatial_compose53i_dy(&cs, buffer, width, height, stride);
1116 }
1117
1118
1119 void ff_snow_horizontal_compose97i(IDWTELEM *b, int width){
1120     IDWTELEM temp[width];
1121     const int w2= (width+1)>>1;
1122
1123     inv_lift (temp   , b      , b   +w2, 1, 1, 1, width,  W_DM, W_DO, W_DS, 0, 1);
1124     inv_lift (temp+w2, b   +w2, temp   , 1, 1, 1, width,  W_CM, W_CO, W_CS, 1, 1);
1125     inv_liftS(b      , temp   , temp+w2, 2, 1, 1, width,  W_BM, W_BO, W_BS, 0, 1);
1126     inv_lift (b+1    , temp+w2, b      , 2, 1, 2, width,  W_AM, W_AO, W_AS, 1, 0);
1127 }
1128
1129 static void vertical_compose97iH0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
1130     int i;
1131
1132     for(i=0; i<width; i++){
1133         b1[i] += (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
1134     }
1135 }
1136
1137 static void vertical_compose97iH1(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
1138     int i;
1139
1140     for(i=0; i<width; i++){
1141         b1[i] -= (W_CM*(b0[i] + b2[i])+W_CO)>>W_CS;
1142     }
1143 }
1144
1145 static void vertical_compose97iL0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
1146     int i;
1147
1148     for(i=0; i<width; i++){
1149 #ifdef liftS
1150         b1[i] += (W_BM*(b0[i] + b2[i])+W_BO)>>W_BS;
1151 #else
1152         b1[i] += (W_BM*(b0[i] + b2[i])+4*b1[i]+W_BO)>>W_BS;
1153 #endif
1154     }
1155 }
1156
1157 static void vertical_compose97iL1(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
1158     int i;
1159
1160     for(i=0; i<width; i++){
1161         b1[i] -= (W_DM*(b0[i] + b2[i])+W_DO)>>W_DS;
1162     }
1163 }
1164
1165 void ff_snow_vertical_compose97i(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, IDWTELEM *b3, IDWTELEM *b4, IDWTELEM *b5, int width){
1166     int i;
1167
1168     for(i=0; i<width; i++){
1169         b4[i] -= (W_DM*(b3[i] + b5[i])+W_DO)>>W_DS;
1170         b3[i] -= (W_CM*(b2[i] + b4[i])+W_CO)>>W_CS;
1171 #ifdef liftS
1172         b2[i] += (W_BM*(b1[i] + b3[i])+W_BO)>>W_BS;
1173 #else
1174         b2[i] += (W_BM*(b1[i] + b3[i])+4*b2[i]+W_BO)>>W_BS;
1175 #endif
1176         b1[i] += (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
1177     }
1178 }
1179
1180 static void spatial_compose97i_buffered_init(DWTCompose *cs, slice_buffer * sb, int height, int stride_line){
1181     cs->b0 = slice_buffer_get_line(sb, mirror(-3-1, height-1) * stride_line);
1182     cs->b1 = slice_buffer_get_line(sb, mirror(-3  , height-1) * stride_line);
1183     cs->b2 = slice_buffer_get_line(sb, mirror(-3+1, height-1) * stride_line);
1184     cs->b3 = slice_buffer_get_line(sb, mirror(-3+2, height-1) * stride_line);
1185     cs->y = -3;
1186 }
1187
1188 static void spatial_compose97i_init(DWTCompose *cs, IDWTELEM *buffer, int height, int stride){
1189     cs->b0 = buffer + mirror(-3-1, height-1)*stride;
1190     cs->b1 = buffer + mirror(-3  , height-1)*stride;
1191     cs->b2 = buffer + mirror(-3+1, height-1)*stride;
1192     cs->b3 = buffer + mirror(-3+2, height-1)*stride;
1193     cs->y = -3;
1194 }
1195
1196 static void spatial_compose97i_dy_buffered(DSPContext *dsp, DWTCompose *cs, slice_buffer * sb, int width, int height, int stride_line){
1197     int y = cs->y;
1198
1199     IDWTELEM *b0= cs->b0;
1200     IDWTELEM *b1= cs->b1;
1201     IDWTELEM *b2= cs->b2;
1202     IDWTELEM *b3= cs->b3;
1203     IDWTELEM *b4= slice_buffer_get_line(sb, mirror(y + 3, height - 1) * stride_line);
1204     IDWTELEM *b5= slice_buffer_get_line(sb, mirror(y + 4, height - 1) * stride_line);
1205
1206     if(y>0 && y+4<height){
1207         dsp->vertical_compose97i(b0, b1, b2, b3, b4, b5, width);
1208     }else{
1209         if(y+3<(unsigned)height) vertical_compose97iL1(b3, b4, b5, width);
1210         if(y+2<(unsigned)height) vertical_compose97iH1(b2, b3, b4, width);
1211         if(y+1<(unsigned)height) vertical_compose97iL0(b1, b2, b3, width);
1212         if(y+0<(unsigned)height) vertical_compose97iH0(b0, b1, b2, width);
1213     }
1214
1215     if(y-1<(unsigned)height) dsp->horizontal_compose97i(b0, width);
1216     if(y+0<(unsigned)height) dsp->horizontal_compose97i(b1, width);
1217
1218     cs->b0=b2;
1219     cs->b1=b3;
1220     cs->b2=b4;
1221     cs->b3=b5;
1222     cs->y += 2;
1223 }
1224
1225 static void spatial_compose97i_dy(DWTCompose *cs, IDWTELEM *buffer, int width, int height, int stride){
1226     int y = cs->y;
1227     IDWTELEM *b0= cs->b0;
1228     IDWTELEM *b1= cs->b1;
1229     IDWTELEM *b2= cs->b2;
1230     IDWTELEM *b3= cs->b3;
1231     IDWTELEM *b4= buffer + mirror(y+3, height-1)*stride;
1232     IDWTELEM *b5= buffer + mirror(y+4, height-1)*stride;
1233
1234     if(y+3<(unsigned)height) vertical_compose97iL1(b3, b4, b5, width);
1235     if(y+2<(unsigned)height) vertical_compose97iH1(b2, b3, b4, width);
1236     if(y+1<(unsigned)height) vertical_compose97iL0(b1, b2, b3, width);
1237     if(y+0<(unsigned)height) vertical_compose97iH0(b0, b1, b2, width);
1238
1239     if(y-1<(unsigned)height) ff_snow_horizontal_compose97i(b0, width);
1240     if(y+0<(unsigned)height) ff_snow_horizontal_compose97i(b1, width);
1241
1242     cs->b0=b2;
1243     cs->b1=b3;
1244     cs->b2=b4;
1245     cs->b3=b5;
1246     cs->y += 2;
1247 }
1248
1249 static void av_unused spatial_compose97i(IDWTELEM *buffer, int width, int height, int stride){
1250     DWTCompose cs;
1251     spatial_compose97i_init(&cs, buffer, height, stride);
1252     while(cs.y <= height)
1253         spatial_compose97i_dy(&cs, buffer, width, height, stride);
1254 }
1255
1256 static void ff_spatial_idwt_buffered_init(DWTCompose *cs, slice_buffer * sb, int width, int height, int stride_line, int type, int decomposition_count){
1257     int level;
1258     for(level=decomposition_count-1; level>=0; level--){
1259         switch(type){
1260         case DWT_97: spatial_compose97i_buffered_init(cs+level, sb, height>>level, stride_line<<level); break;
1261         case DWT_53: spatial_compose53i_buffered_init(cs+level, sb, height>>level, stride_line<<level); break;
1262         }
1263     }
1264 }
1265
1266 static void ff_spatial_idwt_buffered_slice(DSPContext *dsp, DWTCompose *cs, slice_buffer * slice_buf, int width, int height, int stride_line, int type, int decomposition_count, int y){
1267     const int support = type==1 ? 3 : 5;
1268     int level;
1269     if(type==2) return;
1270
1271     for(level=decomposition_count-1; level>=0; level--){
1272         while(cs[level].y <= FFMIN((y>>level)+support, height>>level)){
1273             switch(type){
1274             case DWT_97: spatial_compose97i_dy_buffered(dsp, cs+level, slice_buf, width>>level, height>>level, stride_line<<level);
1275                 break;
1276             case DWT_53: spatial_compose53i_dy_buffered(cs+level, slice_buf, width>>level, height>>level, stride_line<<level);
1277                 break;
1278             }
1279         }
1280     }
1281 }
1282
1283 static inline void unpack_coeffs(SnowContext *s, SubBand *b, SubBand * parent, int orientation){
1284     const int w= b->width;
1285     const int h= b->height;
1286     int x,y;
1287
1288     int run, runs;
1289     x_and_coeff *xc= b->x_coeff;
1290     x_and_coeff *prev_xc= NULL;
1291     x_and_coeff *prev2_xc= xc;
1292     x_and_coeff *parent_xc= parent ? parent->x_coeff : NULL;
1293     x_and_coeff *prev_parent_xc= parent_xc;
1294
1295     runs= get_symbol2(&s->c, b->state[30], 0);
1296     if(runs-- > 0) run= get_symbol2(&s->c, b->state[1], 3);
1297     else           run= INT_MAX;
1298
1299     for(y=0; y<h; y++){
1300         int v=0;
1301         int lt=0, t=0, rt=0;
1302
1303         if(y && prev_xc->x == 0){
1304             rt= prev_xc->coeff;
1305         }
1306         for(x=0; x<w; x++){
1307             int p=0;
1308             const int l= v;
1309
1310             lt= t; t= rt;
1311
1312             if(y){
1313                 if(prev_xc->x <= x)
1314                     prev_xc++;
1315                 if(prev_xc->x == x + 1)
1316                     rt= prev_xc->coeff;
1317                 else
1318                     rt=0;
1319             }
1320             if(parent_xc){
1321                 if(x>>1 > parent_xc->x){
1322                     parent_xc++;
1323                 }
1324                 if(x>>1 == parent_xc->x){
1325                     p= parent_xc->coeff;
1326                 }
1327             }
1328             if(/*ll|*/l|lt|t|rt|p){
1329                 int context= av_log2(/*FFABS(ll) + */3*(l>>1) + (lt>>1) + (t&~1) + (rt>>1) + (p>>1));
1330
1331                 v=get_rac(&s->c, &b->state[0][context]);
1332                 if(v){
1333                     v= 2*(get_symbol2(&s->c, b->state[context + 2], context-4) + 1);
1334                     v+=get_rac(&s->c, &b->state[0][16 + 1 + 3 + quant3bA[l&0xFF] + 3*quant3bA[t&0xFF]]);
1335
1336                     xc->x=x;
1337                     (xc++)->coeff= v;
1338                 }
1339             }else{
1340                 if(!run){
1341                     if(runs-- > 0) run= get_symbol2(&s->c, b->state[1], 3);
1342                     else           run= INT_MAX;
1343                     v= 2*(get_symbol2(&s->c, b->state[0 + 2], 0-4) + 1);
1344                     v+=get_rac(&s->c, &b->state[0][16 + 1 + 3]);
1345
1346                     xc->x=x;
1347                     (xc++)->coeff= v;
1348                 }else{
1349                     int max_run;
1350                     run--;
1351                     v=0;
1352
1353                     if(y) max_run= FFMIN(run, prev_xc->x - x - 2);
1354                     else  max_run= FFMIN(run, w-x-1);
1355                     if(parent_xc)
1356                         max_run= FFMIN(max_run, 2*parent_xc->x - x - 1);
1357                     x+= max_run;
1358                     run-= max_run;
1359                 }
1360             }
1361         }
1362         (xc++)->x= w+1; //end marker
1363         prev_xc= prev2_xc;
1364         prev2_xc= xc;
1365
1366         if(parent_xc){
1367             if(y&1){
1368                 while(parent_xc->x != parent->width+1)
1369                     parent_xc++;
1370                 parent_xc++;
1371                 prev_parent_xc= parent_xc;
1372             }else{
1373                 parent_xc= prev_parent_xc;
1374             }
1375         }
1376     }
1377
1378     (xc++)->x= w+1; //end marker
1379 }
1380
1381 static inline void decode_subband_slice_buffered(SnowContext *s, SubBand *b, slice_buffer * sb, int start_y, int h, int save_state[1]){
1382     const int w= b->width;
1383     int y;
1384     const int qlog= av_clip(s->qlog + b->qlog, 0, QROOT*16);
1385     int qmul= qexp[qlog&(QROOT-1)]<<(qlog>>QSHIFT);
1386     int qadd= (s->qbias*qmul)>>QBIAS_SHIFT;
1387     int new_index = 0;
1388
1389     if(b->ibuf == s->spatial_idwt_buffer || s->qlog == LOSSLESS_QLOG){
1390         qadd= 0;
1391         qmul= 1<<QEXPSHIFT;
1392     }
1393
1394     /* If we are on the second or later slice, restore our index. */
1395     if (start_y != 0)
1396         new_index = save_state[0];
1397
1398
1399     for(y=start_y; y<h; y++){
1400         int x = 0;
1401         int v;
1402         IDWTELEM * line = slice_buffer_get_line(sb, y * b->stride_line + b->buf_y_offset) + b->buf_x_offset;
1403         memset(line, 0, b->width*sizeof(IDWTELEM));
1404         v = b->x_coeff[new_index].coeff;
1405         x = b->x_coeff[new_index++].x;
1406         while(x < w){
1407             register int t= ( (v>>1)*qmul + qadd)>>QEXPSHIFT;
1408             register int u= -(v&1);
1409             line[x] = (t^u) - u;
1410
1411             v = b->x_coeff[new_index].coeff;
1412             x = b->x_coeff[new_index++].x;
1413         }
1414     }
1415
1416     /* Save our variables for the next slice. */
1417     save_state[0] = new_index;
1418
1419     return;
1420 }
1421
1422 static void reset_contexts(SnowContext *s){ //FIXME better initial contexts
1423     int plane_index, level, orientation;
1424
1425     for(plane_index=0; plane_index<3; plane_index++){
1426         for(level=0; level<MAX_DECOMPOSITIONS; level++){
1427             for(orientation=level ? 1:0; orientation<4; orientation++){
1428                 memset(s->plane[plane_index].band[level][orientation].state, MID_STATE, sizeof(s->plane[plane_index].band[level][orientation].state));
1429             }
1430         }
1431     }
1432     memset(s->header_state, MID_STATE, sizeof(s->header_state));
1433     memset(s->block_state, MID_STATE, sizeof(s->block_state));
1434 }
1435
1436 static int alloc_blocks(SnowContext *s){
1437     int w= -((-s->avctx->width )>>LOG2_MB_SIZE);
1438     int h= -((-s->avctx->height)>>LOG2_MB_SIZE);
1439
1440     s->b_width = w;
1441     s->b_height= h;
1442
1443     av_free(s->block);
1444     s->block= av_mallocz(w * h * sizeof(BlockNode) << (s->block_max_depth*2));
1445     return 0;
1446 }
1447
1448 static inline void copy_rac_state(RangeCoder *d, RangeCoder *s){
1449     uint8_t *bytestream= d->bytestream;
1450     uint8_t *bytestream_start= d->bytestream_start;
1451     *d= *s;
1452     d->bytestream= bytestream;
1453     d->bytestream_start= bytestream_start;
1454 }
1455
1456 static inline void set_blocks(SnowContext *s, int level, int x, int y, int l, int cb, int cr, int mx, int my, int ref, int type){
1457     const int w= s->b_width << s->block_max_depth;
1458     const int rem_depth= s->block_max_depth - level;
1459     const int index= (x + y*w) << rem_depth;
1460     const int block_w= 1<<rem_depth;
1461     BlockNode block;
1462     int i,j;
1463
1464     block.color[0]= l;
1465     block.color[1]= cb;
1466     block.color[2]= cr;
1467     block.mx= mx;
1468     block.my= my;
1469     block.ref= ref;
1470     block.type= type;
1471     block.level= level;
1472
1473     for(j=0; j<block_w; j++){
1474         for(i=0; i<block_w; i++){
1475             s->block[index + i + j*w]= block;
1476         }
1477     }
1478 }
1479
1480 static inline void init_ref(MotionEstContext *c, uint8_t *src[3], uint8_t *ref[3], uint8_t *ref2[3], int x, int y, int ref_index){
1481     const int offset[3]= {
1482           y*c->  stride + x,
1483         ((y*c->uvstride + x)>>1),
1484         ((y*c->uvstride + x)>>1),
1485     };
1486     int i;
1487     for(i=0; i<3; i++){
1488         c->src[0][i]= src [i];
1489         c->ref[0][i]= ref [i] + offset[i];
1490     }
1491     assert(!ref_index);
1492 }
1493
1494 static inline void pred_mv(SnowContext *s, int *mx, int *my, int ref,
1495                            const BlockNode *left, const BlockNode *top, const BlockNode *tr){
1496     if(s->ref_frames == 1){
1497         *mx = mid_pred(left->mx, top->mx, tr->mx);
1498         *my = mid_pred(left->my, top->my, tr->my);
1499     }else{
1500         const int *scale = scale_mv_ref[ref];
1501         *mx = mid_pred((left->mx * scale[left->ref] + 128) >>8,
1502                        (top ->mx * scale[top ->ref] + 128) >>8,
1503                        (tr  ->mx * scale[tr  ->ref] + 128) >>8);
1504         *my = mid_pred((left->my * scale[left->ref] + 128) >>8,
1505                        (top ->my * scale[top ->ref] + 128) >>8,
1506                        (tr  ->my * scale[tr  ->ref] + 128) >>8);
1507     }
1508 }
1509
1510 static av_always_inline int same_block(BlockNode *a, BlockNode *b){
1511     if((a->type&BLOCK_INTRA) && (b->type&BLOCK_INTRA)){
1512         return !((a->color[0] - b->color[0]) | (a->color[1] - b->color[1]) | (a->color[2] - b->color[2]));
1513     }else{
1514         return !((a->mx - b->mx) | (a->my - b->my) | (a->ref - b->ref) | ((a->type ^ b->type)&BLOCK_INTRA));
1515     }
1516 }
1517
1518 static void decode_q_branch(SnowContext *s, int level, int x, int y){
1519     const int w= s->b_width << s->block_max_depth;
1520     const int rem_depth= s->block_max_depth - level;
1521     const int index= (x + y*w) << rem_depth;
1522     int trx= (x+1)<<rem_depth;
1523     const BlockNode *left  = x ? &s->block[index-1] : &null_block;
1524     const BlockNode *top   = y ? &s->block[index-w] : &null_block;
1525     const BlockNode *tl    = y && x ? &s->block[index-w-1] : left;
1526     const BlockNode *tr    = y && trx<w && ((x&1)==0 || level==0) ? &s->block[index-w+(1<<rem_depth)] : tl; //FIXME use lt
1527     int s_context= 2*left->level + 2*top->level + tl->level + tr->level;
1528
1529     if(s->keyframe){
1530         set_blocks(s, level, x, y, null_block.color[0], null_block.color[1], null_block.color[2], null_block.mx, null_block.my, null_block.ref, BLOCK_INTRA);
1531         return;
1532     }
1533
1534     if(level==s->block_max_depth || get_rac(&s->c, &s->block_state[4 + s_context])){
1535         int type, mx, my;
1536         int l = left->color[0];
1537         int cb= left->color[1];
1538         int cr= left->color[2];
1539         int ref = 0;
1540         int ref_context= av_log2(2*left->ref) + av_log2(2*top->ref);
1541         int mx_context= av_log2(2*FFABS(left->mx - top->mx)) + 0*av_log2(2*FFABS(tr->mx - top->mx));
1542         int my_context= av_log2(2*FFABS(left->my - top->my)) + 0*av_log2(2*FFABS(tr->my - top->my));
1543
1544         type= get_rac(&s->c, &s->block_state[1 + left->type + top->type]) ? BLOCK_INTRA : 0;
1545
1546         if(type){
1547             pred_mv(s, &mx, &my, 0, left, top, tr);
1548             l += get_symbol(&s->c, &s->block_state[32], 1);
1549             cb+= get_symbol(&s->c, &s->block_state[64], 1);
1550             cr+= get_symbol(&s->c, &s->block_state[96], 1);
1551         }else{
1552             if(s->ref_frames > 1)
1553                 ref= get_symbol(&s->c, &s->block_state[128 + 1024 + 32*ref_context], 0);
1554             pred_mv(s, &mx, &my, ref, left, top, tr);
1555             mx+= get_symbol(&s->c, &s->block_state[128 + 32*(mx_context + 16*!!ref)], 1);
1556             my+= get_symbol(&s->c, &s->block_state[128 + 32*(my_context + 16*!!ref)], 1);
1557         }
1558         set_blocks(s, level, x, y, l, cb, cr, mx, my, ref, type);
1559     }else{
1560         decode_q_branch(s, level+1, 2*x+0, 2*y+0);
1561         decode_q_branch(s, level+1, 2*x+1, 2*y+0);
1562         decode_q_branch(s, level+1, 2*x+0, 2*y+1);
1563         decode_q_branch(s, level+1, 2*x+1, 2*y+1);
1564     }
1565 }
1566
1567 static void decode_blocks(SnowContext *s){
1568     int x, y;
1569     int w= s->b_width;
1570     int h= s->b_height;
1571
1572     for(y=0; y<h; y++){
1573         for(x=0; x<w; x++){
1574             decode_q_branch(s, 0, x, y);
1575         }
1576     }
1577 }
1578
1579 static void mc_block(Plane *p, uint8_t *dst, const uint8_t *src, uint8_t *tmp, int stride, int b_w, int b_h, int dx, int dy){
1580     static const uint8_t weight[64]={
1581     8,7,6,5,4,3,2,1,
1582     7,7,0,0,0,0,0,1,
1583     6,0,6,0,0,0,2,0,
1584     5,0,0,5,0,3,0,0,
1585     4,0,0,0,4,0,0,0,
1586     3,0,0,5,0,3,0,0,
1587     2,0,6,0,0,0,2,0,
1588     1,7,0,0,0,0,0,1,
1589     };
1590
1591     static const uint8_t brane[256]={
1592     0x00,0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x11,0x12,0x12,0x12,0x12,0x12,0x12,0x12,
1593     0x04,0x05,0xcc,0xcc,0xcc,0xcc,0xcc,0x41,0x15,0x16,0xcc,0xcc,0xcc,0xcc,0xcc,0x52,
1594     0x04,0xcc,0x05,0xcc,0xcc,0xcc,0x41,0xcc,0x15,0xcc,0x16,0xcc,0xcc,0xcc,0x52,0xcc,
1595     0x04,0xcc,0xcc,0x05,0xcc,0x41,0xcc,0xcc,0x15,0xcc,0xcc,0x16,0xcc,0x52,0xcc,0xcc,
1596     0x04,0xcc,0xcc,0xcc,0x41,0xcc,0xcc,0xcc,0x15,0xcc,0xcc,0xcc,0x16,0xcc,0xcc,0xcc,
1597     0x04,0xcc,0xcc,0x41,0xcc,0x05,0xcc,0xcc,0x15,0xcc,0xcc,0x52,0xcc,0x16,0xcc,0xcc,
1598     0x04,0xcc,0x41,0xcc,0xcc,0xcc,0x05,0xcc,0x15,0xcc,0x52,0xcc,0xcc,0xcc,0x16,0xcc,
1599     0x04,0x41,0xcc,0xcc,0xcc,0xcc,0xcc,0x05,0x15,0x52,0xcc,0xcc,0xcc,0xcc,0xcc,0x16,
1600     0x44,0x45,0x45,0x45,0x45,0x45,0x45,0x45,0x55,0x56,0x56,0x56,0x56,0x56,0x56,0x56,
1601     0x48,0x49,0xcc,0xcc,0xcc,0xcc,0xcc,0x85,0x59,0x5A,0xcc,0xcc,0xcc,0xcc,0xcc,0x96,
1602     0x48,0xcc,0x49,0xcc,0xcc,0xcc,0x85,0xcc,0x59,0xcc,0x5A,0xcc,0xcc,0xcc,0x96,0xcc,
1603     0x48,0xcc,0xcc,0x49,0xcc,0x85,0xcc,0xcc,0x59,0xcc,0xcc,0x5A,0xcc,0x96,0xcc,0xcc,
1604     0x48,0xcc,0xcc,0xcc,0x49,0xcc,0xcc,0xcc,0x59,0xcc,0xcc,0xcc,0x96,0xcc,0xcc,0xcc,
1605     0x48,0xcc,0xcc,0x85,0xcc,0x49,0xcc,0xcc,0x59,0xcc,0xcc,0x96,0xcc,0x5A,0xcc,0xcc,
1606     0x48,0xcc,0x85,0xcc,0xcc,0xcc,0x49,0xcc,0x59,0xcc,0x96,0xcc,0xcc,0xcc,0x5A,0xcc,
1607     0x48,0x85,0xcc,0xcc,0xcc,0xcc,0xcc,0x49,0x59,0x96,0xcc,0xcc,0xcc,0xcc,0xcc,0x5A,
1608     };
1609
1610     static const uint8_t needs[16]={
1611     0,1,0,0,
1612     2,4,2,0,
1613     0,1,0,0,
1614     15
1615     };
1616
1617     int x, y, b, r, l;
1618     int16_t tmpIt   [64*(32+HTAPS_MAX)];
1619     uint8_t tmp2t[3][stride*(32+HTAPS_MAX)];
1620     int16_t *tmpI= tmpIt;
1621     uint8_t *tmp2= tmp2t[0];
1622     const uint8_t *hpel[11];
1623     assert(dx<16 && dy<16);
1624     r= brane[dx + 16*dy]&15;
1625     l= brane[dx + 16*dy]>>4;
1626
1627     b= needs[l] | needs[r];
1628     if(p && !p->diag_mc)
1629         b= 15;
1630
1631     if(b&5){
1632         for(y=0; y < b_h+HTAPS_MAX-1; y++){
1633             for(x=0; x < b_w; x++){
1634                 int a_1=src[x + HTAPS_MAX/2-4];
1635                 int a0= src[x + HTAPS_MAX/2-3];
1636                 int a1= src[x + HTAPS_MAX/2-2];
1637                 int a2= src[x + HTAPS_MAX/2-1];
1638                 int a3= src[x + HTAPS_MAX/2+0];
1639                 int a4= src[x + HTAPS_MAX/2+1];
1640                 int a5= src[x + HTAPS_MAX/2+2];
1641                 int a6= src[x + HTAPS_MAX/2+3];
1642                 int am=0;
1643                 if(!p || p->fast_mc){
1644                     am= 20*(a2+a3) - 5*(a1+a4) + (a0+a5);
1645                     tmpI[x]= am;
1646                     am= (am+16)>>5;
1647                 }else{
1648                     am= p->hcoeff[0]*(a2+a3) + p->hcoeff[1]*(a1+a4) + p->hcoeff[2]*(a0+a5) + p->hcoeff[3]*(a_1+a6);
1649                     tmpI[x]= am;
1650                     am= (am+32)>>6;
1651                 }
1652
1653                 if(am&(~255)) am= ~(am>>31);
1654                 tmp2[x]= am;
1655             }
1656             tmpI+= 64;
1657             tmp2+= stride;
1658             src += stride;
1659         }
1660         src -= stride*y;
1661     }
1662     src += HTAPS_MAX/2 - 1;
1663     tmp2= tmp2t[1];
1664
1665     if(b&2){
1666         for(y=0; y < b_h; y++){
1667             for(x=0; x < b_w+1; x++){
1668                 int a_1=src[x + (HTAPS_MAX/2-4)*stride];
1669                 int a0= src[x + (HTAPS_MAX/2-3)*stride];
1670                 int a1= src[x + (HTAPS_MAX/2-2)*stride];
1671                 int a2= src[x + (HTAPS_MAX/2-1)*stride];
1672                 int a3= src[x + (HTAPS_MAX/2+0)*stride];
1673                 int a4= src[x + (HTAPS_MAX/2+1)*stride];
1674                 int a5= src[x + (HTAPS_MAX/2+2)*stride];
1675                 int a6= src[x + (HTAPS_MAX/2+3)*stride];
1676                 int am=0;
1677                 if(!p || p->fast_mc)
1678                     am= (20*(a2+a3) - 5*(a1+a4) + (a0+a5) + 16)>>5;
1679                 else
1680                     am= (p->hcoeff[0]*(a2+a3) + p->hcoeff[1]*(a1+a4) + p->hcoeff[2]*(a0+a5) + p->hcoeff[3]*(a_1+a6) + 32)>>6;
1681
1682                 if(am&(~255)) am= ~(am>>31);
1683                 tmp2[x]= am;
1684             }
1685             src += stride;
1686             tmp2+= stride;
1687         }
1688         src -= stride*y;
1689     }
1690     src += stride*(HTAPS_MAX/2 - 1);
1691     tmp2= tmp2t[2];
1692     tmpI= tmpIt;
1693     if(b&4){
1694         for(y=0; y < b_h; y++){
1695             for(x=0; x < b_w; x++){
1696                 int a_1=tmpI[x + (HTAPS_MAX/2-4)*64];
1697                 int a0= tmpI[x + (HTAPS_MAX/2-3)*64];
1698                 int a1= tmpI[x + (HTAPS_MAX/2-2)*64];
1699                 int a2= tmpI[x + (HTAPS_MAX/2-1)*64];
1700                 int a3= tmpI[x + (HTAPS_MAX/2+0)*64];
1701                 int a4= tmpI[x + (HTAPS_MAX/2+1)*64];
1702                 int a5= tmpI[x + (HTAPS_MAX/2+2)*64];
1703                 int a6= tmpI[x + (HTAPS_MAX/2+3)*64];
1704                 int am=0;
1705                 if(!p || p->fast_mc)
1706                     am= (20*(a2+a3) - 5*(a1+a4) + (a0+a5) + 512)>>10;
1707                 else
1708                     am= (p->hcoeff[0]*(a2+a3) + p->hcoeff[1]*(a1+a4) + p->hcoeff[2]*(a0+a5) + p->hcoeff[3]*(a_1+a6) + 2048)>>12;
1709                 if(am&(~255)) am= ~(am>>31);
1710                 tmp2[x]= am;
1711             }
1712             tmpI+= 64;
1713             tmp2+= stride;
1714         }
1715     }
1716
1717     hpel[ 0]= src;
1718     hpel[ 1]= tmp2t[0] + stride*(HTAPS_MAX/2-1);
1719     hpel[ 2]= src + 1;
1720
1721     hpel[ 4]= tmp2t[1];
1722     hpel[ 5]= tmp2t[2];
1723     hpel[ 6]= tmp2t[1] + 1;
1724
1725     hpel[ 8]= src + stride;
1726     hpel[ 9]= hpel[1] + stride;
1727     hpel[10]= hpel[8] + 1;
1728
1729     if(b==15){
1730         const uint8_t *src1= hpel[dx/8 + dy/8*4  ];
1731         const uint8_t *src2= hpel[dx/8 + dy/8*4+1];
1732         const uint8_t *src3= hpel[dx/8 + dy/8*4+4];
1733         const uint8_t *src4= hpel[dx/8 + dy/8*4+5];
1734         dx&=7;
1735         dy&=7;
1736         for(y=0; y < b_h; y++){
1737             for(x=0; x < b_w; x++){
1738                 dst[x]= ((8-dx)*(8-dy)*src1[x] + dx*(8-dy)*src2[x]+
1739                          (8-dx)*   dy *src3[x] + dx*   dy *src4[x]+32)>>6;
1740             }
1741             src1+=stride;
1742             src2+=stride;
1743             src3+=stride;
1744             src4+=stride;
1745             dst +=stride;
1746         }
1747     }else{
1748         const uint8_t *src1= hpel[l];
1749         const uint8_t *src2= hpel[r];
1750         int a= weight[((dx&7) + (8*(dy&7)))];
1751         int b= 8-a;
1752         for(y=0; y < b_h; y++){
1753             for(x=0; x < b_w; x++){
1754                 dst[x]= (a*src1[x] + b*src2[x] + 4)>>3;
1755             }
1756             src1+=stride;
1757             src2+=stride;
1758             dst +=stride;
1759         }
1760     }
1761 }
1762
1763 #define mca(dx,dy,b_w)\
1764 static void mc_block_hpel ## dx ## dy ## b_w(uint8_t *dst, const uint8_t *src, int stride, int h){\
1765     uint8_t tmp[stride*(b_w+HTAPS_MAX-1)];\
1766     assert(h==b_w);\
1767     mc_block(NULL, dst, src-(HTAPS_MAX/2-1)-(HTAPS_MAX/2-1)*stride, tmp, stride, b_w, b_w, dx, dy);\
1768 }
1769
1770 mca( 0, 0,16)
1771 mca( 8, 0,16)
1772 mca( 0, 8,16)
1773 mca( 8, 8,16)
1774 mca( 0, 0,8)
1775 mca( 8, 0,8)
1776 mca( 0, 8,8)
1777 mca( 8, 8,8)
1778
1779 static void pred_block(SnowContext *s, uint8_t *dst, uint8_t *tmp, int stride, int sx, int sy, int b_w, int b_h, BlockNode *block, int plane_index, int w, int h){
1780     if(block->type & BLOCK_INTRA){
1781         int x, y;
1782         const int color = block->color[plane_index];
1783         const int color4= color*0x01010101;
1784         if(b_w==32){
1785             for(y=0; y < b_h; y++){
1786                 *(uint32_t*)&dst[0 + y*stride]= color4;
1787                 *(uint32_t*)&dst[4 + y*stride]= color4;
1788                 *(uint32_t*)&dst[8 + y*stride]= color4;
1789                 *(uint32_t*)&dst[12+ y*stride]= color4;
1790                 *(uint32_t*)&dst[16+ y*stride]= color4;
1791                 *(uint32_t*)&dst[20+ y*stride]= color4;
1792                 *(uint32_t*)&dst[24+ y*stride]= color4;
1793                 *(uint32_t*)&dst[28+ y*stride]= color4;
1794             }
1795         }else if(b_w==16){
1796             for(y=0; y < b_h; y++){
1797                 *(uint32_t*)&dst[0 + y*stride]= color4;
1798                 *(uint32_t*)&dst[4 + y*stride]= color4;
1799                 *(uint32_t*)&dst[8 + y*stride]= color4;
1800                 *(uint32_t*)&dst[12+ y*stride]= color4;
1801             }
1802         }else if(b_w==8){
1803             for(y=0; y < b_h; y++){
1804                 *(uint32_t*)&dst[0 + y*stride]= color4;
1805                 *(uint32_t*)&dst[4 + y*stride]= color4;
1806             }
1807         }else if(b_w==4){
1808             for(y=0; y < b_h; y++){
1809                 *(uint32_t*)&dst[0 + y*stride]= color4;
1810             }
1811         }else{
1812             for(y=0; y < b_h; y++){
1813                 for(x=0; x < b_w; x++){
1814                     dst[x + y*stride]= color;
1815                 }
1816             }
1817         }
1818     }else{
1819         uint8_t *src= s->last_picture[block->ref].data[plane_index];
1820         const int scale= plane_index ?  s->mv_scale : 2*s->mv_scale;
1821         int mx= block->mx*scale;
1822         int my= block->my*scale;
1823         const int dx= mx&15;
1824         const int dy= my&15;
1825         const int tab_index= 3 - (b_w>>2) + (b_w>>4);
1826         sx += (mx>>4) - (HTAPS_MAX/2-1);
1827         sy += (my>>4) - (HTAPS_MAX/2-1);
1828         src += sx + sy*stride;
1829         if(   (unsigned)sx >= w - b_w - (HTAPS_MAX-2)
1830            || (unsigned)sy >= h - b_h - (HTAPS_MAX-2)){
1831             ff_emulated_edge_mc(tmp + MB_SIZE, src, stride, b_w+HTAPS_MAX-1, b_h+HTAPS_MAX-1, sx, sy, w, h);
1832             src= tmp + MB_SIZE;
1833         }
1834 //        assert(b_w == b_h || 2*b_w == b_h || b_w == 2*b_h);
1835 //        assert(!(b_w&(b_w-1)));
1836         assert(b_w>1 && b_h>1);
1837         assert((tab_index>=0 && tab_index<4) || b_w==32);
1838         if((dx&3) || (dy&3) || !(b_w == b_h || 2*b_w == b_h || b_w == 2*b_h) || (b_w&(b_w-1)) || !s->plane[plane_index].fast_mc )
1839             mc_block(&s->plane[plane_index], dst, src, tmp, stride, b_w, b_h, dx, dy);
1840         else if(b_w==32){
1841             int y;
1842             for(y=0; y<b_h; y+=16){
1843                 s->dsp.put_h264_qpel_pixels_tab[0][dy+(dx>>2)](dst + y*stride, src + 3 + (y+3)*stride,stride);
1844                 s->dsp.put_h264_qpel_pixels_tab[0][dy+(dx>>2)](dst + 16 + y*stride, src + 19 + (y+3)*stride,stride);
1845             }
1846         }else if(b_w==b_h)
1847             s->dsp.put_h264_qpel_pixels_tab[tab_index  ][dy+(dx>>2)](dst,src + 3 + 3*stride,stride);
1848         else if(b_w==2*b_h){
1849             s->dsp.put_h264_qpel_pixels_tab[tab_index+1][dy+(dx>>2)](dst    ,src + 3       + 3*stride,stride);
1850             s->dsp.put_h264_qpel_pixels_tab[tab_index+1][dy+(dx>>2)](dst+b_h,src + 3 + b_h + 3*stride,stride);
1851         }else{
1852             assert(2*b_w==b_h);
1853             s->dsp.put_h264_qpel_pixels_tab[tab_index  ][dy+(dx>>2)](dst           ,src + 3 + 3*stride           ,stride);
1854             s->dsp.put_h264_qpel_pixels_tab[tab_index  ][dy+(dx>>2)](dst+b_w*stride,src + 3 + 3*stride+b_w*stride,stride);
1855         }
1856     }
1857 }
1858
1859 void ff_snow_inner_add_yblock(const uint8_t *obmc, const int obmc_stride, uint8_t * * block, int b_w, int b_h,
1860                               int src_x, int src_y, int src_stride, slice_buffer * sb, int add, uint8_t * dst8){
1861     int y, x;
1862     IDWTELEM * dst;
1863     for(y=0; y<b_h; y++){
1864         //FIXME ugly misuse of obmc_stride
1865         const uint8_t *obmc1= obmc + y*obmc_stride;
1866         const uint8_t *obmc2= obmc1+ (obmc_stride>>1);
1867         const uint8_t *obmc3= obmc1+ obmc_stride*(obmc_stride>>1);
1868         const uint8_t *obmc4= obmc3+ (obmc_stride>>1);
1869         dst = slice_buffer_get_line(sb, src_y + y);
1870         for(x=0; x<b_w; x++){
1871             int v=   obmc1[x] * block[3][x + y*src_stride]
1872                     +obmc2[x] * block[2][x + y*src_stride]
1873                     +obmc3[x] * block[1][x + y*src_stride]
1874                     +obmc4[x] * block[0][x + y*src_stride];
1875
1876             v <<= 8 - LOG2_OBMC_MAX;
1877             if(FRAC_BITS != 8){
1878                 v >>= 8 - FRAC_BITS;
1879             }
1880             if(add){
1881                 v += dst[x + src_x];
1882                 v = (v + (1<<(FRAC_BITS-1))) >> FRAC_BITS;
1883                 if(v&(~255)) v= ~(v>>31);
1884                 dst8[x + y*src_stride] = v;
1885             }else{
1886                 dst[x + src_x] -= v;
1887             }
1888         }
1889     }
1890 }
1891
1892 //FIXME name cleanup (b_w, block_w, b_width stuff)
1893 static av_always_inline void add_yblock(SnowContext *s, int sliced, slice_buffer *sb, IDWTELEM *dst, uint8_t *dst8, const uint8_t *obmc, int src_x, int src_y, int b_w, int b_h, int w, int h, int dst_stride, int src_stride, int obmc_stride, int b_x, int b_y, int add, int offset_dst, int plane_index){
1894     const int b_width = s->b_width  << s->block_max_depth;
1895     const int b_height= s->b_height << s->block_max_depth;
1896     const int b_stride= b_width;
1897     BlockNode *lt= &s->block[b_x + b_y*b_stride];
1898     BlockNode *rt= lt+1;
1899     BlockNode *lb= lt+b_stride;
1900     BlockNode *rb= lb+1;
1901     uint8_t *block[4];
1902     int tmp_step= src_stride >= 7*MB_SIZE ? MB_SIZE : MB_SIZE*src_stride;
1903     uint8_t *tmp = s->scratchbuf;
1904     uint8_t *ptmp;
1905     int x,y;
1906
1907     if(b_x<0){
1908         lt= rt;
1909         lb= rb;
1910     }else if(b_x + 1 >= b_width){
1911         rt= lt;
1912         rb= lb;
1913     }
1914     if(b_y<0){
1915         lt= lb;
1916         rt= rb;
1917     }else if(b_y + 1 >= b_height){
1918         lb= lt;
1919         rb= rt;
1920     }
1921
1922     if(src_x<0){ //FIXME merge with prev & always round internal width up to *16
1923         obmc -= src_x;
1924         b_w += src_x;
1925         if(!sliced && !offset_dst)
1926             dst -= src_x;
1927         src_x=0;
1928     }else if(src_x + b_w > w){
1929         b_w = w - src_x;
1930     }
1931     if(src_y<0){
1932         obmc -= src_y*obmc_stride;
1933         b_h += src_y;
1934         if(!sliced && !offset_dst)
1935             dst -= src_y*dst_stride;
1936         src_y=0;
1937     }else if(src_y + b_h> h){
1938         b_h = h - src_y;
1939     }
1940
1941     if(b_w<=0 || b_h<=0) return;
1942
1943     assert(src_stride > 2*MB_SIZE + 5);
1944
1945     if(!sliced && offset_dst)
1946         dst += src_x + src_y*dst_stride;
1947     dst8+= src_x + src_y*src_stride;
1948 //    src += src_x + src_y*src_stride;
1949
1950     ptmp= tmp + 3*tmp_step;
1951     block[0]= ptmp;
1952     ptmp+=tmp_step;
1953     pred_block(s, block[0], tmp, src_stride, src_x, src_y, b_w, b_h, lt, plane_index, w, h);
1954
1955     if(same_block(lt, rt)){
1956         block[1]= block[0];
1957     }else{
1958         block[1]= ptmp;
1959         ptmp+=tmp_step;
1960         pred_block(s, block[1], tmp, src_stride, src_x, src_y, b_w, b_h, rt, plane_index, w, h);
1961     }
1962
1963     if(same_block(lt, lb)){
1964         block[2]= block[0];
1965     }else if(same_block(rt, lb)){
1966         block[2]= block[1];
1967     }else{
1968         block[2]= ptmp;
1969         ptmp+=tmp_step;
1970         pred_block(s, block[2], tmp, src_stride, src_x, src_y, b_w, b_h, lb, plane_index, w, h);
1971     }
1972
1973     if(same_block(lt, rb) ){
1974         block[3]= block[0];
1975     }else if(same_block(rt, rb)){
1976         block[3]= block[1];
1977     }else if(same_block(lb, rb)){
1978         block[3]= block[2];
1979     }else{
1980         block[3]= ptmp;
1981         pred_block(s, block[3], tmp, src_stride, src_x, src_y, b_w, b_h, rb, plane_index, w, h);
1982     }
1983 #if 0
1984     for(y=0; y<b_h; y++){
1985         for(x=0; x<b_w; x++){
1986             int v=   obmc [x + y*obmc_stride] * block[3][x + y*src_stride] * (256/OBMC_MAX);
1987             if(add) dst[x + y*dst_stride] += v;
1988             else    dst[x + y*dst_stride] -= v;
1989         }
1990     }
1991     for(y=0; y<b_h; y++){
1992         uint8_t *obmc2= obmc + (obmc_stride>>1);
1993         for(x=0; x<b_w; x++){
1994             int v=   obmc2[x + y*obmc_stride] * block[2][x + y*src_stride] * (256/OBMC_MAX);
1995             if(add) dst[x + y*dst_stride] += v;
1996             else    dst[x + y*dst_stride] -= v;
1997         }
1998     }
1999     for(y=0; y<b_h; y++){
2000         uint8_t *obmc3= obmc + obmc_stride*(obmc_stride>>1);
2001         for(x=0; x<b_w; x++){
2002             int v=   obmc3[x + y*obmc_stride] * block[1][x + y*src_stride] * (256/OBMC_MAX);
2003             if(add) dst[x + y*dst_stride] += v;
2004             else    dst[x + y*dst_stride] -= v;
2005         }
2006     }
2007     for(y=0; y<b_h; y++){
2008         uint8_t *obmc3= obmc + obmc_stride*(obmc_stride>>1);
2009         uint8_t *obmc4= obmc3+ (obmc_stride>>1);
2010         for(x=0; x<b_w; x++){
2011             int v=   obmc4[x + y*obmc_stride] * block[0][x + y*src_stride] * (256/OBMC_MAX);
2012             if(add) dst[x + y*dst_stride] += v;
2013             else    dst[x + y*dst_stride] -= v;
2014         }
2015     }
2016 #else
2017     if(sliced){
2018         s->dsp.inner_add_yblock(obmc, obmc_stride, block, b_w, b_h, src_x,src_y, src_stride, sb, add, dst8);
2019     }else{
2020         for(y=0; y<b_h; y++){
2021             //FIXME ugly misuse of obmc_stride
2022             const uint8_t *obmc1= obmc + y*obmc_stride;
2023             const uint8_t *obmc2= obmc1+ (obmc_stride>>1);
2024             const uint8_t *obmc3= obmc1+ obmc_stride*(obmc_stride>>1);
2025             const uint8_t *obmc4= obmc3+ (obmc_stride>>1);
2026             for(x=0; x<b_w; x++){
2027                 int v=   obmc1[x] * block[3][x + y*src_stride]
2028                         +obmc2[x] * block[2][x + y*src_stride]
2029                         +obmc3[x] * block[1][x + y*src_stride]
2030                         +obmc4[x] * block[0][x + y*src_stride];
2031
2032                 v <<= 8 - LOG2_OBMC_MAX;
2033                 if(FRAC_BITS != 8){
2034                     v >>= 8 - FRAC_BITS;
2035                 }
2036                 if(add){
2037                     v += dst[x + y*dst_stride];
2038                     v = (v + (1<<(FRAC_BITS-1))) >> FRAC_BITS;
2039                     if(v&(~255)) v= ~(v>>31);
2040                     dst8[x + y*src_stride] = v;
2041                 }else{
2042                     dst[x + y*dst_stride] -= v;
2043                 }
2044             }
2045         }
2046     }
2047 #endif /* 0 */
2048 }
2049
2050 static av_always_inline void predict_slice_buffered(SnowContext *s, slice_buffer * sb, IDWTELEM * old_buffer, int plane_index, int add, int mb_y){
2051     Plane *p= &s->plane[plane_index];
2052     const int mb_w= s->b_width  << s->block_max_depth;
2053     const int mb_h= s->b_height << s->block_max_depth;
2054     int x, y, mb_x;
2055     int block_size = MB_SIZE >> s->block_max_depth;
2056     int block_w    = plane_index ? block_size/2 : block_size;
2057     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];
2058     int obmc_stride= plane_index ? block_size : 2*block_size;
2059     int ref_stride= s->current_picture.linesize[plane_index];
2060     uint8_t *dst8= s->current_picture.data[plane_index];
2061     int w= p->width;
2062     int h= p->height;
2063
2064     if(s->keyframe || (s->avctx->debug&512)){
2065         if(mb_y==mb_h)
2066             return;
2067
2068         if(add){
2069             for(y=block_w*mb_y; y<FFMIN(h,block_w*(mb_y+1)); y++){
2070 //                DWTELEM * line = slice_buffer_get_line(sb, y);
2071                 IDWTELEM * line = sb->line[y];
2072                 for(x=0; x<w; x++){
2073 //                    int v= buf[x + y*w] + (128<<FRAC_BITS) + (1<<(FRAC_BITS-1));
2074                     int v= line[x] + (128<<FRAC_BITS) + (1<<(FRAC_BITS-1));
2075                     v >>= FRAC_BITS;
2076                     if(v&(~255)) v= ~(v>>31);
2077                     dst8[x + y*ref_stride]= v;
2078                 }
2079             }
2080         }else{
2081             for(y=block_w*mb_y; y<FFMIN(h,block_w*(mb_y+1)); y++){
2082 //                DWTELEM * line = slice_buffer_get_line(sb, y);
2083                 IDWTELEM * line = sb->line[y];
2084                 for(x=0; x<w; x++){
2085                     line[x] -= 128 << FRAC_BITS;
2086 //                    buf[x + y*w]-= 128<<FRAC_BITS;
2087                 }
2088             }
2089         }
2090
2091         return;
2092     }
2093
2094     for(mb_x=0; mb_x<=mb_w; mb_x++){
2095         add_yblock(s, 1, sb, old_buffer, dst8, obmc,
2096                    block_w*mb_x - block_w/2,
2097                    block_w*mb_y - block_w/2,
2098                    block_w, block_w,
2099                    w, h,
2100                    w, ref_stride, obmc_stride,
2101                    mb_x - 1, mb_y - 1,
2102                    add, 0, plane_index);
2103     }
2104 }
2105
2106 static av_always_inline void predict_slice(SnowContext *s, IDWTELEM *buf, int plane_index, int add, int mb_y){
2107     Plane *p= &s->plane[plane_index];
2108     const int mb_w= s->b_width  << s->block_max_depth;
2109     const int mb_h= s->b_height << s->block_max_depth;
2110     int x, y, mb_x;
2111     int block_size = MB_SIZE >> s->block_max_depth;
2112     int block_w    = plane_index ? block_size/2 : block_size;
2113     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];
2114     const int obmc_stride= plane_index ? block_size : 2*block_size;
2115     int ref_stride= s->current_picture.linesize[plane_index];
2116     uint8_t *dst8= s->current_picture.data[plane_index];
2117     int w= p->width;
2118     int h= p->height;
2119
2120     if(s->keyframe || (s->avctx->debug&512)){
2121         if(mb_y==mb_h)
2122             return;
2123
2124         if(add){
2125             for(y=block_w*mb_y; y<FFMIN(h,block_w*(mb_y+1)); y++){
2126                 for(x=0; x<w; x++){
2127                     int v= buf[x + y*w] + (128<<FRAC_BITS) + (1<<(FRAC_BITS-1));
2128                     v >>= FRAC_BITS;
2129                     if(v&(~255)) v= ~(v>>31);
2130                     dst8[x + y*ref_stride]= v;
2131                 }
2132             }
2133         }else{
2134             for(y=block_w*mb_y; y<FFMIN(h,block_w*(mb_y+1)); y++){
2135                 for(x=0; x<w; x++){
2136                     buf[x + y*w]-= 128<<FRAC_BITS;
2137                 }
2138             }
2139         }
2140
2141         return;
2142     }
2143
2144     for(mb_x=0; mb_x<=mb_w; mb_x++){
2145         add_yblock(s, 0, NULL, buf, dst8, obmc,
2146                    block_w*mb_x - block_w/2,
2147                    block_w*mb_y - block_w/2,
2148                    block_w, block_w,
2149                    w, h,
2150                    w, ref_stride, obmc_stride,
2151                    mb_x - 1, mb_y - 1,
2152                    add, 1, plane_index);
2153     }
2154 }
2155
2156 static av_always_inline void predict_plane(SnowContext *s, IDWTELEM *buf, int plane_index, int add){
2157     const int mb_h= s->b_height << s->block_max_depth;
2158     int mb_y;
2159     for(mb_y=0; mb_y<=mb_h; mb_y++)
2160         predict_slice(s, buf, plane_index, add, mb_y);
2161 }
2162
2163 static void dequantize_slice_buffered(SnowContext *s, slice_buffer * sb, SubBand *b, IDWTELEM *src, int stride, int start_y, int end_y){
2164     const int w= b->width;
2165     const int qlog= av_clip(s->qlog + b->qlog, 0, QROOT*16);
2166     const int qmul= qexp[qlog&(QROOT-1)]<<(qlog>>QSHIFT);
2167     const int qadd= (s->qbias*qmul)>>QBIAS_SHIFT;
2168     int x,y;
2169
2170     if(s->qlog == LOSSLESS_QLOG) return;
2171
2172     for(y=start_y; y<end_y; y++){
2173 //        DWTELEM * line = slice_buffer_get_line_from_address(sb, src + (y * stride));
2174         IDWTELEM * line = slice_buffer_get_line(sb, (y * b->stride_line) + b->buf_y_offset) + b->buf_x_offset;
2175         for(x=0; x<w; x++){
2176             int i= line[x];
2177             if(i<0){
2178                 line[x]= -((-i*qmul + qadd)>>(QEXPSHIFT)); //FIXME try different bias
2179             }else if(i>0){
2180                 line[x]=  (( i*qmul + qadd)>>(QEXPSHIFT));
2181             }
2182         }
2183     }
2184 }
2185
2186 static void correlate_slice_buffered(SnowContext *s, slice_buffer * sb, SubBand *b, IDWTELEM *src, int stride, int inverse, int use_median, int start_y, int end_y){
2187     const int w= b->width;
2188     int x,y;
2189
2190     IDWTELEM * line=0; // silence silly "could be used without having been initialized" warning
2191     IDWTELEM * prev;
2192
2193     if (start_y != 0)
2194         line = slice_buffer_get_line(sb, ((start_y - 1) * b->stride_line) + b->buf_y_offset) + b->buf_x_offset;
2195
2196     for(y=start_y; y<end_y; y++){
2197         prev = line;
2198 //        line = slice_buffer_get_line_from_address(sb, src + (y * stride));
2199         line = slice_buffer_get_line(sb, (y * b->stride_line) + b->buf_y_offset) + b->buf_x_offset;
2200         for(x=0; x<w; x++){
2201             if(x){
2202                 if(use_median){
2203                     if(y && x+1<w) line[x] += mid_pred(line[x - 1], prev[x], prev[x + 1]);
2204                     else  line[x] += line[x - 1];
2205                 }else{
2206                     if(y) line[x] += mid_pred(line[x - 1], prev[x], line[x - 1] + prev[x] - prev[x - 1]);
2207                     else  line[x] += line[x - 1];
2208                 }
2209             }else{
2210                 if(y) line[x] += prev[x];
2211             }
2212         }
2213     }
2214 }
2215
2216 static void decode_qlogs(SnowContext *s){
2217     int plane_index, level, orientation;
2218
2219     for(plane_index=0; plane_index<3; plane_index++){
2220         for(level=0; level<s->spatial_decomposition_count; level++){
2221             for(orientation=level ? 1:0; orientation<4; orientation++){
2222                 int q;
2223                 if     (plane_index==2) q= s->plane[1].band[level][orientation].qlog;
2224                 else if(orientation==2) q= s->plane[plane_index].band[level][1].qlog;
2225                 else                    q= get_symbol(&s->c, s->header_state, 1);
2226                 s->plane[plane_index].band[level][orientation].qlog= q;
2227             }
2228         }
2229     }
2230 }
2231
2232 #define GET_S(dst, check) \
2233     tmp= get_symbol(&s->c, s->header_state, 0);\
2234     if(!(check)){\
2235         av_log(s->avctx, AV_LOG_ERROR, "Error " #dst " is %d\n", tmp);\
2236         return -1;\
2237     }\
2238     dst= tmp;
2239
2240 static int decode_header(SnowContext *s){
2241     int plane_index, tmp;
2242     uint8_t kstate[32];
2243
2244     memset(kstate, MID_STATE, sizeof(kstate));
2245
2246     s->keyframe= get_rac(&s->c, kstate);
2247     if(s->keyframe || s->always_reset){
2248         reset_contexts(s);
2249         s->spatial_decomposition_type=
2250         s->qlog=
2251         s->qbias=
2252         s->mv_scale=
2253         s->block_max_depth= 0;
2254     }
2255     if(s->keyframe){
2256         GET_S(s->version, tmp <= 0U)
2257         s->always_reset= get_rac(&s->c, s->header_state);
2258         s->temporal_decomposition_type= get_symbol(&s->c, s->header_state, 0);
2259         s->temporal_decomposition_count= get_symbol(&s->c, s->header_state, 0);
2260         GET_S(s->spatial_decomposition_count, 0 < tmp && tmp <= MAX_DECOMPOSITIONS)
2261         s->colorspace_type= get_symbol(&s->c, s->header_state, 0);
2262         s->chroma_h_shift= get_symbol(&s->c, s->header_state, 0);
2263         s->chroma_v_shift= get_symbol(&s->c, s->header_state, 0);
2264         s->spatial_scalability= get_rac(&s->c, s->header_state);
2265 //        s->rate_scalability= get_rac(&s->c, s->header_state);
2266         GET_S(s->max_ref_frames, tmp < (unsigned)MAX_REF_FRAMES)
2267         s->max_ref_frames++;
2268
2269         decode_qlogs(s);
2270     }
2271
2272     if(!s->keyframe){
2273         if(get_rac(&s->c, s->header_state)){
2274             for(plane_index=0; plane_index<2; plane_index++){
2275                 int htaps, i, sum=0;
2276                 Plane *p= &s->plane[plane_index];
2277                 p->diag_mc= get_rac(&s->c, s->header_state);
2278                 htaps= get_symbol(&s->c, s->header_state, 0)*2 + 2;
2279                 if((unsigned)htaps > HTAPS_MAX || htaps==0)
2280                     return -1;
2281                 p->htaps= htaps;
2282                 for(i= htaps/2; i; i--){
2283                     p->hcoeff[i]= get_symbol(&s->c, s->header_state, 0) * (1-2*(i&1));
2284                     sum += p->hcoeff[i];
2285                 }
2286                 p->hcoeff[0]= 32-sum;
2287             }
2288             s->plane[2].diag_mc= s->plane[1].diag_mc;
2289             s->plane[2].htaps  = s->plane[1].htaps;
2290             memcpy(s->plane[2].hcoeff, s->plane[1].hcoeff, sizeof(s->plane[1].hcoeff));
2291         }
2292         if(get_rac(&s->c, s->header_state)){
2293             GET_S(s->spatial_decomposition_count, 0 < tmp && tmp <= MAX_DECOMPOSITIONS)
2294             decode_qlogs(s);
2295         }
2296     }
2297
2298     s->spatial_decomposition_type+= get_symbol(&s->c, s->header_state, 1);
2299     if(s->spatial_decomposition_type > 1U){
2300         av_log(s->avctx, AV_LOG_ERROR, "spatial_decomposition_type %d not supported", s->spatial_decomposition_type);
2301         return -1;
2302     }
2303     if(FFMIN(s->avctx-> width>>s->chroma_h_shift,
2304              s->avctx->height>>s->chroma_v_shift) >> (s->spatial_decomposition_count-1) <= 0){
2305         av_log(s->avctx, AV_LOG_ERROR, "spatial_decomposition_count %d too large for size", s->spatial_decomposition_count);
2306         return -1;
2307     }
2308
2309     s->qlog           += get_symbol(&s->c, s->header_state, 1);
2310     s->mv_scale       += get_symbol(&s->c, s->header_state, 1);
2311     s->qbias          += get_symbol(&s->c, s->header_state, 1);
2312     s->block_max_depth+= get_symbol(&s->c, s->header_state, 1);
2313     if(s->block_max_depth > 1 || s->block_max_depth < 0){
2314         av_log(s->avctx, AV_LOG_ERROR, "block_max_depth= %d is too large", s->block_max_depth);
2315         s->block_max_depth= 0;
2316         return -1;
2317     }
2318
2319     return 0;
2320 }
2321
2322 static void init_qexp(void){
2323     int i;
2324     double v=128;
2325
2326     for(i=0; i<QROOT; i++){
2327         qexp[i]= lrintf(v);
2328         v *= pow(2, 1.0 / QROOT);
2329     }
2330 }
2331
2332 static av_cold int common_init(AVCodecContext *avctx){
2333     SnowContext *s = avctx->priv_data;
2334     int width, height;
2335     int i, j;
2336
2337     s->avctx= avctx;
2338     s->max_ref_frames=1; //just make sure its not an invalid value in case of no initial keyframe
2339
2340     dsputil_init(&s->dsp, avctx);
2341
2342 #define mcf(dx,dy)\
2343     s->dsp.put_qpel_pixels_tab       [0][dy+dx/4]=\
2344     s->dsp.put_no_rnd_qpel_pixels_tab[0][dy+dx/4]=\
2345         s->dsp.put_h264_qpel_pixels_tab[0][dy+dx/4];\
2346     s->dsp.put_qpel_pixels_tab       [1][dy+dx/4]=\
2347     s->dsp.put_no_rnd_qpel_pixels_tab[1][dy+dx/4]=\
2348         s->dsp.put_h264_qpel_pixels_tab[1][dy+dx/4];
2349
2350     mcf( 0, 0)
2351     mcf( 4, 0)
2352     mcf( 8, 0)
2353     mcf(12, 0)
2354     mcf( 0, 4)
2355     mcf( 4, 4)
2356     mcf( 8, 4)
2357     mcf(12, 4)
2358     mcf( 0, 8)
2359     mcf( 4, 8)
2360     mcf( 8, 8)
2361     mcf(12, 8)
2362     mcf( 0,12)
2363     mcf( 4,12)
2364     mcf( 8,12)
2365     mcf(12,12)
2366
2367 #define mcfh(dx,dy)\
2368     s->dsp.put_pixels_tab       [0][dy/4+dx/8]=\
2369     s->dsp.put_no_rnd_pixels_tab[0][dy/4+dx/8]=\
2370         mc_block_hpel ## dx ## dy ## 16;\
2371     s->dsp.put_pixels_tab       [1][dy/4+dx/8]=\
2372     s->dsp.put_no_rnd_pixels_tab[1][dy/4+dx/8]=\
2373         mc_block_hpel ## dx ## dy ## 8;
2374
2375     mcfh(0, 0)
2376     mcfh(8, 0)
2377     mcfh(0, 8)
2378     mcfh(8, 8)
2379
2380     if(!qexp[0])
2381         init_qexp();
2382
2383 //    dec += FFMAX(s->chroma_h_shift, s->chroma_v_shift);
2384
2385     width= s->avctx->width;
2386     height= s->avctx->height;
2387
2388     s->spatial_idwt_buffer= av_mallocz(width*height*sizeof(IDWTELEM));
2389     s->spatial_dwt_buffer= av_mallocz(width*height*sizeof(DWTELEM)); //FIXME this does not belong here
2390
2391     for(i=0; i<MAX_REF_FRAMES; i++)
2392         for(j=0; j<MAX_REF_FRAMES; j++)
2393             scale_mv_ref[i][j] = 256*(i+1)/(j+1);
2394
2395     s->avctx->get_buffer(s->avctx, &s->mconly_picture);
2396     s->scratchbuf = av_malloc(s->mconly_picture.linesize[0]*7*MB_SIZE);
2397
2398     return 0;
2399 }
2400
2401 static int common_init_after_header(AVCodecContext *avctx){
2402     SnowContext *s = avctx->priv_data;
2403     int plane_index, level, orientation;
2404
2405     for(plane_index=0; plane_index<3; plane_index++){
2406         int w= s->avctx->width;
2407         int h= s->avctx->height;
2408
2409         if(plane_index){
2410             w>>= s->chroma_h_shift;
2411             h>>= s->chroma_v_shift;
2412         }
2413         s->plane[plane_index].width = w;
2414         s->plane[plane_index].height= h;
2415
2416         for(level=s->spatial_decomposition_count-1; level>=0; level--){
2417             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2418                 SubBand *b= &s->plane[plane_index].band[level][orientation];
2419
2420                 b->buf= s->spatial_dwt_buffer;
2421                 b->level= level;
2422                 b->stride= s->plane[plane_index].width << (s->spatial_decomposition_count - level);
2423                 b->width = (w + !(orientation&1))>>1;
2424                 b->height= (h + !(orientation>1))>>1;
2425
2426                 b->stride_line = 1 << (s->spatial_decomposition_count - level);
2427                 b->buf_x_offset = 0;
2428                 b->buf_y_offset = 0;
2429
2430                 if(orientation&1){
2431                     b->buf += (w+1)>>1;
2432                     b->buf_x_offset = (w+1)>>1;
2433                 }
2434                 if(orientation>1){
2435                     b->buf += b->stride>>1;
2436                     b->buf_y_offset = b->stride_line >> 1;
2437                 }
2438                 b->ibuf= s->spatial_idwt_buffer + (b->buf - s->spatial_dwt_buffer);
2439
2440                 if(level)
2441                     b->parent= &s->plane[plane_index].band[level-1][orientation];
2442                 //FIXME avoid this realloc
2443                 av_freep(&b->x_coeff);
2444                 b->x_coeff=av_mallocz(((b->width+1) * b->height+1)*sizeof(x_and_coeff));
2445             }
2446             w= (w+1)>>1;
2447             h= (h+1)>>1;
2448         }
2449     }
2450
2451     return 0;
2452 }
2453
2454 #define QUANTIZE2 0
2455
2456 #if QUANTIZE2==1
2457 #define Q2_STEP 8
2458
2459 static void find_sse(SnowContext *s, Plane *p, int *score, int score_stride, IDWTELEM *r0, IDWTELEM *r1, int level, int orientation){
2460     SubBand *b= &p->band[level][orientation];
2461     int x, y;
2462     int xo=0;
2463     int yo=0;
2464     int step= 1 << (s->spatial_decomposition_count - level);
2465
2466     if(orientation&1)
2467         xo= step>>1;
2468     if(orientation&2)
2469         yo= step>>1;
2470
2471     //FIXME bias for nonzero ?
2472     //FIXME optimize
2473     memset(score, 0, sizeof(*score)*score_stride*((p->height + Q2_STEP-1)/Q2_STEP));
2474     for(y=0; y<p->height; y++){
2475         for(x=0; x<p->width; x++){
2476             int sx= (x-xo + step/2) / step / Q2_STEP;
2477             int sy= (y-yo + step/2) / step / Q2_STEP;
2478             int v= r0[x + y*p->width] - r1[x + y*p->width];
2479             assert(sx>=0 && sy>=0 && sx < score_stride);
2480             v= ((v+8)>>4)<<4;
2481             score[sx + sy*score_stride] += v*v;
2482             assert(score[sx + sy*score_stride] >= 0);
2483         }
2484     }
2485 }
2486
2487 static void dequantize_all(SnowContext *s, Plane *p, IDWTELEM *buffer, int width, int height){
2488     int level, orientation;
2489
2490     for(level=0; level<s->spatial_decomposition_count; level++){
2491         for(orientation=level ? 1 : 0; orientation<4; orientation++){
2492             SubBand *b= &p->band[level][orientation];
2493             IDWTELEM *dst= buffer + (b->ibuf - s->spatial_idwt_buffer);
2494
2495             dequantize(s, b, dst, b->stride);
2496         }
2497     }
2498 }
2499
2500 static void dwt_quantize(SnowContext *s, Plane *p, DWTELEM *buffer, int width, int height, int stride, int type){
2501     int level, orientation, ys, xs, x, y, pass;
2502     IDWTELEM best_dequant[height * stride];
2503     IDWTELEM idwt2_buffer[height * stride];
2504     const int score_stride= (width + 10)/Q2_STEP;
2505     int best_score[(width + 10)/Q2_STEP * (height + 10)/Q2_STEP]; //FIXME size
2506     int score[(width + 10)/Q2_STEP * (height + 10)/Q2_STEP]; //FIXME size
2507     int threshold= (s->m.lambda * s->m.lambda) >> 6;
2508
2509     //FIXME pass the copy cleanly ?
2510
2511 //    memcpy(dwt_buffer, buffer, height * stride * sizeof(DWTELEM));
2512     ff_spatial_dwt(buffer, width, height, stride, type, s->spatial_decomposition_count);
2513
2514     for(level=0; level<s->spatial_decomposition_count; level++){
2515         for(orientation=level ? 1 : 0; orientation<4; orientation++){
2516             SubBand *b= &p->band[level][orientation];
2517             IDWTELEM *dst= best_dequant + (b->ibuf - s->spatial_idwt_buffer);
2518              DWTELEM *src=       buffer + (b-> buf - s->spatial_dwt_buffer);
2519             assert(src == b->buf); // code does not depend on this but it is true currently
2520
2521             quantize(s, b, dst, src, b->stride, s->qbias);
2522         }
2523     }
2524     for(pass=0; pass<1; pass++){
2525         if(s->qbias == 0) //keyframe
2526             continue;
2527         for(level=0; level<s->spatial_decomposition_count; level++){
2528             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2529                 SubBand *b= &p->band[level][orientation];
2530                 IDWTELEM *dst= idwt2_buffer + (b->ibuf - s->spatial_idwt_buffer);
2531                 IDWTELEM *best_dst= best_dequant + (b->ibuf - s->spatial_idwt_buffer);
2532
2533                 for(ys= 0; ys<Q2_STEP; ys++){
2534                     for(xs= 0; xs<Q2_STEP; xs++){
2535                         memcpy(idwt2_buffer, best_dequant, height * stride * sizeof(IDWTELEM));
2536                         dequantize_all(s, p, idwt2_buffer, width, height);
2537                         ff_spatial_idwt(idwt2_buffer, width, height, stride, type, s->spatial_decomposition_count);
2538                         find_sse(s, p, best_score, score_stride, idwt2_buffer, s->spatial_idwt_buffer, level, orientation);
2539                         memcpy(idwt2_buffer, best_dequant, height * stride * sizeof(IDWTELEM));
2540                         for(y=ys; y<b->height; y+= Q2_STEP){
2541                             for(x=xs; x<b->width; x+= Q2_STEP){
2542                                 if(dst[x + y*b->stride]<0) dst[x + y*b->stride]++;
2543                                 if(dst[x + y*b->stride]>0) dst[x + y*b->stride]--;
2544                                 //FIXME try more than just --
2545                             }
2546                         }
2547                         dequantize_all(s, p, idwt2_buffer, width, height);
2548                         ff_spatial_idwt(idwt2_buffer, width, height, stride, type, s->spatial_decomposition_count);
2549                         find_sse(s, p, score, score_stride, idwt2_buffer, s->spatial_idwt_buffer, level, orientation);
2550                         for(y=ys; y<b->height; y+= Q2_STEP){
2551                             for(x=xs; x<b->width; x+= Q2_STEP){
2552                                 int score_idx= x/Q2_STEP + (y/Q2_STEP)*score_stride;
2553                                 if(score[score_idx] <= best_score[score_idx] + threshold){
2554                                     best_score[score_idx]= score[score_idx];
2555                                     if(best_dst[x + y*b->stride]<0) best_dst[x + y*b->stride]++;
2556                                     if(best_dst[x + y*b->stride]>0) best_dst[x + y*b->stride]--;
2557                                     //FIXME copy instead
2558                                 }
2559                             }
2560                         }
2561                     }
2562                 }
2563             }
2564         }
2565     }
2566     memcpy(s->spatial_idwt_buffer, best_dequant, height * stride * sizeof(IDWTELEM)); //FIXME work with that directly instead of copy at the end
2567 }
2568
2569 #endif /* QUANTIZE2==1 */
2570
2571 #define USE_HALFPEL_PLANE 0
2572
2573 static void halfpel_interpol(SnowContext *s, uint8_t *halfpel[4][4], AVFrame *frame){
2574     int p,x,y;
2575
2576     assert(!(s->avctx->flags & CODEC_FLAG_EMU_EDGE));
2577
2578     for(p=0; p<3; p++){
2579         int is_chroma= !!p;
2580         int w= s->avctx->width  >>is_chroma;
2581         int h= s->avctx->height >>is_chroma;
2582         int ls= frame->linesize[p];
2583         uint8_t *src= frame->data[p];
2584
2585         halfpel[1][p]= (uint8_t*)av_malloc(ls * (h+2*EDGE_WIDTH)) + EDGE_WIDTH*(1+ls);
2586         halfpel[2][p]= (uint8_t*)av_malloc(ls * (h+2*EDGE_WIDTH)) + EDGE_WIDTH*(1+ls);
2587         halfpel[3][p]= (uint8_t*)av_malloc(ls * (h+2*EDGE_WIDTH)) + EDGE_WIDTH*(1+ls);
2588
2589         halfpel[0][p]= src;
2590         for(y=0; y<h; y++){
2591             for(x=0; x<w; x++){
2592                 int i= y*ls + x;
2593
2594                 halfpel[1][p][i]= (20*(src[i] + src[i+1]) - 5*(src[i-1] + src[i+2]) + (src[i-2] + src[i+3]) + 16 )>>5;
2595             }
2596         }
2597         for(y=0; y<h; y++){
2598             for(x=0; x<w; x++){
2599                 int i= y*ls + x;
2600
2601                 halfpel[2][p][i]= (20*(src[i] + src[i+ls]) - 5*(src[i-ls] + src[i+2*ls]) + (src[i-2*ls] + src[i+3*ls]) + 16 )>>5;
2602             }
2603         }
2604         src= halfpel[1][p];
2605         for(y=0; y<h; y++){
2606             for(x=0; x<w; x++){
2607                 int i= y*ls + x;
2608
2609                 halfpel[3][p][i]= (20*(src[i] + src[i+ls]) - 5*(src[i-ls] + src[i+2*ls]) + (src[i-2*ls] + src[i+3*ls]) + 16 )>>5;
2610             }
2611         }
2612
2613 //FIXME border!
2614     }
2615 }
2616
2617 static void release_buffer(AVCodecContext *avctx){
2618     SnowContext *s = avctx->priv_data;
2619     int i;
2620
2621     if(s->last_picture[s->max_ref_frames-1].data[0]){
2622         avctx->release_buffer(avctx, &s->last_picture[s->max_ref_frames-1]);
2623         for(i=0; i<9; i++)
2624             if(s->halfpel_plane[s->max_ref_frames-1][1+i/3][i%3])
2625                 av_free(s->halfpel_plane[s->max_ref_frames-1][1+i/3][i%3] - EDGE_WIDTH*(1+s->current_picture.linesize[i%3]));
2626     }
2627 }
2628
2629 static int frame_start(SnowContext *s){
2630    AVFrame tmp;
2631    int w= s->avctx->width; //FIXME round up to x16 ?
2632    int h= s->avctx->height;
2633
2634     if(s->current_picture.data[0]){
2635         s->dsp.draw_edges(s->current_picture.data[0], s->current_picture.linesize[0], w   , h   , EDGE_WIDTH  );
2636         s->dsp.draw_edges(s->current_picture.data[1], s->current_picture.linesize[1], w>>1, h>>1, EDGE_WIDTH/2);
2637         s->dsp.draw_edges(s->current_picture.data[2], s->current_picture.linesize[2], w>>1, h>>1, EDGE_WIDTH/2);
2638     }
2639
2640     release_buffer(s->avctx);
2641
2642     tmp= s->last_picture[s->max_ref_frames-1];
2643     memmove(s->last_picture+1, s->last_picture, (s->max_ref_frames-1)*sizeof(AVFrame));
2644     memmove(s->halfpel_plane+1, s->halfpel_plane, (s->max_ref_frames-1)*sizeof(void*)*4*4);
2645     if(USE_HALFPEL_PLANE && s->current_picture.data[0])
2646         halfpel_interpol(s, s->halfpel_plane[0], &s->current_picture);
2647     s->last_picture[0]= s->current_picture;
2648     s->current_picture= tmp;
2649
2650     if(s->keyframe){
2651         s->ref_frames= 0;
2652     }else{
2653         int i;
2654         for(i=0; i<s->max_ref_frames && s->last_picture[i].data[0]; i++)
2655             if(i && s->last_picture[i-1].key_frame)
2656                 break;
2657         s->ref_frames= i;
2658         if(s->ref_frames==0){
2659             av_log(s->avctx,AV_LOG_ERROR, "No reference frames\n");
2660             return -1;
2661         }
2662     }
2663
2664     s->current_picture.reference= 1;
2665     if(s->avctx->get_buffer(s->avctx, &s->current_picture) < 0){
2666         av_log(s->avctx, AV_LOG_ERROR, "get_buffer() failed\n");
2667         return -1;
2668     }
2669
2670     s->current_picture.key_frame= s->keyframe;
2671
2672     return 0;
2673 }
2674
2675 static av_cold void common_end(SnowContext *s){
2676     int plane_index, level, orientation, i;
2677
2678     av_freep(&s->spatial_dwt_buffer);
2679     av_freep(&s->spatial_idwt_buffer);
2680
2681     s->m.me.temp= NULL;
2682     av_freep(&s->m.me.scratchpad);
2683     av_freep(&s->m.me.map);
2684     av_freep(&s->m.me.score_map);
2685     av_freep(&s->m.obmc_scratchpad);
2686
2687     av_freep(&s->block);
2688     av_freep(&s->scratchbuf);
2689
2690     for(i=0; i<MAX_REF_FRAMES; i++){
2691         av_freep(&s->ref_mvs[i]);
2692         av_freep(&s->ref_scores[i]);
2693         if(s->last_picture[i].data[0])
2694             s->avctx->release_buffer(s->avctx, &s->last_picture[i]);
2695     }
2696
2697     for(plane_index=0; plane_index<3; plane_index++){
2698         for(level=s->spatial_decomposition_count-1; level>=0; level--){
2699             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2700                 SubBand *b= &s->plane[plane_index].band[level][orientation];
2701
2702                 av_freep(&b->x_coeff);
2703             }
2704         }
2705     }
2706     if (s->mconly_picture.data[0])
2707         s->avctx->release_buffer(s->avctx, &s->mconly_picture);
2708     if (s->current_picture.data[0])
2709         s->avctx->release_buffer(s->avctx, &s->current_picture);
2710 }
2711
2712 static av_cold int decode_init(AVCodecContext *avctx)
2713 {
2714     avctx->pix_fmt= PIX_FMT_YUV420P;
2715
2716     common_init(avctx);
2717
2718     return 0;
2719 }
2720
2721 static int decode_frame(AVCodecContext *avctx, void *data, int *data_size, AVPacket *avpkt){
2722     const uint8_t *buf = avpkt->data;
2723     int buf_size = avpkt->size;
2724     SnowContext *s = avctx->priv_data;
2725     RangeCoder * const c= &s->c;
2726     int bytes_read;
2727     AVFrame *picture = data;
2728     int level, orientation, plane_index;
2729
2730     ff_init_range_decoder(c, buf, buf_size);
2731     ff_build_rac_states(c, 0.05*(1LL<<32), 256-8);
2732
2733     s->current_picture.pict_type= FF_I_TYPE; //FIXME I vs. P
2734     if(decode_header(s)<0)
2735         return -1;
2736     common_init_after_header(avctx);
2737
2738     // realloc slice buffer for the case that spatial_decomposition_count changed
2739     slice_buffer_destroy(&s->sb);
2740     slice_buffer_init(&s->sb, s->plane[0].height, (MB_SIZE >> s->block_max_depth) + s->spatial_decomposition_count * 8 + 1, s->plane[0].width, s->spatial_idwt_buffer);
2741
2742     for(plane_index=0; plane_index<3; plane_index++){
2743         Plane *p= &s->plane[plane_index];
2744         p->fast_mc= p->diag_mc && p->htaps==6 && p->hcoeff[0]==40
2745                                               && p->hcoeff[1]==-10
2746                                               && p->hcoeff[2]==2;
2747     }
2748
2749     alloc_blocks(s);
2750
2751     if(frame_start(s) < 0)
2752         return -1;
2753     //keyframe flag duplication mess FIXME
2754     if(avctx->debug&FF_DEBUG_PICT_INFO)
2755         av_log(avctx, AV_LOG_ERROR, "keyframe:%d qlog:%d\n", s->keyframe, s->qlog);
2756
2757     decode_blocks(s);
2758
2759     for(plane_index=0; plane_index<3; plane_index++){
2760         Plane *p= &s->plane[plane_index];
2761         int w= p->width;
2762         int h= p->height;
2763         int x, y;
2764         int decode_state[MAX_DECOMPOSITIONS][4][1]; /* Stored state info for unpack_coeffs. 1 variable per instance. */
2765
2766         if(s->avctx->debug&2048){
2767             memset(s->spatial_dwt_buffer, 0, sizeof(DWTELEM)*w*h);
2768             predict_plane(s, s->spatial_idwt_buffer, plane_index, 1);
2769
2770             for(y=0; y<h; y++){
2771                 for(x=0; x<w; x++){
2772                     int v= s->current_picture.data[plane_index][y*s->current_picture.linesize[plane_index] + x];
2773                     s->mconly_picture.data[plane_index][y*s->mconly_picture.linesize[plane_index] + x]= v;
2774                 }
2775             }
2776         }
2777
2778         {
2779         for(level=0; level<s->spatial_decomposition_count; level++){
2780             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2781                 SubBand *b= &p->band[level][orientation];
2782                 unpack_coeffs(s, b, b->parent, orientation);
2783             }
2784         }
2785         }
2786
2787         {
2788         const int mb_h= s->b_height << s->block_max_depth;
2789         const int block_size = MB_SIZE >> s->block_max_depth;
2790         const int block_w    = plane_index ? block_size/2 : block_size;
2791         int mb_y;
2792         DWTCompose cs[MAX_DECOMPOSITIONS];
2793         int yd=0, yq=0;
2794         int y;
2795         int end_y;
2796
2797         ff_spatial_idwt_buffered_init(cs, &s->sb, w, h, 1, s->spatial_decomposition_type, s->spatial_decomposition_count);
2798         for(mb_y=0; mb_y<=mb_h; mb_y++){
2799
2800             int slice_starty = block_w*mb_y;
2801             int slice_h = block_w*(mb_y+1);
2802             if (!(s->keyframe || s->avctx->debug&512)){
2803                 slice_starty = FFMAX(0, slice_starty - (block_w >> 1));
2804                 slice_h -= (block_w >> 1);
2805             }
2806
2807             for(level=0; level<s->spatial_decomposition_count; level++){
2808                 for(orientation=level ? 1 : 0; orientation<4; orientation++){
2809                     SubBand *b= &p->band[level][orientation];
2810                     int start_y;
2811                     int end_y;
2812                     int our_mb_start = mb_y;
2813                     int our_mb_end = (mb_y + 1);
2814                     const int extra= 3;
2815                     start_y = (mb_y ? ((block_w * our_mb_start) >> (s->spatial_decomposition_count - level)) + s->spatial_decomposition_count - level + extra: 0);
2816                     end_y = (((block_w * our_mb_end) >> (s->spatial_decomposition_count - level)) + s->spatial_decomposition_count - level + extra);
2817                     if (!(s->keyframe || s->avctx->debug&512)){
2818                         start_y = FFMAX(0, start_y - (block_w >> (1+s->spatial_decomposition_count - level)));
2819                         end_y = FFMAX(0, end_y - (block_w >> (1+s->spatial_decomposition_count - level)));
2820                     }
2821                     start_y = FFMIN(b->height, start_y);
2822                     end_y = FFMIN(b->height, end_y);
2823
2824                     if (start_y != end_y){
2825                         if (orientation == 0){
2826                             SubBand * correlate_band = &p->band[0][0];
2827                             int correlate_end_y = FFMIN(b->height, end_y + 1);
2828                             int correlate_start_y = FFMIN(b->height, (start_y ? start_y + 1 : 0));
2829                             decode_subband_slice_buffered(s, correlate_band, &s->sb, correlate_start_y, correlate_end_y, decode_state[0][0]);
2830                             correlate_slice_buffered(s, &s->sb, correlate_band, correlate_band->ibuf, correlate_band->stride, 1, 0, correlate_start_y, correlate_end_y);
2831                             dequantize_slice_buffered(s, &s->sb, correlate_band, correlate_band->ibuf, correlate_band->stride, start_y, end_y);
2832                         }
2833                         else
2834                             decode_subband_slice_buffered(s, b, &s->sb, start_y, end_y, decode_state[level][orientation]);
2835                     }
2836                 }
2837             }
2838
2839             for(; yd<slice_h; yd+=4){
2840                 ff_spatial_idwt_buffered_slice(&s->dsp, cs, &s->sb, w, h, 1, s->spatial_decomposition_type, s->spatial_decomposition_count, yd);
2841             }
2842
2843             if(s->qlog == LOSSLESS_QLOG){
2844                 for(; yq<slice_h && yq<h; yq++){
2845                     IDWTELEM * line = slice_buffer_get_line(&s->sb, yq);
2846                     for(x=0; x<w; x++){
2847                         line[x] <<= FRAC_BITS;
2848                     }
2849                 }
2850             }
2851
2852             predict_slice_buffered(s, &s->sb, s->spatial_idwt_buffer, plane_index, 1, mb_y);
2853
2854             y = FFMIN(p->height, slice_starty);
2855             end_y = FFMIN(p->height, slice_h);
2856             while(y < end_y)
2857                 slice_buffer_release(&s->sb, y++);
2858         }
2859
2860         slice_buffer_flush(&s->sb);
2861         }
2862
2863     }
2864
2865     emms_c();
2866
2867     release_buffer(avctx);
2868
2869     if(!(s->avctx->debug&2048))
2870         *picture= s->current_picture;
2871     else
2872         *picture= s->mconly_picture;
2873
2874     *data_size = sizeof(AVFrame);
2875
2876     bytes_read= c->bytestream - c->bytestream_start;
2877     if(bytes_read ==0) av_log(s->avctx, AV_LOG_ERROR, "error at end of frame\n"); //FIXME
2878
2879     return bytes_read;
2880 }
2881
2882 static av_cold int decode_end(AVCodecContext *avctx)
2883 {
2884     SnowContext *s = avctx->priv_data;
2885
2886     slice_buffer_destroy(&s->sb);
2887
2888     common_end(s);
2889
2890     return 0;
2891 }
2892
2893 AVCodec snow_decoder = {
2894     "snow",
2895     CODEC_TYPE_VIDEO,
2896     CODEC_ID_SNOW,
2897     sizeof(SnowContext),
2898     decode_init,
2899     NULL,
2900     decode_end,
2901     decode_frame,
2902     CODEC_CAP_DR1 /*| CODEC_CAP_DRAW_HORIZ_BAND*/,
2903     NULL,
2904     .long_name = NULL_IF_CONFIG_SMALL("Snow"),
2905 };
2906
2907 #if CONFIG_SNOW_ENCODER
2908 static av_cold int encode_init(AVCodecContext *avctx)
2909 {
2910     SnowContext *s = avctx->priv_data;
2911     int plane_index;
2912
2913     if(avctx->strict_std_compliance > FF_COMPLIANCE_EXPERIMENTAL){
2914         av_log(avctx, AV_LOG_ERROR, "This codec is under development, files encoded with it may not be decodable with future versions!!!\n"
2915                "Use vstrict=-2 / -strict -2 to use it anyway.\n");
2916         return -1;
2917     }
2918
2919     if(avctx->prediction_method == DWT_97
2920        && (avctx->flags & CODEC_FLAG_QSCALE)
2921        && avctx->global_quality == 0){
2922         av_log(avctx, AV_LOG_ERROR, "The 9/7 wavelet is incompatible with lossless mode.\n");
2923         return -1;
2924     }
2925
2926     s->spatial_decomposition_type= avctx->prediction_method; //FIXME add decorrelator type r transform_type
2927
2928     s->mv_scale       = (avctx->flags & CODEC_FLAG_QPEL) ? 2 : 4;
2929     s->block_max_depth= (avctx->flags & CODEC_FLAG_4MV ) ? 1 : 0;
2930
2931     for(plane_index=0; plane_index<3; plane_index++){
2932         s->plane[plane_index].diag_mc= 1;
2933         s->plane[plane_index].htaps= 6;
2934         s->plane[plane_index].hcoeff[0]=  40;
2935         s->plane[plane_index].hcoeff[1]= -10;
2936         s->plane[plane_index].hcoeff[2]=   2;
2937         s->plane[plane_index].fast_mc= 1;
2938     }
2939
2940     common_init(avctx);
2941     alloc_blocks(s);
2942
2943     s->version=0;
2944
2945     s->m.avctx   = avctx;
2946     s->m.flags   = avctx->flags;
2947     s->m.bit_rate= avctx->bit_rate;
2948
2949     s->m.me.temp      =
2950     s->m.me.scratchpad= av_mallocz((avctx->width+64)*2*16*2*sizeof(uint8_t));
2951     s->m.me.map       = av_mallocz(ME_MAP_SIZE*sizeof(uint32_t));
2952     s->m.me.score_map = av_mallocz(ME_MAP_SIZE*sizeof(uint32_t));
2953     s->m.obmc_scratchpad= av_mallocz(MB_SIZE*MB_SIZE*12*sizeof(uint32_t));
2954     h263_encode_init(&s->m); //mv_penalty
2955
2956     s->max_ref_frames = FFMAX(FFMIN(avctx->refs, MAX_REF_FRAMES), 1);
2957
2958     if(avctx->flags&CODEC_FLAG_PASS1){
2959         if(!avctx->stats_out)
2960             avctx->stats_out = av_mallocz(256);
2961     }
2962     if((avctx->flags&CODEC_FLAG_PASS2) || !(avctx->flags&CODEC_FLAG_QSCALE)){
2963         if(ff_rate_control_init(&s->m) < 0)
2964             return -1;
2965     }
2966     s->pass1_rc= !(avctx->flags & (CODEC_FLAG_QSCALE|CODEC_FLAG_PASS2));
2967
2968     avctx->coded_frame= &s->current_picture;
2969     switch(avctx->pix_fmt){
2970 //    case PIX_FMT_YUV444P:
2971 //    case PIX_FMT_YUV422P:
2972     case PIX_FMT_YUV420P:
2973     case PIX_FMT_GRAY8:
2974 //    case PIX_FMT_YUV411P:
2975 //    case PIX_FMT_YUV410P:
2976         s->colorspace_type= 0;
2977         break;
2978 /*    case PIX_FMT_RGB32:
2979         s->colorspace= 1;
2980         break;*/
2981     default:
2982         av_log(avctx, AV_LOG_ERROR, "pixel format not supported\n");
2983         return -1;
2984     }
2985 //    avcodec_get_chroma_sub_sample(avctx->pix_fmt, &s->chroma_h_shift, &s->chroma_v_shift);
2986     s->chroma_h_shift= 1;
2987     s->chroma_v_shift= 1;
2988
2989     ff_set_cmp(&s->dsp, s->dsp.me_cmp, s->avctx->me_cmp);
2990     ff_set_cmp(&s->dsp, s->dsp.me_sub_cmp, s->avctx->me_sub_cmp);
2991
2992     s->avctx->get_buffer(s->avctx, &s->input_picture);
2993
2994     if(s->avctx->me_method == ME_ITER){
2995         int i;
2996         int size= s->b_width * s->b_height << 2*s->block_max_depth;
2997         for(i=0; i<s->max_ref_frames; i++){
2998             s->ref_mvs[i]= av_mallocz(size*sizeof(int16_t[2]));
2999             s->ref_scores[i]= av_mallocz(size*sizeof(uint32_t));
3000         }
3001     }
3002
3003     return 0;
3004 }
3005
3006 //near copy & paste from dsputil, FIXME
3007 static int pix_sum(uint8_t * pix, int line_size, int w)
3008 {
3009     int s, i, j;
3010
3011     s = 0;
3012     for (i = 0; i < w; i++) {
3013         for (j = 0; j < w; j++) {
3014             s += pix[0];
3015             pix ++;
3016         }
3017         pix += line_size - w;
3018     }
3019     return s;
3020 }
3021
3022 //near copy & paste from dsputil, FIXME
3023 static int pix_norm1(uint8_t * pix, int line_size, int w)
3024 {
3025     int s, i, j;
3026     uint32_t *sq = ff_squareTbl + 256;
3027
3028     s = 0;
3029     for (i = 0; i < w; i++) {
3030         for (j = 0; j < w; j ++) {
3031             s += sq[pix[0]];
3032             pix ++;
3033         }
3034         pix += line_size - w;
3035     }
3036     return s;
3037 }
3038
3039 //FIXME copy&paste
3040 #define P_LEFT P[1]
3041 #define P_TOP P[2]
3042 #define P_TOPRIGHT P[3]
3043 #define P_MEDIAN P[4]
3044 #define P_MV1 P[9]
3045 #define FLAG_QPEL   1 //must be 1
3046
3047 static int encode_q_branch(SnowContext *s, int level, int x, int y){
3048     uint8_t p_buffer[1024];
3049     uint8_t i_buffer[1024];
3050     uint8_t p_state[sizeof(s->block_state)];
3051     uint8_t i_state[sizeof(s->block_state)];
3052     RangeCoder pc, ic;
3053     uint8_t *pbbak= s->c.bytestream;
3054     uint8_t *pbbak_start= s->c.bytestream_start;
3055     int score, score2, iscore, i_len, p_len, block_s, sum, base_bits;
3056     const int w= s->b_width  << s->block_max_depth;
3057     const int h= s->b_height << s->block_max_depth;
3058     const int rem_depth= s->block_max_depth - level;
3059     const int index= (x + y*w) << rem_depth;
3060     const int block_w= 1<<(LOG2_MB_SIZE - level);
3061     int trx= (x+1)<<rem_depth;
3062     int try= (y+1)<<rem_depth;
3063     const BlockNode *left  = x ? &s->block[index-1] : &null_block;
3064     const BlockNode *top   = y ? &s->block[index-w] : &null_block;
3065     const BlockNode *right = trx<w ? &s->block[index+1] : &null_block;
3066     const BlockNode *bottom= try<h ? &s->block[index+w] : &null_block;
3067     const BlockNode *tl    = y && x ? &s->block[index-w-1] : left;
3068     const BlockNode *tr    = y && trx<w && ((x&1)==0 || level==0) ? &s->block[index-w+(1<<rem_depth)] : tl; //FIXME use lt
3069     int pl = left->color[0];
3070     int pcb= left->color[1];
3071     int pcr= left->color[2];
3072     int pmx, pmy;
3073     int mx=0, my=0;
3074     int l,cr,cb;
3075     const int stride= s->current_picture.linesize[0];
3076     const int uvstride= s->current_picture.linesize[1];
3077     uint8_t *current_data[3]= { s->input_picture.data[0] + (x + y*  stride)*block_w,
3078                                 s->input_picture.data[1] + (x + y*uvstride)*block_w/2,
3079                                 s->input_picture.data[2] + (x + y*uvstride)*block_w/2};
3080     int P[10][2];
3081     int16_t last_mv[3][2];
3082     int qpel= !!(s->avctx->flags & CODEC_FLAG_QPEL); //unused
3083     const int shift= 1+qpel;
3084     MotionEstContext *c= &s->m.me;
3085     int ref_context= av_log2(2*left->ref) + av_log2(2*top->ref);
3086     int mx_context= av_log2(2*FFABS(left->mx - top->mx));
3087     int my_context= av_log2(2*FFABS(left->my - top->my));
3088     int s_context= 2*left->level + 2*top->level + tl->level + tr->level;
3089     int ref, best_ref, ref_score, ref_mx, ref_my;
3090
3091     assert(sizeof(s->block_state) >= 256);
3092     if(s->keyframe){
3093         set_blocks(s, level, x, y, pl, pcb, pcr, 0, 0, 0, BLOCK_INTRA);
3094         return 0;
3095     }
3096
3097 //    clip predictors / edge ?
3098
3099     P_LEFT[0]= left->mx;
3100     P_LEFT[1]= left->my;
3101     P_TOP [0]= top->mx;
3102     P_TOP [1]= top->my;
3103     P_TOPRIGHT[0]= tr->mx;
3104     P_TOPRIGHT[1]= tr->my;
3105
3106     last_mv[0][0]= s->block[index].mx;
3107     last_mv[0][1]= s->block[index].my;
3108     last_mv[1][0]= right->mx;
3109     last_mv[1][1]= right->my;
3110     last_mv[2][0]= bottom->mx;
3111     last_mv[2][1]= bottom->my;
3112
3113     s->m.mb_stride=2;
3114     s->m.mb_x=
3115     s->m.mb_y= 0;
3116     c->skip= 0;
3117
3118     assert(c->  stride ==   stride);
3119     assert(c->uvstride == uvstride);
3120
3121     c->penalty_factor    = get_penalty_factor(s->lambda, s->lambda2, c->avctx->me_cmp);
3122     c->sub_penalty_factor= get_penalty_factor(s->lambda, s->lambda2, c->avctx->me_sub_cmp);
3123     c->mb_penalty_factor = get_penalty_factor(s->lambda, s->lambda2, c->avctx->mb_cmp);
3124     c->current_mv_penalty= c->mv_penalty[s->m.f_code=1] + MAX_MV;
3125
3126     c->xmin = - x*block_w - 16+3;
3127     c->ymin = - y*block_w - 16+3;
3128     c->xmax = - (x+1)*block_w + (w<<(LOG2_MB_SIZE - s->block_max_depth)) + 16-3;
3129     c->ymax = - (y+1)*block_w + (h<<(LOG2_MB_SIZE - s->block_max_depth)) + 16-3;
3130
3131     if(P_LEFT[0]     > (c->xmax<<shift)) P_LEFT[0]    = (c->xmax<<shift);
3132     if(P_LEFT[1]     > (c->ymax<<shift)) P_LEFT[1]    = (c->ymax<<shift);
3133     if(P_TOP[0]      > (c->xmax<<shift)) P_TOP[0]     = (c->xmax<<shift);
3134     if(P_TOP[1]      > (c->ymax<<shift)) P_TOP[1]     = (c->ymax<<shift);
3135     if(P_TOPRIGHT[0] < (c->xmin<<shift)) P_TOPRIGHT[0]= (c->xmin<<shift);
3136     if(P_TOPRIGHT[0] > (c->xmax<<shift)) P_TOPRIGHT[0]= (c->xmax<<shift); //due to pmx no clip
3137     if(P_TOPRIGHT[1] > (c->ymax<<shift)) P_TOPRIGHT[1]= (c->ymax<<shift);
3138
3139     P_MEDIAN[0]= mid_pred(P_LEFT[0], P_TOP[0], P_TOPRIGHT[0]);
3140     P_MEDIAN[1]= mid_pred(P_LEFT[1], P_TOP[1], P_TOPRIGHT[1]);
3141
3142     if (!y) {
3143         c->pred_x= P_LEFT[0];
3144         c->pred_y= P_LEFT[1];
3145     } else {
3146         c->pred_x = P_MEDIAN[0];
3147         c->pred_y = P_MEDIAN[1];
3148     }
3149
3150     score= INT_MAX;
3151     best_ref= 0;
3152     for(ref=0; ref<s->ref_frames; ref++){
3153         init_ref(c, current_data, s->last_picture[ref].data, NULL, block_w*x, block_w*y, 0);
3154
3155         ref_score= ff_epzs_motion_search(&s->m, &ref_mx, &ref_my, P, 0, /*ref_index*/ 0, last_mv,
3156                                          (1<<16)>>shift, level-LOG2_MB_SIZE+4, block_w);
3157
3158         assert(ref_mx >= c->xmin);
3159         assert(ref_mx <= c->xmax);
3160         assert(ref_my >= c->ymin);
3161         assert(ref_my <= c->ymax);
3162
3163         ref_score= c->sub_motion_search(&s->m, &ref_mx, &ref_my, ref_score, 0, 0, level-LOG2_MB_SIZE+4, block_w);
3164         ref_score= ff_get_mb_score(&s->m, ref_mx, ref_my, 0, 0, level-LOG2_MB_SIZE+4, block_w, 0);
3165         ref_score+= 2*av_log2(2*ref)*c->penalty_factor;
3166         if(s->ref_mvs[ref]){
3167             s->ref_mvs[ref][index][0]= ref_mx;
3168             s->ref_mvs[ref][index][1]= ref_my;
3169             s->ref_scores[ref][index]= ref_score;
3170         }
3171         if(score > ref_score){
3172             score= ref_score;
3173             best_ref= ref;
3174             mx= ref_mx;
3175             my= ref_my;
3176         }
3177     }
3178     //FIXME if mb_cmp != SSE then intra cannot be compared currently and mb_penalty vs. lambda2
3179
3180   //  subpel search
3181     base_bits= get_rac_count(&s->c) - 8*(s->c.bytestream - s->c.bytestream_start);
3182     pc= s->c;
3183     pc.bytestream_start=
3184     pc.bytestream= p_buffer; //FIXME end/start? and at the other stoo
3185     memcpy(p_state, s->block_state, sizeof(s->block_state));
3186
3187     if(level!=s->block_max_depth)
3188         put_rac(&pc, &p_state[4 + s_context], 1);
3189     put_rac(&pc, &p_state[1 + left->type + top->type], 0);
3190     if(s->ref_frames > 1)
3191         put_symbol(&pc, &p_state[128 + 1024 + 32*ref_context], best_ref, 0);
3192     pred_mv(s, &pmx, &pmy, best_ref, left, top, tr);
3193     put_symbol(&pc, &p_state[128 + 32*(mx_context + 16*!!best_ref)], mx - pmx, 1);
3194     put_symbol(&pc, &p_state[128 + 32*(my_context + 16*!!best_ref)], my - pmy, 1);
3195     p_len= pc.bytestream - pc.bytestream_start;
3196     score += (s->lambda2*(get_rac_count(&pc)-base_bits))>>FF_LAMBDA_SHIFT;
3197
3198     block_s= block_w*block_w;
3199     sum = pix_sum(current_data[0], stride, block_w);
3200     l= (sum + block_s/2)/block_s;
3201     iscore = pix_norm1(current_data[0], stride, block_w) - 2*l*sum + l*l*block_s;
3202
3203     block_s= block_w*block_w>>2;
3204     sum = pix_sum(current_data[1], uvstride, block_w>>1);
3205     cb= (sum + block_s/2)/block_s;
3206 //    iscore += pix_norm1(&current_mb[1][0], uvstride, block_w>>1) - 2*cb*sum + cb*cb*block_s;
3207     sum = pix_sum(current_data[2], uvstride, block_w>>1);
3208     cr= (sum + block_s/2)/block_s;
3209 //    iscore += pix_norm1(&current_mb[2][0], uvstride, block_w>>1) - 2*cr*sum + cr*cr*block_s;
3210
3211     ic= s->c;
3212     ic.bytestream_start=
3213     ic.bytestream= i_buffer; //FIXME end/start? and at the other stoo
3214     memcpy(i_state, s->block_state, sizeof(s->block_state));
3215     if(level!=s->block_max_depth)
3216         put_rac(&ic, &i_state[4 + s_context], 1);
3217     put_rac(&ic, &i_state[1 + left->type + top->type], 1);
3218     put_symbol(&ic, &i_state[32],  l-pl , 1);
3219     put_symbol(&ic, &i_state[64], cb-pcb, 1);
3220     put_symbol(&ic, &i_state[96], cr-pcr, 1);
3221     i_len= ic.bytestream - ic.bytestream_start;
3222     iscore += (s->lambda2*(get_rac_count(&ic)-base_bits))>>FF_LAMBDA_SHIFT;
3223
3224 //    assert(score==256*256*256*64-1);
3225     assert(iscore < 255*255*256 + s->lambda2*10);
3226     assert(iscore >= 0);
3227     assert(l>=0 && l<=255);
3228     assert(pl>=0 && pl<=255);
3229
3230     if(level==0){
3231         int varc= iscore >> 8;
3232         int vard= score >> 8;
3233         if (vard <= 64 || vard < varc)
3234             c->scene_change_score+= ff_sqrt(vard) - ff_sqrt(varc);
3235         else
3236             c->scene_change_score+= s->m.qscale;
3237     }
3238
3239     if(level!=s->block_max_depth){
3240         put_rac(&s->c, &s->block_state[4 + s_context], 0);
3241         score2 = encode_q_branch(s, level+1, 2*x+0, 2*y+0);
3242         score2+= encode_q_branch(s, level+1, 2*x+1, 2*y+0);
3243         score2+= encode_q_branch(s, level+1, 2*x+0, 2*y+1);
3244         score2+= encode_q_branch(s, level+1, 2*x+1, 2*y+1);
3245         score2+= s->lambda2>>FF_LAMBDA_SHIFT; //FIXME exact split overhead
3246
3247         if(score2 < score && score2 < iscore)
3248             return score2;
3249     }
3250
3251     if(iscore < score){
3252         pred_mv(s, &pmx, &pmy, 0, left, top, tr);
3253         memcpy(pbbak, i_buffer, i_len);
3254         s->c= ic;
3255         s->c.bytestream_start= pbbak_start;
3256         s->c.bytestream= pbbak + i_len;
3257         set_blocks(s, level, x, y, l, cb, cr, pmx, pmy, 0, BLOCK_INTRA);
3258         memcpy(s->block_state, i_state, sizeof(s->block_state));
3259         return iscore;
3260     }else{
3261         memcpy(pbbak, p_buffer, p_len);
3262         s->c= pc;
3263         s->c.bytestream_start= pbbak_start;
3264         s->c.bytestream= pbbak + p_len;
3265         set_blocks(s, level, x, y, pl, pcb, pcr, mx, my, best_ref, 0);
3266         memcpy(s->block_state, p_state, sizeof(s->block_state));
3267         return score;
3268     }
3269 }
3270
3271 static void encode_q_branch2(SnowContext *s, int level, int x, int y){
3272     const int w= s->b_width  << s->block_max_depth;
3273     const int rem_depth= s->block_max_depth - level;
3274     const int index= (x + y*w) << rem_depth;
3275     int trx= (x+1)<<rem_depth;
3276     BlockNode *b= &s->block[index];
3277     const BlockNode *left  = x ? &s->block[index-1] : &null_block;
3278     const BlockNode *top   = y ? &s->block[index-w] : &null_block;
3279     const BlockNode *tl    = y && x ? &s->block[index-w-1] : left;
3280     const BlockNode *tr    = y && trx<w && ((x&1)==0 || level==0) ? &s->block[index-w+(1<<rem_depth)] : tl; //FIXME use lt
3281     int pl = left->color[0];
3282     int pcb= left->color[1];
3283     int pcr= left->color[2];
3284     int pmx, pmy;
3285     int ref_context= av_log2(2*left->ref) + av_log2(2*top->ref);
3286     int mx_context= av_log2(2*FFABS(left->mx - top->mx)) + 16*!!b->ref;
3287     int my_context= av_log2(2*FFABS(left->my - top->my)) + 16*!!b->ref;
3288     int s_context= 2*left->level + 2*top->level + tl->level + tr->level;
3289
3290     if(s->keyframe){
3291         set_blocks(s, level, x, y, pl, pcb, pcr, 0, 0, 0, BLOCK_INTRA);
3292         return;
3293     }
3294
3295     if(level!=s->block_max_depth){
3296         if(same_block(b,b+1) && same_block(b,b+w) && same_block(b,b+w+1)){
3297             put_rac(&s->c, &s->block_state[4 + s_context], 1);
3298         }else{
3299             put_rac(&s->c, &s->block_state[4 + s_context], 0);
3300             encode_q_branch2(s, level+1, 2*x+0, 2*y+0);
3301             encode_q_branch2(s, level+1, 2*x+1, 2*y+0);
3302             encode_q_branch2(s, level+1, 2*x+0, 2*y+1);
3303             encode_q_branch2(s, level+1, 2*x+1, 2*y+1);
3304             return;
3305         }
3306     }
3307     if(b->type & BLOCK_INTRA){
3308         pred_mv(s, &pmx, &pmy, 0, left, top, tr);
3309         put_rac(&s->c, &s->block_state[1 + (left->type&1) + (top->type&1)], 1);
3310         put_symbol(&s->c, &s->block_state[32], b->color[0]-pl , 1);
3311         put_symbol(&s->c, &s->block_state[64], b->color[1]-pcb, 1);
3312         put_symbol(&s->c, &s->block_state[96], b->color[2]-pcr, 1);
3313         set_blocks(s, level, x, y, b->color[0], b->color[1], b->color[2], pmx, pmy, 0, BLOCK_INTRA);
3314     }else{
3315         pred_mv(s, &pmx, &pmy, b->ref, left, top, tr);
3316         put_rac(&s->c, &s->block_state[1 + (left->type&1) + (top->type&1)], 0);
3317         if(s->ref_frames > 1)
3318             put_symbol(&s->c, &s->block_state[128 + 1024 + 32*ref_context], b->ref, 0);
3319         put_symbol(&s->c, &s->block_state[128 + 32*mx_context], b->mx - pmx, 1);
3320         put_symbol(&s->c, &s->block_state[128 + 32*my_context], b->my - pmy, 1);
3321         set_blocks(s, level, x, y, pl, pcb, pcr, b->mx, b->my, b->ref, 0);
3322     }
3323 }
3324
3325 static int get_dc(SnowContext *s, int mb_x, int mb_y, int plane_index){
3326     int i, x2, y2;
3327     Plane *p= &s->plane[plane_index];
3328     const int block_size = MB_SIZE >> s->block_max_depth;
3329     const int block_w    = plane_index ? block_size/2 : block_size;
3330     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];
3331     const int obmc_stride= plane_index ? block_size : 2*block_size;
3332     const int ref_stride= s->current_picture.linesize[plane_index];
3333     uint8_t *src= s-> input_picture.data[plane_index];
3334     IDWTELEM *dst= (IDWTELEM*)s->m.obmc_scratchpad + plane_index*block_size*block_size*4; //FIXME change to unsigned
3335     const int b_stride = s->b_width << s->block_max_depth;
3336     const int w= p->width;
3337     const int h= p->height;
3338     int index= mb_x + mb_y*b_stride;
3339     BlockNode *b= &s->block[index];
3340     BlockNode backup= *b;
3341     int ab=0;
3342     int aa=0;
3343
3344     b->type|= BLOCK_INTRA;
3345     b->color[plane_index]= 0;
3346     memset(dst, 0, obmc_stride*obmc_stride*sizeof(IDWTELEM));
3347
3348     for(i=0; i<4; i++){
3349         int mb_x2= mb_x + (i &1) - 1;
3350         int mb_y2= mb_y + (i>>1) - 1;
3351         int x= block_w*mb_x2 + block_w/2;
3352         int y= block_w*mb_y2 + block_w/2;
3353
3354         add_yblock(s, 0, NULL, dst + ((i&1)+(i>>1)*obmc_stride)*block_w, NULL, obmc,
3355                     x, y, block_w, block_w, w, h, obmc_stride, ref_stride, obmc_stride, mb_x2, mb_y2, 0, 0, plane_index);
3356
3357         for(y2= FFMAX(y, 0); y2<FFMIN(h, y+block_w); y2++){
3358             for(x2= FFMAX(x, 0); x2<FFMIN(w, x+block_w); x2++){
3359                 int index= x2-(block_w*mb_x - block_w/2) + (y2-(block_w*mb_y - block_w/2))*obmc_stride;
3360                 int obmc_v= obmc[index];
3361                 int d;
3362                 if(y<0) obmc_v += obmc[index + block_w*obmc_stride];
3363                 if(x<0) obmc_v += obmc[index + block_w];
3364                 if(y+block_w>h) obmc_v += obmc[index - block_w*obmc_stride];
3365                 if(x+block_w>w) obmc_v += obmc[index - block_w];
3366                 //FIXME precalculate this or simplify it somehow else
3367
3368                 d = -dst[index] + (1<<(FRAC_BITS-1));
3369                 dst[index] = d;
3370                 ab += (src[x2 + y2*ref_stride] - (d>>FRAC_BITS)) * obmc_v;
3371                 aa += obmc_v * obmc_v; //FIXME precalculate this
3372             }
3373         }
3374     }
3375     *b= backup;
3376
3377     return av_clip(((ab<<LOG2_OBMC_MAX) + aa/2)/aa, 0, 255); //FIXME we should not need clipping
3378 }
3379
3380 static inline int get_block_bits(SnowContext *s, int x, int y, int w){
3381     const int b_stride = s->b_width << s->block_max_depth;
3382     const int b_height = s->b_height<< s->block_max_depth;
3383     int index= x + y*b_stride;
3384     const BlockNode *b     = &s->block[index];
3385     const BlockNode *left  = x ? &s->block[index-1] : &null_block;
3386     const BlockNode *top   = y ? &s->block[index-b_stride] : &null_block;
3387     const BlockNode *tl    = y && x ? &s->block[index-b_stride-1] : left;
3388     const BlockNode *tr    = y && x+w<b_stride ? &s->block[index-b_stride+w] : tl;
3389     int dmx, dmy;
3390 //  int mx_context= av_log2(2*FFABS(left->mx - top->mx));
3391 //  int my_context= av_log2(2*FFABS(left->my - top->my));
3392
3393     if(x<0 || x>=b_stride || y>=b_height)
3394         return 0;
3395 /*
3396 1            0      0
3397 01X          1-2    1
3398 001XX        3-6    2-3
3399 0001XXX      7-14   4-7
3400 00001XXXX   15-30   8-15
3401 */
3402 //FIXME try accurate rate
3403 //FIXME intra and inter predictors if surrounding blocks are not the same type
3404     if(b->type & BLOCK_INTRA){
3405         return 3+2*( av_log2(2*FFABS(left->color[0] - b->color[0]))
3406                    + av_log2(2*FFABS(left->color[1] - b->color[1]))
3407                    + av_log2(2*FFABS(left->color[2] - b->color[2])));
3408     }else{
3409         pred_mv(s, &dmx, &dmy, b->ref, left, top, tr);
3410         dmx-= b->mx;
3411         dmy-= b->my;
3412         return 2*(1 + av_log2(2*FFABS(dmx)) //FIXME kill the 2* can be merged in lambda
3413                     + av_log2(2*FFABS(dmy))
3414                     + av_log2(2*b->ref));
3415     }
3416 }
3417
3418 static int get_block_rd(SnowContext *s, int mb_x, int mb_y, int plane_index, const uint8_t *obmc_edged){
3419     Plane *p= &s->plane[plane_index];
3420     const int block_size = MB_SIZE >> s->block_max_depth;
3421     const int block_w    = plane_index ? block_size/2 : block_size;
3422     const int obmc_stride= plane_index ? block_size : 2*block_size;
3423     const int ref_stride= s->current_picture.linesize[plane_index];
3424     uint8_t *dst= s->current_picture.data[plane_index];
3425     uint8_t *src= s->  input_picture.data[plane_index];
3426     IDWTELEM *pred= (IDWTELEM*)s->m.obmc_scratchpad + plane_index*block_size*block_size*4;
3427     uint8_t *cur = s->scratchbuf;
3428     uint8_t tmp[ref_stride*(2*MB_SIZE+HTAPS_MAX-1)];
3429     const int b_stride = s->b_width << s->block_max_depth;
3430     const int b_height = s->b_height<< s->block_max_depth;
3431     const int w= p->width;
3432     const int h= p->height;
3433     int distortion;
3434     int rate= 0;
3435     const int penalty_factor= get_penalty_factor(s->lambda, s->lambda2, s->avctx->me_cmp);
3436     int sx= block_w*mb_x - block_w/2;
3437     int sy= block_w*mb_y - block_w/2;
3438     int x0= FFMAX(0,-sx);
3439     int y0= FFMAX(0,-sy);
3440     int x1= FFMIN(block_w*2, w-sx);
3441     int y1= FFMIN(block_w*2, h-sy);
3442     int i,x,y;
3443
3444     pred_block(s, cur, tmp, ref_stride, sx, sy, block_w*2, block_w*2, &s->block[mb_x + mb_y*b_stride], plane_index, w, h);
3445
3446     for(y=y0; y<y1; y++){
3447         const uint8_t *obmc1= obmc_edged + y*obmc_stride;
3448         const IDWTELEM *pred1 = pred + y*obmc_stride;
3449         uint8_t *cur1 = cur + y*ref_stride;
3450         uint8_t *dst1 = dst + sx + (sy+y)*ref_stride;
3451         for(x=x0; x<x1; x++){
3452 #if FRAC_BITS >= LOG2_OBMC_MAX
3453             int v = (cur1[x] * obmc1[x]) << (FRAC_BITS - LOG2_OBMC_MAX);
3454 #else
3455             int v = (cur1[x] * obmc1[x] + (1<<(LOG2_OBMC_MAX - FRAC_BITS-1))) >> (LOG2_OBMC_MAX - FRAC_BITS);
3456 #endif
3457             v = (v + pred1[x]) >> FRAC_BITS;
3458             if(v&(~255)) v= ~(v>>31);
3459             dst1[x] = v;
3460         }
3461     }
3462
3463     /* copy the regions where obmc[] = (uint8_t)256 */
3464     if(LOG2_OBMC_MAX == 8
3465         && (mb_x == 0 || mb_x == b_stride-1)
3466         && (mb_y == 0 || mb_y == b_height-1)){
3467         if(mb_x == 0)
3468             x1 = block_w;
3469         else
3470             x0 = block_w;
3471         if(mb_y == 0)
3472             y1 = block_w;
3473         else
3474             y0 = block_w;
3475         for(y=y0; y<y1; y++)
3476             memcpy(dst + sx+x0 + (sy+y)*ref_stride, cur + x0 + y*ref_stride, x1-x0);
3477     }
3478
3479     if(block_w==16){
3480         /* FIXME rearrange dsputil to fit 32x32 cmp functions */
3481         /* FIXME check alignment of the cmp wavelet vs the encoding wavelet */
3482         /* FIXME cmps overlap but do not cover the wavelet's whole support.
3483          * So improving the score of one block is not strictly guaranteed
3484          * to improve the score of the whole frame, thus iterative motion
3485          * estimation does not always converge. */
3486         if(s->avctx->me_cmp == FF_CMP_W97)
3487             distortion = w97_32_c(&s->m, src + sx + sy*ref_stride, dst + sx + sy*ref_stride, ref_stride, 32);
3488         else if(s->avctx->me_cmp == FF_CMP_W53)
3489             distortion = w53_32_c(&s->m, src + sx + sy*ref_stride, dst + sx + sy*ref_stride, ref_stride, 32);
3490         else{
3491             distortion = 0;
3492             for(i=0; i<4; i++){
3493                 int off = sx+16*(i&1) + (sy+16*(i>>1))*ref_stride;
3494                 distortion += s->dsp.me_cmp[0](&s->m, src + off, dst + off, ref_stride, 16);
3495             }
3496         }
3497     }else{
3498         assert(block_w==8);
3499         distortion = s->dsp.me_cmp[0](&s->m, src + sx + sy*ref_stride, dst + sx + sy*ref_stride, ref_stride, block_w*2);
3500     }
3501
3502     if(plane_index==0){
3503         for(i=0; i<4; i++){
3504 /* ..RRr
3505  * .RXx.
3506  * rxx..
3507  */
3508             rate += get_block_bits(s, mb_x + (i&1) - (i>>1), mb_y + (i>>1), 1);
3509         }
3510         if(mb_x == b_stride-2)
3511             rate += get_block_bits(s, mb_x + 1, mb_y + 1, 1);
3512     }
3513     return distortion + rate*penalty_factor;
3514 }
3515
3516 static int get_4block_rd(SnowContext *s, int mb_x, int mb_y, int plane_index){
3517     int i, y2;
3518     Plane *p= &s->plane[plane_index];
3519     const int block_size = MB_SIZE >> s->block_max_depth;
3520     const int block_w    = plane_index ? block_size/2 : block_size;
3521     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];
3522     const int obmc_stride= plane_index ? block_size : 2*block_size;
3523     const int ref_stride= s->current_picture.linesize[plane_index];
3524     uint8_t *dst= s->current_picture.data[plane_index];
3525     uint8_t *src= s-> input_picture.data[plane_index];
3526     //FIXME zero_dst is const but add_yblock changes dst if add is 0 (this is never the case for dst=zero_dst
3527     // const has only been removed from zero_dst to suppress a warning
3528     static IDWTELEM zero_dst[4096]; //FIXME
3529     const int b_stride = s->b_width << s->block_max_depth;
3530     const int w= p->width;
3531     const int h= p->height;
3532     int distortion= 0;
3533     int rate= 0;
3534     const int penalty_factor= get_penalty_factor(s->lambda, s->lambda2, s->avctx->me_cmp);
3535
3536     for(i=0; i<9; i++){
3537         int mb_x2= mb_x + (i%3) - 1;
3538         int mb_y2= mb_y + (i/3) - 1;
3539         int x= block_w*mb_x2 + block_w/2;
3540         int y= block_w*mb_y2 + block_w/2;
3541
3542         add_yblock(s, 0, NULL, zero_dst, dst, obmc,
3543                    x, y, block_w, block_w, w, h, /*dst_stride*/0, ref_stride, obmc_stride, mb_x2, mb_y2, 1, 1, plane_index);
3544
3545         //FIXME find a cleaner/simpler way to skip the outside stuff
3546         for(y2= y; y2<0; y2++)
3547             memcpy(dst + x + y2*ref_stride, src + x + y2*ref_stride, block_w);
3548         for(y2= h; y2<y+block_w; y2++)
3549             memcpy(dst + x + y2*ref_stride, src + x + y2*ref_stride, block_w);
3550         if(x<0){
3551             for(y2= y; y2<y+block_w; y2++)
3552                 memcpy(dst + x + y2*ref_stride, src + x + y2*ref_stride, -x);
3553         }
3554         if(x+block_w > w){
3555             for(y2= y; y2<y+block_w; y2++)
3556                 memcpy(dst + w + y2*ref_stride, src + w + y2*ref_stride, x+block_w - w);
3557         }
3558
3559         assert(block_w== 8 || block_w==16);
3560         distortion += s->dsp.me_cmp[block_w==8](&s->m, src + x + y*ref_stride, dst + x + y*ref_stride, ref_stride, block_w);
3561     }
3562
3563     if(plane_index==0){
3564         BlockNode *b= &s->block[mb_x+mb_y*b_stride];
3565         int merged= same_block(b,b+1) && same_block(b,b+b_stride) && same_block(b,b+b_stride+1);
3566
3567 /* ..RRRr
3568  * .RXXx.
3569  * .RXXx.
3570  * rxxx.
3571  */
3572         if(merged)
3573             rate = get_block_bits(s, mb_x, mb_y, 2);
3574         for(i=merged?4:0; i<9; i++){
3575             static const int dxy[9][2] = {{0,0},{1,0},{0,1},{1,1},{2,0},{2,1},{-1,2},{0,2},{1,2}};
3576             rate += get_block_bits(s, mb_x + dxy[i][0], mb_y + dxy[i][1], 1);
3577         }
3578     }
3579     return distortion + rate*penalty_factor;
3580 }
3581
3582 static void ff_spatial_idwt_init(DWTCompose *cs, IDWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count){
3583     int level;
3584     for(level=decomposition_count-1; level>=0; level--){
3585         switch(type){
3586         case DWT_97: spatial_compose97i_init(cs+level, buffer, height>>level, stride<<level); break;
3587         case DWT_53: spatial_compose53i_init(cs+level, buffer, height>>level, stride<<level); break;
3588         }
3589     }
3590 }
3591
3592 static void ff_spatial_idwt_slice(DWTCompose *cs, IDWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count, int y){
3593     const int support = type==1 ? 3 : 5;
3594     int level;
3595     if(type==2) return;
3596
3597     for(level=decomposition_count-1; level>=0; level--){
3598         while(cs[level].y <= FFMIN((y>>level)+support, height>>level)){
3599             switch(type){
3600             case DWT_97: spatial_compose97i_dy(cs+level, buffer, width>>level, height>>level, stride<<level);
3601                 break;