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