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