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