Optimize horizontal_compose53i() by inlining inv_lift() and optimizing the
[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) vertical_compose53iL0(b1, b2, b3, width);
1074         if(y+0<(unsigned)height) vertical_compose53iH0(b0, b1, b2, width);
1075
1076         if(y-1<(unsigned)height) horizontal_compose53i(b0, width);
1077         if(y+0<(unsigned)height) horizontal_compose53i(b1, width);
1078
1079     cs->b0 = b2;
1080     cs->b1 = b3;
1081     cs->y += 2;
1082 }
1083
1084 static void spatial_compose53i_dy(DWTCompose *cs, IDWTELEM *buffer, int width, int height, int stride){
1085     int y= cs->y;
1086     IDWTELEM *b0= cs->b0;
1087     IDWTELEM *b1= cs->b1;
1088     IDWTELEM *b2= buffer + mirror(y+1, height-1)*stride;
1089     IDWTELEM *b3= buffer + mirror(y+2, height-1)*stride;
1090
1091         if(y+1<(unsigned)height) vertical_compose53iL0(b1, b2, b3, width);
1092         if(y+0<(unsigned)height) vertical_compose53iH0(b0, b1, b2, width);
1093
1094         if(y-1<(unsigned)height) horizontal_compose53i(b0, width);
1095         if(y+0<(unsigned)height) horizontal_compose53i(b1, width);
1096
1097     cs->b0 = b2;
1098     cs->b1 = b3;
1099     cs->y += 2;
1100 }
1101
1102 static void av_unused spatial_compose53i(IDWTELEM *buffer, int width, int height, int stride){
1103     DWTCompose cs;
1104     spatial_compose53i_init(&cs, buffer, height, stride);
1105     while(cs.y <= height)
1106         spatial_compose53i_dy(&cs, buffer, width, height, stride);
1107 }
1108
1109
1110 void ff_snow_horizontal_compose97i(IDWTELEM *b, int width){
1111     IDWTELEM temp[width];
1112     const int w2= (width+1)>>1;
1113
1114     inv_lift (temp   , b      , b   +w2, 1, 1, 1, width,  W_DM, W_DO, W_DS, 0, 1);
1115     inv_lift (temp+w2, b   +w2, temp   , 1, 1, 1, width,  W_CM, W_CO, W_CS, 1, 1);
1116     inv_liftS(b      , temp   , temp+w2, 2, 1, 1, width,  W_BM, W_BO, W_BS, 0, 1);
1117     inv_lift (b+1    , temp+w2, b      , 2, 1, 2, width,  W_AM, W_AO, W_AS, 1, 0);
1118 }
1119
1120 static void vertical_compose97iH0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
1121     int i;
1122
1123     for(i=0; i<width; i++){
1124         b1[i] += (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
1125     }
1126 }
1127
1128 static void vertical_compose97iH1(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
1129     int i;
1130
1131     for(i=0; i<width; i++){
1132         b1[i] -= (W_CM*(b0[i] + b2[i])+W_CO)>>W_CS;
1133     }
1134 }
1135
1136 static void vertical_compose97iL0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
1137     int i;
1138
1139     for(i=0; i<width; i++){
1140 #ifdef liftS
1141         b1[i] += (W_BM*(b0[i] + b2[i])+W_BO)>>W_BS;
1142 #else
1143         b1[i] += (W_BM*(b0[i] + b2[i])+4*b1[i]+W_BO)>>W_BS;
1144 #endif
1145     }
1146 }
1147
1148 static void vertical_compose97iL1(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
1149     int i;
1150
1151     for(i=0; i<width; i++){
1152         b1[i] -= (W_DM*(b0[i] + b2[i])+W_DO)>>W_DS;
1153     }
1154 }
1155
1156 void ff_snow_vertical_compose97i(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, IDWTELEM *b3, IDWTELEM *b4, IDWTELEM *b5, int width){
1157     int i;
1158
1159     for(i=0; i<width; i++){
1160         b4[i] -= (W_DM*(b3[i] + b5[i])+W_DO)>>W_DS;
1161         b3[i] -= (W_CM*(b2[i] + b4[i])+W_CO)>>W_CS;
1162 #ifdef liftS
1163         b2[i] += (W_BM*(b1[i] + b3[i])+W_BO)>>W_BS;
1164 #else
1165         b2[i] += (W_BM*(b1[i] + b3[i])+4*b2[i]+W_BO)>>W_BS;
1166 #endif
1167         b1[i] += (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
1168     }
1169 }
1170
1171 static void spatial_compose97i_buffered_init(DWTCompose *cs, slice_buffer * sb, int height, int stride_line){
1172     cs->b0 = slice_buffer_get_line(sb, mirror(-3-1, height-1) * stride_line);
1173     cs->b1 = slice_buffer_get_line(sb, mirror(-3  , height-1) * stride_line);
1174     cs->b2 = slice_buffer_get_line(sb, mirror(-3+1, height-1) * stride_line);
1175     cs->b3 = slice_buffer_get_line(sb, mirror(-3+2, height-1) * stride_line);
1176     cs->y = -3;
1177 }
1178
1179 static void spatial_compose97i_init(DWTCompose *cs, IDWTELEM *buffer, int height, int stride){
1180     cs->b0 = buffer + mirror(-3-1, height-1)*stride;
1181     cs->b1 = buffer + mirror(-3  , height-1)*stride;
1182     cs->b2 = buffer + mirror(-3+1, height-1)*stride;
1183     cs->b3 = buffer + mirror(-3+2, height-1)*stride;
1184     cs->y = -3;
1185 }
1186
1187 static void spatial_compose97i_dy_buffered(DSPContext *dsp, DWTCompose *cs, slice_buffer * sb, int width, int height, int stride_line){
1188     int y = cs->y;
1189
1190     IDWTELEM *b0= cs->b0;
1191     IDWTELEM *b1= cs->b1;
1192     IDWTELEM *b2= cs->b2;
1193     IDWTELEM *b3= cs->b3;
1194     IDWTELEM *b4= slice_buffer_get_line(sb, mirror(y + 3, height - 1) * stride_line);
1195     IDWTELEM *b5= slice_buffer_get_line(sb, mirror(y + 4, height - 1) * stride_line);
1196
1197     if(y>0 && y+4<height){
1198         dsp->vertical_compose97i(b0, b1, b2, b3, b4, b5, width);
1199     }else{
1200         if(y+3<(unsigned)height) vertical_compose97iL1(b3, b4, b5, width);
1201         if(y+2<(unsigned)height) vertical_compose97iH1(b2, b3, b4, width);
1202         if(y+1<(unsigned)height) vertical_compose97iL0(b1, b2, b3, width);
1203         if(y+0<(unsigned)height) vertical_compose97iH0(b0, b1, b2, width);
1204     }
1205
1206     if(y-1<(unsigned)height) dsp->horizontal_compose97i(b0, width);
1207     if(y+0<(unsigned)height) dsp->horizontal_compose97i(b1, width);
1208
1209     cs->b0=b2;
1210     cs->b1=b3;
1211     cs->b2=b4;
1212     cs->b3=b5;
1213     cs->y += 2;
1214 }
1215
1216 static void spatial_compose97i_dy(DWTCompose *cs, IDWTELEM *buffer, int width, int height, int stride){
1217     int y = cs->y;
1218     IDWTELEM *b0= cs->b0;
1219     IDWTELEM *b1= cs->b1;
1220     IDWTELEM *b2= cs->b2;
1221     IDWTELEM *b3= cs->b3;
1222     IDWTELEM *b4= buffer + mirror(y+3, height-1)*stride;
1223     IDWTELEM *b5= buffer + mirror(y+4, height-1)*stride;
1224
1225     if(y+3<(unsigned)height) vertical_compose97iL1(b3, b4, b5, width);
1226     if(y+2<(unsigned)height) vertical_compose97iH1(b2, b3, b4, width);
1227     if(y+1<(unsigned)height) vertical_compose97iL0(b1, b2, b3, width);
1228     if(y+0<(unsigned)height) vertical_compose97iH0(b0, b1, b2, width);
1229
1230     if(y-1<(unsigned)height) ff_snow_horizontal_compose97i(b0, width);
1231     if(y+0<(unsigned)height) ff_snow_horizontal_compose97i(b1, width);
1232
1233     cs->b0=b2;
1234     cs->b1=b3;
1235     cs->b2=b4;
1236     cs->b3=b5;
1237     cs->y += 2;
1238 }
1239
1240 static void av_unused spatial_compose97i(IDWTELEM *buffer, int width, int height, int stride){
1241     DWTCompose cs;
1242     spatial_compose97i_init(&cs, buffer, height, stride);
1243     while(cs.y <= height)
1244         spatial_compose97i_dy(&cs, buffer, width, height, stride);
1245 }
1246
1247 static void ff_spatial_idwt_buffered_init(DWTCompose *cs, slice_buffer * sb, int width, int height, int stride_line, int type, int decomposition_count){
1248     int level;
1249     for(level=decomposition_count-1; level>=0; level--){
1250         switch(type){
1251         case DWT_97: spatial_compose97i_buffered_init(cs+level, sb, height>>level, stride_line<<level); break;
1252         case DWT_53: spatial_compose53i_buffered_init(cs+level, sb, height>>level, stride_line<<level); break;
1253         }
1254     }
1255 }
1256
1257 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){
1258     const int support = type==1 ? 3 : 5;
1259     int level;
1260     if(type==2) return;
1261
1262     for(level=decomposition_count-1; level>=0; level--){
1263         while(cs[level].y <= FFMIN((y>>level)+support, height>>level)){
1264             switch(type){
1265             case DWT_97: spatial_compose97i_dy_buffered(dsp, cs+level, slice_buf, width>>level, height>>level, stride_line<<level);
1266                 break;
1267             case DWT_53: spatial_compose53i_dy_buffered(cs+level, slice_buf, width>>level, height>>level, stride_line<<level);
1268                 break;
1269             }
1270         }
1271     }
1272 }
1273
1274 static inline void unpack_coeffs(SnowContext *s, SubBand *b, SubBand * parent, int orientation){
1275     const int w= b->width;
1276     const int h= b->height;
1277     int x,y;
1278
1279     int run, runs;
1280     x_and_coeff *xc= b->x_coeff;
1281     x_and_coeff *prev_xc= NULL;
1282     x_and_coeff *prev2_xc= xc;
1283     x_and_coeff *parent_xc= parent ? parent->x_coeff : NULL;
1284     x_and_coeff *prev_parent_xc= parent_xc;
1285
1286     runs= get_symbol2(&s->c, b->state[30], 0);
1287     if(runs-- > 0) run= get_symbol2(&s->c, b->state[1], 3);
1288     else           run= INT_MAX;
1289
1290     for(y=0; y<h; y++){
1291         int v=0;
1292         int lt=0, t=0, rt=0;
1293
1294         if(y && prev_xc->x == 0){
1295             rt= prev_xc->coeff;
1296         }
1297         for(x=0; x<w; x++){
1298             int p=0;
1299             const int l= v;
1300
1301             lt= t; t= rt;
1302
1303             if(y){
1304                 if(prev_xc->x <= x)
1305                     prev_xc++;
1306                 if(prev_xc->x == x + 1)
1307                     rt= prev_xc->coeff;
1308                 else
1309                     rt=0;
1310             }
1311             if(parent_xc){
1312                 if(x>>1 > parent_xc->x){
1313                     parent_xc++;
1314                 }
1315                 if(x>>1 == parent_xc->x){
1316                     p= parent_xc->coeff;
1317                 }
1318             }
1319             if(/*ll|*/l|lt|t|rt|p){
1320                 int context= av_log2(/*FFABS(ll) + */3*(l>>1) + (lt>>1) + (t&~1) + (rt>>1) + (p>>1));
1321
1322                 v=get_rac(&s->c, &b->state[0][context]);
1323                 if(v){
1324                     v= 2*(get_symbol2(&s->c, b->state[context + 2], context-4) + 1);
1325                     v+=get_rac(&s->c, &b->state[0][16 + 1 + 3 + quant3bA[l&0xFF] + 3*quant3bA[t&0xFF]]);
1326
1327                     xc->x=x;
1328                     (xc++)->coeff= v;
1329                 }
1330             }else{
1331                 if(!run){
1332                     if(runs-- > 0) run= get_symbol2(&s->c, b->state[1], 3);
1333                     else           run= INT_MAX;
1334                     v= 2*(get_symbol2(&s->c, b->state[0 + 2], 0-4) + 1);
1335                     v+=get_rac(&s->c, &b->state[0][16 + 1 + 3]);
1336
1337                     xc->x=x;
1338                     (xc++)->coeff= v;
1339                 }else{
1340                     int max_run;
1341                     run--;
1342                     v=0;
1343
1344                     if(y) max_run= FFMIN(run, prev_xc->x - x - 2);
1345                     else  max_run= FFMIN(run, w-x-1);
1346                     if(parent_xc)
1347                         max_run= FFMIN(max_run, 2*parent_xc->x - x - 1);
1348                     x+= max_run;
1349                     run-= max_run;
1350                 }
1351             }
1352         }
1353         (xc++)->x= w+1; //end marker
1354         prev_xc= prev2_xc;
1355         prev2_xc= xc;
1356
1357         if(parent_xc){
1358             if(y&1){
1359                 while(parent_xc->x != parent->width+1)
1360                     parent_xc++;
1361                 parent_xc++;
1362                 prev_parent_xc= parent_xc;
1363             }else{
1364                 parent_xc= prev_parent_xc;
1365             }
1366         }
1367     }
1368
1369     (xc++)->x= w+1; //end marker
1370 }
1371
1372 static inline void decode_subband_slice_buffered(SnowContext *s, SubBand *b, slice_buffer * sb, int start_y, int h, int save_state[1]){
1373     const int w= b->width;
1374     int y;
1375     const int qlog= av_clip(s->qlog + b->qlog, 0, QROOT*16);
1376     int qmul= qexp[qlog&(QROOT-1)]<<(qlog>>QSHIFT);
1377     int qadd= (s->qbias*qmul)>>QBIAS_SHIFT;
1378     int new_index = 0;
1379
1380     if(b->ibuf == s->spatial_idwt_buffer || s->qlog == LOSSLESS_QLOG){
1381         qadd= 0;
1382         qmul= 1<<QEXPSHIFT;
1383     }
1384
1385     /* If we are on the second or later slice, restore our index. */
1386     if (start_y != 0)
1387         new_index = save_state[0];
1388
1389
1390     for(y=start_y; y<h; y++){
1391         int x = 0;
1392         int v;
1393         IDWTELEM * line = slice_buffer_get_line(sb, y * b->stride_line + b->buf_y_offset) + b->buf_x_offset;
1394         memset(line, 0, b->width*sizeof(IDWTELEM));
1395         v = b->x_coeff[new_index].coeff;
1396         x = b->x_coeff[new_index++].x;
1397         while(x < w){
1398             register int t= ( (v>>1)*qmul + qadd)>>QEXPSHIFT;
1399             register int u= -(v&1);
1400             line[x] = (t^u) - u;
1401
1402             v = b->x_coeff[new_index].coeff;
1403             x = b->x_coeff[new_index++].x;
1404         }
1405     }
1406
1407     /* Save our variables for the next slice. */
1408     save_state[0] = new_index;
1409
1410     return;
1411 }
1412
1413 static void reset_contexts(SnowContext *s){ //FIXME better initial contexts
1414     int plane_index, level, orientation;
1415
1416     for(plane_index=0; plane_index<3; plane_index++){
1417         for(level=0; level<MAX_DECOMPOSITIONS; level++){
1418             for(orientation=level ? 1:0; orientation<4; orientation++){
1419                 memset(s->plane[plane_index].band[level][orientation].state, MID_STATE, sizeof(s->plane[plane_index].band[level][orientation].state));
1420             }
1421         }
1422     }
1423     memset(s->header_state, MID_STATE, sizeof(s->header_state));
1424     memset(s->block_state, MID_STATE, sizeof(s->block_state));
1425 }
1426
1427 static int alloc_blocks(SnowContext *s){
1428     int w= -((-s->avctx->width )>>LOG2_MB_SIZE);
1429     int h= -((-s->avctx->height)>>LOG2_MB_SIZE);
1430
1431     s->b_width = w;
1432     s->b_height= h;
1433
1434     av_free(s->block);
1435     s->block= av_mallocz(w * h * sizeof(BlockNode) << (s->block_max_depth*2));
1436     return 0;
1437 }
1438
1439 static inline void copy_rac_state(RangeCoder *d, RangeCoder *s){
1440     uint8_t *bytestream= d->bytestream;
1441     uint8_t *bytestream_start= d->bytestream_start;
1442     *d= *s;
1443     d->bytestream= bytestream;
1444     d->bytestream_start= bytestream_start;
1445 }
1446
1447 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){
1448     const int w= s->b_width << s->block_max_depth;
1449     const int rem_depth= s->block_max_depth - level;
1450     const int index= (x + y*w) << rem_depth;
1451     const int block_w= 1<<rem_depth;
1452     BlockNode block;
1453     int i,j;
1454
1455     block.color[0]= l;
1456     block.color[1]= cb;
1457     block.color[2]= cr;
1458     block.mx= mx;
1459     block.my= my;
1460     block.ref= ref;
1461     block.type= type;
1462     block.level= level;
1463
1464     for(j=0; j<block_w; j++){
1465         for(i=0; i<block_w; i++){
1466             s->block[index + i + j*w]= block;
1467         }
1468     }
1469 }
1470
1471 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){
1472     const int offset[3]= {
1473           y*c->  stride + x,
1474         ((y*c->uvstride + x)>>1),
1475         ((y*c->uvstride + x)>>1),
1476     };
1477     int i;
1478     for(i=0; i<3; i++){
1479         c->src[0][i]= src [i];
1480         c->ref[0][i]= ref [i] + offset[i];
1481     }
1482     assert(!ref_index);
1483 }
1484
1485 static inline void pred_mv(SnowContext *s, int *mx, int *my, int ref,
1486                            const BlockNode *left, const BlockNode *top, const BlockNode *tr){
1487     if(s->ref_frames == 1){
1488         *mx = mid_pred(left->mx, top->mx, tr->mx);
1489         *my = mid_pred(left->my, top->my, tr->my);
1490     }else{
1491         const int *scale = scale_mv_ref[ref];
1492         *mx = mid_pred((left->mx * scale[left->ref] + 128) >>8,
1493                        (top ->mx * scale[top ->ref] + 128) >>8,
1494                        (tr  ->mx * scale[tr  ->ref] + 128) >>8);
1495         *my = mid_pred((left->my * scale[left->ref] + 128) >>8,
1496                        (top ->my * scale[top ->ref] + 128) >>8,
1497                        (tr  ->my * scale[tr  ->ref] + 128) >>8);
1498     }
1499 }
1500
1501 static av_always_inline int same_block(BlockNode *a, BlockNode *b){
1502     if((a->type&BLOCK_INTRA) && (b->type&BLOCK_INTRA)){
1503         return !((a->color[0] - b->color[0]) | (a->color[1] - b->color[1]) | (a->color[2] - b->color[2]));
1504     }else{
1505         return !((a->mx - b->mx) | (a->my - b->my) | (a->ref - b->ref) | ((a->type ^ b->type)&BLOCK_INTRA));
1506     }
1507 }
1508
1509 static void decode_q_branch(SnowContext *s, int level, int x, int y){
1510     const int w= s->b_width << s->block_max_depth;
1511     const int rem_depth= s->block_max_depth - level;
1512     const int index= (x + y*w) << rem_depth;
1513     int trx= (x+1)<<rem_depth;
1514     const BlockNode *left  = x ? &s->block[index-1] : &null_block;
1515     const BlockNode *top   = y ? &s->block[index-w] : &null_block;
1516     const BlockNode *tl    = y && x ? &s->block[index-w-1] : left;
1517     const BlockNode *tr    = y && trx<w && ((x&1)==0 || level==0) ? &s->block[index-w+(1<<rem_depth)] : tl; //FIXME use lt
1518     int s_context= 2*left->level + 2*top->level + tl->level + tr->level;
1519
1520     if(s->keyframe){
1521         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);
1522         return;
1523     }
1524
1525     if(level==s->block_max_depth || get_rac(&s->c, &s->block_state[4 + s_context])){
1526         int type, mx, my;
1527         int l = left->color[0];
1528         int cb= left->color[1];
1529         int cr= left->color[2];
1530         int ref = 0;
1531         int ref_context= av_log2(2*left->ref) + av_log2(2*top->ref);
1532         int mx_context= av_log2(2*FFABS(left->mx - top->mx)) + 0*av_log2(2*FFABS(tr->mx - top->mx));
1533         int my_context= av_log2(2*FFABS(left->my - top->my)) + 0*av_log2(2*FFABS(tr->my - top->my));
1534
1535         type= get_rac(&s->c, &s->block_state[1 + left->type + top->type]) ? BLOCK_INTRA : 0;
1536
1537         if(type){
1538             pred_mv(s, &mx, &my, 0, left, top, tr);
1539             l += get_symbol(&s->c, &s->block_state[32], 1);
1540             cb+= get_symbol(&s->c, &s->block_state[64], 1);
1541             cr+= get_symbol(&s->c, &s->block_state[96], 1);
1542         }else{
1543             if(s->ref_frames > 1)
1544                 ref= get_symbol(&s->c, &s->block_state[128 + 1024 + 32*ref_context], 0);
1545             pred_mv(s, &mx, &my, ref, left, top, tr);
1546             mx+= get_symbol(&s->c, &s->block_state[128 + 32*(mx_context + 16*!!ref)], 1);
1547             my+= get_symbol(&s->c, &s->block_state[128 + 32*(my_context + 16*!!ref)], 1);
1548         }
1549         set_blocks(s, level, x, y, l, cb, cr, mx, my, ref, type);
1550     }else{
1551         decode_q_branch(s, level+1, 2*x+0, 2*y+0);
1552         decode_q_branch(s, level+1, 2*x+1, 2*y+0);
1553         decode_q_branch(s, level+1, 2*x+0, 2*y+1);
1554         decode_q_branch(s, level+1, 2*x+1, 2*y+1);
1555     }
1556 }
1557
1558 static void decode_blocks(SnowContext *s){
1559     int x, y;
1560     int w= s->b_width;
1561     int h= s->b_height;
1562
1563     for(y=0; y<h; y++){
1564         for(x=0; x<w; x++){
1565             decode_q_branch(s, 0, x, y);
1566         }
1567     }
1568 }
1569
1570 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){
1571     static const uint8_t weight[64]={
1572     8,7,6,5,4,3,2,1,
1573     7,7,0,0,0,0,0,1,
1574     6,0,6,0,0,0,2,0,
1575     5,0,0,5,0,3,0,0,
1576     4,0,0,0,4,0,0,0,
1577     3,0,0,5,0,3,0,0,
1578     2,0,6,0,0,0,2,0,
1579     1,7,0,0,0,0,0,1,
1580     };
1581
1582     static const uint8_t brane[256]={
1583     0x00,0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x11,0x12,0x12,0x12,0x12,0x12,0x12,0x12,
1584     0x04,0x05,0xcc,0xcc,0xcc,0xcc,0xcc,0x41,0x15,0x16,0xcc,0xcc,0xcc,0xcc,0xcc,0x52,
1585     0x04,0xcc,0x05,0xcc,0xcc,0xcc,0x41,0xcc,0x15,0xcc,0x16,0xcc,0xcc,0xcc,0x52,0xcc,
1586     0x04,0xcc,0xcc,0x05,0xcc,0x41,0xcc,0xcc,0x15,0xcc,0xcc,0x16,0xcc,0x52,0xcc,0xcc,
1587     0x04,0xcc,0xcc,0xcc,0x41,0xcc,0xcc,0xcc,0x15,0xcc,0xcc,0xcc,0x16,0xcc,0xcc,0xcc,
1588     0x04,0xcc,0xcc,0x41,0xcc,0x05,0xcc,0xcc,0x15,0xcc,0xcc,0x52,0xcc,0x16,0xcc,0xcc,
1589     0x04,0xcc,0x41,0xcc,0xcc,0xcc,0x05,0xcc,0x15,0xcc,0x52,0xcc,0xcc,0xcc,0x16,0xcc,
1590     0x04,0x41,0xcc,0xcc,0xcc,0xcc,0xcc,0x05,0x15,0x52,0xcc,0xcc,0xcc,0xcc,0xcc,0x16,
1591     0x44,0x45,0x45,0x45,0x45,0x45,0x45,0x45,0x55,0x56,0x56,0x56,0x56,0x56,0x56,0x56,
1592     0x48,0x49,0xcc,0xcc,0xcc,0xcc,0xcc,0x85,0x59,0x5A,0xcc,0xcc,0xcc,0xcc,0xcc,0x96,
1593     0x48,0xcc,0x49,0xcc,0xcc,0xcc,0x85,0xcc,0x59,0xcc,0x5A,0xcc,0xcc,0xcc,0x96,0xcc,
1594     0x48,0xcc,0xcc,0x49,0xcc,0x85,0xcc,0xcc,0x59,0xcc,0xcc,0x5A,0xcc,0x96,0xcc,0xcc,
1595     0x48,0xcc,0xcc,0xcc,0x49,0xcc,0xcc,0xcc,0x59,0xcc,0xcc,0xcc,0x96,0xcc,0xcc,0xcc,
1596     0x48,0xcc,0xcc,0x85,0xcc,0x49,0xcc,0xcc,0x59,0xcc,0xcc,0x96,0xcc,0x5A,0xcc,0xcc,
1597     0x48,0xcc,0x85,0xcc,0xcc,0xcc,0x49,0xcc,0x59,0xcc,0x96,0xcc,0xcc,0xcc,0x5A,0xcc,
1598     0x48,0x85,0xcc,0xcc,0xcc,0xcc,0xcc,0x49,0x59,0x96,0xcc,0xcc,0xcc,0xcc,0xcc,0x5A,
1599     };
1600
1601     static const uint8_t needs[16]={
1602     0,1,0,0,
1603     2,4,2,0,
1604     0,1,0,0,
1605     15
1606     };
1607
1608     int x, y, b, r, l;
1609     int16_t tmpIt   [64*(32+HTAPS_MAX)];
1610     uint8_t tmp2t[3][stride*(32+HTAPS_MAX)];
1611     int16_t *tmpI= tmpIt;
1612     uint8_t *tmp2= tmp2t[0];
1613     const uint8_t *hpel[11];
1614     assert(dx<16 && dy<16);
1615     r= brane[dx + 16*dy]&15;
1616     l= brane[dx + 16*dy]>>4;
1617
1618     b= needs[l] | needs[r];
1619     if(p && !p->diag_mc)
1620         b= 15;
1621
1622     if(b&5){
1623         for(y=0; y < b_h+HTAPS_MAX-1; y++){
1624             for(x=0; x < b_w; x++){
1625                 int a_1=src[x + HTAPS_MAX/2-4];
1626                 int a0= src[x + HTAPS_MAX/2-3];
1627                 int a1= src[x + HTAPS_MAX/2-2];
1628                 int a2= src[x + HTAPS_MAX/2-1];
1629                 int a3= src[x + HTAPS_MAX/2+0];
1630                 int a4= src[x + HTAPS_MAX/2+1];
1631                 int a5= src[x + HTAPS_MAX/2+2];
1632                 int a6= src[x + HTAPS_MAX/2+3];
1633                 int am=0;
1634                 if(!p || p->fast_mc){
1635                     am= 20*(a2+a3) - 5*(a1+a4) + (a0+a5);
1636                     tmpI[x]= am;
1637                     am= (am+16)>>5;
1638                 }else{
1639                     am= p->hcoeff[0]*(a2+a3) + p->hcoeff[1]*(a1+a4) + p->hcoeff[2]*(a0+a5) + p->hcoeff[3]*(a_1+a6);
1640                     tmpI[x]= am;
1641                     am= (am+32)>>6;
1642                 }
1643
1644                 if(am&(~255)) am= ~(am>>31);
1645                 tmp2[x]= am;
1646             }
1647             tmpI+= 64;
1648             tmp2+= stride;
1649             src += stride;
1650         }
1651         src -= stride*y;
1652     }
1653     src += HTAPS_MAX/2 - 1;
1654     tmp2= tmp2t[1];
1655
1656     if(b&2){
1657         for(y=0; y < b_h; y++){
1658             for(x=0; x < b_w+1; x++){
1659                 int a_1=src[x + (HTAPS_MAX/2-4)*stride];
1660                 int a0= src[x + (HTAPS_MAX/2-3)*stride];
1661                 int a1= src[x + (HTAPS_MAX/2-2)*stride];
1662                 int a2= src[x + (HTAPS_MAX/2-1)*stride];
1663                 int a3= src[x + (HTAPS_MAX/2+0)*stride];
1664                 int a4= src[x + (HTAPS_MAX/2+1)*stride];
1665                 int a5= src[x + (HTAPS_MAX/2+2)*stride];
1666                 int a6= src[x + (HTAPS_MAX/2+3)*stride];
1667                 int am=0;
1668                 if(!p || p->fast_mc)
1669                     am= (20*(a2+a3) - 5*(a1+a4) + (a0+a5) + 16)>>5;
1670                 else
1671                     am= (p->hcoeff[0]*(a2+a3) + p->hcoeff[1]*(a1+a4) + p->hcoeff[2]*(a0+a5) + p->hcoeff[3]*(a_1+a6) + 32)>>6;
1672
1673                 if(am&(~255)) am= ~(am>>31);
1674                 tmp2[x]= am;
1675             }
1676             src += stride;
1677             tmp2+= stride;
1678         }
1679         src -= stride*y;
1680     }
1681     src += stride*(HTAPS_MAX/2 - 1);
1682     tmp2= tmp2t[2];
1683     tmpI= tmpIt;
1684     if(b&4){
1685         for(y=0; y < b_h; y++){
1686             for(x=0; x < b_w; x++){
1687                 int a_1=tmpI[x + (HTAPS_MAX/2-4)*64];
1688                 int a0= tmpI[x + (HTAPS_MAX/2-3)*64];
1689                 int a1= tmpI[x + (HTAPS_MAX/2-2)*64];
1690                 int a2= tmpI[x + (HTAPS_MAX/2-1)*64];
1691                 int a3= tmpI[x + (HTAPS_MAX/2+0)*64];
1692                 int a4= tmpI[x + (HTAPS_MAX/2+1)*64];
1693                 int a5= tmpI[x + (HTAPS_MAX/2+2)*64];
1694                 int a6= tmpI[x + (HTAPS_MAX/2+3)*64];
1695                 int am=0;
1696                 if(!p || p->fast_mc)
1697                     am= (20*(a2+a3) - 5*(a1+a4) + (a0+a5) + 512)>>10;
1698                 else
1699                     am= (p->hcoeff[0]*(a2+a3) + p->hcoeff[1]*(a1+a4) + p->hcoeff[2]*(a0+a5) + p->hcoeff[3]*(a_1+a6) + 2048)>>12;
1700                 if(am&(~255)) am= ~(am>>31);
1701                 tmp2[x]= am;
1702             }
1703             tmpI+= 64;
1704             tmp2+= stride;
1705         }
1706     }
1707
1708     hpel[ 0]= src;
1709     hpel[ 1]= tmp2t[0] + stride*(HTAPS_MAX/2-1);
1710     hpel[ 2]= src + 1;
1711
1712     hpel[ 4]= tmp2t[1];
1713     hpel[ 5]= tmp2t[2];
1714     hpel[ 6]= tmp2t[1] + 1;
1715
1716     hpel[ 8]= src + stride;
1717     hpel[ 9]= hpel[1] + stride;
1718     hpel[10]= hpel[8] + 1;
1719
1720     if(b==15){
1721         const uint8_t *src1= hpel[dx/8 + dy/8*4  ];
1722         const uint8_t *src2= hpel[dx/8 + dy/8*4+1];
1723         const uint8_t *src3= hpel[dx/8 + dy/8*4+4];
1724         const uint8_t *src4= hpel[dx/8 + dy/8*4+5];
1725         dx&=7;
1726         dy&=7;
1727         for(y=0; y < b_h; y++){
1728             for(x=0; x < b_w; x++){
1729                 dst[x]= ((8-dx)*(8-dy)*src1[x] + dx*(8-dy)*src2[x]+
1730                          (8-dx)*   dy *src3[x] + dx*   dy *src4[x]+32)>>6;
1731             }
1732             src1+=stride;
1733             src2+=stride;
1734             src3+=stride;
1735             src4+=stride;
1736             dst +=stride;
1737         }
1738     }else{
1739         const uint8_t *src1= hpel[l];
1740         const uint8_t *src2= hpel[r];
1741         int a= weight[((dx&7) + (8*(dy&7)))];
1742         int b= 8-a;
1743         for(y=0; y < b_h; y++){
1744             for(x=0; x < b_w; x++){
1745                 dst[x]= (a*src1[x] + b*src2[x] + 4)>>3;
1746             }
1747             src1+=stride;
1748             src2+=stride;
1749             dst +=stride;
1750         }
1751     }
1752 }
1753
1754 #define mca(dx,dy,b_w)\
1755 static void mc_block_hpel ## dx ## dy ## b_w(uint8_t *dst, const uint8_t *src, int stride, int h){\
1756     uint8_t tmp[stride*(b_w+HTAPS_MAX-1)];\
1757     assert(h==b_w);\
1758     mc_block(NULL, dst, src-(HTAPS_MAX/2-1)-(HTAPS_MAX/2-1)*stride, tmp, stride, b_w, b_w, dx, dy);\
1759 }
1760
1761 mca( 0, 0,16)
1762 mca( 8, 0,16)
1763 mca( 0, 8,16)
1764 mca( 8, 8,16)
1765 mca( 0, 0,8)
1766 mca( 8, 0,8)
1767 mca( 0, 8,8)
1768 mca( 8, 8,8)
1769
1770 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){
1771     if(block->type & BLOCK_INTRA){
1772         int x, y;
1773         const int color = block->color[plane_index];
1774         const int color4= color*0x01010101;
1775         if(b_w==32){
1776             for(y=0; y < b_h; y++){
1777                 *(uint32_t*)&dst[0 + y*stride]= color4;
1778                 *(uint32_t*)&dst[4 + y*stride]= color4;
1779                 *(uint32_t*)&dst[8 + y*stride]= color4;
1780                 *(uint32_t*)&dst[12+ y*stride]= color4;
1781                 *(uint32_t*)&dst[16+ y*stride]= color4;
1782                 *(uint32_t*)&dst[20+ y*stride]= color4;
1783                 *(uint32_t*)&dst[24+ y*stride]= color4;
1784                 *(uint32_t*)&dst[28+ y*stride]= color4;
1785             }
1786         }else if(b_w==16){
1787             for(y=0; y < b_h; y++){
1788                 *(uint32_t*)&dst[0 + y*stride]= color4;
1789                 *(uint32_t*)&dst[4 + y*stride]= color4;
1790                 *(uint32_t*)&dst[8 + y*stride]= color4;
1791                 *(uint32_t*)&dst[12+ y*stride]= color4;
1792             }
1793         }else if(b_w==8){
1794             for(y=0; y < b_h; y++){
1795                 *(uint32_t*)&dst[0 + y*stride]= color4;
1796                 *(uint32_t*)&dst[4 + y*stride]= color4;
1797             }
1798         }else if(b_w==4){
1799             for(y=0; y < b_h; y++){
1800                 *(uint32_t*)&dst[0 + y*stride]= color4;
1801             }
1802         }else{
1803             for(y=0; y < b_h; y++){
1804                 for(x=0; x < b_w; x++){
1805                     dst[x + y*stride]= color;
1806                 }
1807             }
1808         }
1809     }else{
1810         uint8_t *src= s->last_picture[block->ref].data[plane_index];
1811         const int scale= plane_index ?  s->mv_scale : 2*s->mv_scale;
1812         int mx= block->mx*scale;
1813         int my= block->my*scale;
1814         const int dx= mx&15;
1815         const int dy= my&15;
1816         const int tab_index= 3 - (b_w>>2) + (b_w>>4);
1817         sx += (mx>>4) - (HTAPS_MAX/2-1);
1818         sy += (my>>4) - (HTAPS_MAX/2-1);
1819         src += sx + sy*stride;
1820         if(   (unsigned)sx >= w - b_w - (HTAPS_MAX-2)
1821            || (unsigned)sy >= h - b_h - (HTAPS_MAX-2)){
1822             ff_emulated_edge_mc(tmp + MB_SIZE, src, stride, b_w+HTAPS_MAX-1, b_h+HTAPS_MAX-1, sx, sy, w, h);
1823             src= tmp + MB_SIZE;
1824         }
1825 //        assert(b_w == b_h || 2*b_w == b_h || b_w == 2*b_h);
1826 //        assert(!(b_w&(b_w-1)));
1827         assert(b_w>1 && b_h>1);
1828         assert((tab_index>=0 && tab_index<4) || b_w==32);
1829         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 )
1830             mc_block(&s->plane[plane_index], dst, src, tmp, stride, b_w, b_h, dx, dy);
1831         else if(b_w==32){
1832             int y;
1833             for(y=0; y<b_h; y+=16){
1834                 s->dsp.put_h264_qpel_pixels_tab[0][dy+(dx>>2)](dst + y*stride, src + 3 + (y+3)*stride,stride);
1835                 s->dsp.put_h264_qpel_pixels_tab[0][dy+(dx>>2)](dst + 16 + y*stride, src + 19 + (y+3)*stride,stride);
1836             }
1837         }else if(b_w==b_h)
1838             s->dsp.put_h264_qpel_pixels_tab[tab_index  ][dy+(dx>>2)](dst,src + 3 + 3*stride,stride);
1839         else if(b_w==2*b_h){
1840             s->dsp.put_h264_qpel_pixels_tab[tab_index+1][dy+(dx>>2)](dst    ,src + 3       + 3*stride,stride);
1841             s->dsp.put_h264_qpel_pixels_tab[tab_index+1][dy+(dx>>2)](dst+b_h,src + 3 + b_h + 3*stride,stride);
1842         }else{
1843             assert(2*b_w==b_h);
1844             s->dsp.put_h264_qpel_pixels_tab[tab_index  ][dy+(dx>>2)](dst           ,src + 3 + 3*stride           ,stride);
1845             s->dsp.put_h264_qpel_pixels_tab[tab_index  ][dy+(dx>>2)](dst+b_w*stride,src + 3 + 3*stride+b_w*stride,stride);
1846         }
1847     }
1848 }
1849
1850 void ff_snow_inner_add_yblock(const uint8_t *obmc, const int obmc_stride, uint8_t * * block, int b_w, int b_h,
1851                               int src_x, int src_y, int src_stride, slice_buffer * sb, int add, uint8_t * dst8){
1852     int y, x;
1853     IDWTELEM * dst;
1854     for(y=0; y<b_h; y++){
1855         //FIXME ugly misuse of obmc_stride
1856         const uint8_t *obmc1= obmc + y*obmc_stride;
1857         const uint8_t *obmc2= obmc1+ (obmc_stride>>1);
1858         const uint8_t *obmc3= obmc1+ obmc_stride*(obmc_stride>>1);
1859         const uint8_t *obmc4= obmc3+ (obmc_stride>>1);
1860         dst = slice_buffer_get_line(sb, src_y + y);
1861         for(x=0; x<b_w; x++){
1862             int v=   obmc1[x] * block[3][x + y*src_stride]
1863                     +obmc2[x] * block[2][x + y*src_stride]
1864                     +obmc3[x] * block[1][x + y*src_stride]
1865                     +obmc4[x] * block[0][x + y*src_stride];
1866
1867             v <<= 8 - LOG2_OBMC_MAX;
1868             if(FRAC_BITS != 8){
1869                 v >>= 8 - FRAC_BITS;
1870             }
1871             if(add){
1872                 v += dst[x + src_x];
1873                 v = (v + (1<<(FRAC_BITS-1))) >> FRAC_BITS;
1874                 if(v&(~255)) v= ~(v>>31);
1875                 dst8[x + y*src_stride] = v;
1876             }else{
1877                 dst[x + src_x] -= v;
1878             }
1879         }
1880     }
1881 }
1882
1883 //FIXME name cleanup (b_w, block_w, b_width stuff)
1884 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){
1885     const int b_width = s->b_width  << s->block_max_depth;
1886     const int b_height= s->b_height << s->block_max_depth;
1887     const int b_stride= b_width;
1888     BlockNode *lt= &s->block[b_x + b_y*b_stride];
1889     BlockNode *rt= lt+1;
1890     BlockNode *lb= lt+b_stride;
1891     BlockNode *rb= lb+1;
1892     uint8_t *block[4];
1893     int tmp_step= src_stride >= 7*MB_SIZE ? MB_SIZE : MB_SIZE*src_stride;
1894     uint8_t *tmp = s->scratchbuf;
1895     uint8_t *ptmp;
1896     int x,y;
1897
1898     if(b_x<0){
1899         lt= rt;
1900         lb= rb;
1901     }else if(b_x + 1 >= b_width){
1902         rt= lt;
1903         rb= lb;
1904     }
1905     if(b_y<0){
1906         lt= lb;
1907         rt= rb;
1908     }else if(b_y + 1 >= b_height){
1909         lb= lt;
1910         rb= rt;
1911     }
1912
1913     if(src_x<0){ //FIXME merge with prev & always round internal width up to *16
1914         obmc -= src_x;
1915         b_w += src_x;
1916         if(!sliced && !offset_dst)
1917             dst -= src_x;
1918         src_x=0;
1919     }else if(src_x + b_w > w){
1920         b_w = w - src_x;
1921     }
1922     if(src_y<0){
1923         obmc -= src_y*obmc_stride;
1924         b_h += src_y;
1925         if(!sliced && !offset_dst)
1926             dst -= src_y*dst_stride;
1927         src_y=0;
1928     }else if(src_y + b_h> h){
1929         b_h = h - src_y;
1930     }
1931
1932     if(b_w<=0 || b_h<=0) return;
1933
1934     assert(src_stride > 2*MB_SIZE + 5);
1935
1936     if(!sliced && offset_dst)
1937         dst += src_x + src_y*dst_stride;
1938     dst8+= src_x + src_y*src_stride;
1939 //    src += src_x + src_y*src_stride;
1940
1941     ptmp= tmp + 3*tmp_step;
1942     block[0]= ptmp;
1943     ptmp+=tmp_step;
1944     pred_block(s, block[0], tmp, src_stride, src_x, src_y, b_w, b_h, lt, plane_index, w, h);
1945
1946     if(same_block(lt, rt)){
1947         block[1]= block[0];
1948     }else{
1949         block[1]= ptmp;
1950         ptmp+=tmp_step;
1951         pred_block(s, block[1], tmp, src_stride, src_x, src_y, b_w, b_h, rt, plane_index, w, h);
1952     }
1953
1954     if(same_block(lt, lb)){
1955         block[2]= block[0];
1956     }else if(same_block(rt, lb)){
1957         block[2]= block[1];
1958     }else{
1959         block[2]= ptmp;
1960         ptmp+=tmp_step;
1961         pred_block(s, block[2], tmp, src_stride, src_x, src_y, b_w, b_h, lb, plane_index, w, h);
1962     }
1963
1964     if(same_block(lt, rb) ){
1965         block[3]= block[0];
1966     }else if(same_block(rt, rb)){
1967         block[3]= block[1];
1968     }else if(same_block(lb, rb)){
1969         block[3]= block[2];
1970     }else{
1971         block[3]= ptmp;
1972         pred_block(s, block[3], tmp, src_stride, src_x, src_y, b_w, b_h, rb, plane_index, w, h);
1973     }
1974 #if 0
1975     for(y=0; y<b_h; y++){
1976         for(x=0; x<b_w; x++){
1977             int v=   obmc [x + y*obmc_stride] * block[3][x + y*src_stride] * (256/OBMC_MAX);
1978             if(add) dst[x + y*dst_stride] += v;
1979             else    dst[x + y*dst_stride] -= v;
1980         }
1981     }
1982     for(y=0; y<b_h; y++){
1983         uint8_t *obmc2= obmc + (obmc_stride>>1);
1984         for(x=0; x<b_w; x++){
1985             int v=   obmc2[x + y*obmc_stride] * block[2][x + y*src_stride] * (256/OBMC_MAX);
1986             if(add) dst[x + y*dst_stride] += v;
1987             else    dst[x + y*dst_stride] -= v;
1988         }
1989     }
1990     for(y=0; y<b_h; y++){
1991         uint8_t *obmc3= obmc + obmc_stride*(obmc_stride>>1);
1992         for(x=0; x<b_w; x++){
1993             int v=   obmc3[x + y*obmc_stride] * block[1][x + y*src_stride] * (256/OBMC_MAX);
1994             if(add) dst[x + y*dst_stride] += v;
1995             else    dst[x + y*dst_stride] -= v;
1996         }
1997     }
1998     for(y=0; y<b_h; y++){
1999         uint8_t *obmc3= obmc + obmc_stride*(obmc_stride>>1);
2000         uint8_t *obmc4= obmc3+ (obmc_stride>>1);
2001         for(x=0; x<b_w; x++){
2002             int v=   obmc4[x + y*obmc_stride] * block[0][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 #else
2008     if(sliced){
2009         s->dsp.inner_add_yblock(obmc, obmc_stride, block, b_w, b_h, src_x,src_y, src_stride, sb, add, dst8);
2010     }else{
2011         for(y=0; y<b_h; y++){
2012             //FIXME ugly misuse of obmc_stride
2013             const uint8_t *obmc1= obmc + y*obmc_stride;
2014             const uint8_t *obmc2= obmc1+ (obmc_stride>>1);
2015             const uint8_t *obmc3= obmc1+ obmc_stride*(obmc_stride>>1);
2016             const uint8_t *obmc4= obmc3+ (obmc_stride>>1);
2017             for(x=0; x<b_w; x++){
2018                 int v=   obmc1[x] * block[3][x + y*src_stride]
2019                         +obmc2[x] * block[2][x + y*src_stride]
2020                         +obmc3[x] * block[1][x + y*src_stride]
2021                         +obmc4[x] * block[0][x + y*src_stride];
2022
2023                 v <<= 8 - LOG2_OBMC_MAX;
2024                 if(FRAC_BITS != 8){
2025                     v >>= 8 - FRAC_BITS;
2026                 }
2027                 if(add){
2028                     v += dst[x + y*dst_stride];
2029                     v = (v + (1<<(FRAC_BITS-1))) >> FRAC_BITS;
2030                     if(v&(~255)) v= ~(v>>31);
2031                     dst8[x + y*src_stride] = v;
2032                 }else{
2033                     dst[x + y*dst_stride] -= v;
2034                 }
2035             }
2036         }
2037     }
2038 #endif /* 0 */
2039 }
2040
2041 static av_always_inline void predict_slice_buffered(SnowContext *s, slice_buffer * sb, IDWTELEM * old_buffer, int plane_index, int add, int mb_y){
2042     Plane *p= &s->plane[plane_index];
2043     const int mb_w= s->b_width  << s->block_max_depth;
2044     const int mb_h= s->b_height << s->block_max_depth;
2045     int x, y, mb_x;
2046     int block_size = MB_SIZE >> s->block_max_depth;
2047     int block_w    = plane_index ? block_size/2 : block_size;
2048     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];
2049     int obmc_stride= plane_index ? block_size : 2*block_size;
2050     int ref_stride= s->current_picture.linesize[plane_index];
2051     uint8_t *dst8= s->current_picture.data[plane_index];
2052     int w= p->width;
2053     int h= p->height;
2054
2055     if(s->keyframe || (s->avctx->debug&512)){
2056         if(mb_y==mb_h)
2057             return;
2058
2059         if(add){
2060             for(y=block_w*mb_y; y<FFMIN(h,block_w*(mb_y+1)); y++){
2061 //                DWTELEM * line = slice_buffer_get_line(sb, y);
2062                 IDWTELEM * line = sb->line[y];
2063                 for(x=0; x<w; x++){
2064 //                    int v= buf[x + y*w] + (128<<FRAC_BITS) + (1<<(FRAC_BITS-1));
2065                     int v= line[x] + (128<<FRAC_BITS) + (1<<(FRAC_BITS-1));
2066                     v >>= FRAC_BITS;
2067                     if(v&(~255)) v= ~(v>>31);
2068                     dst8[x + y*ref_stride]= v;
2069                 }
2070             }
2071         }else{
2072             for(y=block_w*mb_y; y<FFMIN(h,block_w*(mb_y+1)); y++){
2073 //                DWTELEM * line = slice_buffer_get_line(sb, y);
2074                 IDWTELEM * line = sb->line[y];
2075                 for(x=0; x<w; x++){
2076                     line[x] -= 128 << FRAC_BITS;
2077 //                    buf[x + y*w]-= 128<<FRAC_BITS;
2078                 }
2079             }
2080         }
2081
2082         return;
2083     }
2084
2085     for(mb_x=0; mb_x<=mb_w; mb_x++){
2086         add_yblock(s, 1, sb, old_buffer, dst8, obmc,
2087                    block_w*mb_x - block_w/2,
2088                    block_w*mb_y - block_w/2,
2089                    block_w, block_w,
2090                    w, h,
2091                    w, ref_stride, obmc_stride,
2092                    mb_x - 1, mb_y - 1,
2093                    add, 0, plane_index);
2094     }
2095 }
2096
2097 static av_always_inline void predict_slice(SnowContext *s, IDWTELEM *buf, int plane_index, int add, int mb_y){
2098     Plane *p= &s->plane[plane_index];
2099     const int mb_w= s->b_width  << s->block_max_depth;
2100     const int mb_h= s->b_height << s->block_max_depth;
2101     int x, y, mb_x;
2102     int block_size = MB_SIZE >> s->block_max_depth;
2103     int block_w    = plane_index ? block_size/2 : block_size;
2104     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];
2105     const int obmc_stride= plane_index ? block_size : 2*block_size;
2106     int ref_stride= s->current_picture.linesize[plane_index];
2107     uint8_t *dst8= s->current_picture.data[plane_index];
2108     int w= p->width;
2109     int h= p->height;
2110
2111     if(s->keyframe || (s->avctx->debug&512)){
2112         if(mb_y==mb_h)
2113             return;
2114
2115         if(add){
2116             for(y=block_w*mb_y; y<FFMIN(h,block_w*(mb_y+1)); y++){
2117                 for(x=0; x<w; x++){
2118                     int v= buf[x + y*w] + (128<<FRAC_BITS) + (1<<(FRAC_BITS-1));
2119                     v >>= FRAC_BITS;
2120                     if(v&(~255)) v= ~(v>>31);
2121                     dst8[x + y*ref_stride]= v;
2122                 }
2123             }
2124         }else{
2125             for(y=block_w*mb_y; y<FFMIN(h,block_w*(mb_y+1)); y++){
2126                 for(x=0; x<w; x++){
2127                     buf[x + y*w]-= 128<<FRAC_BITS;
2128                 }
2129             }
2130         }
2131
2132         return;
2133     }
2134
2135     for(mb_x=0; mb_x<=mb_w; mb_x++){
2136         add_yblock(s, 0, NULL, buf, dst8, obmc,
2137                    block_w*mb_x - block_w/2,
2138                    block_w*mb_y - block_w/2,
2139                    block_w, block_w,
2140                    w, h,
2141                    w, ref_stride, obmc_stride,
2142                    mb_x - 1, mb_y - 1,
2143                    add, 1, plane_index);
2144     }
2145 }
2146
2147 static av_always_inline void predict_plane(SnowContext *s, IDWTELEM *buf, int plane_index, int add){
2148     const int mb_h= s->b_height << s->block_max_depth;
2149     int mb_y;
2150     for(mb_y=0; mb_y<=mb_h; mb_y++)
2151         predict_slice(s, buf, plane_index, add, mb_y);
2152 }
2153
2154 static void dequantize_slice_buffered(SnowContext *s, slice_buffer * sb, SubBand *b, IDWTELEM *src, int stride, int start_y, int end_y){
2155     const int w= b->width;
2156     const int qlog= av_clip(s->qlog + b->qlog, 0, QROOT*16);
2157     const int qmul= qexp[qlog&(QROOT-1)]<<(qlog>>QSHIFT);
2158     const int qadd= (s->qbias*qmul)>>QBIAS_SHIFT;
2159     int x,y;
2160
2161     if(s->qlog == LOSSLESS_QLOG) return;
2162
2163     for(y=start_y; y<end_y; y++){
2164 //        DWTELEM * line = slice_buffer_get_line_from_address(sb, src + (y * stride));
2165         IDWTELEM * line = slice_buffer_get_line(sb, (y * b->stride_line) + b->buf_y_offset) + b->buf_x_offset;
2166         for(x=0; x<w; x++){
2167             int i= line[x];
2168             if(i<0){
2169                 line[x]= -((-i*qmul + qadd)>>(QEXPSHIFT)); //FIXME try different bias
2170             }else if(i>0){
2171                 line[x]=  (( i*qmul + qadd)>>(QEXPSHIFT));
2172             }
2173         }
2174     }
2175 }
2176
2177 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){
2178     const int w= b->width;
2179     int x,y;
2180
2181     IDWTELEM * line=0; // silence silly "could be used without having been initialized" warning
2182     IDWTELEM * prev;
2183
2184     if (start_y != 0)
2185         line = slice_buffer_get_line(sb, ((start_y - 1) * b->stride_line) + b->buf_y_offset) + b->buf_x_offset;
2186
2187     for(y=start_y; y<end_y; y++){
2188         prev = line;
2189 //        line = slice_buffer_get_line_from_address(sb, src + (y * stride));
2190         line = slice_buffer_get_line(sb, (y * b->stride_line) + b->buf_y_offset) + b->buf_x_offset;
2191         for(x=0; x<w; x++){
2192             if(x){
2193                 if(use_median){
2194                     if(y && x+1<w) line[x] += mid_pred(line[x - 1], prev[x], prev[x + 1]);
2195                     else  line[x] += line[x - 1];
2196                 }else{
2197                     if(y) line[x] += mid_pred(line[x - 1], prev[x], line[x - 1] + prev[x] - prev[x - 1]);
2198                     else  line[x] += line[x - 1];
2199                 }
2200             }else{
2201                 if(y) line[x] += prev[x];
2202             }
2203         }
2204     }
2205 }
2206
2207 static void decode_qlogs(SnowContext *s){
2208     int plane_index, level, orientation;
2209
2210     for(plane_index=0; plane_index<3; plane_index++){
2211         for(level=0; level<s->spatial_decomposition_count; level++){
2212             for(orientation=level ? 1:0; orientation<4; orientation++){
2213                 int q;
2214                 if     (plane_index==2) q= s->plane[1].band[level][orientation].qlog;
2215                 else if(orientation==2) q= s->plane[plane_index].band[level][1].qlog;
2216                 else                    q= get_symbol(&s->c, s->header_state, 1);
2217                 s->plane[plane_index].band[level][orientation].qlog= q;
2218             }
2219         }
2220     }
2221 }
2222
2223 #define GET_S(dst, check) \
2224     tmp= get_symbol(&s->c, s->header_state, 0);\
2225     if(!(check)){\
2226         av_log(s->avctx, AV_LOG_ERROR, "Error " #dst " is %d\n", tmp);\
2227         return -1;\
2228     }\
2229     dst= tmp;
2230
2231 static int decode_header(SnowContext *s){
2232     int plane_index, tmp;
2233     uint8_t kstate[32];
2234
2235     memset(kstate, MID_STATE, sizeof(kstate));
2236
2237     s->keyframe= get_rac(&s->c, kstate);
2238     if(s->keyframe || s->always_reset){
2239         reset_contexts(s);
2240         s->spatial_decomposition_type=
2241         s->qlog=
2242         s->qbias=
2243         s->mv_scale=
2244         s->block_max_depth= 0;
2245     }
2246     if(s->keyframe){
2247         GET_S(s->version, tmp <= 0U)
2248         s->always_reset= get_rac(&s->c, s->header_state);
2249         s->temporal_decomposition_type= get_symbol(&s->c, s->header_state, 0);
2250         s->temporal_decomposition_count= get_symbol(&s->c, s->header_state, 0);
2251         GET_S(s->spatial_decomposition_count, 0 < tmp && tmp <= MAX_DECOMPOSITIONS)
2252         s->colorspace_type= get_symbol(&s->c, s->header_state, 0);
2253         s->chroma_h_shift= get_symbol(&s->c, s->header_state, 0);
2254         s->chroma_v_shift= get_symbol(&s->c, s->header_state, 0);
2255         s->spatial_scalability= get_rac(&s->c, s->header_state);
2256 //        s->rate_scalability= get_rac(&s->c, s->header_state);
2257         GET_S(s->max_ref_frames, tmp < (unsigned)MAX_REF_FRAMES)
2258         s->max_ref_frames++;
2259
2260         decode_qlogs(s);
2261     }
2262
2263     if(!s->keyframe){
2264         if(get_rac(&s->c, s->header_state)){
2265             for(plane_index=0; plane_index<2; plane_index++){
2266                 int htaps, i, sum=0;
2267                 Plane *p= &s->plane[plane_index];
2268                 p->diag_mc= get_rac(&s->c, s->header_state);
2269                 htaps= get_symbol(&s->c, s->header_state, 0)*2 + 2;
2270                 if((unsigned)htaps > HTAPS_MAX || htaps==0)
2271                     return -1;
2272                 p->htaps= htaps;
2273                 for(i= htaps/2; i; i--){
2274                     p->hcoeff[i]= get_symbol(&s->c, s->header_state, 0) * (1-2*(i&1));
2275                     sum += p->hcoeff[i];
2276                 }
2277                 p->hcoeff[0]= 32-sum;
2278             }
2279             s->plane[2].diag_mc= s->plane[1].diag_mc;
2280             s->plane[2].htaps  = s->plane[1].htaps;
2281             memcpy(s->plane[2].hcoeff, s->plane[1].hcoeff, sizeof(s->plane[1].hcoeff));
2282         }
2283         if(get_rac(&s->c, s->header_state)){
2284             GET_S(s->spatial_decomposition_count, 0 < tmp && tmp <= MAX_DECOMPOSITIONS)
2285             decode_qlogs(s);
2286         }
2287     }
2288
2289     s->spatial_decomposition_type+= get_symbol(&s->c, s->header_state, 1);
2290     if(s->spatial_decomposition_type > 1U){
2291         av_log(s->avctx, AV_LOG_ERROR, "spatial_decomposition_type %d not supported", s->spatial_decomposition_type);
2292         return -1;
2293     }
2294     if(FFMIN(s->avctx-> width>>s->chroma_h_shift,
2295              s->avctx->height>>s->chroma_v_shift) >> (s->spatial_decomposition_count-1) <= 0){
2296         av_log(s->avctx, AV_LOG_ERROR, "spatial_decomposition_count %d too large for size", s->spatial_decomposition_count);
2297         return -1;
2298     }
2299
2300     s->qlog           += get_symbol(&s->c, s->header_state, 1);
2301     s->mv_scale       += get_symbol(&s->c, s->header_state, 1);
2302     s->qbias          += get_symbol(&s->c, s->header_state, 1);
2303     s->block_max_depth+= get_symbol(&s->c, s->header_state, 1);
2304     if(s->block_max_depth > 1 || s->block_max_depth < 0){
2305         av_log(s->avctx, AV_LOG_ERROR, "block_max_depth= %d is too large", s->block_max_depth);
2306         s->block_max_depth= 0;
2307         return -1;
2308     }
2309
2310     return 0;
2311 }
2312
2313 static void init_qexp(void){
2314     int i;
2315     double v=128;
2316
2317     for(i=0; i<QROOT; i++){
2318         qexp[i]= lrintf(v);
2319         v *= pow(2, 1.0 / QROOT);
2320     }
2321 }
2322
2323 static av_cold int common_init(AVCodecContext *avctx){
2324     SnowContext *s = avctx->priv_data;
2325     int width, height;
2326     int i, j;
2327
2328     s->avctx= avctx;
2329     s->max_ref_frames=1; //just make sure its not an invalid value in case of no initial keyframe
2330
2331     dsputil_init(&s->dsp, avctx);
2332
2333 #define mcf(dx,dy)\
2334     s->dsp.put_qpel_pixels_tab       [0][dy+dx/4]=\
2335     s->dsp.put_no_rnd_qpel_pixels_tab[0][dy+dx/4]=\
2336         s->dsp.put_h264_qpel_pixels_tab[0][dy+dx/4];\
2337     s->dsp.put_qpel_pixels_tab       [1][dy+dx/4]=\
2338     s->dsp.put_no_rnd_qpel_pixels_tab[1][dy+dx/4]=\
2339         s->dsp.put_h264_qpel_pixels_tab[1][dy+dx/4];
2340
2341     mcf( 0, 0)
2342     mcf( 4, 0)
2343     mcf( 8, 0)
2344     mcf(12, 0)
2345     mcf( 0, 4)
2346     mcf( 4, 4)
2347     mcf( 8, 4)
2348     mcf(12, 4)
2349     mcf( 0, 8)
2350     mcf( 4, 8)
2351     mcf( 8, 8)
2352     mcf(12, 8)
2353     mcf( 0,12)
2354     mcf( 4,12)
2355     mcf( 8,12)
2356     mcf(12,12)
2357
2358 #define mcfh(dx,dy)\
2359     s->dsp.put_pixels_tab       [0][dy/4+dx/8]=\
2360     s->dsp.put_no_rnd_pixels_tab[0][dy/4+dx/8]=\
2361         mc_block_hpel ## dx ## dy ## 16;\
2362     s->dsp.put_pixels_tab       [1][dy/4+dx/8]=\
2363     s->dsp.put_no_rnd_pixels_tab[1][dy/4+dx/8]=\
2364         mc_block_hpel ## dx ## dy ## 8;
2365
2366     mcfh(0, 0)
2367     mcfh(8, 0)
2368     mcfh(0, 8)
2369     mcfh(8, 8)
2370
2371     if(!qexp[0])
2372         init_qexp();
2373
2374 //    dec += FFMAX(s->chroma_h_shift, s->chroma_v_shift);
2375
2376     width= s->avctx->width;
2377     height= s->avctx->height;
2378
2379     s->spatial_idwt_buffer= av_mallocz(width*height*sizeof(IDWTELEM));
2380     s->spatial_dwt_buffer= av_mallocz(width*height*sizeof(DWTELEM)); //FIXME this does not belong here
2381
2382     for(i=0; i<MAX_REF_FRAMES; i++)
2383         for(j=0; j<MAX_REF_FRAMES; j++)
2384             scale_mv_ref[i][j] = 256*(i+1)/(j+1);
2385
2386     s->avctx->get_buffer(s->avctx, &s->mconly_picture);
2387     s->scratchbuf = av_malloc(s->mconly_picture.linesize[0]*7*MB_SIZE);
2388
2389     return 0;
2390 }
2391
2392 static int common_init_after_header(AVCodecContext *avctx){
2393     SnowContext *s = avctx->priv_data;
2394     int plane_index, level, orientation;
2395
2396     for(plane_index=0; plane_index<3; plane_index++){
2397         int w= s->avctx->width;
2398         int h= s->avctx->height;
2399
2400         if(plane_index){
2401             w>>= s->chroma_h_shift;
2402             h>>= s->chroma_v_shift;
2403         }
2404         s->plane[plane_index].width = w;
2405         s->plane[plane_index].height= h;
2406
2407         for(level=s->spatial_decomposition_count-1; level>=0; level--){
2408             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2409                 SubBand *b= &s->plane[plane_index].band[level][orientation];
2410
2411                 b->buf= s->spatial_dwt_buffer;
2412                 b->level= level;
2413                 b->stride= s->plane[plane_index].width << (s->spatial_decomposition_count - level);
2414                 b->width = (w + !(orientation&1))>>1;
2415                 b->height= (h + !(orientation>1))>>1;
2416
2417                 b->stride_line = 1 << (s->spatial_decomposition_count - level);
2418                 b->buf_x_offset = 0;
2419                 b->buf_y_offset = 0;
2420
2421                 if(orientation&1){
2422                     b->buf += (w+1)>>1;
2423                     b->buf_x_offset = (w+1)>>1;
2424                 }
2425                 if(orientation>1){
2426                     b->buf += b->stride>>1;
2427                     b->buf_y_offset = b->stride_line >> 1;
2428                 }
2429                 b->ibuf= s->spatial_idwt_buffer + (b->buf - s->spatial_dwt_buffer);
2430
2431                 if(level)
2432                     b->parent= &s->plane[plane_index].band[level-1][orientation];
2433                 //FIXME avoid this realloc
2434                 av_freep(&b->x_coeff);
2435                 b->x_coeff=av_mallocz(((b->width+1) * b->height+1)*sizeof(x_and_coeff));
2436             }
2437             w= (w+1)>>1;
2438             h= (h+1)>>1;
2439         }
2440     }
2441
2442     return 0;
2443 }
2444
2445 #define QUANTIZE2 0
2446
2447 #if QUANTIZE2==1
2448 #define Q2_STEP 8
2449
2450 static void find_sse(SnowContext *s, Plane *p, int *score, int score_stride, IDWTELEM *r0, IDWTELEM *r1, int level, int orientation){
2451     SubBand *b= &p->band[level][orientation];
2452     int x, y;
2453     int xo=0;
2454     int yo=0;
2455     int step= 1 << (s->spatial_decomposition_count - level);
2456
2457     if(orientation&1)
2458         xo= step>>1;
2459     if(orientation&2)
2460         yo= step>>1;
2461
2462     //FIXME bias for nonzero ?
2463     //FIXME optimize
2464     memset(score, 0, sizeof(*score)*score_stride*((p->height + Q2_STEP-1)/Q2_STEP));
2465     for(y=0; y<p->height; y++){
2466         for(x=0; x<p->width; x++){
2467             int sx= (x-xo + step/2) / step / Q2_STEP;
2468             int sy= (y-yo + step/2) / step / Q2_STEP;
2469             int v= r0[x + y*p->width] - r1[x + y*p->width];
2470             assert(sx>=0 && sy>=0 && sx < score_stride);
2471             v= ((v+8)>>4)<<4;
2472             score[sx + sy*score_stride] += v*v;
2473             assert(score[sx + sy*score_stride] >= 0);
2474         }
2475     }
2476 }
2477
2478 static void dequantize_all(SnowContext *s, Plane *p, IDWTELEM *buffer, int width, int height){
2479     int level, orientation;
2480
2481     for(level=0; level<s->spatial_decomposition_count; level++){
2482         for(orientation=level ? 1 : 0; orientation<4; orientation++){
2483             SubBand *b= &p->band[level][orientation];
2484             IDWTELEM *dst= buffer + (b->ibuf - s->spatial_idwt_buffer);
2485
2486             dequantize(s, b, dst, b->stride);
2487         }
2488     }
2489 }
2490
2491 static void dwt_quantize(SnowContext *s, Plane *p, DWTELEM *buffer, int width, int height, int stride, int type){
2492     int level, orientation, ys, xs, x, y, pass;
2493     IDWTELEM best_dequant[height * stride];
2494     IDWTELEM idwt2_buffer[height * stride];
2495     const int score_stride= (width + 10)/Q2_STEP;
2496     int best_score[(width + 10)/Q2_STEP * (height + 10)/Q2_STEP]; //FIXME size
2497     int score[(width + 10)/Q2_STEP * (height + 10)/Q2_STEP]; //FIXME size
2498     int threshold= (s->m.lambda * s->m.lambda) >> 6;
2499
2500     //FIXME pass the copy cleanly ?
2501
2502 //    memcpy(dwt_buffer, buffer, height * stride * sizeof(DWTELEM));
2503     ff_spatial_dwt(buffer, width, height, stride, type, s->spatial_decomposition_count);
2504
2505     for(level=0; level<s->spatial_decomposition_count; level++){
2506         for(orientation=level ? 1 : 0; orientation<4; orientation++){
2507             SubBand *b= &p->band[level][orientation];
2508             IDWTELEM *dst= best_dequant + (b->ibuf - s->spatial_idwt_buffer);
2509              DWTELEM *src=       buffer + (b-> buf - s->spatial_dwt_buffer);
2510             assert(src == b->buf); // code does not depend on this but it is true currently
2511
2512             quantize(s, b, dst, src, b->stride, s->qbias);
2513         }
2514     }
2515     for(pass=0; pass<1; pass++){
2516         if(s->qbias == 0) //keyframe
2517             continue;
2518         for(level=0; level<s->spatial_decomposition_count; level++){
2519             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2520                 SubBand *b= &p->band[level][orientation];
2521                 IDWTELEM *dst= idwt2_buffer + (b->ibuf - s->spatial_idwt_buffer);
2522                 IDWTELEM *best_dst= best_dequant + (b->ibuf - s->spatial_idwt_buffer);
2523
2524                 for(ys= 0; ys<Q2_STEP; ys++){
2525                     for(xs= 0; xs<Q2_STEP; xs++){
2526                         memcpy(idwt2_buffer, best_dequant, height * stride * sizeof(IDWTELEM));
2527                         dequantize_all(s, p, idwt2_buffer, width, height);
2528                         ff_spatial_idwt(idwt2_buffer, width, height, stride, type, s->spatial_decomposition_count);
2529                         find_sse(s, p, best_score, score_stride, idwt2_buffer, s->spatial_idwt_buffer, level, orientation);
2530                         memcpy(idwt2_buffer, best_dequant, height * stride * sizeof(IDWTELEM));
2531                         for(y=ys; y<b->height; y+= Q2_STEP){
2532                             for(x=xs; x<b->width; x+= Q2_STEP){
2533                                 if(dst[x + y*b->stride]<0) dst[x + y*b->stride]++;
2534                                 if(dst[x + y*b->stride]>0) dst[x + y*b->stride]--;
2535                                 //FIXME try more than just --
2536                             }
2537                         }
2538                         dequantize_all(s, p, idwt2_buffer, width, height);
2539                         ff_spatial_idwt(idwt2_buffer, width, height, stride, type, s->spatial_decomposition_count);
2540                         find_sse(s, p, score, score_stride, idwt2_buffer, s->spatial_idwt_buffer, level, orientation);
2541                         for(y=ys; y<b->height; y+= Q2_STEP){
2542                             for(x=xs; x<b->width; x+= Q2_STEP){
2543                                 int score_idx= x/Q2_STEP + (y/Q2_STEP)*score_stride;
2544                                 if(score[score_idx] <= best_score[score_idx] + threshold){
2545                                     best_score[score_idx]= score[score_idx];
2546                                     if(best_dst[x + y*b->stride]<0) best_dst[x + y*b->stride]++;
2547                                     if(best_dst[x + y*b->stride]>0) best_dst[x + y*b->stride]--;
2548                                     //FIXME copy instead
2549                                 }
2550                             }
2551                         }
2552                     }
2553                 }
2554             }
2555         }
2556     }
2557     memcpy(s->spatial_idwt_buffer, best_dequant, height * stride * sizeof(IDWTELEM)); //FIXME work with that directly instead of copy at the end
2558 }
2559
2560 #endif /* QUANTIZE2==1 */
2561
2562 #define USE_HALFPEL_PLANE 0
2563
2564 static void halfpel_interpol(SnowContext *s, uint8_t *halfpel[4][4], AVFrame *frame){
2565     int p,x,y;
2566
2567     assert(!(s->avctx->flags & CODEC_FLAG_EMU_EDGE));
2568
2569     for(p=0; p<3; p++){
2570         int is_chroma= !!p;
2571         int w= s->avctx->width  >>is_chroma;
2572         int h= s->avctx->height >>is_chroma;
2573         int ls= frame->linesize[p];
2574         uint8_t *src= frame->data[p];
2575
2576         halfpel[1][p]= (uint8_t*)av_malloc(ls * (h+2*EDGE_WIDTH)) + EDGE_WIDTH*(1+ls);
2577         halfpel[2][p]= (uint8_t*)av_malloc(ls * (h+2*EDGE_WIDTH)) + EDGE_WIDTH*(1+ls);
2578         halfpel[3][p]= (uint8_t*)av_malloc(ls * (h+2*EDGE_WIDTH)) + EDGE_WIDTH*(1+ls);
2579
2580         halfpel[0][p]= src;
2581         for(y=0; y<h; y++){
2582             for(x=0; x<w; x++){
2583                 int i= y*ls + x;
2584
2585                 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;
2586             }
2587         }
2588         for(y=0; y<h; y++){
2589             for(x=0; x<w; x++){
2590                 int i= y*ls + x;
2591
2592                 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;
2593             }
2594         }
2595         src= halfpel[1][p];
2596         for(y=0; y<h; y++){
2597             for(x=0; x<w; x++){
2598                 int i= y*ls + x;
2599
2600                 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;
2601             }
2602         }
2603
2604 //FIXME border!
2605     }
2606 }
2607
2608 static void release_buffer(AVCodecContext *avctx){
2609     SnowContext *s = avctx->priv_data;
2610     int i;
2611
2612     if(s->last_picture[s->max_ref_frames-1].data[0]){
2613         avctx->release_buffer(avctx, &s->last_picture[s->max_ref_frames-1]);
2614         for(i=0; i<9; i++)
2615             if(s->halfpel_plane[s->max_ref_frames-1][1+i/3][i%3])
2616                 av_free(s->halfpel_plane[s->max_ref_frames-1][1+i/3][i%3] - EDGE_WIDTH*(1+s->current_picture.linesize[i%3]));
2617     }
2618 }
2619
2620 static int frame_start(SnowContext *s){
2621    AVFrame tmp;
2622    int w= s->avctx->width; //FIXME round up to x16 ?
2623    int h= s->avctx->height;
2624
2625     if(s->current_picture.data[0]){
2626         s->dsp.draw_edges(s->current_picture.data[0], s->current_picture.linesize[0], w   , h   , EDGE_WIDTH  );
2627         s->dsp.draw_edges(s->current_picture.data[1], s->current_picture.linesize[1], w>>1, h>>1, EDGE_WIDTH/2);
2628         s->dsp.draw_edges(s->current_picture.data[2], s->current_picture.linesize[2], w>>1, h>>1, EDGE_WIDTH/2);
2629     }
2630
2631     release_buffer(s->avctx);
2632
2633     tmp= s->last_picture[s->max_ref_frames-1];
2634     memmove(s->last_picture+1, s->last_picture, (s->max_ref_frames-1)*sizeof(AVFrame));
2635     memmove(s->halfpel_plane+1, s->halfpel_plane, (s->max_ref_frames-1)*sizeof(void*)*4*4);
2636     if(USE_HALFPEL_PLANE && s->current_picture.data[0])
2637         halfpel_interpol(s, s->halfpel_plane[0], &s->current_picture);
2638     s->last_picture[0]= s->current_picture;
2639     s->current_picture= tmp;
2640
2641     if(s->keyframe){
2642         s->ref_frames= 0;
2643     }else{
2644         int i;
2645         for(i=0; i<s->max_ref_frames && s->last_picture[i].data[0]; i++)
2646             if(i && s->last_picture[i-1].key_frame)
2647                 break;
2648         s->ref_frames= i;
2649         if(s->ref_frames==0){
2650             av_log(s->avctx,AV_LOG_ERROR, "No reference frames\n");
2651             return -1;
2652         }
2653     }
2654
2655     s->current_picture.reference= 1;
2656     if(s->avctx->get_buffer(s->avctx, &s->current_picture) < 0){
2657         av_log(s->avctx, AV_LOG_ERROR, "get_buffer() failed\n");
2658         return -1;
2659     }
2660
2661     s->current_picture.key_frame= s->keyframe;
2662
2663     return 0;
2664 }
2665
2666 static av_cold void common_end(SnowContext *s){
2667     int plane_index, level, orientation, i;
2668
2669     av_freep(&s->spatial_dwt_buffer);
2670     av_freep(&s->spatial_idwt_buffer);
2671
2672     s->m.me.temp= NULL;
2673     av_freep(&s->m.me.scratchpad);
2674     av_freep(&s->m.me.map);
2675     av_freep(&s->m.me.score_map);
2676     av_freep(&s->m.obmc_scratchpad);
2677
2678     av_freep(&s->block);
2679     av_freep(&s->scratchbuf);
2680
2681     for(i=0; i<MAX_REF_FRAMES; i++){
2682         av_freep(&s->ref_mvs[i]);
2683         av_freep(&s->ref_scores[i]);
2684         if(s->last_picture[i].data[0])
2685             s->avctx->release_buffer(s->avctx, &s->last_picture[i]);
2686     }
2687
2688     for(plane_index=0; plane_index<3; plane_index++){
2689         for(level=s->spatial_decomposition_count-1; level>=0; level--){
2690             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2691                 SubBand *b= &s->plane[plane_index].band[level][orientation];
2692
2693                 av_freep(&b->x_coeff);
2694             }
2695         }
2696     }
2697     if (s->mconly_picture.data[0])
2698         s->avctx->release_buffer(s->avctx, &s->mconly_picture);
2699     if (s->current_picture.data[0])
2700         s->avctx->release_buffer(s->avctx, &s->current_picture);
2701 }
2702
2703 static av_cold int decode_init(AVCodecContext *avctx)
2704 {
2705     avctx->pix_fmt= PIX_FMT_YUV420P;
2706
2707     common_init(avctx);
2708
2709     return 0;
2710 }
2711
2712 static int decode_frame(AVCodecContext *avctx, void *data, int *data_size, AVPacket *avpkt){
2713     const uint8_t *buf = avpkt->data;
2714     int buf_size = avpkt->size;
2715     SnowContext *s = avctx->priv_data;
2716     RangeCoder * const c= &s->c;
2717     int bytes_read;
2718     AVFrame *picture = data;
2719     int level, orientation, plane_index;
2720
2721     ff_init_range_decoder(c, buf, buf_size);
2722     ff_build_rac_states(c, 0.05*(1LL<<32), 256-8);
2723
2724     s->current_picture.pict_type= FF_I_TYPE; //FIXME I vs. P
2725     if(decode_header(s)<0)
2726         return -1;
2727     common_init_after_header(avctx);
2728
2729     // realloc slice buffer for the case that spatial_decomposition_count changed
2730     slice_buffer_destroy(&s->sb);
2731     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);
2732
2733     for(plane_index=0; plane_index<3; plane_index++){
2734         Plane *p= &s->plane[plane_index];
2735         p->fast_mc= p->diag_mc && p->htaps==6 && p->hcoeff[0]==40
2736                                               && p->hcoeff[1]==-10
2737                                               && p->hcoeff[2]==2;
2738     }
2739
2740     alloc_blocks(s);
2741
2742     if(frame_start(s) < 0)
2743         return -1;
2744     //keyframe flag duplication mess FIXME
2745     if(avctx->debug&FF_DEBUG_PICT_INFO)
2746         av_log(avctx, AV_LOG_ERROR, "keyframe:%d qlog:%d\n", s->keyframe, s->qlog);
2747
2748     decode_blocks(s);
2749
2750     for(plane_index=0; plane_index<3; plane_index++){
2751         Plane *p= &s->plane[plane_index];
2752         int w= p->width;
2753         int h= p->height;
2754         int x, y;
2755         int decode_state[MAX_DECOMPOSITIONS][4][1]; /* Stored state info for unpack_coeffs. 1 variable per instance. */
2756
2757         if(s->avctx->debug&2048){
2758             memset(s->spatial_dwt_buffer, 0, sizeof(DWTELEM)*w*h);
2759             predict_plane(s, s->spatial_idwt_buffer, plane_index, 1);
2760
2761             for(y=0; y<h; y++){
2762                 for(x=0; x<w; x++){
2763                     int v= s->current_picture.data[plane_index][y*s->current_picture.linesize[plane_index] + x];
2764                     s->mconly_picture.data[plane_index][y*s->mconly_picture.linesize[plane_index] + x]= v;
2765                 }
2766             }
2767         }
2768
2769         {
2770         for(level=0; level<s->spatial_decomposition_count; level++){
2771             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2772                 SubBand *b= &p->band[level][orientation];
2773                 unpack_coeffs(s, b, b->parent, orientation);
2774             }
2775         }
2776         }
2777
2778         {
2779         const int mb_h= s->b_height << s->block_max_depth;
2780         const int block_size = MB_SIZE >> s->block_max_depth;
2781         const int block_w    = plane_index ? block_size/2 : block_size;
2782         int mb_y;
2783         DWTCompose cs[MAX_DECOMPOSITIONS];
2784         int yd=0, yq=0;
2785         int y;
2786         int end_y;
2787
2788         ff_spatial_idwt_buffered_init(cs, &s->sb, w, h, 1, s->spatial_decomposition_type, s->spatial_decomposition_count);
2789         for(mb_y=0; mb_y<=mb_h; mb_y++){
2790
2791             int slice_starty = block_w*mb_y;
2792             int slice_h = block_w*(mb_y+1);
2793             if (!(s->keyframe || s->avctx->debug&512)){
2794                 slice_starty = FFMAX(0, slice_starty - (block_w >> 1));
2795                 slice_h -= (block_w >> 1);
2796             }
2797
2798             for(level=0; level<s->spatial_decomposition_count; level++){
2799                 for(orientation=level ? 1 : 0; orientation<4; orientation++){
2800                     SubBand *b= &p->band[level][orientation];
2801                     int start_y;
2802                     int end_y;
2803                     int our_mb_start = mb_y;
2804                     int our_mb_end = (mb_y + 1);
2805                     const int extra= 3;
2806                     start_y = (mb_y ? ((block_w * our_mb_start) >> (s->spatial_decomposition_count - level)) + s->spatial_decomposition_count - level + extra: 0);
2807                     end_y = (((block_w * our_mb_end) >> (s->spatial_decomposition_count - level)) + s->spatial_decomposition_count - level + extra);
2808                     if (!(s->keyframe || s->avctx->debug&512)){
2809                         start_y = FFMAX(0, start_y - (block_w >> (1+s->spatial_decomposition_count - level)));
2810                         end_y = FFMAX(0, end_y - (block_w >> (1+s->spatial_decomposition_count - level)));
2811                     }
2812                     start_y = FFMIN(b->height, start_y);
2813                     end_y = FFMIN(b->height, end_y);
2814
2815                     if (start_y != end_y){
2816                         if (orientation == 0){
2817                             SubBand * correlate_band = &p->band[0][0];
2818                             int correlate_end_y = FFMIN(b->height, end_y + 1);
2819                             int correlate_start_y = FFMIN(b->height, (start_y ? start_y + 1 : 0));
2820                             decode_subband_slice_buffered(s, correlate_band, &s->sb, correlate_start_y, correlate_end_y, decode_state[0][0]);
2821                             correlate_slice_buffered(s, &s->sb, correlate_band, correlate_band->ibuf, correlate_band->stride, 1, 0, correlate_start_y, correlate_end_y);
2822                             dequantize_slice_buffered(s, &s->sb, correlate_band, correlate_band->ibuf, correlate_band->stride, start_y, end_y);
2823                         }
2824                         else
2825                             decode_subband_slice_buffered(s, b, &s->sb, start_y, end_y, decode_state[level][orientation]);
2826                     }
2827                 }
2828             }
2829
2830             for(; yd<slice_h; yd+=4){
2831                 ff_spatial_idwt_buffered_slice(&s->dsp, cs, &s->sb, w, h, 1, s->spatial_decomposition_type, s->spatial_decomposition_count, yd);
2832             }
2833
2834             if(s->qlog == LOSSLESS_QLOG){
2835                 for(; yq<slice_h && yq<h; yq++){
2836                     IDWTELEM * line = slice_buffer_get_line(&s->sb, yq);
2837                     for(x=0; x<w; x++){
2838                         line[x] <<= FRAC_BITS;
2839                     }
2840                 }
2841             }
2842
2843             predict_slice_buffered(s, &s->sb, s->spatial_idwt_buffer, plane_index, 1, mb_y);
2844
2845             y = FFMIN(p->height, slice_starty);
2846             end_y = FFMIN(p->height, slice_h);
2847             while(y < end_y)
2848                 slice_buffer_release(&s->sb, y++);
2849         }
2850
2851         slice_buffer_flush(&s->sb);
2852         }
2853
2854     }
2855
2856     emms_c();
2857
2858     release_buffer(avctx);
2859
2860     if(!(s->avctx->debug&2048))
2861         *picture= s->current_picture;
2862     else
2863         *picture= s->mconly_picture;
2864
2865     *data_size = sizeof(AVFrame);
2866
2867     bytes_read= c->bytestream - c->bytestream_start;
2868     if(bytes_read ==0) av_log(s->avctx, AV_LOG_ERROR, "error at end of frame\n"); //FIXME
2869
2870     return bytes_read;
2871 }
2872
2873 static av_cold int decode_end(AVCodecContext *avctx)
2874 {
2875     SnowContext *s = avctx->priv_data;
2876
2877     slice_buffer_destroy(&s->sb);
2878
2879     common_end(s);
2880
2881     return 0;
2882 }
2883
2884 AVCodec snow_decoder = {
2885     "snow",
2886     CODEC_TYPE_VIDEO,
2887     CODEC_ID_SNOW,
2888     sizeof(SnowContext),
2889     decode_init,
2890     NULL,
2891     decode_end,
2892     decode_frame,
2893     CODEC_CAP_DR1 /*| CODEC_CAP_DRAW_HORIZ_BAND*/,
2894     NULL,
2895     .long_name = NULL_IF_CONFIG_SMALL("Snow"),
2896 };
2897
2898 #if CONFIG_SNOW_ENCODER
2899 static av_cold int encode_init(AVCodecContext *avctx)
2900 {
2901     SnowContext *s = avctx->priv_data;
2902     int plane_index;
2903
2904     if(avctx->strict_std_compliance > FF_COMPLIANCE_EXPERIMENTAL){
2905         av_log(avctx, AV_LOG_ERROR, "This codec is under development, files encoded with it may not be decodable with future versions!!!\n"
2906                "Use vstrict=-2 / -strict -2 to use it anyway.\n");
2907         return -1;
2908     }
2909
2910     if(avctx->prediction_method == DWT_97
2911        && (avctx->flags & CODEC_FLAG_QSCALE)
2912        && avctx->global_quality == 0){
2913         av_log(avctx, AV_LOG_ERROR, "The 9/7 wavelet is incompatible with lossless mode.\n");
2914         return -1;
2915     }
2916
2917     s->spatial_decomposition_type= avctx->prediction_method; //FIXME add decorrelator type r transform_type
2918
2919     s->mv_scale       = (avctx->flags & CODEC_FLAG_QPEL) ? 2 : 4;
2920     s->block_max_depth= (avctx->flags & CODEC_FLAG_4MV ) ? 1 : 0;
2921
2922     for(plane_index=0; plane_index<3; plane_index++){
2923         s->plane[plane_index].diag_mc= 1;
2924         s->plane[plane_index].htaps= 6;
2925         s->plane[plane_index].hcoeff[0]=  40;
2926         s->plane[plane_index].hcoeff[1]= -10;
2927         s->plane[plane_index].hcoeff[2]=   2;
2928         s->plane[plane_index].fast_mc= 1;
2929     }
2930
2931     common_init(avctx);
2932     alloc_blocks(s);
2933
2934     s->version=0;
2935
2936     s->m.avctx   = avctx;
2937     s->m.flags   = avctx->flags;
2938     s->m.bit_rate= avctx->bit_rate;
2939
2940     s->m.me.temp      =
2941     s->m.me.scratchpad= av_mallocz((avctx->width+64)*2*16*2*sizeof(uint8_t));
2942     s->m.me.map       = av_mallocz(ME_MAP_SIZE*sizeof(uint32_t));
2943     s->m.me.score_map = av_mallocz(ME_MAP_SIZE*sizeof(uint32_t));
2944     s->m.obmc_scratchpad= av_mallocz(MB_SIZE*MB_SIZE*12*sizeof(uint32_t));
2945     h263_encode_init(&s->m); //mv_penalty
2946
2947     s->max_ref_frames = FFMAX(FFMIN(avctx->refs, MAX_REF_FRAMES), 1);
2948
2949     if(avctx->flags&CODEC_FLAG_PASS1){
2950         if(!avctx->stats_out)
2951             avctx->stats_out = av_mallocz(256);
2952     }
2953     if((avctx->flags&CODEC_FLAG_PASS2) || !(avctx->flags&CODEC_FLAG_QSCALE)){
2954         if(ff_rate_control_init(&s->m) < 0)
2955             return -1;
2956     }
2957     s->pass1_rc= !(avctx->flags & (CODEC_FLAG_QSCALE|CODEC_FLAG_PASS2));
2958
2959     avctx->coded_frame= &s->current_picture;
2960     switch(avctx->pix_fmt){
2961 //    case PIX_FMT_YUV444P:
2962 //    case PIX_FMT_YUV422P:
2963     case PIX_FMT_YUV420P:
2964     case PIX_FMT_GRAY8:
2965 //    case PIX_FMT_YUV411P:
2966 //    case PIX_FMT_YUV410P:
2967         s->colorspace_type= 0;
2968         break;
2969 /*    case PIX_FMT_RGB32:
2970         s->colorspace= 1;
2971         break;*/
2972     default:
2973         av_log(avctx, AV_LOG_ERROR, "pixel format not supported\n");
2974         return -1;
2975     }
2976 //    avcodec_get_chroma_sub_sample(avctx->pix_fmt, &s->chroma_h_shift, &s->chroma_v_shift);
2977     s->chroma_h_shift= 1;
2978     s->chroma_v_shift= 1;
2979
2980     ff_set_cmp(&s->dsp, s->dsp.me_cmp, s->avctx->me_cmp);
2981     ff_set_cmp(&s->dsp, s->dsp.me_sub_cmp, s->avctx->me_sub_cmp);
2982
2983     s->avctx->get_buffer(s->avctx, &s->input_picture);
2984
2985     if(s->avctx->me_method == ME_ITER){
2986         int i;
2987         int size= s->b_width * s->b_height << 2*s->block_max_depth;
2988         for(i=0; i<s->max_ref_frames; i++){
2989             s->ref_mvs[i]= av_mallocz(size*sizeof(int16_t[2]));
2990             s->ref_scores[i]= av_mallocz(size*sizeof(uint32_t));
2991         }
2992     }
2993
2994     return 0;
2995 }
2996
2997 //near copy & paste from dsputil, FIXME
2998 static int pix_sum(uint8_t * pix, int line_size, int w)
2999 {
3000     int s, i, j;
3001
3002     s = 0;
3003     for (i = 0; i < w; i++) {
3004         for (j = 0; j < w; j++) {
3005             s += pix[0];
3006             pix ++;
3007         }
3008         pix += line_size - w;
3009     }
3010     return s;
3011 }
3012
3013 //near copy & paste from dsputil, FIXME
3014 static int pix_norm1(uint8_t * pix, int line_size, int w)
3015 {
3016     int s, i, j;
3017     uint32_t *sq = ff_squareTbl + 256;
3018
3019     s = 0;
3020     for (i = 0; i < w; i++) {
3021         for (j = 0; j < w; j ++) {
3022             s += sq[pix[0]];
3023             pix ++;
3024         }
3025         pix += line_size - w;
3026     }
3027     return s;
3028 }
3029
3030 //FIXME copy&paste
3031 #define P_LEFT P[1]
3032 #define P_TOP P[2]
3033 #define P_TOPRIGHT P[3]
3034 #define P_MEDIAN P[4]
3035 #define P_MV1 P[9]
3036 #define FLAG_QPEL   1 //must be 1
3037
3038 static int encode_q_branch(SnowContext *s, int level, int x, int y){
3039     uint8_t p_buffer[1024];
3040     uint8_t i_buffer[1024];
3041     uint8_t p_state[sizeof(s->block_state)];
3042     uint8_t i_state[sizeof(s->block_state)];
3043     RangeCoder pc, ic;
3044     uint8_t *pbbak= s->c.bytestream;
3045     uint8_t *pbbak_start= s->c.bytestream_start;
3046     int score, score2, iscore, i_len, p_len, block_s, sum, base_bits;
3047     const int w= s->b_width  << s->block_max_depth;
3048     const int h= s->b_height << s->block_max_depth;
3049     const int rem_depth= s->block_max_depth - level;
3050     const int index= (x + y*w) << rem_depth;
3051     const int block_w= 1<<(LOG2_MB_SIZE - level);
3052     int trx= (x+1)<<rem_depth;
3053     int try= (y+1)<<rem_depth;
3054     const BlockNode *left  = x ? &s->block[index-1] : &null_block;
3055     const BlockNode *top   = y ? &s->block[index-w] : &null_block;
3056     const BlockNode *right = trx<w ? &s->block[index+1] : &null_block;
3057     const BlockNode *bottom= try<h ? &s->block[index+w] : &null_block;
3058     const BlockNode *tl    = y && x ? &s->block[index-w-1] : left;
3059     const BlockNode *tr    = y && trx<w && ((x&1)==0 || level==0) ? &s->block[index-w+(1<<rem_depth)] : tl; //FIXME use lt
3060     int pl = left->color[0];
3061     int pcb= left->color[1];
3062     int pcr= left->color[2];
3063     int pmx, pmy;
3064     int mx=0, my=0;
3065     int l,cr,cb;
3066     const int stride= s->current_picture.linesize[0];
3067     const int uvstride= s->current_picture.linesize[1];
3068     uint8_t *current_data[3]= { s->input_picture.data[0] + (x + y*  stride)*block_w,
3069                                 s->input_picture.data[1] + (x + y*uvstride)*block_w/2,
3070                                 s->input_picture.data[2] + (x + y*uvstride)*block_w/2};
3071     int P[10][2];
3072     int16_t last_mv[3][2];
3073     int qpel= !!(s->avctx->flags & CODEC_FLAG_QPEL); //unused
3074     const int shift= 1+qpel;
3075     MotionEstContext *c= &s->m.me;
3076     int ref_context= av_log2(2*left->ref) + av_log2(2*top->ref);
3077     int mx_context= av_log2(2*FFABS(left->mx - top->mx));
3078     int my_context= av_log2(2*FFABS(left->my - top->my));
3079     int s_context= 2*left->level + 2*top->level + tl->level + tr->level;
3080     int ref, best_ref, ref_score, ref_mx, ref_my;
3081
3082     assert(sizeof(s->block_state) >= 256);
3083     if(s->keyframe){
3084         set_blocks(s, level, x, y, pl, pcb, pcr, 0, 0, 0, BLOCK_INTRA);
3085         return 0;
3086     }
3087
3088 //    clip predictors / edge ?
3089
3090     P_LEFT[0]= left->mx;
3091     P_LEFT[1]= left->my;
3092     P_TOP [0]= top->mx;
3093     P_TOP [1]= top->my;
3094     P_TOPRIGHT[0]= tr->mx;
3095     P_TOPRIGHT[1]= tr->my;
3096
3097     last_mv[0][0]= s->block[index].mx;
3098     last_mv[0][1]= s->block[index].my;
3099     last_mv[1][0]= right->mx;
3100     last_mv[1][1]= right->my;
3101     last_mv[2][0]= bottom->mx;
3102     last_mv[2][1]= bottom->my;
3103
3104     s->m.mb_stride=2;
3105     s->m.mb_x=
3106     s->m.mb_y= 0;
3107     c->skip= 0;
3108
3109     assert(c->  stride ==   stride);
3110     assert(c->uvstride == uvstride);
3111
3112     c->penalty_factor    = get_penalty_factor(s->lambda, s->lambda2, c->avctx->me_cmp);
3113     c->sub_penalty_factor= get_penalty_factor(s->lambda, s->lambda2, c->avctx->me_sub_cmp);
3114     c->mb_penalty_factor = get_penalty_factor(s->lambda, s->lambda2, c->avctx->mb_cmp);
3115     c->current_mv_penalty= c->mv_penalty[s->m.f_code=1] + MAX_MV;
3116
3117     c->xmin = - x*block_w - 16+3;
3118     c->ymin = - y*block_w - 16+3;
3119     c->xmax = - (x+1)*block_w + (w<<(LOG2_MB_SIZE - s->block_max_depth)) + 16-3;
3120     c->ymax = - (y+1)*block_w + (h<<(LOG2_MB_SIZE - s->block_max_depth)) + 16-3;
3121
3122     if(P_LEFT[0]     > (c->xmax<<shift)) P_LEFT[0]    = (c->xmax<<shift);
3123     if(P_LEFT[1]     > (c->ymax<<shift)) P_LEFT[1]    = (c->ymax<<shift);
3124     if(P_TOP[0]      > (c->xmax<<shift)) P_TOP[0]     = (c->xmax<<shift);
3125     if(P_TOP[1]      > (c->ymax<<shift)) P_TOP[1]     = (c->ymax<<shift);
3126     if(P_TOPRIGHT[0] < (c->xmin<<shift)) P_TOPRIGHT[0]= (c->xmin<<shift);
3127     if(P_TOPRIGHT[0] > (c->xmax<<shift)) P_TOPRIGHT[0]= (c->xmax<<shift); //due to pmx no clip
3128     if(P_TOPRIGHT[1] > (c->ymax<<shift)) P_TOPRIGHT[1]= (c->ymax<<shift);
3129
3130     P_MEDIAN[0]= mid_pred(P_LEFT[0], P_TOP[0], P_TOPRIGHT[0]);
3131     P_MEDIAN[1]= mid_pred(P_LEFT[1], P_TOP[1], P_TOPRIGHT[1]);
3132
3133     if (!y) {
3134         c->pred_x= P_LEFT[0];
3135         c->pred_y= P_LEFT[1];
3136     } else {
3137         c->pred_x = P_MEDIAN[0];
3138         c->pred_y = P_MEDIAN[1];
3139     }
3140
3141     score= INT_MAX;
3142     best_ref= 0;
3143     for(ref=0; ref<s->ref_frames; ref++){
3144         init_ref(c, current_data, s->last_picture[ref].data, NULL, block_w*x, block_w*y, 0);
3145
3146         ref_score= ff_epzs_motion_search(&s->m, &ref_mx, &ref_my, P, 0, /*ref_index*/ 0, last_mv,
3147                                          (1<<16)>>shift, level-LOG2_MB_SIZE+4, block_w);
3148
3149         assert(ref_mx >= c->xmin);
3150         assert(ref_mx <= c->xmax);
3151         assert(ref_my >= c->ymin);
3152         assert(ref_my <= c->ymax);
3153
3154         ref_score= c->sub_motion_search(&s->m, &ref_mx, &ref_my, ref_score, 0, 0, level-LOG2_MB_SIZE+4, block_w);
3155         ref_score= ff_get_mb_score(&s->m, ref_mx, ref_my, 0, 0, level-LOG2_MB_SIZE+4, block_w, 0);
3156         ref_score+= 2*av_log2(2*ref)*c->penalty_factor;
3157         if(s->ref_mvs[ref]){
3158             s->ref_mvs[ref][index][0]= ref_mx;
3159             s->ref_mvs[ref][index][1]= ref_my;
3160             s->ref_scores[ref][index]= ref_score;
3161         }
3162         if(score > ref_score){
3163             score= ref_score;
3164             best_ref= ref;
3165             mx= ref_mx;
3166             my= ref_my;
3167         }
3168     }
3169     //FIXME if mb_cmp != SSE then intra cannot be compared currently and mb_penalty vs. lambda2
3170
3171   //  subpel search
3172     base_bits= get_rac_count(&s->c) - 8*(s->c.bytestream - s->c.bytestream_start);
3173     pc= s->c;
3174     pc.bytestream_start=
3175     pc.bytestream= p_buffer; //FIXME end/start? and at the other stoo
3176     memcpy(p_state, s->block_state, sizeof(s->block_state));
3177
3178     if(level!=s->block_max_depth)
3179         put_rac(&pc, &p_state[4 + s_context], 1);
3180     put_rac(&pc, &p_state[1 + left->type + top->type], 0);
3181     if(s->ref_frames > 1)
3182         put_symbol(&pc, &p_state[128 + 1024 + 32*ref_context], best_ref, 0);
3183     pred_mv(s, &pmx, &pmy, best_ref, left, top, tr);
3184     put_symbol(&pc, &p_state[128 + 32*(mx_context + 16*!!best_ref)], mx - pmx, 1);
3185     put_symbol(&pc, &p_state[128 + 32*(my_context + 16*!!best_ref)], my - pmy, 1);
3186     p_len= pc.bytestream - pc.bytestream_start;
3187     score += (s->lambda2*(get_rac_count(&pc)-base_bits))>>FF_LAMBDA_SHIFT;
3188
3189     block_s= block_w*block_w;
3190     sum = pix_sum(current_data[0], stride, block_w);
3191     l= (sum + block_s/2)/block_s;
3192     iscore = pix_norm1(current_data[0], stride, block_w) - 2*l*sum + l*l*block_s;
3193
3194     block_s= block_w*block_w>>2;
3195     sum = pix_sum(current_data[1], uvstride, block_w>>1);
3196     cb= (sum + block_s/2)/block_s;
3197 //    iscore += pix_norm1(&current_mb[1][0], uvstride, block_w>>1) - 2*cb*sum + cb*cb*block_s;
3198     sum = pix_sum(current_data[2], uvstride, block_w>>1);
3199     cr= (sum + block_s/2)/block_s;
3200 //    iscore += pix_norm1(&current_mb[2][0], uvstride, block_w>>1) - 2*cr*sum + cr*cr*block_s;
3201
3202     ic= s->c;
3203     ic.bytestream_start=
3204     ic.bytestream= i_buffer; //FIXME end/start? and at the other stoo
3205     memcpy(i_state, s->block_state, sizeof(s->block_state));
3206     if(level!=s->block_max_depth)
3207         put_rac(&ic, &i_state[4 + s_context], 1);
3208     put_rac(&ic, &i_state[1 + left->type + top->type], 1);
3209     put_symbol(&ic, &i_state[32],  l-pl , 1);
3210     put_symbol(&ic, &i_state[64], cb-pcb, 1);
3211     put_symbol(&ic, &i_state[96], cr-pcr, 1);
3212     i_len= ic.bytestream - ic.bytestream_start;
3213     iscore += (s->lambda2*(get_rac_count(&ic)-base_bits))>>FF_LAMBDA_SHIFT;
3214
3215 //    assert(score==256*256*256*64-1);
3216     assert(iscore < 255*255*256 + s->lambda2*10);
3217     assert(iscore >= 0);
3218     assert(l>=0 && l<=255);
3219     assert(pl>=0 && pl<=255);
3220
3221     if(level==0){
3222         int varc= iscore >> 8;
3223         int vard= score >> 8;
3224         if (vard <= 64 || vard < varc)
3225             c->scene_change_score+= ff_sqrt(vard) - ff_sqrt(varc);
3226         else
3227             c->scene_change_score+= s->m.qscale;
3228     }
3229
3230     if(level!=s->block_max_depth){
3231         put_rac(&s->c, &s->block_state[4 + s_context], 0);
3232         score2 = encode_q_branch(s, level+1, 2*x+0, 2*y+0);
3233         score2+= encode_q_branch(s, level+1, 2*x+1, 2*y+0);
3234         score2+= encode_q_branch(s, level+1, 2*x+0, 2*y+1);
3235         score2+= encode_q_branch(s, level+1, 2*x+1, 2*y+1);
3236         score2+= s->lambda2>>FF_LAMBDA_SHIFT; //FIXME exact split overhead
3237
3238         if(score2 < score && score2 < iscore)
3239             return score2;
3240     }
3241
3242     if(iscore < score){
3243         pred_mv(s, &pmx, &pmy, 0, left, top, tr);
3244         memcpy(pbbak, i_buffer, i_len);
3245         s->c= ic;
3246         s->c.bytestream_start= pbbak_start;
3247         s->c.bytestream= pbbak + i_len;
3248         set_blocks(s, level, x, y, l, cb, cr, pmx, pmy, 0, BLOCK_INTRA);
3249         memcpy(s->block_state, i_state, sizeof(s->block_state));
3250         return iscore;
3251     }else{
3252         memcpy(pbbak, p_buffer, p_len);
3253         s->c= pc;
3254         s->c.bytestream_start= pbbak_start;
3255         s->c.bytestream= pbbak + p_len;
3256         set_blocks(s, level, x, y, pl, pcb, pcr, mx, my, best_ref, 0);
3257         memcpy(s->block_state, p_state, sizeof(s->block_state));
3258         return score;
3259     }
3260 }
3261
3262 static void encode_q_branch2(SnowContext *s, int level, int x, int y){
3263     const int w= s->b_width  << s->block_max_depth;
3264     const int rem_depth= s->block_max_depth - level;
3265     const int index= (x + y*w) << rem_depth;
3266     int trx= (x+1)<<rem_depth;
3267     BlockNode *b= &s->block[index];
3268     const BlockNode *left  = x ? &s->block[index-1] : &null_block;
3269     const BlockNode *top   = y ? &s->block[index-w] : &null_block;
3270     const BlockNode *tl    = y && x ? &s->block[index-w-1] : left;
3271     const BlockNode *tr    = y && trx<w && ((x&1)==0 || level==0) ? &s->block[index-w+(1<<rem_depth)] : tl; //FIXME use lt
3272     int pl = left->color[0];
3273     int pcb= left->color[1];
3274     int pcr= left->color[2];
3275     int pmx, pmy;
3276     int ref_context= av_log2(2*left->ref) + av_log2(2*top->ref);
3277     int mx_context= av_log2(2*FFABS(left->mx - top->mx)) + 16*!!b->ref;
3278     int my_context= av_log2(2*FFABS(left->my - top->my)) + 16*!!b->ref;
3279     int s_context= 2*left->level + 2*top->level + tl->level + tr->level;
3280
3281     if(s->keyframe){
3282         set_blocks(s, level, x, y, pl, pcb, pcr, 0, 0, 0, BLOCK_INTRA);
3283         return;
3284     }
3285
3286     if(level!=s->block_max_depth){
3287         if(same_block(b,b+1) && same_block(b,b+w) && same_block(b,b+w+1)){
3288             put_rac(&s->c, &s->block_state[4 + s_context], 1);
3289         }else{
3290             put_rac(&s->c, &s->block_state[4 + s_context], 0);
3291             encode_q_branch2(s, level+1, 2*x+0, 2*y+0);
3292             encode_q_branch2(s, level+1, 2*x+1, 2*y+0);
3293             encode_q_branch2(s, level+1, 2*x+0, 2*y+1);
3294             encode_q_branch2(s, level+1, 2*x+1, 2*y+1);
3295             return;
3296         }
3297     }
3298     if(b->type & BLOCK_INTRA){
3299         pred_mv(s, &pmx, &pmy, 0, left, top, tr);
3300         put_rac(&s->c, &s->block_state[1 + (left->type&1) + (top->type&1)], 1);
3301         put_symbol(&s->c, &s->block_state[32], b->color[0]-pl , 1);
3302         put_symbol(&s->c, &s->block_state[64], b->color[1]-pcb, 1);
3303         put_symbol(&s->c, &s->block_state[96], b->color[2]-pcr, 1);
3304         set_blocks(s, level, x, y, b->color[0], b->color[1], b->color[2], pmx, pmy, 0, BLOCK_INTRA);
3305     }else{
3306         pred_mv(s, &pmx, &pmy, b->ref, left, top, tr);
3307         put_rac(&s->c, &s->block_state[1 + (left->type&1) + (top->type&1)], 0);
3308         if(s->ref_frames > 1)
3309             put_symbol(&s->c, &s->block_state[128 + 1024 + 32*ref_context], b->ref, 0);
3310         put_symbol(&s->c, &s->block_state[128 + 32*mx_context], b->mx - pmx, 1);
3311         put_symbol(&s->c, &s->block_state[128 + 32*my_context], b->my - pmy, 1);
3312         set_blocks(s, level, x, y, pl, pcb, pcr, b->mx, b->my, b->ref, 0);
3313     }
3314 }
3315
3316 static int get_dc(SnowContext *s, int mb_x, int mb_y, int plane_index){
3317     int i, x2, y2;
3318     Plane *p= &s->plane[plane_index];
3319     const int block_size = MB_SIZE >> s->block_max_depth;
3320     const int block_w    = plane_index ? block_size/2 : block_size;
3321     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];
3322     const int obmc_stride= plane_index ? block_size : 2*block_size;
3323     const int ref_stride= s->current_picture.linesize[plane_index];
3324     uint8_t *src= s-> input_picture.data[plane_index];
3325     IDWTELEM *dst= (IDWTELEM*)s->m.obmc_scratchpad + plane_index*block_size*block_size*4; //FIXME change to unsigned
3326     const int b_stride = s->b_width << s->block_max_depth;
3327     const int w= p->width;
3328     const int h= p->height;
3329     int index= mb_x + mb_y*b_stride;
3330     BlockNode *b= &s->block[index];
3331     BlockNode backup= *b;
3332     int ab=0;
3333     int aa=0;
3334
3335     b->type|= BLOCK_INTRA;
3336     b->color[plane_index]= 0;
3337     memset(dst, 0, obmc_stride*obmc_stride*sizeof(IDWTELEM));
3338
3339     for(i=0; i<4; i++){
3340         int mb_x2= mb_x + (i &1) - 1;
3341         int mb_y2= mb_y + (i>>1) - 1;
3342         int x= block_w*mb_x2 + block_w/2;
3343         int y= block_w*mb_y2 + block_w/2;
3344
3345         add_yblock(s, 0, NULL, dst + ((i&1)+(i>>1)*obmc_stride)*block_w, NULL, obmc,
3346                     x, y, block_w, block_w, w, h, obmc_stride, ref_stride, obmc_stride, mb_x2, mb_y2, 0, 0, plane_index);
3347
3348         for(y2= FFMAX(y, 0); y2<FFMIN(h, y+block_w); y2++){
3349             for(x2= FFMAX(x, 0); x2<FFMIN(w, x+block_w); x2++){
3350                 int index= x2-(block_w*mb_x - block_w/2) + (y2-(block_w*mb_y - block_w/2))*obmc_stride;
3351                 int obmc_v= obmc[index];
3352                 int d;
3353                 if(y<0) obmc_v += obmc[index + block_w*obmc_stride];
3354                 if(x<0) obmc_v += obmc[index + block_w];
3355                 if(y+block_w>h) obmc_v += obmc[index - block_w*obmc_stride];
3356                 if(x+block_w>w) obmc_v += obmc[index - block_w];
3357                 //FIXME precalculate this or simplify it somehow else
3358
3359                 d = -dst[index] + (1<<(FRAC_BITS-1));
3360                 dst[index] = d;
3361                 ab += (src[x2 + y2*ref_stride] - (d>>FRAC_BITS)) * obmc_v;
3362                 aa += obmc_v * obmc_v; //FIXME precalculate this
3363             }
3364         }
3365     }
3366     *b= backup;
3367
3368     return av_clip(((ab<<LOG2_OBMC_MAX) + aa/2)/aa, 0, 255); //FIXME we should not need clipping
3369 }
3370
3371 static inline int get_block_bits(SnowContext *s, int x, int y, int w){
3372     const int b_stride = s->b_width << s->block_max_depth;
3373     const int b_height = s->b_height<< s->block_max_depth;
3374     int index= x + y*b_stride;
3375     const BlockNode *b     = &s->block[index];
3376     const BlockNode *left  = x ? &s->block[index-1] : &null_block;
3377     const BlockNode *top   = y ? &s->block[index-b_stride] : &null_block;
3378     const BlockNode *tl    = y && x ? &s->block[index-b_stride-1] : left;
3379     const BlockNode *tr    = y && x+w<b_stride ? &s->block[index-b_stride+w] : tl;
3380     int dmx, dmy;
3381 //  int mx_context= av_log2(2*FFABS(left->mx - top->mx));
3382 //  int my_context= av_log2(2*FFABS(left->my - top->my));
3383
3384     if(x<0 || x>=b_stride || y>=b_height)
3385         return 0;
3386 /*
3387 1            0      0
3388 01X          1-2    1
3389 001XX        3-6    2-3
3390 0001XXX      7-14   4-7
3391 00001XXXX   15-30   8-15
3392 */
3393 //FIXME try accurate rate
3394 //FIXME intra and inter predictors if surrounding blocks are not the same type
3395     if(b->type & BLOCK_INTRA){
3396         return 3+2*( av_log2(2*FFABS(left->color[0] - b->color[0]))
3397                    + av_log2(2*FFABS(left->color[1] - b->color[1]))
3398                    + av_log2(2*FFABS(left->color[2] - b->color[2])));
3399     }else{
3400         pred_mv(s, &dmx, &dmy, b->ref, left, top, tr);
3401         dmx-= b->mx;
3402         dmy-= b->my;
3403         return 2*(1 + av_log2(2*FFABS(dmx)) //FIXME kill the 2* can be merged in lambda
3404                     + av_log2(2*FFABS(dmy))
3405                     + av_log2(2*b->ref));
3406     }
3407 }
3408
3409 static int get_block_rd(SnowContext *s, int mb_x, int mb_y, int plane_index, const uint8_t *obmc_edged){
3410     Plane *p= &s->plane[plane_index];
3411     const int block_size = MB_SIZE >> s->block_max_depth;
3412     const int block_w    = plane_index ? block_size/2 : block_size;
3413     const int obmc_stride= plane_index ? block_size : 2*block_size;
3414     const int ref_stride= s->current_picture.linesize[plane_index];
3415     uint8_t *dst= s->current_picture.data[plane_index];
3416     uint8_t *src= s->  input_picture.data[plane_index];
3417     IDWTELEM *pred= (IDWTELEM*)s->m.obmc_scratchpad + plane_index*block_size*block_size*4;
3418     uint8_t *cur = s->scratchbuf;
3419     uint8_t tmp[ref_stride*(2*MB_SIZE+HTAPS_MAX-1)];
3420     const int b_stride = s->b_width << s->block_max_depth;
3421     const int b_height = s->b_height<< s->block_max_depth;
3422     const int w= p->width;
3423     const int h= p->height;
3424     int distortion;
3425     int rate= 0;
3426     const int penalty_factor= get_penalty_factor(s->lambda, s->lambda2, s->avctx->me_cmp);
3427     int sx= block_w*mb_x - block_w/2;
3428     int sy= block_w*mb_y - block_w/2;
3429     int x0= FFMAX(0,-sx);
3430     int y0= FFMAX(0,-sy);
3431     int x1= FFMIN(block_w*2, w-sx);
3432     int y1= FFMIN(block_w*2, h-sy);
3433     int i,x,y;
3434
3435     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);
3436
3437     for(y=y0; y<y1; y++){
3438         const uint8_t *obmc1= obmc_edged + y*obmc_stride;
3439         const IDWTELEM *pred1 = pred + y*obmc_stride;
3440         uint8_t *cur1 = cur + y*ref_stride;
3441         uint8_t *dst1 = dst + sx + (sy+y)*ref_stride;
3442         for(x=x0; x<x1; x++){
3443 #if FRAC_BITS >= LOG2_OBMC_MAX
3444             int v = (cur1[x] * obmc1[x]) << (FRAC_BITS - LOG2_OBMC_MAX);
3445 #else
3446             int v = (cur1[x] * obmc1[x] + (1<<(LOG2_OBMC_MAX - FRAC_BITS-1))) >> (LOG2_OBMC_MAX - FRAC_BITS);
3447 #endif
3448             v = (v + pred1[x]) >> FRAC_BITS;
3449             if(v&(~255)) v= ~(v>>31);
3450             dst1[x] = v;
3451         }
3452     }
3453
3454     /* copy the regions where obmc[] = (uint8_t)256 */
3455     if(LOG2_OBMC_MAX == 8
3456         && (mb_x == 0 || mb_x == b_stride-1)
3457         && (mb_y == 0 || mb_y == b_height-1)){
3458         if(mb_x == 0)
3459             x1 = block_w;
3460         else
3461             x0 = block_w;
3462         if(mb_y == 0)
3463             y1 = block_w;
3464         else
3465             y0 = block_w;
3466         for(y=y0; y<y1; y++)
3467             memcpy(dst + sx+x0 + (sy+y)*ref_stride, cur + x0 + y*ref_stride, x1-x0);
3468     }
3469
3470     if(block_w==16){
3471         /* FIXME rearrange dsputil to fit 32x32 cmp functions */
3472         /* FIXME check alignment of the cmp wavelet vs the encoding wavelet */
3473         /* FIXME cmps overlap but do not cover the wavelet's whole support.
3474          * So improving the score of one block is not strictly guaranteed
3475          * to improve the score of the whole frame, thus iterative motion
3476          * estimation does not always converge. */
3477         if(s->avctx->me_cmp == FF_CMP_W97)
3478             distortion = w97_32_c(&s->m, src + sx + sy*ref_stride, dst + sx + sy*ref_stride, ref_stride, 32);
3479         else if(s->avctx->me_cmp == FF_CMP_W53)
3480             distortion = w53_32_c(&s->m, src + sx + sy*ref_stride, dst + sx + sy*ref_stride, ref_stride, 32);
3481         else{
3482             distortion = 0;
3483             for(i=0; i<4; i++){
3484                 int off = sx+16*(i&1) + (sy+16*(i>>1))*ref_stride;
3485                 distortion += s->dsp.me_cmp[0](&s->m, src + off, dst + off, ref_stride, 16);
3486             }
3487         }
3488     }else{
3489         assert(block_w==8);
3490         distortion = s->dsp.me_cmp[0](&s->m, src + sx + sy*ref_stride, dst + sx + sy*ref_stride, ref_stride, block_w*2);
3491     }
3492
3493     if(plane_index==0){
3494         for(i=0; i<4; i++){
3495 /* ..RRr
3496  * .RXx.
3497  * rxx..
3498  */
3499             rate += get_block_bits(s, mb_x + (i&1) - (i>>1), mb_y + (i>>1), 1);
3500         }
3501         if(mb_x == b_stride-2)
3502             rate += get_block_bits(s, mb_x + 1, mb_y + 1, 1);
3503     }
3504     return distortion + rate*penalty_factor;
3505 }
3506
3507 static int get_4block_rd(SnowContext *s, int mb_x, int mb_y, int plane_index){
3508     int i, y2;
3509     Plane *p= &s->plane[plane_index];
3510     const int block_size = MB_SIZE >> s->block_max_depth;
3511     const int block_w    = plane_index ? block_size/2 : block_size;
3512     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];
3513     const int obmc_stride= plane_index ? block_size : 2*block_size;
3514     const int ref_stride= s->current_picture.linesize[plane_index];
3515     uint8_t *dst= s->current_picture.data[plane_index];
3516     uint8_t *src= s-> input_picture.data[plane_index];
3517     //FIXME zero_dst is const but add_yblock changes dst if add is 0 (this is never the case for dst=zero_dst
3518     // const has only been removed from zero_dst to suppress a warning
3519     static IDWTELEM zero_dst[4096]; //FIXME
3520     const int b_stride = s->b_width << s->block_max_depth;
3521     const int w= p->width;
3522     const int h= p->height;
3523     int distortion= 0;
3524     int rate= 0;
3525     const int penalty_factor= get_penalty_factor(s->lambda, s->lambda2, s->avctx->me_cmp);
3526
3527     for(i=0; i<9; i++){
3528         int mb_x2= mb_x + (i%3) - 1;
3529         int mb_y2= mb_y + (i/3) - 1;
3530         int x= block_w*mb_x2 + block_w/2;
3531         int y= block_w*mb_y2 + block_w/2;
3532
3533         add_yblock(s, 0, NULL, zero_dst, dst, obmc,
3534                    x, y, block_w, block_w, w, h, /*dst_stride*/0, ref_stride, obmc_stride, mb_x2, mb_y2, 1, 1, plane_index);
3535
3536         //FIXME find a cleaner/simpler way to skip the outside stuff
3537         for(y2= y; y2<0; y2++)
3538             memcpy(dst + x + y2*ref_stride, src + x + y2*ref_stride, block_w);
3539         for(y2= h; y2<y+block_w; y2++)
3540             memcpy(dst + x + y2*ref_stride, src + x + y2*ref_stride, block_w);
3541         if(x<0){
3542             for(y2= y; y2<y+block_w; y2++)
3543                 memcpy(dst + x + y2*ref_stride, src + x + y2*ref_stride, -x);
3544         }
3545         if(x+block_w > w){
3546             for(y2= y; y2<y+block_w; y2++)
3547                 memcpy(dst + w + y2*ref_stride, src + w + y2*ref_stride, x+block_w - w);
3548         }
3549
3550         assert(block_w== 8 || block_w==16);
3551         distortion += s->dsp.me_cmp[block_w==8](&s->m, src + x + y*ref_stride, dst + x + y*ref_stride, ref_stride, block_w);
3552     }
3553
3554     if(plane_index==0){
3555         BlockNode *b= &s->block[mb_x+mb_y*b_stride];
3556         int merged= same_block(b,b+1) && same_block(b,b+b_stride) && same_block(b,b+b_stride+1);
3557
3558 /* ..RRRr
3559  * .RXXx.
3560  * .RXXx.
3561  * rxxx.
3562  */
3563         if(merged)
3564             rate = get_block_bits(s, mb_x, mb_y, 2);
3565         for(i=merged?4:0; i<9; i++){
3566             static const int dxy[9][2] = {{0,0},{1,0},{0,1},{1,1},{2,0},{2,1},{-1,2},{0,2},{1,2}};
3567             rate += get_block_bits(s, mb_x + dxy[i][0], mb_y + dxy[i][1], 1);
3568         }
3569     }
3570     return distortion + rate*penalty_factor;
3571 }
3572
3573 static void ff_spatial_idwt_init(DWTCompose *cs, IDWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count){
3574     int level;
3575     for(level=decomposition_count-1; level>=0; level--){
3576         switch(type){
3577         case DWT_97: spatial_compose97i_init(cs+level, buffer, height>>level, stride<<level); break;
3578         case DWT_53: spatial_compose53i_init(cs+level, buffer, height>>level, stride<<level); break;
3579         }
3580     }
3581 }
3582
3583 static void ff_spatial_idwt_slice(DWTCompose *cs, IDWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count, int y){
3584     const int support = type==1 ? 3 : 5;
3585     int level;
3586     if(type==2) return;
3587
3588     for(level=decomposition_count-1; level>=0; level--){
3589         while(cs[level].y <= FFMIN((y>>level)+support, height>>level)){
3590             switch(type){
3591             case DWT_97: spatial_compose97i_dy(cs+level, buffer, width>>level, height>>level, stride<<level);
3592                 break;
3593             case DWT_53: spatial_compose53i_dy(cs+level, buffer, width>>level, height>>level, stride<<level);
3594                 break;
3595             }
3596         }
3597     }
3598 }
3599
3600 static void ff_spatial_idwt(IDWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count){
3601         DWTCompose cs[MAX_DECOMPOSITIONS];
3602         int y;
3603         ff_spatial_idwt_init(cs, buffer, width, height, stride, type, decomposition_count);
3604         for(y=0; y<height; y+=4)
3605             ff_spatial_idwt_slice(cs, buffer, width, height, stride, type, decomposition_count, y);
3606 }
3607
3608 static int encode_subband_c0run(SnowContext *s, SubBand *b, IDWTELEM *src, IDWTELEM *parent, int stride, int orientation){
3609     const int w= b->width;
3610     const int h= b->height;
3611     int x, y;
3612
3613     if(1){
3614         int run=0;
3615         int runs[w*h];
3616         int run_index=0;
3617         int max_index;
3618
3619         for(y=0; y<h; y++){
3620             for(x=0; x<w; x++){
3621                 int v, p=0;
3622                 int /*ll=0, */l=0, lt=0, t=0, rt=0;
3623                 v= src[x + y*stride];
3624
3625                 if(y){
3626                     t= src[x + (y-1)*stride];
3627                     if(x){
3628                         lt= src[x - 1 + (y-1)*stride];
3629                     }
3630                     if(x + 1 < w){
3631                         rt= src[x + 1 + (y-1)*stride];
3632                     }
3633                 }
3634                 if(x){
3635                     l= src[x - 1 + y*stride];
3636                     /*if(x > 1){
3637                         if(orientation==1) ll= src[y + (x-2)*stride];
3638                         else               ll= src[x - 2 + y*stride];
3639                     }*/
3640                 }
3641                 if(parent){
3642                     int px= x>>1;
3643                     int py= y>>1;
3644                     if(px<b->parent->width && py<b->parent->height)
3645                         p= parent[px + py*2*stride];
3646                 }
3647                 if(!(/*ll|*/l|lt|t|rt|p)){
3648                     if(v){
3649                         runs[run_index++]= run;
3650                         run=0;
3651                     }else{
3652                         run++;
3653                     }
3654                 }
3655             }
3656         }
3657         max_index= run_index;
3658         runs[run_index++]= run;
3659         run_index=0;
3660         run= runs[run_index++];
3661
3662         put_symbol2(&s->c, b->state[30], max_index, 0);
3663         if(run_index <= max_index)
3664             put_symbol2(&s->c, b->state[1], run, 3);
3665
3666         for(y=0; y<h; y++){
3667             if(s->c.bytestream_end - s->c.bytestream < w*40){
3668                 av_log(s->avctx, AV_LOG_ERROR, "encoded frame too large\n");
3669                 return -1;
3670             }
3671             for(x=0; x<w; x++){
3672                 int v, p=0;
3673                 int /*ll=0, */l=0, lt=0, t=0, rt=0;
3674                 v= src[x + y*stride];
3675
3676                 if(y){
3677                     t= src[x + (y-1)*stride];
3678                     if(x){
3679                         lt= src[x - 1 + (y-1)*stride];
3680                     }
3681                     if(x + 1 < w){
3682                         rt= src[x + 1 + (y-1)*stride];
3683                     }
3684                 }
3685                 if(x){
3686                     l= src[x - 1 + y*stride];
3687                     /*if(x > 1){
3688                         if(orientation==1) ll= src[y + (x-2)*stride];
3689                         else               ll= src[x - 2 + y*stride];
3690                     }*/
3691                 }
3692                 if(parent){
3693                     int px= x>>1;
3694                     int py= y>>1;
3695                     if(px<b->parent->width && py<b->parent->height)
3696                         p= parent[px + py*2*stride];
3697                 }
3698                 if(/*ll|*/l|lt|t|rt|p){
3699                     int context= av_log2(/*FFABS(ll) + */3*FFABS(l) + FFABS(lt) + 2*FFABS(t) + FFABS(rt) + FFABS(p));
3700
3701                     put_rac(&s->c, &b->state[0][context], !!v);
3702                 }else{
3703                     if(!run){
3704                         run= runs[run_index++];
3705
3706                         if(run_index <= max_index)
3707                             put_symbol2(&s->c, b->state[1], run, 3);
3708                         assert(v);
3709                     }else{
3710                         run--;
3711                         assert(!v);
3712                     }
3713                 }
3714                 if(v){
3715                     int context= av_log2(/*FFABS(ll) + */3*FFABS(l) + FFABS(lt) + 2*FFABS(t) + FFABS(rt) + FFABS(p));
3716                     int l2= 2*FFABS(l) + (l<0);
3717                     int t2= 2*FFABS(t) + (t<0);
3718
3719                     put_symbol2(&s->c, b->state[context + 2], FFABS(v)-1, context-4);
3720                     put_rac(&s->c, &b->state[0][16 + 1 + 3 + quant3bA[l2&0xFF] + 3*quant3bA[t2&0xFF]], v<0);
3721                 }
3722             }
3723         }
3724     }
3725     return 0;
3726 }
3727
3728 static int encode_subband(SnowContext *s, SubBand *b, IDWTELEM *src, IDWTELEM *parent, int stride, int orientation){
3729 //    encode_subband_qtree(s, b, src, parent, stride, orientation);
3730 //    encode_subband_z0run(s, b, src, parent, stride, orientation);
3731     return encode_subband_c0run(s, b, src, parent, stride, orientation);
3732 //    encode_subband_dzr(s, b, src, parent, stride, orientation);
3733 }
3734
3735 static av_always_inline int check_block(SnowContext *s, int mb_x, int mb_y, int p[3], int intra, const uint8_t *obmc_edged, int *best_rd){
3736     const int b_stride= s->b_width << s->block_max_depth;
3737     BlockNode *block= &s->block[mb_x + mb_y * b_stride];
3738     BlockNode backup= *block;
3739     int rd, index, value;
3740
3741     assert(mb_x>=0 && mb_y>=0);
3742     assert(mb_x<b_stride);
3743
3744     if(intra){
3745         block->color[0] = p[0];
3746         block->color[1] = p[1];
3747         block->color[2] = p[2];
3748         block->type |= BLOCK_INTRA;
3749     }else{
3750         index= (p[0] + 31*p[1]) & (ME_CACHE_SIZE-1);
3751         value= s->me_cache_generation + (p[0]>>10) + (p[1]<<6) + (block->ref<<12);
3752         if(s->me_cache[index] == value)
3753             return 0;
3754         s->me_cache[index]= value;
3755
3756         block->mx= p[0];
3757         block->my= p[1];
3758         block->type &= ~BLOCK_INTRA;
3759     }
3760
3761     rd= get_block_rd(s, mb_x, mb_y, 0, obmc_edged);
3762
3763 //FIXME chroma
3764     if(rd < *best_rd){
3765         *best_rd= rd;
3766         return 1;
3767     }else{
3768         *block= backup;
3769         return 0;
3770     }
3771 }
3772
3773 /* special case for int[2] args we discard afterwards,
3774  * fixes compilation problem with gcc 2.95 */
3775 static av_always_inline int check_block_inter(SnowContext *s, int mb_x, int mb_y, int p0, int p1, const uint8_t *obmc_edged, int *best_rd){
3776     int p[2] = {p0, p1};
3777     return check_block(s, mb_x, mb_y, p, 0, obmc_edged, best_rd);
3778 }
3779
3780 static av_always_inline int check_4block_inter(SnowContext *s, int mb_x, int mb_y, int p0, int p1, int ref, int *best_rd){
3781     const int b_stride= s->b_width << s->block_max_depth;
3782     BlockNode *block= &s->block[mb_x + mb_y * b_stride];
3783     BlockNode backup[4]= {block[0], block[1], block[b_stride], block[b_stride+1]};
3784     int rd, index, value;
3785
3786     assert(mb_x>=0 && mb_y>=0);
3787     assert(mb_x<b_stride);
3788     assert(((mb_x|mb_y)&1) == 0);
3789
3790     index= (p0 + 31*p1) & (ME_CACHE_SIZE-1);
3791     value= s->me_cache_generation + (p0>>10) + (p1<<6) + (block->ref<<12);
3792     if(s->me_cache[index] == value)
3793         return 0;
3794     s->me_cache[index]= value;
3795
3796     block->mx= p0;
3797     block->my= p1;
3798     block->ref= ref;
3799     block->type &= ~BLOCK_INTRA;
3800     block[1]= block[b_stride]= block[b_stride+1]= *block;
3801
3802     rd= get_4block_rd(s, mb_x, mb_y, 0);
3803
3804 //FIXME chroma
3805     if(rd < *best_rd){
3806         *best_rd= rd;
3807         return 1;
3808     }else{
3809         block[0]= backup[0];
3810         block[1]= backup[1];
3811         block[b_stride]= backup[2];
3812         block[b_stride+1]= backup[3];
3813         return 0;
3814     }
3815 }
3816
3817 static void iterative_me(SnowContext *s){
3818     int pass, mb_x, mb_y;
3819     const int b_width = s->b_width  << s->block_max_depth;
3820     const int b_height= s->b_height << s->block_max_depth;
3821     const int b_stride= b_width;
3822     int color[3];
3823
3824     {
3825         RangeCoder r = s->c;
3826         uint8_t state[sizeof(s->block_state)];
3827         memcpy(state, s->block_state, sizeof(s->block_state));
3828         for(mb_y= 0; mb_y<s->b_height; mb_y++)
3829             for(mb_x= 0; mb_x<s->b_width; mb_x++)
3830                 encode_q_branch(s, 0, mb_x, mb_y);
3831         s->c = r;
3832         memcpy(s->block_state, state, sizeof(s->block_state));
3833     }
3834
3835     for(pass=0; pass<25; pass++){
3836         int change= 0;
3837
3838         for(mb_y= 0; mb_y<b_height; mb_y++){
3839             for(mb_x= 0; mb_x<b_width; mb_x++){
3840                 int dia_change, i, j, ref;
3841                 int best_rd= INT_MAX, ref_rd;
3842                 BlockNode backup, ref_b;
3843                 const int index= mb_x + mb_y * b_stride;
3844                 BlockNode *block= &s->block[index];
3845                 BlockNode *tb =                   mb_y            ? &s->block[index-b_stride  ] : NULL;
3846                 BlockNode *lb = mb_x                              ? &s->block[index         -1] : NULL;
3847                 BlockNode *rb = mb_x+1<b_width                    ? &s->block[index         +1] : NULL;
3848                 BlockNode *bb =                   mb_y+1<b_height ? &s->block[index+b_stride  ] : NULL;
3849                 BlockNode *tlb= mb_x           && mb_y            ? &s->block[index-b_stride-1] : NULL;
3850                 BlockNode *trb= mb_x+1<b_width && mb_y            ? &s->block[index-b_stride+1] : NULL;
3851                 BlockNode *blb= mb_x           && mb_y+1<b_height ? &s->block[index+b_stride-1] : NULL;
3852                 BlockNode *brb= mb_x+1<b_width && mb_y+1<b_height ? &s->block[index+b_stride+1] : NULL;
3853                 const int b_w= (MB_SIZE >> s->block_max_depth);
3854                 uint8_t obmc_edged[b_w*2][b_w*2];
3855
3856                 if(pass && (block->type & BLOCK_OPT))
3857                     continue;
3858                 block->type |= BLOCK_OPT;
3859
3860                 backup= *block;
3861
3862                 if(!s->me_cache_generation)
3863                     memset(s->me_cache, 0, sizeof(s->me_cache));
3864                 s->me_cache_generation += 1<<22;
3865
3866                 //FIXME precalculate
3867                 {
3868                     int x, y;
3869                     memcpy(obmc_edged, obmc_tab[s->block_max_depth], b_w*b_w*4);
3870                     if(mb_x==0)
3871                         for(y=0; y<b_w*2; y++)
3872                             memset(obmc_edged[y], obmc_edged[y][0] + obmc_edged[y][b_w-1], b_w);
3873                     if(mb_x==b_stride-1)
3874                         for(y=0; y<b_w*2; y++)
3875                             memset(obmc_edged[y]+b_w, obmc_edged[y][b_w] + obmc_edged[y][b_w*2-1], b_w);
3876                     if(mb_y==0){
3877                         for(x=0; x<b_w*2; x++)
3878                             obmc_edged[0][x] += obmc_edged[b_w-1][x];
3879                         for(y=1; y<b_w; y++)
3880                             memcpy(obmc_edged[y], obmc_edged[0], b_w*2);
3881                     }
3882                     if(mb_y==b_height-1){
3883                         for(x=0; x<b_w*2; x++)
3884                             obmc_edged[b_w*2-1][x] += obmc_edged[b_w][x];
3885                         for(y=b_w; y<b_w*2-1; y++)
3886                             memcpy(obmc_edged[y], obmc_edged[b_w*2-1], b_w*2);
3887                     }
3888                 }
3889
3890                 //skip stuff outside the picture
3891                 if(mb_x==0 || mb_y==0 || mb_x==b_width-1 || mb_y==b_height-1){
3892                     uint8_t *src= s->  input_picture.data[0];
3893                     uint8_t *dst= s->current_picture.data[0];
3894                     const int stride= s->current_picture.linesize[0];
3895                     const int block_w= MB_SIZE >> s->block_max_depth;
3896                     const int sx= block_w*mb_x - block_w/2;
3897                     const int sy= block_w*mb_y - block_w/2;
3898                     const int w= s->plane[0].width;
3899                     const int h= s->plane[0].height;
3900                     int y;
3901
3902                     for(y=sy; y<0; y++)
3903                         memcpy(dst + sx + y*stride, src + sx + y*stride, block_w*2);
3904                     for(y=h; y<sy+block_w*2; y++)
3905                         memcpy(dst + sx + y*stride, src + sx + y*stride, block_w*2);
3906                     if(sx<0){
3907                         for(y=sy; y<sy+block_w*2; y++)
3908                             memcpy(dst + sx + y*stride, src + sx + y*stride, -sx);
3909                     }
3910                     if(sx+block_w*2 > w){
3911                         for(y=sy; y<sy+block_w*2; y++)
3912                             memcpy(dst + w + y*stride, src + w + y*stride, sx+block_w*2 - w);
3913                     }
3914                 }
3915
3916                 // intra(black) = neighbors' contribution to the current block
3917                 for(i=0; i<3; i++)
3918                     color[i]= get_dc(s, mb_x, mb_y, i);
3919
3920                 // get previous score (cannot be cached due to OBMC)
3921                 if(pass > 0 && (block->type&BLOCK_INTRA)){
3922                     int color0[3]= {block->color[0], block->color[1], block->color[2]};
3923                     check_block(s, mb_x, mb_y, color0, 1, *obmc_edged, &best_rd);
3924                 }else
3925                     check_block_inter(s, mb_x, mb_y, block->mx, block->my, *obmc_edged, &best_rd);
3926
3927                 ref_b= *block;
3928                 ref_rd= best_rd;
3929                 for(ref=0; ref < s->ref_frames; ref++){
3930                     int16_t (*mvr)[2]= &s->ref_mvs[ref][index];
3931                     if(s->ref_scores[ref][index] > s->ref_scores[ref_b.ref][index]*3/2) //FIXME tune threshold
3932                         continue;
3933                     block->ref= ref;
3934                     best_rd= INT_MAX;
3935
3936                     check_block_inter(s, mb_x, mb_y, mvr[0][0], mvr[0][1], *obmc_edged, &best_rd);
3937                     check_block_inter(s, mb_x, mb_y, 0, 0, *obmc_edged, &best_rd);
3938                     if(tb)
3939                         check_block_inter(s, mb_x, mb_y, mvr[-b_stride][0], mvr[-b_stride][1], *obmc_edged, &best_rd);
3940                     if(lb)
3941                         check_block_inter(s, mb_x, mb_y, mvr[-1][0], mvr[-1][1], *obmc_edged, &best_rd);
3942                     if(rb)
3943                         check_block_inter(s, mb_x, mb_y, mvr[1][0], mvr[1][1], *obmc_edged, &best_rd);
3944                     if(bb)
3945                         check_block_inter(s, mb_x, mb_y, mvr[b_stride][0], mvr[b_stride][1], *obmc_edged, &best_rd);
3946
3947                     /* fullpel ME */
3948                     //FIXME avoid subpel interpolation / round to nearest integer
3949                     do{
3950                         dia_change=0;
3951                         for(i=0; i<FFMAX(s->avctx->dia_size, 1); i++){
3952                             for(j=0; j<i; j++){
3953                                 dia_change |= check_block_inter(s, mb_x, mb_y, block->mx+4*(i-j), block->my+(4*j), *obmc_edged, &best_rd);
3954                                 dia_change |= check_block_inter(s, mb_x, mb_y, block->mx-4*(i-j), block->my-(4*j), *obmc_edged, &best_rd);
3955                                 dia_change |= check_block_inter(s, mb_x, mb_y, block->mx+4*(i-j), block->my-(4*j), *obmc_edged, &best_rd);
3956                                 dia_change |= check_block_inter(s, mb_x, mb_y, block->mx-4*(i-j), block->my+(4*j), *obmc_edged, &best_rd);
3957                             }
3958                         }
3959                     }while(dia_change);
3960                     /* subpel ME */
3961                     do{
3962                         static const int square[8][2]= {{+1, 0},{-1, 0},{ 0,+1},{ 0,-1},{+1,+1},{-1,-1},{+1,-1},{-1,+1},};
3963                         dia_change=0;
3964                         for(i=0; i<8; i++)
3965                             dia_change |= check_block_inter(s, mb_x, mb_y, block->mx+square[i][0], block->my+square[i][1], *obmc_edged, &best_rd);
3966                     }while(dia_change);
3967                     //FIXME or try the standard 2 pass qpel or similar
3968
3969                     mvr[0][0]= block->mx;
3970                     mvr[0][1]= block->my;
3971                     if(ref_rd > best_rd){
3972                         ref_rd= best_rd;
3973                         ref_b= *block;
3974                     }
3975                 }
3976                 best_rd= ref_rd;
3977                 *block= ref_b;
3978 #if 1
3979                 check_block(s, mb_x, mb_y, color, 1, *obmc_edged, &best_rd);
3980                 //FIXME RD style color selection
3981 #endif
3982                 if(!same_block(block, &backup)){
3983                     if(tb ) tb ->type &= ~BLOCK_OPT;
3984                     if(lb ) lb ->type &= ~BLOCK_OPT;
3985                     if(rb ) rb ->type &= ~BLOCK_OPT;
3986                     if(bb ) bb ->type &= ~BLOCK_OPT;
3987                     if(tlb) tlb->type &= ~BLOCK_OPT;
3988                     if(trb) trb->type &= ~BLOCK_OPT;
3989                     if(blb) blb->type &= ~BLOCK_OPT;
3990                     if(brb) brb->type &= ~BLOCK_OPT;
3991                     change ++;
3992                 }
3993             }
3994         }
3995         av_log(s->avctx, AV_LOG_ERROR, "pass:%d changed:%d\n", pass, change);
3996         if(!change)
3997             break;
3998     }
3999
4000     if(s->block_max_depth == 1){
4001         int change= 0;
4002         for(mb_y= 0; mb_y<b_height; mb_y+=2){
4003             for(mb_x= 0; mb_x<b_width; mb_x+=2){
4004                 int i;
4005                 int best_rd, init_rd;
4006                 const int index= mb_x + mb_y * b_stride;
4007                 BlockNode *b[4];
4008
4009                 b[0]= &s->block[index];
4010                 b[1]= b[0]+1;
4011                 b[2]= b[0]+b_stride;
4012                 b[3]= b[2]+1;
4013                 if(same_block(b[0], b[1]) &&
4014                    same_block(b[0], b[2]) &&
4015                    same_block(b[0], b[3]))
4016                     continue;
4017
4018                 if(!s->me_cache_generation)
4019                     memset(s->me_cache, 0, sizeof(s->me_cache));
4020                 s->me_cache_generation += 1<<22;
4021
4022                 init_rd= best_rd= get_4block_rd(s, mb_x, mb_y, 0);
4023
4024                 //FIXME more multiref search?
4025                 check_4block_inter(s, mb_x, mb_y,
4026                                    (b[0]->mx + b[1]->mx + b[2]->mx + b[3]->mx + 2) >> 2,
4027                                    (b[0]->my + b[1]->my + b[2]->my + b[3]->my + 2) >> 2, 0, &best_rd);
4028
4029                 for(i=0; i<4; i++)
4030                     if(!(b[i]->type&BLOCK_INTRA))
4031                         check_4block_inter(s, mb_x, mb_y, b[i]->mx, b[i]->my, b[i]->ref, &best_rd);
4032
4033                 if(init_rd != best_rd)
4034                     change++;
4035             }
4036         }
4037         av_log(s->avctx, AV_LOG_ERROR, "pass:4mv changed:%d\n", change*4);
4038     }
4039 }
4040
4041 static void encode_blocks(SnowContext *s, int search){
4042     int x, y;
4043     int w= s->b_width;
4044     int h= s->b_height;
4045
4046     if(s->avctx->me_method == ME_ITER && !s->keyframe && search)
4047         iterative_me(s);
4048
4049     for(y=0; y<h; y++){
4050         if(s->c.bytestream_end - s->c.bytestream < w*MB_SIZE*MB_SIZE*3){ //FIXME nicer limit
4051             av_log(s->avctx, AV_LOG_ERROR, "encoded frame too large\n");
4052             return;
4053         }
4054         for(x=0; x<w; x++){
4055             if(s->avctx->me_method == ME_ITER || !search)
4056                 encode_q_branch2(s, 0, x, y);
4057             else
4058                 encode_q_branch (s, 0, x, y);
4059         }
4060     }
4061 }
4062
4063 static void quantize(SnowContext *s, SubBand *b, IDWTELEM *dst, DWTELEM *src, int stride, int bias){
4064     const int w= b->width;
4065     const int h= b->height;
4066     const int qlog= av_clip(s->qlog + b->qlog, 0, QROOT*16);
4067     const int qmul= qexp[qlog&(QROOT-1)]<<((qlog>>QSHIFT) + ENCODER_EXTRA_BITS);
4068     int x,y, thres1, thres2;
4069
4070     if(s->qlog == LOSSLESS_QLOG){
4071         for(y=0; y<h; y++)
4072             for(x=0; x<w; x++)
4073                 dst[x + y*stride]= src[x + y*stride];
4074         return;
4075     }
4076
4077     bias= bias ? 0 : (3*qmul)>>3;
4078     thres1= ((qmul - bias)>>QEXPSHIFT) - 1;
4079     thres2= 2*thres1;
4080
4081     if(!bias){
4082         for(y=0; y<h; y++){
4083             for(x=0; x<w; x++){
4084                 int i= src[x + y*stride];
4085
4086                 if((unsigned)(i+thres1) > thres2){
4087                     if(i>=0){
4088                         i<<= QEXPSHIFT;
4089                         i/= qmul; //FIXME optimize
4090                         dst[x + y*stride]=  i;
4091                     }else{
4092                         i= -i;
4093                         i<<= QEXPSHIFT;
4094                         i/= qmul; //FIXME optimize
4095                         dst[x + y*stride]= -i;
4096                     }
4097                 }else
4098                     dst[x + y*stride]= 0;
4099             }
4100         }
4101     }else{
4102         for(y=0; y<h; y++){
4103             for(x=0; x<w; x++){
4104                 int i= src[x + y*stride];
4105
4106                 if((unsigned)(i+thres1) > thres2){
4107                     if(i>=0){
4108                         i<<= QEXPSHIFT;
4109                         i= (i + bias) / qmul; //FIXME optimize
4110                         dst[x + y*stride]=  i;
4111                     }else{
4112                         i= -i;
4113                         i<<= QEXPSHIFT;
4114                         i= (i + bias) / qmul; //FIXME optimize
4115                         dst[x + y*stride]= -i;
4116                     }
4117                 }else
4118                     dst[x + y*stride]= 0;
4119             }
4120         }
4121     }
4122 }
4123
4124 static void dequantize(SnowContext *s, SubBand *b, IDWTELEM *src, int stride){
4125     const int w= b->width;
4126     const int h= b->height;
4127     const int qlog= av_clip(s->qlog + b->qlog, 0, QROOT*16);
4128     const int qmul= qexp[qlog&(QROOT-1)]<<(qlog>>QSHIFT);
4129     const int qadd= (s->qbias*qmul)>>QBIAS_SHIFT;
4130     int x,y;
4131
4132     if(s->qlog == LOSSLESS_QLOG) return;
4133
4134     for(y=0; y<h; y++){
4135         for(x=0; x<w; x++){
4136             int i= src[x + y*stride];
4137             if(i<0){
4138                 src[x + y*stride]= -((-i*qmul + qadd)>>(QEXPSHIFT)); //FIXME try different bias
4139             }else if(i>0){
4140                 src[x + y*stride]=  (( i*qmul + qadd)>>(QEXPSHIFT));
4141             }
4142         }
4143     }
4144 }
4145
4146 static void decorrelate(SnowContext *s, SubBand *b, IDWTELEM *src, int stride, int inverse, int use_median){
4147     const int w= b->width;
4148     const int h= b->height;
4149     int x,y;
4150
4151     for(y=h-1; y>=0; y--){
4152         for(x=w-1; x>=0; x--){
4153             int i= x + y*stride;
4154
4155             if(x){
4156                 if(use_median){
4157                     if(y && x+1<w) src[i] -= mid_pred(src[i - 1], src[i - stride], src[i - stride + 1]);
4158                     else  src[i] -= src[i - 1];
4159                 }else{
4160                     if(y) src[i] -= mid_pred(src[i - 1], src[i - stride], src[i - 1] + src[i - stride] - src[i - 1 - stride]);