Snow 1pass ratecontrol
[ffmpeg.git] / libavcodec / snow.c
1 /*
2  * Copyright (C) 2004 Michael Niedermayer <michaelni@gmx.at>
3  *
4  * This library is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU Lesser General Public
6  * License as published by the Free Software Foundation; either
7  * version 2 of the License, or (at your option) any later version.
8  *
9  * This library is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12  * Lesser General Public License for more details.
13  *
14  * You should have received a copy of the GNU Lesser General Public
15  * License along with this library; if not, write to the Free Software
16  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17  */
18
19 #include "avcodec.h"
20 #include "common.h"
21 #include "dsputil.h"
22 #include "snow.h"
23
24 #include "rangecoder.h"
25
26 #include "mpegvideo.h"
27
28 #undef NDEBUG
29 #include <assert.h>
30
31 static const int8_t quant3[256]={
32  0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
33  1, 1, 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, 0,
48 };
49 static const int8_t quant3b[256]={
50  0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
51  1, 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 };
67 static const int8_t quant3bA[256]={
68  0, 0, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
69  1,-1, 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 };
85 static const int8_t quant5[256]={
86  0, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
87  2, 2, 2, 2, 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,-1,-1,-1,
102 };
103 static const int8_t quant7[256]={
104  0, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
105  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
106  2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
107  3, 3, 3, 3, 3, 3, 3, 3, 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,-2,-2,-2,-2,-2,-2,-2,
118 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
119 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-1,-1,
120 };
121 static const int8_t quant9[256]={
122  0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3,
123  3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
124  4, 4, 4, 4, 4, 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,-3,-3,-3,-3,
137 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-2,-2,-2,-2,-1,-1,
138 };
139 static const int8_t quant11[256]={
140  0, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4,
141  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
142  4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
143  5, 5, 5, 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,-4,-4,
154 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
155 -4,-4,-4,-4,-4,-3,-3,-3,-3,-3,-3,-3,-2,-2,-2,-1,
156 };
157 static const int8_t quant13[256]={
158  0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4,
159  4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
160  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
161  5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
162  6, 6, 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,-5,
171 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
172 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
173 -4,-4,-4,-4,-4,-4,-4,-4,-4,-3,-3,-3,-3,-2,-2,-1,
174 };
175
176 #if 0 //64*cubic
177 static const uint8_t obmc32[1024]={
178  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,
179  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,
180  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,
181  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,
182  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,
183  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,
184  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,
185  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,
186  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,
187  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,
188  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,
189  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,
190  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,
191  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,
192  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,
193  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,
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, 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,
196  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,
197  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,
198  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,
199  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,
200  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,
201  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,
202  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,
203  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,
204  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,
205  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,
206  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,
207  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,
208  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,
209  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,
210 //error:0.000022
211 };
212 static const uint8_t obmc16[256]={
213  0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
214  0, 1, 1, 2, 4, 5, 5, 6, 6, 5, 5, 4, 2, 1, 1, 0,
215  0, 1, 4, 6, 9,11,13,15,15,13,11, 9, 6, 4, 1, 0,
216  0, 2, 6,11,15,20,24,26,26,24,20,15,11, 6, 2, 0,
217  0, 4, 9,15,23,29,34,38,38,34,29,23,15, 9, 4, 0,
218  0, 5,11,20,29,38,45,49,49,45,38,29,20,11, 5, 0,
219  1, 5,13,24,34,45,53,57,57,53,45,34,24,13, 5, 1,
220  1, 6,15,26,38,49,57,62,62,57,49,38,26,15, 6, 1,
221  1, 6,15,26,38,49,57,62,62,57,49,38,26,15, 6, 1,
222  1, 5,13,24,34,45,53,57,57,53,45,34,24,13, 5, 1,
223  0, 5,11,20,29,38,45,49,49,45,38,29,20,11, 5, 0,
224  0, 4, 9,15,23,29,34,38,38,34,29,23,15, 9, 4, 0,
225  0, 2, 6,11,15,20,24,26,26,24,20,15,11, 6, 2, 0,
226  0, 1, 4, 6, 9,11,13,15,15,13,11, 9, 6, 4, 1, 0,
227  0, 1, 1, 2, 4, 5, 5, 6, 6, 5, 5, 4, 2, 1, 1, 0,
228  0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
229 //error:0.000033
230 };
231 #elif 1 // 64*linear
232 static const uint8_t obmc32[1024]={
233   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,
234   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,
235   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,
236   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,
237   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,
238   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,
239   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,
240   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,
241   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,
242   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,
243   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,
244   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,
245   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,
246   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,
247   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,
248   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,
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, 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,
251   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,
252   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,
253   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,
254   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,
255   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,
256   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,
257   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,
258   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,
259   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,
260   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,
261   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,
262   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,
263   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,
264   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,
265  //error:0.000020
266 };
267 static const uint8_t obmc16[256]={
268   0,  4,  4,  8,  8, 12, 12, 16, 16, 12, 12,  8,  8,  4,  4,  0,
269   4,  8, 16, 20, 28, 32, 40, 44, 44, 40, 32, 28, 20, 16,  8,  4,
270   4, 16, 24, 36, 44, 56, 64, 76, 76, 64, 56, 44, 36, 24, 16,  4,
271   8, 20, 36, 48, 64, 76, 92,104,104, 92, 76, 64, 48, 36, 20,  8,
272   8, 28, 44, 64, 80,100,116,136,136,116,100, 80, 64, 44, 28,  8,
273  12, 32, 56, 76,100,120,144,164,164,144,120,100, 76, 56, 32, 12,
274  12, 40, 64, 92,116,144,168,196,196,168,144,116, 92, 64, 40, 12,
275  16, 44, 76,104,136,164,196,224,224,196,164,136,104, 76, 44, 16,
276  16, 44, 76,104,136,164,196,224,224,196,164,136,104, 76, 44, 16,
277  12, 40, 64, 92,116,144,168,196,196,168,144,116, 92, 64, 40, 12,
278  12, 32, 56, 76,100,120,144,164,164,144,120,100, 76, 56, 32, 12,
279   8, 28, 44, 64, 80,100,116,136,136,116,100, 80, 64, 44, 28,  8,
280   8, 20, 36, 48, 64, 76, 92,104,104, 92, 76, 64, 48, 36, 20,  8,
281   4, 16, 24, 36, 44, 56, 64, 76, 76, 64, 56, 44, 36, 24, 16,  4,
282   4,  8, 16, 20, 28, 32, 40, 44, 44, 40, 32, 28, 20, 16,  8,  4,
283   0,  4,  4,  8,  8, 12, 12, 16, 16, 12, 12,  8,  8,  4,  4,  0,
284 //error:0.000015
285 };
286 #else //64*cos
287 static const uint8_t obmc32[1024]={
288  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,
289  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,
290  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,
291  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,
292  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,
293  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,
294  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,
295  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,
296  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,
297  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,
298  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,
299  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,
300  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,
301  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,
302  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,
303  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,
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,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,
306  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,
307  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,
308  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,
309  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,
310  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,
311  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,
312  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,
313  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,
314  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,
315  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,
316  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,
317  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,
318  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,
319  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,
320 //error:0.000022
321 };
322 static const uint8_t obmc16[256]={
323  0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0,
324  0, 0, 1, 2, 3, 4, 5, 5, 5, 5, 4, 3, 2, 1, 0, 0,
325  0, 1, 3, 6, 8,11,13,14,14,13,11, 8, 6, 3, 1, 0,
326  0, 2, 6,10,15,20,24,26,26,24,20,15,10, 6, 2, 0,
327  0, 3, 8,16,23,30,35,38,38,35,30,23,16, 8, 3, 0,
328  1, 4,11,20,30,39,46,49,49,46,39,30,20,11, 4, 1,
329  1, 5,13,24,35,46,54,58,58,54,46,35,24,13, 5, 1,
330  0, 5,14,26,38,49,58,63,63,58,49,38,26,14, 5, 0,
331  0, 5,14,26,38,49,58,63,63,58,49,38,26,14, 5, 0,
332  1, 5,13,24,35,46,54,58,58,54,46,35,24,13, 5, 1,
333  1, 4,11,20,30,39,46,49,49,46,39,30,20,11, 4, 1,
334  0, 3, 8,16,23,30,35,38,38,35,30,23,16, 8, 3, 0,
335  0, 2, 6,10,15,20,24,26,26,24,20,15,10, 6, 2, 0,
336  0, 1, 3, 6, 8,11,13,14,14,13,11, 8, 6, 3, 1, 0,
337  0, 0, 1, 2, 3, 4, 5, 5, 5, 5, 4, 3, 2, 1, 0, 0,
338  0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0,
339 //error:0.000022
340 };
341 #endif
342
343 //linear *64
344 static const uint8_t obmc8[64]={
345   4, 12, 20, 28, 28, 20, 12,  4,
346  12, 36, 60, 84, 84, 60, 36, 12,
347  20, 60,100,140,140,100, 60, 20,
348  28, 84,140,196,196,140, 84, 28,
349  28, 84,140,196,196,140, 84, 28,
350  20, 60,100,140,140,100, 60, 20,
351  12, 36, 60, 84, 84, 60, 36, 12,
352   4, 12, 20, 28, 28, 20, 12,  4,
353 //error:0.000000
354 };
355
356 //linear *64
357 static const uint8_t obmc4[16]={
358  16, 48, 48, 16,
359  48,144,144, 48,
360  48,144,144, 48,
361  16, 48, 48, 16,
362 //error:0.000000
363 };
364
365 static const uint8_t *obmc_tab[4]={
366     obmc32, obmc16, obmc8, obmc4
367 };
368
369 typedef struct BlockNode{
370     int16_t mx;
371     int16_t my;
372     uint8_t color[3];
373     uint8_t type;
374 //#define TYPE_SPLIT    1
375 #define BLOCK_INTRA   1
376 #define BLOCK_OPT     2
377 //#define TYPE_NOCOLOR  4
378     uint8_t level; //FIXME merge into type?
379 }BlockNode;
380
381 static const BlockNode null_block= { //FIXME add border maybe
382     .color= {128,128,128},
383     .mx= 0,
384     .my= 0,
385     .type= 0,
386     .level= 0,
387 };
388
389 #define LOG2_MB_SIZE 4
390 #define MB_SIZE (1<<LOG2_MB_SIZE)
391
392 typedef struct x_and_coeff{
393     int16_t x;
394     uint16_t coeff;
395 } x_and_coeff;
396
397 typedef struct SubBand{
398     int level;
399     int stride;
400     int width;
401     int height;
402     int qlog;                                   ///< log(qscale)/log[2^(1/6)]
403     DWTELEM *buf;
404     int buf_x_offset;
405     int buf_y_offset;
406     int stride_line; ///< Stride measured in lines, not pixels.
407     x_and_coeff * x_coeff;
408     struct SubBand *parent;
409     uint8_t state[/*7*2*/ 7 + 512][32];
410 }SubBand;
411
412 typedef struct Plane{
413     int width;
414     int height;
415     SubBand band[MAX_DECOMPOSITIONS][4];
416 }Plane;
417
418 typedef struct SnowContext{
419 //    MpegEncContext m; // needed for motion estimation, should not be used for anything else, the idea is to make the motion estimation eventually independant of MpegEncContext, so this will be removed then (FIXME/XXX)
420
421     AVCodecContext *avctx;
422     RangeCoder c;
423     DSPContext dsp;
424     AVFrame new_picture;
425     AVFrame input_picture;              ///< new_picture with the internal linesizes
426     AVFrame current_picture;
427     AVFrame last_picture;
428     AVFrame mconly_picture;
429 //     uint8_t q_context[16];
430     uint8_t header_state[32];
431     uint8_t block_state[128 + 32*128];
432     int keyframe;
433     int always_reset;
434     int version;
435     int spatial_decomposition_type;
436     int temporal_decomposition_type;
437     int spatial_decomposition_count;
438     int temporal_decomposition_count;
439     DWTELEM *spatial_dwt_buffer;
440     int colorspace_type;
441     int chroma_h_shift;
442     int chroma_v_shift;
443     int spatial_scalability;
444     int qlog;
445     int lambda;
446     int lambda2;
447     int pass1_rc;
448     int mv_scale;
449     int qbias;
450 #define QBIAS_SHIFT 3
451     int b_width;
452     int b_height;
453     int block_max_depth;
454     Plane plane[MAX_PLANES];
455     BlockNode *block;
456 #define ME_CACHE_SIZE 1024
457     int me_cache[ME_CACHE_SIZE];
458     int me_cache_generation;
459     slice_buffer sb;
460
461     MpegEncContext m; // needed for motion estimation, should not be used for anything else, the idea is to make the motion estimation eventually independant of MpegEncContext, so this will be removed then (FIXME/XXX)
462 }SnowContext;
463
464 typedef struct {
465     DWTELEM *b0;
466     DWTELEM *b1;
467     DWTELEM *b2;
468     DWTELEM *b3;
469     int y;
470 } dwt_compose_t;
471
472 #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)))
473 //#define slice_buffer_get_line(slice_buf, line_num) (slice_buffer_load_line((slice_buf), (line_num)))
474
475 static void iterative_me(SnowContext *s);
476
477 static void slice_buffer_init(slice_buffer * buf, int line_count, int max_allocated_lines, int line_width, DWTELEM * base_buffer)
478 {
479     int i;
480
481     buf->base_buffer = base_buffer;
482     buf->line_count = line_count;
483     buf->line_width = line_width;
484     buf->data_count = max_allocated_lines;
485     buf->line = (DWTELEM * *) av_mallocz (sizeof(DWTELEM *) * line_count);
486     buf->data_stack = (DWTELEM * *) av_malloc (sizeof(DWTELEM *) * max_allocated_lines);
487
488     for (i = 0; i < max_allocated_lines; i++)
489     {
490       buf->data_stack[i] = (DWTELEM *) av_malloc (sizeof(DWTELEM) * line_width);
491     }
492
493     buf->data_stack_top = max_allocated_lines - 1;
494 }
495
496 static DWTELEM * slice_buffer_load_line(slice_buffer * buf, int line)
497 {
498     int offset;
499     DWTELEM * buffer;
500
501 //  av_log(NULL, AV_LOG_DEBUG, "Cache hit: %d\n", line);
502
503     assert(buf->data_stack_top >= 0);
504 //  assert(!buf->line[line]);
505     if (buf->line[line])
506         return buf->line[line];
507
508     offset = buf->line_width * line;
509     buffer = buf->data_stack[buf->data_stack_top];
510     buf->data_stack_top--;
511     buf->line[line] = buffer;
512
513 //  av_log(NULL, AV_LOG_DEBUG, "slice_buffer_load_line: line: %d remaining: %d\n", line, buf->data_stack_top + 1);
514
515     return buffer;
516 }
517
518 static void slice_buffer_release(slice_buffer * buf, int line)
519 {
520     int offset;
521     DWTELEM * buffer;
522
523     assert(line >= 0 && line < buf->line_count);
524     assert(buf->line[line]);
525
526     offset = buf->line_width * line;
527     buffer = buf->line[line];
528     buf->data_stack_top++;
529     buf->data_stack[buf->data_stack_top] = buffer;
530     buf->line[line] = NULL;
531
532 //  av_log(NULL, AV_LOG_DEBUG, "slice_buffer_release: line: %d remaining: %d\n", line, buf->data_stack_top + 1);
533 }
534
535 static void slice_buffer_flush(slice_buffer * buf)
536 {
537     int i;
538     for (i = 0; i < buf->line_count; i++)
539     {
540         if (buf->line[i])
541         {
542 //      av_log(NULL, AV_LOG_DEBUG, "slice_buffer_flush: line: %d \n", i);
543             slice_buffer_release(buf, i);
544         }
545     }
546 }
547
548 static void slice_buffer_destroy(slice_buffer * buf)
549 {
550     int i;
551     slice_buffer_flush(buf);
552
553     for (i = buf->data_count - 1; i >= 0; i--)
554     {
555         assert(buf->data_stack[i]);
556         av_freep(&buf->data_stack[i]);
557     }
558     assert(buf->data_stack);
559     av_freep(&buf->data_stack);
560     assert(buf->line);
561     av_freep(&buf->line);
562 }
563
564 #ifdef __sgi
565 // Avoid a name clash on SGI IRIX
566 #undef qexp
567 #endif
568 #define QEXPSHIFT (7-FRAC_BITS+8) //FIXME try to change this to 0
569 static uint8_t qexp[QROOT];
570
571 static inline int mirror(int v, int m){
572     while((unsigned)v > (unsigned)m){
573         v=-v;
574         if(v<0) v+= 2*m;
575     }
576     return v;
577 }
578
579 static inline void put_symbol(RangeCoder *c, uint8_t *state, int v, int is_signed){
580     int i;
581
582     if(v){
583         const int a= ABS(v);
584         const int e= av_log2(a);
585 #if 1
586         const int el= FFMIN(e, 10);
587         put_rac(c, state+0, 0);
588
589         for(i=0; i<el; i++){
590             put_rac(c, state+1+i, 1);  //1..10
591         }
592         for(; i<e; i++){
593             put_rac(c, state+1+9, 1);  //1..10
594         }
595         put_rac(c, state+1+FFMIN(i,9), 0);
596
597         for(i=e-1; i>=el; i--){
598             put_rac(c, state+22+9, (a>>i)&1); //22..31
599         }
600         for(; i>=0; i--){
601             put_rac(c, state+22+i, (a>>i)&1); //22..31
602         }
603
604         if(is_signed)
605             put_rac(c, state+11 + el, v < 0); //11..21
606 #else
607
608         put_rac(c, state+0, 0);
609         if(e<=9){
610             for(i=0; i<e; i++){
611                 put_rac(c, state+1+i, 1);  //1..10
612             }
613             put_rac(c, state+1+i, 0);
614
615             for(i=e-1; i>=0; i--){
616                 put_rac(c, state+22+i, (a>>i)&1); //22..31
617             }
618
619             if(is_signed)
620                 put_rac(c, state+11 + e, v < 0); //11..21
621         }else{
622             for(i=0; i<e; i++){
623                 put_rac(c, state+1+FFMIN(i,9), 1);  //1..10
624             }
625             put_rac(c, state+1+FFMIN(i,9), 0);
626
627             for(i=e-1; i>=0; i--){
628                 put_rac(c, state+22+FFMIN(i,9), (a>>i)&1); //22..31
629             }
630
631             if(is_signed)
632                 put_rac(c, state+11 + FFMIN(e,10), v < 0); //11..21
633         }
634 #endif
635     }else{
636         put_rac(c, state+0, 1);
637     }
638 }
639
640 static inline int get_symbol(RangeCoder *c, uint8_t *state, int is_signed){
641     if(get_rac(c, state+0))
642         return 0;
643     else{
644         int i, e, a;
645         e= 0;
646         while(get_rac(c, state+1 + FFMIN(e,9))){ //1..10
647             e++;
648         }
649
650         a= 1;
651         for(i=e-1; i>=0; i--){
652             a += a + get_rac(c, state+22 + FFMIN(i,9)); //22..31
653         }
654
655         if(is_signed && get_rac(c, state+11 + FFMIN(e,10))) //11..21
656             return -a;
657         else
658             return a;
659     }
660 }
661
662 static inline void put_symbol2(RangeCoder *c, uint8_t *state, int v, int log2){
663     int i;
664     int r= log2>=0 ? 1<<log2 : 1;
665
666     assert(v>=0);
667     assert(log2>=-4);
668
669     while(v >= r){
670         put_rac(c, state+4+log2, 1);
671         v -= r;
672         log2++;
673         if(log2>0) r+=r;
674     }
675     put_rac(c, state+4+log2, 0);
676
677     for(i=log2-1; i>=0; i--){
678         put_rac(c, state+31-i, (v>>i)&1);
679     }
680 }
681
682 static inline int get_symbol2(RangeCoder *c, uint8_t *state, int log2){
683     int i;
684     int r= log2>=0 ? 1<<log2 : 1;
685     int v=0;
686
687     assert(log2>=-4);
688
689     while(get_rac(c, state+4+log2)){
690         v+= r;
691         log2++;
692         if(log2>0) r+=r;
693     }
694
695     for(i=log2-1; i>=0; i--){
696         v+= get_rac(c, state+31-i)<<i;
697     }
698
699     return v;
700 }
701
702 static always_inline void lift(DWTELEM *dst, DWTELEM *src, DWTELEM *ref, int dst_step, int src_step, int ref_step, int width, int mul, int add, int shift, int highpass, int inverse){
703     const int mirror_left= !highpass;
704     const int mirror_right= (width&1) ^ highpass;
705     const int w= (width>>1) - 1 + (highpass & width);
706     int i;
707
708 #define LIFT(src, ref, inv) ((src) + ((inv) ? - (ref) : + (ref)))
709     if(mirror_left){
710         dst[0] = LIFT(src[0], ((mul*2*ref[0]+add)>>shift), inverse);
711         dst += dst_step;
712         src += src_step;
713     }
714
715     for(i=0; i<w; i++){
716         dst[i*dst_step] = LIFT(src[i*src_step], ((mul*(ref[i*ref_step] + ref[(i+1)*ref_step])+add)>>shift), inverse);
717     }
718
719     if(mirror_right){
720         dst[w*dst_step] = LIFT(src[w*src_step], ((mul*2*ref[w*ref_step]+add)>>shift), inverse);
721     }
722 }
723
724 #ifndef lift5
725 static always_inline void lift5(DWTELEM *dst, DWTELEM *src, DWTELEM *ref, int dst_step, int src_step, int ref_step, int width, int mul, int add, int shift, int highpass, int inverse){
726     const int mirror_left= !highpass;
727     const int mirror_right= (width&1) ^ highpass;
728     const int w= (width>>1) - 1 + (highpass & width);
729     int i;
730
731     if(mirror_left){
732         int r= 3*2*ref[0];
733         r += r>>4;
734         r += r>>8;
735         dst[0] = LIFT(src[0], ((r+add)>>shift), inverse);
736         dst += dst_step;
737         src += src_step;
738     }
739
740     for(i=0; i<w; i++){
741         int r= 3*(ref[i*ref_step] + ref[(i+1)*ref_step]);
742         r += r>>4;
743         r += r>>8;
744         dst[i*dst_step] = LIFT(src[i*src_step], ((r+add)>>shift), inverse);
745     }
746
747     if(mirror_right){
748         int r= 3*2*ref[w*ref_step];
749         r += r>>4;
750         r += r>>8;
751         dst[w*dst_step] = LIFT(src[w*src_step], ((r+add)>>shift), inverse);
752     }
753 }
754 #endif
755
756 #ifndef liftS
757 static always_inline void liftS(DWTELEM *dst, DWTELEM *src, DWTELEM *ref, int dst_step, int src_step, int ref_step, int width, int mul, int add, int shift, int highpass, int inverse){
758     const int mirror_left= !highpass;
759     const int mirror_right= (width&1) ^ highpass;
760     const int w= (width>>1) - 1 + (highpass & width);
761     int i;
762
763     assert(shift == 4);
764 #define LIFTS(src, ref, inv) ((inv) ? (src) - (((ref) - 4*(src))>>shift): (16*4*(src) + 4*(ref) + 8 + (5<<27))/(5*16) - (1<<23))
765     if(mirror_left){
766         dst[0] = LIFTS(src[0], mul*2*ref[0]+add, inverse);
767         dst += dst_step;
768         src += src_step;
769     }
770
771     for(i=0; i<w; i++){
772         dst[i*dst_step] = LIFTS(src[i*src_step], mul*(ref[i*ref_step] + ref[(i+1)*ref_step])+add, inverse);
773     }
774
775     if(mirror_right){
776         dst[w*dst_step] = LIFTS(src[w*src_step], mul*2*ref[w*ref_step]+add, inverse);
777     }
778 }
779 #endif
780
781
782 static void inplace_lift(DWTELEM *dst, int width, int *coeffs, int n, int shift, int start, int inverse){
783     int x, i;
784
785     for(x=start; x<width; x+=2){
786         int64_t sum=0;
787
788         for(i=0; i<n; i++){
789             int x2= x + 2*i - n + 1;
790             if     (x2<     0) x2= -x2;
791             else if(x2>=width) x2= 2*width-x2-2;
792             sum += coeffs[i]*(int64_t)dst[x2];
793         }
794         if(inverse) dst[x] -= (sum + (1<<shift)/2)>>shift;
795         else        dst[x] += (sum + (1<<shift)/2)>>shift;
796     }
797 }
798
799 static void inplace_liftV(DWTELEM *dst, int width, int height, int stride, int *coeffs, int n, int shift, int start, int inverse){
800     int x, y, i;
801     for(y=start; y<height; y+=2){
802         for(x=0; x<width; x++){
803             int64_t sum=0;
804
805             for(i=0; i<n; i++){
806                 int y2= y + 2*i - n + 1;
807                 if     (y2<      0) y2= -y2;
808                 else if(y2>=height) y2= 2*height-y2-2;
809                 sum += coeffs[i]*(int64_t)dst[x + y2*stride];
810             }
811             if(inverse) dst[x + y*stride] -= (sum + (1<<shift)/2)>>shift;
812             else        dst[x + y*stride] += (sum + (1<<shift)/2)>>shift;
813         }
814     }
815 }
816
817 #define SCALEX 1
818 #define LX0 0
819 #define LX1 1
820
821 #if 0 // more accurate 9/7
822 #define N1 2
823 #define SHIFT1 14
824 #define COEFFS1 (int[]){-25987,-25987}
825 #define N2 2
826 #define SHIFT2 19
827 #define COEFFS2 (int[]){-27777,-27777}
828 #define N3 2
829 #define SHIFT3 15
830 #define COEFFS3 (int[]){28931,28931}
831 #define N4 2
832 #define SHIFT4 15
833 #define COEFFS4 (int[]){14533,14533}
834 #elif 1 // 13/7 CRF
835 #define N1 4
836 #define SHIFT1 4
837 #define COEFFS1 (int[]){1,-9,-9,1}
838 #define N2 4
839 #define SHIFT2 4
840 #define COEFFS2 (int[]){-1,5,5,-1}
841 #define N3 0
842 #define SHIFT3 1
843 #define COEFFS3 NULL
844 #define N4 0
845 #define SHIFT4 1
846 #define COEFFS4 NULL
847 #elif 1 // 3/5
848 #define LX0 1
849 #define LX1 0
850 #define SCALEX 0.5
851 #define N1 2
852 #define SHIFT1 1
853 #define COEFFS1 (int[]){1,1}
854 #define N2 2
855 #define SHIFT2 2
856 #define COEFFS2 (int[]){-1,-1}
857 #define N3 0
858 #define SHIFT3 0
859 #define COEFFS3 NULL
860 #define N4 0
861 #define SHIFT4 0
862 #define COEFFS4 NULL
863 #elif 1 // 11/5
864 #define N1 0
865 #define SHIFT1 1
866 #define COEFFS1 NULL
867 #define N2 2
868 #define SHIFT2 2
869 #define COEFFS2 (int[]){-1,-1}
870 #define N3 2
871 #define SHIFT3 0
872 #define COEFFS3 (int[]){-1,-1}
873 #define N4 4
874 #define SHIFT4 7
875 #define COEFFS4 (int[]){-5,29,29,-5}
876 #define SCALEX 4
877 #elif 1 // 9/7 CDF
878 #define N1 2
879 #define SHIFT1 7
880 #define COEFFS1 (int[]){-203,-203}
881 #define N2 2
882 #define SHIFT2 12
883 #define COEFFS2 (int[]){-217,-217}
884 #define N3 2
885 #define SHIFT3 7
886 #define COEFFS3 (int[]){113,113}
887 #define N4 2
888 #define SHIFT4 9
889 #define COEFFS4 (int[]){227,227}
890 #define SCALEX 1
891 #elif 1 // 7/5 CDF
892 #define N1 0
893 #define SHIFT1 1
894 #define COEFFS1 NULL
895 #define N2 2
896 #define SHIFT2 2
897 #define COEFFS2 (int[]){-1,-1}
898 #define N3 2
899 #define SHIFT3 0
900 #define COEFFS3 (int[]){-1,-1}
901 #define N4 2
902 #define SHIFT4 4
903 #define COEFFS4 (int[]){3,3}
904 #elif 1 // 9/7 MN
905 #define N1 4
906 #define SHIFT1 4
907 #define COEFFS1 (int[]){1,-9,-9,1}
908 #define N2 2
909 #define SHIFT2 2
910 #define COEFFS2 (int[]){1,1}
911 #define N3 0
912 #define SHIFT3 1
913 #define COEFFS3 NULL
914 #define N4 0
915 #define SHIFT4 1
916 #define COEFFS4 NULL
917 #else // 13/7 CRF
918 #define N1 4
919 #define SHIFT1 4
920 #define COEFFS1 (int[]){1,-9,-9,1}
921 #define N2 4
922 #define SHIFT2 4
923 #define COEFFS2 (int[]){-1,5,5,-1}
924 #define N3 0
925 #define SHIFT3 1
926 #define COEFFS3 NULL
927 #define N4 0
928 #define SHIFT4 1
929 #define COEFFS4 NULL
930 #endif
931 static void horizontal_decomposeX(DWTELEM *b, int width){
932     DWTELEM temp[width];
933     const int width2= width>>1;
934     const int w2= (width+1)>>1;
935     int x;
936
937     inplace_lift(b, width, COEFFS1, N1, SHIFT1, LX1, 0);
938     inplace_lift(b, width, COEFFS2, N2, SHIFT2, LX0, 0);
939     inplace_lift(b, width, COEFFS3, N3, SHIFT3, LX1, 0);
940     inplace_lift(b, width, COEFFS4, N4, SHIFT4, LX0, 0);
941
942     for(x=0; x<width2; x++){
943         temp[x   ]= b[2*x    ];
944         temp[x+w2]= b[2*x + 1];
945     }
946     if(width&1)
947         temp[x   ]= b[2*x    ];
948     memcpy(b, temp, width*sizeof(int));
949 }
950
951 static void horizontal_composeX(DWTELEM *b, int width){
952     DWTELEM temp[width];
953     const int width2= width>>1;
954     int x;
955     const int w2= (width+1)>>1;
956
957     memcpy(temp, b, width*sizeof(int));
958     for(x=0; x<width2; x++){
959         b[2*x    ]= temp[x   ];
960         b[2*x + 1]= temp[x+w2];
961     }
962     if(width&1)
963         b[2*x    ]= temp[x   ];
964
965     inplace_lift(b, width, COEFFS4, N4, SHIFT4, LX0, 1);
966     inplace_lift(b, width, COEFFS3, N3, SHIFT3, LX1, 1);
967     inplace_lift(b, width, COEFFS2, N2, SHIFT2, LX0, 1);
968     inplace_lift(b, width, COEFFS1, N1, SHIFT1, LX1, 1);
969 }
970
971 static void spatial_decomposeX(DWTELEM *buffer, int width, int height, int stride){
972     int x, y;
973
974     for(y=0; y<height; y++){
975         for(x=0; x<width; x++){
976             buffer[y*stride + x] *= SCALEX;
977         }
978     }
979
980     for(y=0; y<height; y++){
981         horizontal_decomposeX(buffer + y*stride, width);
982     }
983
984     inplace_liftV(buffer, width, height, stride, COEFFS1, N1, SHIFT1, LX1, 0);
985     inplace_liftV(buffer, width, height, stride, COEFFS2, N2, SHIFT2, LX0, 0);
986     inplace_liftV(buffer, width, height, stride, COEFFS3, N3, SHIFT3, LX1, 0);
987     inplace_liftV(buffer, width, height, stride, COEFFS4, N4, SHIFT4, LX0, 0);
988 }
989
990 static void spatial_composeX(DWTELEM *buffer, int width, int height, int stride){
991     int x, y;
992
993     inplace_liftV(buffer, width, height, stride, COEFFS4, N4, SHIFT4, LX0, 1);
994     inplace_liftV(buffer, width, height, stride, COEFFS3, N3, SHIFT3, LX1, 1);
995     inplace_liftV(buffer, width, height, stride, COEFFS2, N2, SHIFT2, LX0, 1);
996     inplace_liftV(buffer, width, height, stride, COEFFS1, N1, SHIFT1, LX1, 1);
997
998     for(y=0; y<height; y++){
999         horizontal_composeX(buffer + y*stride, width);
1000     }
1001
1002     for(y=0; y<height; y++){
1003         for(x=0; x<width; x++){
1004             buffer[y*stride + x] /= SCALEX;
1005         }
1006     }
1007 }
1008
1009 static void horizontal_decompose53i(DWTELEM *b, int width){
1010     DWTELEM temp[width];
1011     const int width2= width>>1;
1012     int x;
1013     const int w2= (width+1)>>1;
1014
1015     for(x=0; x<width2; x++){
1016         temp[x   ]= b[2*x    ];
1017         temp[x+w2]= b[2*x + 1];
1018     }
1019     if(width&1)
1020         temp[x   ]= b[2*x    ];
1021 #if 0
1022     {
1023     int A1,A2,A3,A4;
1024     A2= temp[1       ];
1025     A4= temp[0       ];
1026     A1= temp[0+width2];
1027     A1 -= (A2 + A4)>>1;
1028     A4 += (A1 + 1)>>1;
1029     b[0+width2] = A1;
1030     b[0       ] = A4;
1031     for(x=1; x+1<width2; x+=2){
1032         A3= temp[x+width2];
1033         A4= temp[x+1     ];
1034         A3 -= (A2 + A4)>>1;
1035         A2 += (A1 + A3 + 2)>>2;
1036         b[x+width2] = A3;
1037         b[x       ] = A2;
1038
1039         A1= temp[x+1+width2];
1040         A2= temp[x+2       ];
1041         A1 -= (A2 + A4)>>1;
1042         A4 += (A1 + A3 + 2)>>2;
1043         b[x+1+width2] = A1;
1044         b[x+1       ] = A4;
1045     }
1046     A3= temp[width-1];
1047     A3 -= A2;
1048     A2 += (A1 + A3 + 2)>>2;
1049     b[width -1] = A3;
1050     b[width2-1] = A2;
1051     }
1052 #else
1053     lift(b+w2, temp+w2, temp, 1, 1, 1, width, -1, 0, 1, 1, 0);
1054     lift(b   , temp   , b+w2, 1, 1, 1, width,  1, 2, 2, 0, 0);
1055 #endif
1056 }
1057
1058 static void vertical_decompose53iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
1059     int i;
1060
1061     for(i=0; i<width; i++){
1062         b1[i] -= (b0[i] + b2[i])>>1;
1063     }
1064 }
1065
1066 static void vertical_decompose53iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
1067     int i;
1068
1069     for(i=0; i<width; i++){
1070         b1[i] += (b0[i] + b2[i] + 2)>>2;
1071     }
1072 }
1073
1074 static void spatial_decompose53i(DWTELEM *buffer, int width, int height, int stride){
1075     int y;
1076     DWTELEM *b0= buffer + mirror(-2-1, height-1)*stride;
1077     DWTELEM *b1= buffer + mirror(-2  , height-1)*stride;
1078
1079     for(y=-2; y<height; y+=2){
1080         DWTELEM *b2= buffer + mirror(y+1, height-1)*stride;
1081         DWTELEM *b3= buffer + mirror(y+2, height-1)*stride;
1082
1083 {START_TIMER
1084         if(y+1<(unsigned)height) horizontal_decompose53i(b2, width);
1085         if(y+2<(unsigned)height) horizontal_decompose53i(b3, width);
1086 STOP_TIMER("horizontal_decompose53i")}
1087
1088 {START_TIMER
1089         if(y+1<(unsigned)height) vertical_decompose53iH0(b1, b2, b3, width);
1090         if(y+0<(unsigned)height) vertical_decompose53iL0(b0, b1, b2, width);
1091 STOP_TIMER("vertical_decompose53i*")}
1092
1093         b0=b2;
1094         b1=b3;
1095     }
1096 }
1097
1098 static void horizontal_decompose97i(DWTELEM *b, int width){
1099     DWTELEM temp[width];
1100     const int w2= (width+1)>>1;
1101
1102     lift (temp+w2, b    +1, b      , 1, 2, 2, width, -W_AM, W_AO, W_AS, 1, 0);
1103     liftS(temp   , b      , temp+w2, 1, 2, 1, width, -W_BM, W_BO, W_BS, 0, 0);
1104     lift5(b   +w2, temp+w2, temp   , 1, 1, 1, width,  W_CM, W_CO, W_CS, 1, 0);
1105     lift (b      , temp   , b   +w2, 1, 1, 1, width,  W_DM, W_DO, W_DS, 0, 0);
1106 }
1107
1108
1109 static void vertical_decompose97iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
1110     int i;
1111
1112     for(i=0; i<width; i++){
1113         b1[i] -= (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
1114     }
1115 }
1116
1117 static void vertical_decompose97iH1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
1118     int i;
1119
1120     for(i=0; i<width; i++){
1121 #ifdef lift5
1122         b1[i] += (W_CM*(b0[i] + b2[i])+W_CO)>>W_CS;
1123 #else
1124         int r= 3*(b0[i] + b2[i]);
1125         r+= r>>4;
1126         r+= r>>8;
1127         b1[i] += (r+W_CO)>>W_CS;
1128 #endif
1129     }
1130 }
1131
1132 static void vertical_decompose97iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
1133     int i;
1134
1135     for(i=0; i<width; i++){
1136 #ifdef liftS
1137         b1[i] -= (W_BM*(b0[i] + b2[i])+W_BO)>>W_BS;
1138 #else
1139         b1[i] = (16*4*b1[i] - 4*(b0[i] + b2[i]) + 8*5 + (5<<27)) / (5*16) - (1<<23);
1140 #endif
1141     }
1142 }
1143
1144 static void vertical_decompose97iL1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
1145     int i;
1146
1147     for(i=0; i<width; i++){
1148         b1[i] += (W_DM*(b0[i] + b2[i])+W_DO)>>W_DS;
1149     }
1150 }
1151
1152 static void spatial_decompose97i(DWTELEM *buffer, int width, int height, int stride){
1153     int y;
1154     DWTELEM *b0= buffer + mirror(-4-1, height-1)*stride;
1155     DWTELEM *b1= buffer + mirror(-4  , height-1)*stride;
1156     DWTELEM *b2= buffer + mirror(-4+1, height-1)*stride;
1157     DWTELEM *b3= buffer + mirror(-4+2, height-1)*stride;
1158
1159     for(y=-4; y<height; y+=2){
1160         DWTELEM *b4= buffer + mirror(y+3, height-1)*stride;
1161         DWTELEM *b5= buffer + mirror(y+4, height-1)*stride;
1162
1163 {START_TIMER
1164         if(y+3<(unsigned)height) horizontal_decompose97i(b4, width);
1165         if(y+4<(unsigned)height) horizontal_decompose97i(b5, width);
1166 if(width>400){
1167 STOP_TIMER("horizontal_decompose97i")
1168 }}
1169
1170 {START_TIMER
1171         if(y+3<(unsigned)height) vertical_decompose97iH0(b3, b4, b5, width);
1172         if(y+2<(unsigned)height) vertical_decompose97iL0(b2, b3, b4, width);
1173         if(y+1<(unsigned)height) vertical_decompose97iH1(b1, b2, b3, width);
1174         if(y+0<(unsigned)height) vertical_decompose97iL1(b0, b1, b2, width);
1175
1176 if(width>400){
1177 STOP_TIMER("vertical_decompose97i")
1178 }}
1179
1180         b0=b2;
1181         b1=b3;
1182         b2=b4;
1183         b3=b5;
1184     }
1185 }
1186
1187 void ff_spatial_dwt(DWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count){
1188     int level;
1189
1190     for(level=0; level<decomposition_count; level++){
1191         switch(type){
1192         case 0: spatial_decompose97i(buffer, width>>level, height>>level, stride<<level); break;
1193         case 1: spatial_decompose53i(buffer, width>>level, height>>level, stride<<level); break;
1194         case 2: spatial_decomposeX  (buffer, width>>level, height>>level, stride<<level); break;
1195         }
1196     }
1197 }
1198
1199 static void horizontal_compose53i(DWTELEM *b, int width){
1200     DWTELEM temp[width];
1201     const int width2= width>>1;
1202     const int w2= (width+1)>>1;
1203     int x;
1204
1205 #if 0
1206     int A1,A2,A3,A4;
1207     A2= temp[1       ];
1208     A4= temp[0       ];
1209     A1= temp[0+width2];
1210     A1 -= (A2 + A4)>>1;
1211     A4 += (A1 + 1)>>1;
1212     b[0+width2] = A1;
1213     b[0       ] = A4;
1214     for(x=1; x+1<width2; x+=2){
1215         A3= temp[x+width2];
1216         A4= temp[x+1     ];
1217         A3 -= (A2 + A4)>>1;
1218         A2 += (A1 + A3 + 2)>>2;
1219         b[x+width2] = A3;
1220         b[x       ] = A2;
1221
1222         A1= temp[x+1+width2];
1223         A2= temp[x+2       ];
1224         A1 -= (A2 + A4)>>1;
1225         A4 += (A1 + A3 + 2)>>2;
1226         b[x+1+width2] = A1;
1227         b[x+1       ] = A4;
1228     }
1229     A3= temp[width-1];
1230     A3 -= A2;
1231     A2 += (A1 + A3 + 2)>>2;
1232     b[width -1] = A3;
1233     b[width2-1] = A2;
1234 #else
1235     lift(temp   , b   , b+w2, 1, 1, 1, width,  1, 2, 2, 0, 1);
1236     lift(temp+w2, b+w2, temp, 1, 1, 1, width, -1, 0, 1, 1, 1);
1237 #endif
1238     for(x=0; x<width2; x++){
1239         b[2*x    ]= temp[x   ];
1240         b[2*x + 1]= temp[x+w2];
1241     }
1242     if(width&1)
1243         b[2*x    ]= temp[x   ];
1244 }
1245
1246 static void vertical_compose53iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
1247     int i;
1248
1249     for(i=0; i<width; i++){
1250         b1[i] += (b0[i] + b2[i])>>1;
1251     }
1252 }
1253
1254 static void vertical_compose53iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
1255     int i;
1256
1257     for(i=0; i<width; i++){
1258         b1[i] -= (b0[i] + b2[i] + 2)>>2;
1259     }
1260 }
1261
1262 static void spatial_compose53i_buffered_init(dwt_compose_t *cs, slice_buffer * sb, int height, int stride_line){
1263     cs->b0 = slice_buffer_get_line(sb, mirror(-1-1, height-1) * stride_line);
1264     cs->b1 = slice_buffer_get_line(sb, mirror(-1  , height-1) * stride_line);
1265     cs->y = -1;
1266 }
1267
1268 static void spatial_compose53i_init(dwt_compose_t *cs, DWTELEM *buffer, int height, int stride){
1269     cs->b0 = buffer + mirror(-1-1, height-1)*stride;
1270     cs->b1 = buffer + mirror(-1  , height-1)*stride;
1271     cs->y = -1;
1272 }
1273
1274 static void spatial_compose53i_dy_buffered(dwt_compose_t *cs, slice_buffer * sb, int width, int height, int stride_line){
1275     int y= cs->y;
1276
1277     DWTELEM *b0= cs->b0;
1278     DWTELEM *b1= cs->b1;
1279     DWTELEM *b2= slice_buffer_get_line(sb, mirror(y+1, height-1) * stride_line);
1280     DWTELEM *b3= slice_buffer_get_line(sb, mirror(y+2, height-1) * stride_line);
1281
1282 {START_TIMER
1283         if(y+1<(unsigned)height) vertical_compose53iL0(b1, b2, b3, width);
1284         if(y+0<(unsigned)height) vertical_compose53iH0(b0, b1, b2, width);
1285 STOP_TIMER("vertical_compose53i*")}
1286
1287 {START_TIMER
1288         if(y-1<(unsigned)height) horizontal_compose53i(b0, width);
1289         if(y+0<(unsigned)height) horizontal_compose53i(b1, width);
1290 STOP_TIMER("horizontal_compose53i")}
1291
1292     cs->b0 = b2;
1293     cs->b1 = b3;
1294     cs->y += 2;
1295 }
1296
1297 static void spatial_compose53i_dy(dwt_compose_t *cs, DWTELEM *buffer, int width, int height, int stride){
1298     int y= cs->y;
1299     DWTELEM *b0= cs->b0;
1300     DWTELEM *b1= cs->b1;
1301     DWTELEM *b2= buffer + mirror(y+1, height-1)*stride;
1302     DWTELEM *b3= buffer + mirror(y+2, height-1)*stride;
1303
1304 {START_TIMER
1305         if(y+1<(unsigned)height) vertical_compose53iL0(b1, b2, b3, width);
1306         if(y+0<(unsigned)height) vertical_compose53iH0(b0, b1, b2, width);
1307 STOP_TIMER("vertical_compose53i*")}
1308
1309 {START_TIMER
1310         if(y-1<(unsigned)height) horizontal_compose53i(b0, width);
1311         if(y+0<(unsigned)height) horizontal_compose53i(b1, width);
1312 STOP_TIMER("horizontal_compose53i")}
1313
1314     cs->b0 = b2;
1315     cs->b1 = b3;
1316     cs->y += 2;
1317 }
1318
1319 static void spatial_compose53i(DWTELEM *buffer, int width, int height, int stride){
1320     dwt_compose_t cs;
1321     spatial_compose53i_init(&cs, buffer, height, stride);
1322     while(cs.y <= height)
1323         spatial_compose53i_dy(&cs, buffer, width, height, stride);
1324 }
1325
1326
1327 void ff_snow_horizontal_compose97i(DWTELEM *b, int width){
1328     DWTELEM temp[width];
1329     const int w2= (width+1)>>1;
1330
1331     lift (temp   , b      , b   +w2, 1, 1, 1, width,  W_DM, W_DO, W_DS, 0, 1);
1332     lift5(temp+w2, b   +w2, temp   , 1, 1, 1, width,  W_CM, W_CO, W_CS, 1, 1);
1333     liftS(b      , temp   , temp+w2, 2, 1, 1, width, -W_BM, W_BO, W_BS, 0, 1);
1334     lift (b+1    , temp+w2, b      , 2, 1, 2, width, -W_AM, W_AO, W_AS, 1, 1);
1335 }
1336
1337 static void vertical_compose97iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
1338     int i;
1339
1340     for(i=0; i<width; i++){
1341         b1[i] += (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
1342     }
1343 }
1344
1345 static void vertical_compose97iH1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
1346     int i;
1347
1348     for(i=0; i<width; i++){
1349 #ifdef lift5
1350         b1[i] -= (W_CM*(b0[i] + b2[i])+W_CO)>>W_CS;
1351 #else
1352         int r= 3*(b0[i] + b2[i]);
1353         r+= r>>4;
1354         r+= r>>8;
1355         b1[i] -= (r+W_CO)>>W_CS;
1356 #endif
1357     }
1358 }
1359
1360 static void vertical_compose97iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
1361     int i;
1362
1363     for(i=0; i<width; i++){
1364 #ifdef liftS
1365         b1[i] += (W_BM*(b0[i] + b2[i])+W_BO)>>W_BS;
1366 #else
1367         b1[i] += (W_BM*(b0[i] + b2[i])+4*b1[i]+W_BO)>>W_BS;
1368 #endif
1369     }
1370 }
1371
1372 static void vertical_compose97iL1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
1373     int i;
1374
1375     for(i=0; i<width; i++){
1376         b1[i] -= (W_DM*(b0[i] + b2[i])+W_DO)>>W_DS;
1377     }
1378 }
1379
1380 void ff_snow_vertical_compose97i(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, DWTELEM *b3, DWTELEM *b4, DWTELEM *b5, int width){
1381     int i;
1382
1383     for(i=0; i<width; i++){
1384 #ifndef lift5
1385         int r;
1386 #endif
1387         b4[i] -= (W_DM*(b3[i] + b5[i])+W_DO)>>W_DS;
1388 #ifdef lift5
1389         b3[i] -= (W_CM*(b2[i] + b4[i])+W_CO)>>W_CS;
1390 #else
1391         r= 3*(b2[i] + b4[i]);
1392         r+= r>>4;
1393         r+= r>>8;
1394         b3[i] -= (r+W_CO)>>W_CS;
1395 #endif
1396 #ifdef liftS
1397         b2[i] += (W_BM*(b1[i] + b3[i])+W_BO)>>W_BS;
1398 #else
1399         b2[i] += (W_BM*(b1[i] + b3[i])+4*b2[i]+W_BO)>>W_BS;
1400 #endif
1401         b1[i] += (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
1402     }
1403 }
1404
1405 static void spatial_compose97i_buffered_init(dwt_compose_t *cs, slice_buffer * sb, int height, int stride_line){
1406     cs->b0 = slice_buffer_get_line(sb, mirror(-3-1, height-1) * stride_line);
1407     cs->b1 = slice_buffer_get_line(sb, mirror(-3  , height-1) * stride_line);
1408     cs->b2 = slice_buffer_get_line(sb, mirror(-3+1, height-1) * stride_line);
1409     cs->b3 = slice_buffer_get_line(sb, mirror(-3+2, height-1) * stride_line);
1410     cs->y = -3;
1411 }
1412
1413 static void spatial_compose97i_init(dwt_compose_t *cs, DWTELEM *buffer, int height, int stride){
1414     cs->b0 = buffer + mirror(-3-1, height-1)*stride;
1415     cs->b1 = buffer + mirror(-3  , height-1)*stride;
1416     cs->b2 = buffer + mirror(-3+1, height-1)*stride;
1417     cs->b3 = buffer + mirror(-3+2, height-1)*stride;
1418     cs->y = -3;
1419 }
1420
1421 static void spatial_compose97i_dy_buffered(DSPContext *dsp, dwt_compose_t *cs, slice_buffer * sb, int width, int height, int stride_line){
1422     int y = cs->y;
1423
1424     DWTELEM *b0= cs->b0;
1425     DWTELEM *b1= cs->b1;
1426     DWTELEM *b2= cs->b2;
1427     DWTELEM *b3= cs->b3;
1428     DWTELEM *b4= slice_buffer_get_line(sb, mirror(y + 3, height - 1) * stride_line);
1429     DWTELEM *b5= slice_buffer_get_line(sb, mirror(y + 4, height - 1) * stride_line);
1430
1431 {START_TIMER
1432     if(y>0 && y+4<height){
1433         dsp->vertical_compose97i(b0, b1, b2, b3, b4, b5, width);
1434     }else{
1435         if(y+3<(unsigned)height) vertical_compose97iL1(b3, b4, b5, width);
1436         if(y+2<(unsigned)height) vertical_compose97iH1(b2, b3, b4, width);
1437         if(y+1<(unsigned)height) vertical_compose97iL0(b1, b2, b3, width);
1438         if(y+0<(unsigned)height) vertical_compose97iH0(b0, b1, b2, width);
1439     }
1440 if(width>400){
1441 STOP_TIMER("vertical_compose97i")}}
1442
1443 {START_TIMER
1444         if(y-1<(unsigned)height) dsp->horizontal_compose97i(b0, width);
1445         if(y+0<(unsigned)height) dsp->horizontal_compose97i(b1, width);
1446 if(width>400 && y+0<(unsigned)height){
1447 STOP_TIMER("horizontal_compose97i")}}
1448
1449     cs->b0=b2;
1450     cs->b1=b3;
1451     cs->b2=b4;
1452     cs->b3=b5;
1453     cs->y += 2;
1454 }
1455
1456 static void spatial_compose97i_dy(dwt_compose_t *cs, DWTELEM *buffer, int width, int height, int stride){
1457     int y = cs->y;
1458     DWTELEM *b0= cs->b0;
1459     DWTELEM *b1= cs->b1;
1460     DWTELEM *b2= cs->b2;
1461     DWTELEM *b3= cs->b3;
1462     DWTELEM *b4= buffer + mirror(y+3, height-1)*stride;
1463     DWTELEM *b5= buffer + mirror(y+4, height-1)*stride;
1464
1465 {START_TIMER
1466         if(y+3<(unsigned)height) vertical_compose97iL1(b3, b4, b5, width);
1467         if(y+2<(unsigned)height) vertical_compose97iH1(b2, b3, b4, width);
1468         if(y+1<(unsigned)height) vertical_compose97iL0(b1, b2, b3, width);
1469         if(y+0<(unsigned)height) vertical_compose97iH0(b0, b1, b2, width);
1470 if(width>400){
1471 STOP_TIMER("vertical_compose97i")}}
1472
1473 {START_TIMER
1474         if(y-1<(unsigned)height) ff_snow_horizontal_compose97i(b0, width);
1475         if(y+0<(unsigned)height) ff_snow_horizontal_compose97i(b1, width);
1476 if(width>400 && b0 <= b2){
1477 STOP_TIMER("horizontal_compose97i")}}
1478
1479     cs->b0=b2;
1480     cs->b1=b3;
1481     cs->b2=b4;
1482     cs->b3=b5;
1483     cs->y += 2;
1484 }
1485
1486 static void spatial_compose97i(DWTELEM *buffer, int width, int height, int stride){
1487     dwt_compose_t cs;
1488     spatial_compose97i_init(&cs, buffer, height, stride);
1489     while(cs.y <= height)
1490         spatial_compose97i_dy(&cs, buffer, width, height, stride);
1491 }
1492
1493 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){
1494     int level;
1495     for(level=decomposition_count-1; level>=0; level--){
1496         switch(type){
1497         case 0: spatial_compose97i_buffered_init(cs+level, sb, height>>level, stride_line<<level); break;
1498         case 1: spatial_compose53i_buffered_init(cs+level, sb, height>>level, stride_line<<level); break;
1499         /* not slicified yet */
1500         case 2: /*spatial_composeX(buffer, width>>level, height>>level, stride<<level); break;*/
1501           av_log(NULL, AV_LOG_ERROR, "spatial_composeX neither buffered nor slicified yet.\n"); break;
1502         }
1503     }
1504 }
1505
1506 static void ff_spatial_idwt_init(dwt_compose_t *cs, DWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count){
1507     int level;
1508     for(level=decomposition_count-1; level>=0; level--){
1509         switch(type){
1510         case 0: spatial_compose97i_init(cs+level, buffer, height>>level, stride<<level); break;
1511         case 1: spatial_compose53i_init(cs+level, buffer, height>>level, stride<<level); break;
1512         /* not slicified yet */
1513         case 2: spatial_composeX(buffer, width>>level, height>>level, stride<<level); break;
1514         }
1515     }
1516 }
1517
1518 static void ff_spatial_idwt_slice(dwt_compose_t *cs, DWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count, int y){
1519     const int support = type==1 ? 3 : 5;
1520     int level;
1521     if(type==2) return;
1522
1523     for(level=decomposition_count-1; level>=0; level--){
1524         while(cs[level].y <= FFMIN((y>>level)+support, height>>level)){
1525             switch(type){
1526             case 0: spatial_compose97i_dy(cs+level, buffer, width>>level, height>>level, stride<<level);
1527                     break;
1528             case 1: spatial_compose53i_dy(cs+level, buffer, width>>level, height>>level, stride<<level);
1529                     break;
1530             case 2: break;
1531             }
1532         }
1533     }
1534 }
1535
1536 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){
1537     const int support = type==1 ? 3 : 5;
1538     int level;
1539     if(type==2) return;
1540
1541     for(level=decomposition_count-1; level>=0; level--){
1542         while(cs[level].y <= FFMIN((y>>level)+support, height>>level)){
1543             switch(type){
1544             case 0: spatial_compose97i_dy_buffered(dsp, cs+level, slice_buf, width>>level, height>>level, stride_line<<level);
1545                     break;
1546             case 1: spatial_compose53i_dy_buffered(cs+level, slice_buf, width>>level, height>>level, stride_line<<level);
1547                     break;
1548             case 2: break;
1549             }
1550         }
1551     }
1552 }
1553
1554 static void ff_spatial_idwt(DWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count){
1555     if(type==2){
1556         int level;
1557         for(level=decomposition_count-1; level>=0; level--)
1558             spatial_composeX  (buffer, width>>level, height>>level, stride<<level);
1559     }else{
1560         dwt_compose_t cs[MAX_DECOMPOSITIONS];
1561         int y;
1562         ff_spatial_idwt_init(cs, buffer, width, height, stride, type, decomposition_count);
1563         for(y=0; y<height; y+=4)
1564             ff_spatial_idwt_slice(cs, buffer, width, height, stride, type, decomposition_count, y);
1565     }
1566 }
1567
1568 static int encode_subband_c0run(SnowContext *s, SubBand *b, DWTELEM *src, DWTELEM *parent, int stride, int orientation){
1569     const int w= b->width;
1570     const int h= b->height;
1571     int x, y;
1572
1573     if(1){
1574         int run=0;
1575         int runs[w*h];
1576         int run_index=0;
1577         int max_index;
1578
1579         for(y=0; y<h; y++){
1580             for(x=0; x<w; x++){
1581                 int v, p=0;
1582                 int /*ll=0, */l=0, lt=0, t=0, rt=0;
1583                 v= src[x + y*stride];
1584
1585                 if(y){
1586                     t= src[x + (y-1)*stride];
1587                     if(x){
1588                         lt= src[x - 1 + (y-1)*stride];
1589                     }
1590                     if(x + 1 < w){
1591                         rt= src[x + 1 + (y-1)*stride];
1592                     }
1593                 }
1594                 if(x){
1595                     l= src[x - 1 + y*stride];
1596                     /*if(x > 1){
1597                         if(orientation==1) ll= src[y + (x-2)*stride];
1598                         else               ll= src[x - 2 + y*stride];
1599                     }*/
1600                 }
1601                 if(parent){
1602                     int px= x>>1;
1603                     int py= y>>1;
1604                     if(px<b->parent->width && py<b->parent->height)
1605                         p= parent[px + py*2*stride];
1606                 }
1607                 if(!(/*ll|*/l|lt|t|rt|p)){
1608                     if(v){
1609                         runs[run_index++]= run;
1610                         run=0;
1611                     }else{
1612                         run++;
1613                     }
1614                 }
1615             }
1616         }
1617         max_index= run_index;
1618         runs[run_index++]= run;
1619         run_index=0;
1620         run= runs[run_index++];
1621
1622         put_symbol2(&s->c, b->state[30], max_index, 0);
1623         if(run_index <= max_index)
1624             put_symbol2(&s->c, b->state[1], run, 3);
1625
1626         for(y=0; y<h; y++){
1627             if(s->c.bytestream_end - s->c.bytestream < w*40){
1628                 av_log(s->avctx, AV_LOG_ERROR, "encoded frame too large\n");
1629                 return -1;
1630             }
1631             for(x=0; x<w; x++){
1632                 int v, p=0;
1633                 int /*ll=0, */l=0, lt=0, t=0, rt=0;
1634                 v= src[x + y*stride];
1635
1636                 if(y){
1637                     t= src[x + (y-1)*stride];
1638                     if(x){
1639                         lt= src[x - 1 + (y-1)*stride];
1640                     }
1641                     if(x + 1 < w){
1642                         rt= src[x + 1 + (y-1)*stride];
1643                     }
1644                 }
1645                 if(x){
1646                     l= src[x - 1 + y*stride];
1647                     /*if(x > 1){
1648                         if(orientation==1) ll= src[y + (x-2)*stride];
1649                         else               ll= src[x - 2 + y*stride];
1650                     }*/
1651                 }
1652                 if(parent){
1653                     int px= x>>1;
1654                     int py= y>>1;
1655                     if(px<b->parent->width && py<b->parent->height)
1656                         p= parent[px + py*2*stride];
1657                 }
1658                 if(/*ll|*/l|lt|t|rt|p){
1659                     int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));
1660
1661                     put_rac(&s->c, &b->state[0][context], !!v);
1662                 }else{
1663                     if(!run){
1664                         run= runs[run_index++];
1665
1666                         if(run_index <= max_index)
1667                             put_symbol2(&s->c, b->state[1], run, 3);
1668                         assert(v);
1669                     }else{
1670                         run--;
1671                         assert(!v);
1672                     }
1673                 }
1674                 if(v){
1675                     int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));
1676                     int l2= 2*ABS(l) + (l<0);
1677                     int t2= 2*ABS(t) + (t<0);
1678
1679                     put_symbol2(&s->c, b->state[context + 2], ABS(v)-1, context-4);
1680                     put_rac(&s->c, &b->state[0][16 + 1 + 3 + quant3bA[l2&0xFF] + 3*quant3bA[t2&0xFF]], v<0);
1681                 }
1682             }
1683         }
1684     }
1685     return 0;
1686 }
1687
1688 static int encode_subband(SnowContext *s, SubBand *b, DWTELEM *src, DWTELEM *parent, int stride, int orientation){
1689 //    encode_subband_qtree(s, b, src, parent, stride, orientation);
1690 //    encode_subband_z0run(s, b, src, parent, stride, orientation);
1691     return encode_subband_c0run(s, b, src, parent, stride, orientation);
1692 //    encode_subband_dzr(s, b, src, parent, stride, orientation);
1693 }
1694
1695 static inline void unpack_coeffs(SnowContext *s, SubBand *b, SubBand * parent, int orientation){
1696     const int w= b->width;
1697     const int h= b->height;
1698     int x,y;
1699
1700     if(1){
1701         int run, runs;
1702         x_and_coeff *xc= b->x_coeff;
1703         x_and_coeff *prev_xc= NULL;
1704         x_and_coeff *prev2_xc= xc;
1705         x_and_coeff *parent_xc= parent ? parent->x_coeff : NULL;
1706         x_and_coeff *prev_parent_xc= parent_xc;
1707
1708         runs= get_symbol2(&s->c, b->state[30], 0);
1709         if(runs-- > 0) run= get_symbol2(&s->c, b->state[1], 3);
1710         else           run= INT_MAX;
1711
1712         for(y=0; y<h; y++){
1713             int v=0;
1714             int lt=0, t=0, rt=0;
1715
1716             if(y && prev_xc->x == 0){
1717                 rt= prev_xc->coeff;
1718             }
1719             for(x=0; x<w; x++){
1720                 int p=0;
1721                 const int l= v;
1722
1723                 lt= t; t= rt;
1724
1725                 if(y){
1726                     if(prev_xc->x <= x)
1727                         prev_xc++;
1728                     if(prev_xc->x == x + 1)
1729                         rt= prev_xc->coeff;
1730                     else
1731                         rt=0;
1732                 }
1733                 if(parent_xc){
1734                     if(x>>1 > parent_xc->x){
1735                         parent_xc++;
1736                     }
1737                     if(x>>1 == parent_xc->x){
1738                         p= parent_xc->coeff;
1739                     }
1740                 }
1741                 if(/*ll|*/l|lt|t|rt|p){
1742                     int context= av_log2(/*ABS(ll) + */3*(l>>1) + (lt>>1) + (t&~1) + (rt>>1) + (p>>1));
1743
1744                     v=get_rac(&s->c, &b->state[0][context]);
1745                     if(v){
1746                         v= 2*(get_symbol2(&s->c, b->state[context + 2], context-4) + 1);
1747                         v+=get_rac(&s->c, &b->state[0][16 + 1 + 3 + quant3bA[l&0xFF] + 3*quant3bA[t&0xFF]]);
1748
1749                         xc->x=x;
1750                         (xc++)->coeff= v;
1751                     }
1752                 }else{
1753                     if(!run){
1754                         if(runs-- > 0) run= get_symbol2(&s->c, b->state[1], 3);
1755                         else           run= INT_MAX;
1756                         v= 2*(get_symbol2(&s->c, b->state[0 + 2], 0-4) + 1);
1757                         v+=get_rac(&s->c, &b->state[0][16 + 1 + 3]);
1758
1759                         xc->x=x;
1760                         (xc++)->coeff= v;
1761                     }else{
1762                         int max_run;
1763                         run--;
1764                         v=0;
1765
1766                         if(y) max_run= FFMIN(run, prev_xc->x - x - 2);
1767                         else  max_run= FFMIN(run, w-x-1);
1768                         if(parent_xc)
1769                             max_run= FFMIN(max_run, 2*parent_xc->x - x - 1);
1770                         x+= max_run;
1771                         run-= max_run;
1772                     }
1773                 }
1774             }
1775             (xc++)->x= w+1; //end marker
1776             prev_xc= prev2_xc;
1777             prev2_xc= xc;
1778
1779             if(parent_xc){
1780                 if(y&1){
1781                     while(parent_xc->x != parent->width+1)
1782                         parent_xc++;
1783                     parent_xc++;
1784                     prev_parent_xc= parent_xc;
1785                 }else{
1786                     parent_xc= prev_parent_xc;
1787                 }
1788             }
1789         }
1790
1791         (xc++)->x= w+1; //end marker
1792     }
1793 }
1794
1795 static inline void decode_subband_slice_buffered(SnowContext *s, SubBand *b, slice_buffer * sb, int start_y, int h, int save_state[1]){
1796     const int w= b->width;
1797     int y;
1798     const int qlog= clip(s->qlog + b->qlog, 0, QROOT*16);
1799     int qmul= qexp[qlog&(QROOT-1)]<<(qlog>>QSHIFT);
1800     int qadd= (s->qbias*qmul)>>QBIAS_SHIFT;
1801     int new_index = 0;
1802
1803     START_TIMER
1804
1805     if(b->buf == s->spatial_dwt_buffer || s->qlog == LOSSLESS_QLOG){
1806         qadd= 0;
1807         qmul= 1<<QEXPSHIFT;
1808     }
1809
1810     /* If we are on the second or later slice, restore our index. */
1811     if (start_y != 0)
1812         new_index = save_state[0];
1813
1814
1815     for(y=start_y; y<h; y++){
1816         int x = 0;
1817         int v;
1818         DWTELEM * line = slice_buffer_get_line(sb, y * b->stride_line + b->buf_y_offset) + b->buf_x_offset;
1819         memset(line, 0, b->width*sizeof(DWTELEM));
1820         v = b->x_coeff[new_index].coeff;
1821         x = b->x_coeff[new_index++].x;
1822         while(x < w)
1823         {
1824             register int t= ( (v>>1)*qmul + qadd)>>QEXPSHIFT;
1825             register int u= -(v&1);
1826             line[x] = (t^u) - u;
1827
1828             v = b->x_coeff[new_index].coeff;
1829             x = b->x_coeff[new_index++].x;
1830         }
1831     }
1832     if(w > 200 && start_y != 0/*level+1 == s->spatial_decomposition_count*/){
1833         STOP_TIMER("decode_subband")
1834     }
1835
1836     /* Save our variables for the next slice. */
1837     save_state[0] = new_index;
1838
1839     return;
1840 }
1841
1842 static void reset_contexts(SnowContext *s){
1843     int plane_index, level, orientation;
1844
1845     for(plane_index=0; plane_index<3; plane_index++){
1846         for(level=0; level<s->spatial_decomposition_count; level++){
1847             for(orientation=level ? 1:0; orientation<4; orientation++){
1848                 memset(s->plane[plane_index].band[level][orientation].state, MID_STATE, sizeof(s->plane[plane_index].band[level][orientation].state));
1849             }
1850         }
1851     }
1852     memset(s->header_state, MID_STATE, sizeof(s->header_state));
1853     memset(s->block_state, MID_STATE, sizeof(s->block_state));
1854 }
1855
1856 static int alloc_blocks(SnowContext *s){
1857     int w= -((-s->avctx->width )>>LOG2_MB_SIZE);
1858     int h= -((-s->avctx->height)>>LOG2_MB_SIZE);
1859
1860     s->b_width = w;
1861     s->b_height= h;
1862
1863     s->block= av_mallocz(w * h * sizeof(BlockNode) << (s->block_max_depth*2));
1864     return 0;
1865 }
1866
1867 static inline void copy_rac_state(RangeCoder *d, RangeCoder *s){
1868     uint8_t *bytestream= d->bytestream;
1869     uint8_t *bytestream_start= d->bytestream_start;
1870     *d= *s;
1871     d->bytestream= bytestream;
1872     d->bytestream_start= bytestream_start;
1873 }
1874
1875 //near copy & paste from dsputil, FIXME
1876 static int pix_sum(uint8_t * pix, int line_size, int w)
1877 {
1878     int s, i, j;
1879
1880     s = 0;
1881     for (i = 0; i < w; i++) {
1882         for (j = 0; j < w; j++) {
1883             s += pix[0];
1884             pix ++;
1885         }
1886         pix += line_size - w;
1887     }
1888     return s;
1889 }
1890
1891 //near copy & paste from dsputil, FIXME
1892 static int pix_norm1(uint8_t * pix, int line_size, int w)
1893 {
1894     int s, i, j;
1895     uint32_t *sq = squareTbl + 256;
1896
1897     s = 0;
1898     for (i = 0; i < w; i++) {
1899         for (j = 0; j < w; j ++) {
1900             s += sq[pix[0]];
1901             pix ++;
1902         }
1903         pix += line_size - w;
1904     }
1905     return s;
1906 }
1907
1908 static inline void set_blocks(SnowContext *s, int level, int x, int y, int l, int cb, int cr, int mx, int my, int type){
1909     const int w= s->b_width << s->block_max_depth;
1910     const int rem_depth= s->block_max_depth - level;
1911     const int index= (x + y*w) << rem_depth;
1912     const int block_w= 1<<rem_depth;
1913     BlockNode block;
1914     int i,j;
1915
1916     block.color[0]= l;
1917     block.color[1]= cb;
1918     block.color[2]= cr;
1919     block.mx= mx;
1920     block.my= my;
1921     block.type= type;
1922     block.level= level;
1923
1924     for(j=0; j<block_w; j++){
1925         for(i=0; i<block_w; i++){
1926             s->block[index + i + j*w]= block;
1927         }
1928     }
1929 }
1930
1931 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){
1932     const int offset[3]= {
1933           y*c->  stride + x,
1934         ((y*c->uvstride + x)>>1),
1935         ((y*c->uvstride + x)>>1),
1936     };
1937     int i;
1938     for(i=0; i<3; i++){
1939         c->src[0][i]= src [i];
1940         c->ref[0][i]= ref [i] + offset[i];
1941     }
1942     assert(!ref_index);
1943 }
1944
1945 //FIXME copy&paste
1946 #define P_LEFT P[1]
1947 #define P_TOP P[2]
1948 #define P_TOPRIGHT P[3]
1949 #define P_MEDIAN P[4]
1950 #define P_MV1 P[9]
1951 #define FLAG_QPEL   1 //must be 1
1952
1953 static int encode_q_branch(SnowContext *s, int level, int x, int y){
1954     uint8_t p_buffer[1024];
1955     uint8_t i_buffer[1024];
1956     uint8_t p_state[sizeof(s->block_state)];
1957     uint8_t i_state[sizeof(s->block_state)];
1958     RangeCoder pc, ic;
1959     uint8_t *pbbak= s->c.bytestream;
1960     uint8_t *pbbak_start= s->c.bytestream_start;
1961     int score, score2, iscore, i_len, p_len, block_s, sum;
1962     const int w= s->b_width  << s->block_max_depth;
1963     const int h= s->b_height << s->block_max_depth;
1964     const int rem_depth= s->block_max_depth - level;
1965     const int index= (x + y*w) << rem_depth;
1966     const int block_w= 1<<(LOG2_MB_SIZE - level);
1967     int trx= (x+1)<<rem_depth;
1968     int try= (y+1)<<rem_depth;
1969     BlockNode *left  = x ? &s->block[index-1] : &null_block;
1970     BlockNode *top   = y ? &s->block[index-w] : &null_block;
1971     BlockNode *right = trx<w ? &s->block[index+1] : &null_block;
1972     BlockNode *bottom= try<h ? &s->block[index+w] : &null_block;
1973     BlockNode *tl    = y && x ? &s->block[index-w-1] : left;
1974     BlockNode *tr    = y && trx<w && ((x&1)==0 || level==0) ? &s->block[index-w+(1<<rem_depth)] : tl; //FIXME use lt
1975     int pl = left->color[0];
1976     int pcb= left->color[1];
1977     int pcr= left->color[2];
1978     int pmx= mid_pred(left->mx, top->mx, tr->mx);
1979     int pmy= mid_pred(left->my, top->my, tr->my);
1980     int mx=0, my=0;
1981     int l,cr,cb;
1982     const int stride= s->current_picture.linesize[0];
1983     const int uvstride= s->current_picture.linesize[1];
1984     uint8_t *current_data[3]= { s->input_picture.data[0] + (x + y*  stride)*block_w,
1985                                 s->input_picture.data[1] + (x + y*uvstride)*block_w/2,
1986                                 s->input_picture.data[2] + (x + y*uvstride)*block_w/2};
1987     int P[10][2];
1988     int16_t last_mv[3][2];
1989     int qpel= !!(s->avctx->flags & CODEC_FLAG_QPEL); //unused
1990     const int shift= 1+qpel;
1991     MotionEstContext *c= &s->m.me;
1992     int mx_context= av_log2(2*ABS(left->mx - top->mx));
1993     int my_context= av_log2(2*ABS(left->my - top->my));
1994     int s_context= 2*left->level + 2*top->level + tl->level + tr->level;
1995
1996     assert(sizeof(s->block_state) >= 256);
1997     if(s->keyframe){
1998         set_blocks(s, level, x, y, pl, pcb, pcr, pmx, pmy, BLOCK_INTRA);
1999         return 0;
2000     }
2001
2002 //    clip predictors / edge ?
2003
2004     P_LEFT[0]= left->mx;
2005     P_LEFT[1]= left->my;
2006     P_TOP [0]= top->mx;
2007     P_TOP [1]= top->my;
2008     P_TOPRIGHT[0]= tr->mx;
2009     P_TOPRIGHT[1]= tr->my;
2010
2011     last_mv[0][0]= s->block[index].mx;
2012     last_mv[0][1]= s->block[index].my;
2013     last_mv[1][0]= right->mx;
2014     last_mv[1][1]= right->my;
2015     last_mv[2][0]= bottom->mx;
2016     last_mv[2][1]= bottom->my;
2017
2018     s->m.mb_stride=2;
2019     s->m.mb_x=
2020     s->m.mb_y= 0;
2021     s->m.me.skip= 0;
2022
2023     init_ref(c, current_data, s->last_picture.data, NULL, block_w*x, block_w*y, 0);
2024
2025     assert(s->m.me.  stride ==   stride);
2026     assert(s->m.me.uvstride == uvstride);
2027
2028     c->penalty_factor    = get_penalty_factor(s->lambda, s->lambda2, c->avctx->me_cmp);
2029     c->sub_penalty_factor= get_penalty_factor(s->lambda, s->lambda2, c->avctx->me_sub_cmp);
2030     c->mb_penalty_factor = get_penalty_factor(s->lambda, s->lambda2, c->avctx->mb_cmp);
2031     c->current_mv_penalty= c->mv_penalty[s->m.f_code=1] + MAX_MV;
2032
2033     c->xmin = - x*block_w - 16+2;
2034     c->ymin = - y*block_w - 16+2;
2035     c->xmax = - (x+1)*block_w + (w<<(LOG2_MB_SIZE - s->block_max_depth)) + 16-2;
2036     c->ymax = - (y+1)*block_w + (h<<(LOG2_MB_SIZE - s->block_max_depth)) + 16-2;
2037
2038     if(P_LEFT[0]     > (c->xmax<<shift)) P_LEFT[0]    = (c->xmax<<shift);
2039     if(P_LEFT[1]     > (c->ymax<<shift)) P_LEFT[1]    = (c->ymax<<shift);
2040     if(P_TOP[0]      > (c->xmax<<shift)) P_TOP[0]     = (c->xmax<<shift);
2041     if(P_TOP[1]      > (c->ymax<<shift)) P_TOP[1]     = (c->ymax<<shift);
2042     if(P_TOPRIGHT[0] < (c->xmin<<shift)) P_TOPRIGHT[0]= (c->xmin<<shift);
2043     if(P_TOPRIGHT[0] > (c->xmax<<shift)) P_TOPRIGHT[0]= (c->xmax<<shift); //due to pmx no clip
2044     if(P_TOPRIGHT[1] > (c->ymax<<shift)) P_TOPRIGHT[1]= (c->ymax<<shift);
2045
2046     P_MEDIAN[0]= mid_pred(P_LEFT[0], P_TOP[0], P_TOPRIGHT[0]);
2047     P_MEDIAN[1]= mid_pred(P_LEFT[1], P_TOP[1], P_TOPRIGHT[1]);
2048
2049     if (!y) {
2050         c->pred_x= P_LEFT[0];
2051         c->pred_y= P_LEFT[1];
2052     } else {
2053         c->pred_x = P_MEDIAN[0];
2054         c->pred_y = P_MEDIAN[1];
2055     }
2056
2057     score= ff_epzs_motion_search(&s->m, &mx, &my, P, 0, /*ref_index*/ 0, last_mv,
2058                              (1<<16)>>shift, level-LOG2_MB_SIZE+4, block_w);
2059
2060     assert(mx >= c->xmin);
2061     assert(mx <= c->xmax);
2062     assert(my >= c->ymin);
2063     assert(my <= c->ymax);
2064
2065     score= s->m.me.sub_motion_search(&s->m, &mx, &my, score, 0, 0, level-LOG2_MB_SIZE+4, block_w);
2066     score= ff_get_mb_score(&s->m, mx, my, 0, 0, level-LOG2_MB_SIZE+4, block_w, 0);
2067     //FIXME if mb_cmp != SSE then intra cant be compared currently and mb_penalty vs. lambda2
2068
2069   //  subpel search
2070     pc= s->c;
2071     pc.bytestream_start=
2072     pc.bytestream= p_buffer; //FIXME end/start? and at the other stoo
2073     memcpy(p_state, s->block_state, sizeof(s->block_state));
2074
2075     if(level!=s->block_max_depth)
2076         put_rac(&pc, &p_state[4 + s_context], 1);
2077     put_rac(&pc, &p_state[1 + left->type + top->type], 0);
2078     put_symbol(&pc, &p_state[128 + 32*mx_context], mx - pmx, 1);
2079     put_symbol(&pc, &p_state[128 + 32*my_context], my - pmy, 1);
2080     p_len= pc.bytestream - pc.bytestream_start;
2081     score += (s->lambda2*(p_len*8
2082               + (pc.outstanding_count - s->c.outstanding_count)*8
2083               + (-av_log2(pc.range)    + av_log2(s->c.range))
2084              ))>>FF_LAMBDA_SHIFT;
2085
2086     block_s= block_w*block_w;
2087     sum = pix_sum(current_data[0], stride, block_w);
2088     l= (sum + block_s/2)/block_s;
2089     iscore = pix_norm1(current_data[0], stride, block_w) - 2*l*sum + l*l*block_s;
2090
2091     block_s= block_w*block_w>>2;
2092     sum = pix_sum(current_data[1], uvstride, block_w>>1);
2093     cb= (sum + block_s/2)/block_s;
2094 //    iscore += pix_norm1(&current_mb[1][0], uvstride, block_w>>1) - 2*cb*sum + cb*cb*block_s;
2095     sum = pix_sum(current_data[2], uvstride, block_w>>1);
2096     cr= (sum + block_s/2)/block_s;
2097 //    iscore += pix_norm1(&current_mb[2][0], uvstride, block_w>>1) - 2*cr*sum + cr*cr*block_s;
2098
2099     ic= s->c;
2100     ic.bytestream_start=
2101     ic.bytestream= i_buffer; //FIXME end/start? and at the other stoo
2102     memcpy(i_state, s->block_state, sizeof(s->block_state));
2103     if(level!=s->block_max_depth)
2104         put_rac(&ic, &i_state[4 + s_context], 1);
2105     put_rac(&ic, &i_state[1 + left->type + top->type], 1);
2106     put_symbol(&ic, &i_state[32],  l-pl , 1);
2107     put_symbol(&ic, &i_state[64], cb-pcb, 1);
2108     put_symbol(&ic, &i_state[96], cr-pcr, 1);
2109     i_len= ic.bytestream - ic.bytestream_start;
2110     iscore += (s->lambda2*(i_len*8
2111               + (ic.outstanding_count - s->c.outstanding_count)*8
2112               + (-av_log2(ic.range)    + av_log2(s->c.range))
2113              ))>>FF_LAMBDA_SHIFT;
2114
2115 //    assert(score==256*256*256*64-1);
2116     assert(iscore < 255*255*256 + s->lambda2*10);
2117     assert(iscore >= 0);
2118     assert(l>=0 && l<=255);
2119     assert(pl>=0 && pl<=255);
2120
2121     if(level==0){
2122         int varc= iscore >> 8;
2123         int vard= score >> 8;
2124         if (vard <= 64 || vard < varc)
2125             c->scene_change_score+= ff_sqrt(vard) - ff_sqrt(varc);
2126         else
2127             c->scene_change_score+= s->m.qscale;
2128     }
2129
2130     if(level!=s->block_max_depth){
2131         put_rac(&s->c, &s->block_state[4 + s_context], 0);
2132         score2 = encode_q_branch(s, level+1, 2*x+0, 2*y+0);
2133         score2+= encode_q_branch(s, level+1, 2*x+1, 2*y+0);
2134         score2+= encode_q_branch(s, level+1, 2*x+0, 2*y+1);
2135         score2+= encode_q_branch(s, level+1, 2*x+1, 2*y+1);
2136         score2+= s->lambda2>>FF_LAMBDA_SHIFT; //FIXME exact split overhead
2137
2138         if(score2 < score && score2 < iscore)
2139             return score2;
2140     }
2141
2142     if(iscore < score){
2143         memcpy(pbbak, i_buffer, i_len);
2144         s->c= ic;
2145         s->c.bytestream_start= pbbak_start;
2146         s->c.bytestream= pbbak + i_len;
2147         set_blocks(s, level, x, y, l, cb, cr, pmx, pmy, BLOCK_INTRA);
2148         memcpy(s->block_state, i_state, sizeof(s->block_state));
2149         return iscore;
2150     }else{
2151         memcpy(pbbak, p_buffer, p_len);
2152         s->c= pc;
2153         s->c.bytestream_start= pbbak_start;
2154         s->c.bytestream= pbbak + p_len;
2155         set_blocks(s, level, x, y, pl, pcb, pcr, mx, my, 0);
2156         memcpy(s->block_state, p_state, sizeof(s->block_state));
2157         return score;
2158     }
2159 }
2160
2161 static always_inline int same_block(BlockNode *a, BlockNode *b){
2162     if((a->type&BLOCK_INTRA) && (b->type&BLOCK_INTRA)){
2163         return !((a->color[0] - b->color[0]) | (a->color[1] - b->color[1]) | (a->color[2] - b->color[2]));
2164     }else{
2165         return !((a->mx - b->mx) | (a->my - b->my) | ((a->type ^ b->type)&BLOCK_INTRA));
2166     }
2167 }
2168
2169 static void encode_q_branch2(SnowContext *s, int level, int x, int y){
2170     const int w= s->b_width  << s->block_max_depth;
2171     const int rem_depth= s->block_max_depth - level;
2172     const int index= (x + y*w) << rem_depth;
2173     int trx= (x+1)<<rem_depth;
2174     BlockNode *b= &s->block[index];
2175     BlockNode *left  = x ? &s->block[index-1] : &null_block;
2176     BlockNode *top   = y ? &s->block[index-w] : &null_block;
2177     BlockNode *tl    = y && x ? &s->block[index-w-1] : left;
2178     BlockNode *tr    = y && trx<w && ((x&1)==0 || level==0) ? &s->block[index-w+(1<<rem_depth)] : tl; //FIXME use lt
2179     int pl = left->color[0];
2180     int pcb= left->color[1];
2181     int pcr= left->color[2];
2182     int pmx= mid_pred(left->mx, top->mx, tr->mx);
2183     int pmy= mid_pred(left->my, top->my, tr->my);
2184     int mx_context= av_log2(2*ABS(left->mx - top->mx));
2185     int my_context= av_log2(2*ABS(left->my - top->my));
2186     int s_context= 2*left->level + 2*top->level + tl->level + tr->level;
2187
2188     if(s->keyframe){
2189         set_blocks(s, level, x, y, pl, pcb, pcr, pmx, pmy, BLOCK_INTRA);
2190         return;
2191     }
2192
2193     if(level!=s->block_max_depth){
2194         if(same_block(b,b+1) && same_block(b,b+w) && same_block(b,b+w+1)){
2195             put_rac(&s->c, &s->block_state[4 + s_context], 1);
2196         }else{
2197             put_rac(&s->c, &s->block_state[4 + s_context], 0);
2198             encode_q_branch2(s, level+1, 2*x+0, 2*y+0);
2199             encode_q_branch2(s, level+1, 2*x+1, 2*y+0);
2200             encode_q_branch2(s, level+1, 2*x+0, 2*y+1);
2201             encode_q_branch2(s, level+1, 2*x+1, 2*y+1);
2202             return;
2203         }
2204     }
2205     if(b->type & BLOCK_INTRA){
2206         put_rac(&s->c, &s->block_state[1 + (left->type&1) + (top->type&1)], 1);
2207         put_symbol(&s->c, &s->block_state[32], b->color[0]-pl , 1);
2208         put_symbol(&s->c, &s->block_state[64], b->color[1]-pcb, 1);
2209         put_symbol(&s->c, &s->block_state[96], b->color[2]-pcr, 1);
2210         set_blocks(s, level, x, y, b->color[0], b->color[1], b->color[2], pmx, pmy, BLOCK_INTRA);
2211     }else{
2212         put_rac(&s->c, &s->block_state[1 + (left->type&1) + (top->type&1)], 0);
2213         put_symbol(&s->c, &s->block_state[128 + 32*mx_context], b->mx - pmx, 1);
2214         put_symbol(&s->c, &s->block_state[128 + 32*my_context], b->my - pmy, 1);
2215         set_blocks(s, level, x, y, pl, pcb, pcr, b->mx, b->my, 0);
2216     }
2217 }
2218
2219 static void decode_q_branch(SnowContext *s, int level, int x, int y){
2220     const int w= s->b_width << s->block_max_depth;
2221     const int rem_depth= s->block_max_depth - level;
2222     const int index= (x + y*w) << rem_depth;
2223     int trx= (x+1)<<rem_depth;
2224     BlockNode *left  = x ? &s->block[index-1] : &null_block;
2225     BlockNode *top   = y ? &s->block[index-w] : &null_block;
2226     BlockNode *tl    = y && x ? &s->block[index-w-1] : left;
2227     BlockNode *tr    = y && trx<w && ((x&1)==0 || level==0) ? &s->block[index-w+(1<<rem_depth)] : tl; //FIXME use lt
2228     int s_context= 2*left->level + 2*top->level + tl->level + tr->level;
2229
2230     if(s->keyframe){
2231         set_blocks(s, level, x, y, null_block.color[0], null_block.color[1], null_block.color[2], null_block.mx, null_block.my, BLOCK_INTRA);
2232         return;
2233     }
2234
2235     if(level==s->block_max_depth || get_rac(&s->c, &s->block_state[4 + s_context])){
2236         int type;
2237         int l = left->color[0];
2238         int cb= left->color[1];
2239         int cr= left->color[2];
2240         int mx= mid_pred(left->mx, top->mx, tr->mx);
2241         int my= mid_pred(left->my, top->my, tr->my);
2242         int mx_context= av_log2(2*ABS(left->mx - top->mx)) + 0*av_log2(2*ABS(tr->mx - top->mx));
2243         int my_context= av_log2(2*ABS(left->my - top->my)) + 0*av_log2(2*ABS(tr->my - top->my));
2244
2245         type= get_rac(&s->c, &s->block_state[1 + left->type + top->type]) ? BLOCK_INTRA : 0;
2246
2247         if(type){
2248             l += get_symbol(&s->c, &s->block_state[32], 1);
2249             cb+= get_symbol(&s->c, &s->block_state[64], 1);
2250             cr+= get_symbol(&s->c, &s->block_state[96], 1);
2251         }else{
2252             mx+= get_symbol(&s->c, &s->block_state[128 + 32*mx_context], 1);
2253             my+= get_symbol(&s->c, &s->block_state[128 + 32*my_context], 1);
2254         }
2255         set_blocks(s, level, x, y, l, cb, cr, mx, my, type);
2256     }else{
2257         decode_q_branch(s, level+1, 2*x+0, 2*y+0);
2258         decode_q_branch(s, level+1, 2*x+1, 2*y+0);
2259         decode_q_branch(s, level+1, 2*x+0, 2*y+1);
2260         decode_q_branch(s, level+1, 2*x+1, 2*y+1);
2261     }
2262 }
2263
2264 static void encode_blocks(SnowContext *s){
2265     int x, y;
2266     int w= s->b_width;
2267     int h= s->b_height;
2268
2269     if(s->avctx->me_method == ME_ITER && !s->keyframe)
2270         iterative_me(s);
2271
2272     for(y=0; y<h; y++){
2273         if(s->c.bytestream_end - s->c.bytestream < w*MB_SIZE*MB_SIZE*3){ //FIXME nicer limit
2274             av_log(s->avctx, AV_LOG_ERROR, "encoded frame too large\n");
2275             return;
2276         }
2277         for(x=0; x<w; x++){
2278             if(s->avctx->me_method == ME_ITER)
2279                 encode_q_branch2(s, 0, x, y);
2280             else
2281                 encode_q_branch (s, 0, x, y);
2282         }
2283     }
2284 }
2285
2286 static void decode_blocks(SnowContext *s){
2287     int x, y;
2288     int w= s->b_width;
2289     int h= s->b_height;
2290
2291     for(y=0; y<h; y++){
2292         for(x=0; x<w; x++){
2293             decode_q_branch(s, 0, x, y);
2294         }
2295     }
2296 }
2297
2298 static void mc_block(uint8_t *dst, uint8_t *src, uint8_t *tmp, int stride, int b_w, int b_h, int dx, int dy){
2299     int x, y;
2300 START_TIMER
2301     for(y=0; y < b_h+5; y++){
2302         for(x=0; x < b_w; x++){
2303             int a0= src[x    ];
2304             int a1= src[x + 1];
2305             int a2= src[x + 2];
2306             int a3= src[x + 3];
2307             int a4= src[x + 4];
2308             int a5= src[x + 5];
2309 //            int am= 9*(a1+a2) - (a0+a3);
2310             int am= 20*(a2+a3) - 5*(a1+a4) + (a0+a5);
2311 //            int am= 18*(a2+a3) - 2*(a1+a4);
2312 //             int aL= (-7*a0 + 105*a1 + 35*a2 - 5*a3)>>3;
2313 //             int aR= (-7*a3 + 105*a2 + 35*a1 - 5*a0)>>3;
2314
2315 //            if(b_w==16) am= 8*(a1+a2);
2316
2317             if(dx<8) am = (32*a2*( 8-dx) +    am* dx    + 128)>>8;
2318             else     am = (   am*(16-dx) + 32*a3*(dx-8) + 128)>>8;
2319
2320             /* FIXME Try increasing tmp buffer to 16 bits and not clipping here. Should give marginally better results. - Robert*/
2321             if(am&(~255)) am= ~(am>>31);
2322
2323             tmp[x] = am;
2324
2325 /*            if     (dx< 4) tmp[x + y*stride]= (16*a1*( 4-dx) +    aL* dx     + 32)>>6;
2326             else if(dx< 8) tmp[x + y*stride]= (   aL*( 8-dx) +    am*(dx- 4) + 32)>>6;
2327             else if(dx<12) tmp[x + y*stride]= (   am*(12-dx) +    aR*(dx- 8) + 32)>>6;
2328             else           tmp[x + y*stride]= (   aR*(16-dx) + 16*a2*(dx-12) + 32)>>6;*/
2329         }
2330         tmp += stride;
2331         src += stride;
2332     }
2333     tmp -= (b_h+5)*stride;
2334
2335     for(y=0; y < b_h; y++){
2336         for(x=0; x < b_w; x++){
2337             int a0= tmp[x + 0*stride];
2338             int a1= tmp[x + 1*stride];
2339             int a2= tmp[x + 2*stride];
2340             int a3= tmp[x + 3*stride];
2341             int a4= tmp[x + 4*stride];
2342             int a5= tmp[x + 5*stride];
2343             int am= 20*(a2+a3) - 5*(a1+a4) + (a0+a5);
2344 //            int am= 18*(a2+a3) - 2*(a1+a4);
2345 /*            int aL= (-7*a0 + 105*a1 + 35*a2 - 5*a3)>>3;
2346             int aR= (-7*a3 + 105*a2 + 35*a1 - 5*a0)>>3;*/
2347
2348 //            if(b_w==16) am= 8*(a1+a2);
2349
2350             if(dy<8) am =  (32*a2*( 8-dy) +    am* dy    + 128)>>8;
2351             else     am = (   am*(16-dy) + 32*a3*(dy-8) + 128)>>8;
2352
2353             if(am&(~255)) am= ~(am>>31);
2354
2355             dst[x] = am;
2356 /*            if     (dy< 4) tmp[x + y*stride]= (16*a1*( 4-dy) +    aL* dy     + 32)>>6;
2357             else if(dy< 8) tmp[x + y*stride]= (   aL*( 8-dy) +    am*(dy- 4) + 32)>>6;
2358             else if(dy<12) tmp[x + y*stride]= (   am*(12-dy) +    aR*(dy- 8) + 32)>>6;
2359             else           tmp[x + y*stride]= (   aR*(16-dy) + 16*a2*(dy-12) + 32)>>6;*/
2360         }
2361         dst += stride;
2362         tmp += stride;
2363     }
2364 STOP_TIMER("mc_block")
2365 }
2366
2367 #define mca(dx,dy,b_w)\
2368 static void mc_block_hpel ## dx ## dy ## b_w(uint8_t *dst, uint8_t *src, int stride, int h){\
2369     uint8_t tmp[stride*(b_w+5)];\
2370     assert(h==b_w);\
2371     mc_block(dst, src-2-2*stride, tmp, stride, b_w, b_w, dx, dy);\
2372 }
2373
2374 mca( 0, 0,16)
2375 mca( 8, 0,16)
2376 mca( 0, 8,16)
2377 mca( 8, 8,16)
2378 mca( 0, 0,8)
2379 mca( 8, 0,8)
2380 mca( 0, 8,8)
2381 mca( 8, 8,8)
2382
2383 static void pred_block(SnowContext *s, uint8_t *dst, uint8_t *src, uint8_t *tmp, int stride, int sx, int sy, int b_w, int b_h, BlockNode *block, int plane_index, int w, int h){
2384     if(block->type & BLOCK_INTRA){
2385         int x, y;
2386         const int color = block->color[plane_index];
2387         const int color4= color*0x01010101;
2388         if(b_w==32){
2389             for(y=0; y < b_h; y++){
2390                 *(uint32_t*)&dst[0 + y*stride]= color4;
2391                 *(uint32_t*)&dst[4 + y*stride]= color4;
2392                 *(uint32_t*)&dst[8 + y*stride]= color4;
2393                 *(uint32_t*)&dst[12+ y*stride]= color4;
2394                 *(uint32_t*)&dst[16+ y*stride]= color4;
2395                 *(uint32_t*)&dst[20+ y*stride]= color4;
2396                 *(uint32_t*)&dst[24+ y*stride]= color4;
2397                 *(uint32_t*)&dst[28+ y*stride]= color4;
2398             }
2399         }else if(b_w==16){
2400             for(y=0; y < b_h; y++){
2401                 *(uint32_t*)&dst[0 + y*stride]= color4;
2402                 *(uint32_t*)&dst[4 + y*stride]= color4;
2403                 *(uint32_t*)&dst[8 + y*stride]= color4;
2404                 *(uint32_t*)&dst[12+ y*stride]= color4;
2405             }
2406         }else if(b_w==8){
2407             for(y=0; y < b_h; y++){
2408                 *(uint32_t*)&dst[0 + y*stride]= color4;
2409                 *(uint32_t*)&dst[4 + y*stride]= color4;
2410             }
2411         }else if(b_w==4){
2412             for(y=0; y < b_h; y++){
2413                 *(uint32_t*)&dst[0 + y*stride]= color4;
2414             }
2415         }else{
2416             for(y=0; y < b_h; y++){
2417                 for(x=0; x < b_w; x++){
2418                     dst[x + y*stride]= color;
2419                 }
2420             }
2421         }
2422     }else{
2423         const int scale= plane_index ?  s->mv_scale : 2*s->mv_scale;
2424         int mx= block->mx*scale;
2425         int my= block->my*scale;
2426         const int dx= mx&15;
2427         const int dy= my&15;
2428         const int tab_index= 3 - (b_w>>2) + (b_w>>4);
2429         sx += (mx>>4) - 2;
2430         sy += (my>>4) - 2;
2431         src += sx + sy*stride;
2432         if(   (unsigned)sx >= w - b_w - 4
2433            || (unsigned)sy >= h - b_h - 4){
2434             ff_emulated_edge_mc(tmp + MB_SIZE, src, stride, b_w+5, b_h+5, sx, sy, w, h);
2435             src= tmp + MB_SIZE;
2436         }
2437 //        assert(b_w == b_h || 2*b_w == b_h || b_w == 2*b_h);
2438 //        assert(!(b_w&(b_w-1)));
2439         assert(b_w>1 && b_h>1);
2440         assert(tab_index>=0 && tab_index<4 || b_w==32);
2441         if((dx&3) || (dy&3) || !(b_w == b_h || 2*b_w == b_h || b_w == 2*b_h) || (b_w&(b_w-1)))
2442             mc_block(dst, src, tmp, stride, b_w, b_h, dx, dy);
2443         else if(b_w==32){
2444             int y;
2445             for(y=0; y<b_h; y+=16){
2446                 s->dsp.put_h264_qpel_pixels_tab[0][dy+(dx>>2)](dst + y*stride, src + 2 + (y+2)*stride,stride);
2447                 s->dsp.put_h264_qpel_pixels_tab[0][dy+(dx>>2)](dst + 16 + y*stride, src + 18 + (y+2)*stride,stride);
2448             }
2449         }else if(b_w==b_h)
2450             s->dsp.put_h264_qpel_pixels_tab[tab_index  ][dy+(dx>>2)](dst,src + 2 + 2*stride,stride);
2451         else if(b_w==2*b_h){
2452             s->dsp.put_h264_qpel_pixels_tab[tab_index+1][dy+(dx>>2)](dst    ,src + 2       + 2*stride,stride);
2453             s->dsp.put_h264_qpel_pixels_tab[tab_index+1][dy+(dx>>2)](dst+b_h,src + 2 + b_h + 2*stride,stride);
2454         }else{
2455             assert(2*b_w==b_h);
2456             s->dsp.put_h264_qpel_pixels_tab[tab_index  ][dy+(dx>>2)](dst           ,src + 2 + 2*stride           ,stride);
2457             s->dsp.put_h264_qpel_pixels_tab[tab_index  ][dy+(dx>>2)](dst+b_w*stride,src + 2 + 2*stride+b_w*stride,stride);
2458         }
2459     }
2460 }
2461
2462 void ff_snow_inner_add_yblock(uint8_t *obmc, const int obmc_stride, uint8_t * * block, int b_w, int b_h,
2463                               int src_x, int src_y, int src_stride, slice_buffer * sb, int add, uint8_t * dst8){
2464     int y, x;
2465     DWTELEM * dst;
2466     for(y=0; y<b_h; y++){
2467         //FIXME ugly missue of obmc_stride
2468         uint8_t *obmc1= obmc + y*obmc_stride;
2469         uint8_t *obmc2= obmc1+ (obmc_stride>>1);
2470         uint8_t *obmc3= obmc1+ obmc_stride*(obmc_stride>>1);
2471         uint8_t *obmc4= obmc3+ (obmc_stride>>1);
2472         dst = slice_buffer_get_line(sb, src_y + y);
2473         for(x=0; x<b_w; x++){
2474             int v=   obmc1[x] * block[3][x + y*src_stride]
2475                     +obmc2[x] * block[2][x + y*src_stride]
2476                     +obmc3[x] * block[1][x + y*src_stride]
2477                     +obmc4[x] * block[0][x + y*src_stride];
2478
2479             v <<= 8 - LOG2_OBMC_MAX;
2480             if(FRAC_BITS != 8){
2481                 v += 1<<(7 - FRAC_BITS);
2482                 v >>= 8 - FRAC_BITS;
2483             }
2484             if(add){
2485                 v += dst[x + src_x];
2486                 v = (v + (1<<(FRAC_BITS-1))) >> FRAC_BITS;
2487                 if(v&(~255)) v= ~(v>>31);
2488                 dst8[x + y*src_stride] = v;
2489             }else{
2490                 dst[x + src_x] -= v;
2491             }
2492         }
2493     }
2494 }
2495
2496 //FIXME name clenup (b_w, block_w, b_width stuff)
2497 static always_inline void add_yblock_buffered(SnowContext *s, slice_buffer * sb, DWTELEM *old_dst, uint8_t *dst8, uint8_t *src, 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 plane_index){
2498     DWTELEM * dst = NULL;
2499     const int b_width = s->b_width  << s->block_max_depth;
2500     const int b_height= s->b_height << s->block_max_depth;
2501     const int b_stride= b_width;
2502     BlockNode *lt= &s->block[b_x + b_y*b_stride];
2503     BlockNode *rt= lt+1;
2504     BlockNode *lb= lt+b_stride;
2505     BlockNode *rb= lb+1;
2506     uint8_t *block[4];
2507     int tmp_step= src_stride >= 7*MB_SIZE ? MB_SIZE : MB_SIZE*src_stride;
2508     uint8_t tmp[src_stride*7*MB_SIZE]; //FIXME align
2509     uint8_t *ptmp;
2510     int x,y;
2511
2512     if(b_x<0){
2513         lt= rt;
2514         lb= rb;
2515     }else if(b_x + 1 >= b_width){
2516         rt= lt;
2517         rb= lb;
2518     }
2519     if(b_y<0){
2520         lt= lb;
2521         rt= rb;
2522     }else if(b_y + 1 >= b_height){
2523         lb= lt;
2524         rb= rt;
2525     }
2526
2527     if(src_x<0){ //FIXME merge with prev & always round internal width upto *16
2528         obmc -= src_x;
2529         b_w += src_x;
2530         src_x=0;
2531     }else if(src_x + b_w > w){
2532         b_w = w - src_x;
2533     }
2534     if(src_y<0){
2535         obmc -= src_y*obmc_stride;
2536         b_h += src_y;
2537         src_y=0;
2538     }else if(src_y + b_h> h){
2539         b_h = h - src_y;
2540     }
2541
2542     if(b_w<=0 || b_h<=0) return;
2543
2544 assert(src_stride > 2*MB_SIZE + 5);
2545 //    old_dst += src_x + src_y*dst_stride;
2546     dst8+= src_x + src_y*src_stride;
2547 //    src += src_x + src_y*src_stride;
2548
2549     ptmp= tmp + 3*tmp_step;
2550     block[0]= ptmp;
2551     ptmp+=tmp_step;
2552     pred_block(s, block[0], src, tmp, src_stride, src_x, src_y, b_w, b_h, lt, plane_index, w, h);
2553
2554     if(same_block(lt, rt)){
2555         block[1]= block[0];
2556     }else{
2557         block[1]= ptmp;
2558         ptmp+=tmp_step;
2559         pred_block(s, block[1], src, tmp, src_stride, src_x, src_y, b_w, b_h, rt, plane_index, w, h);
2560     }
2561
2562     if(same_block(lt, lb)){
2563         block[2]= block[0];
2564     }else if(same_block(rt, lb)){
2565         block[2]= block[1];
2566     }else{
2567         block[2]= ptmp;
2568         ptmp+=tmp_step;
2569         pred_block(s, block[2], src, tmp, src_stride, src_x, src_y, b_w, b_h, lb, plane_index, w, h);
2570     }
2571
2572     if(same_block(lt, rb) ){
2573         block[3]= block[0];
2574     }else if(same_block(rt, rb)){
2575         block[3]= block[1];
2576     }else if(same_block(lb, rb)){
2577         block[3]= block[2];
2578     }else{
2579         block[3]= ptmp;
2580         pred_block(s, block[3], src, tmp, src_stride, src_x, src_y, b_w, b_h, rb, plane_index, w, h);
2581     }
2582 #if 0
2583     for(y=0; y<b_h; y++){
2584         for(x=0; x<b_w; x++){
2585             int v=   obmc [x + y*obmc_stride] * block[3][x + y*src_stride] * (256/OBMC_MAX);
2586             if(add) dst[x + y*dst_stride] += v;
2587             else    dst[x + y*dst_stride] -= v;
2588         }
2589     }
2590     for(y=0; y<b_h; y++){
2591         uint8_t *obmc2= obmc + (obmc_stride>>1);
2592         for(x=0; x<b_w; x++){
2593             int v=   obmc2[x + y*obmc_stride] * block[2][x + y*src_stride] * (256/OBMC_MAX);
2594             if(add) dst[x + y*dst_stride] += v;
2595             else    dst[x + y*dst_stride] -= v;
2596         }
2597     }
2598     for(y=0; y<b_h; y++){
2599         uint8_t *obmc3= obmc + obmc_stride*(obmc_stride>>1);
2600         for(x=0; x<b_w; x++){
2601             int v=   obmc3[x + y*obmc_stride] * block[1][x + y*src_stride] * (256/OBMC_MAX);
2602             if(add) dst[x + y*dst_stride] += v;
2603             else    dst[x + y*dst_stride] -= v;
2604         }
2605     }
2606     for(y=0; y<b_h; y++){
2607         uint8_t *obmc3= obmc + obmc_stride*(obmc_stride>>1);
2608         uint8_t *obmc4= obmc3+ (obmc_stride>>1);
2609         for(x=0; x<b_w; x++){
2610             int v=   obmc4[x + y*obmc_stride] * block[0][x + y*src_stride] * (256/OBMC_MAX);
2611             if(add) dst[x + y*dst_stride] += v;
2612             else    dst[x + y*dst_stride] -= v;
2613         }
2614     }
2615 #else
2616 {
2617
2618     START_TIMER
2619
2620     s->dsp.inner_add_yblock(obmc, obmc_stride, block, b_w, b_h, src_x,src_y, src_stride, sb, add, dst8);
2621         STOP_TIMER("Inner add y block")
2622 }
2623 #endif
2624 }
2625
2626 //FIXME name clenup (b_w, block_w, b_width stuff)
2627 static always_inline void add_yblock(SnowContext *s, DWTELEM *dst, uint8_t *dst8, uint8_t *src, 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){
2628     const int b_width = s->b_width  << s->block_max_depth;
2629     const int b_height= s->b_height << s->block_max_depth;
2630     const int b_stride= b_width;
2631     BlockNode *lt= &s->block[b_x + b_y*b_stride];
2632     BlockNode *rt= lt+1;
2633     BlockNode *lb= lt+b_stride;
2634     BlockNode *rb= lb+1;
2635     uint8_t *block[4];
2636     int tmp_step= src_stride >= 7*MB_SIZE ? MB_SIZE : MB_SIZE*src_stride;
2637     uint8_t tmp[src_stride*7*MB_SIZE]; //FIXME align
2638     uint8_t *ptmp;
2639     int x,y;
2640
2641     if(b_x<0){
2642         lt= rt;
2643         lb= rb;
2644     }else if(b_x + 1 >= b_width){
2645         rt= lt;
2646         rb= lb;
2647     }
2648     if(b_y<0){
2649         lt= lb;
2650         rt= rb;
2651     }else if(b_y + 1 >= b_height){
2652         lb= lt;
2653         rb= rt;
2654     }
2655
2656     if(src_x<0){ //FIXME merge with prev & always round internal width upto *16
2657         obmc -= src_x;
2658         b_w += src_x;
2659         if(!offset_dst)
2660             dst -= src_x;
2661         src_x=0;
2662     }else if(src_x + b_w > w){
2663         b_w = w - src_x;
2664     }
2665     if(src_y<0){
2666         obmc -= src_y*obmc_stride;
2667         b_h += src_y;
2668         if(!offset_dst)
2669             dst -= src_y*dst_stride;
2670         src_y=0;
2671     }else if(src_y + b_h> h){
2672         b_h = h - src_y;
2673     }
2674
2675     if(b_w<=0 || b_h<=0) return;
2676
2677 assert(src_stride > 2*MB_SIZE + 5);
2678     if(offset_dst)
2679         dst += src_x + src_y*dst_stride;
2680     dst8+= src_x + src_y*src_stride;
2681 //    src += src_x + src_y*src_stride;
2682
2683     ptmp= tmp + 3*tmp_step;
2684     block[0]= ptmp;
2685     ptmp+=tmp_step;
2686     pred_block(s, block[0], src, tmp, src_stride, src_x, src_y, b_w, b_h, lt, plane_index, w, h);
2687
2688     if(same_block(lt, rt)){
2689         block[1]= block[0];
2690     }else{
2691         block[1]= ptmp;
2692         ptmp+=tmp_step;
2693         pred_block(s, block[1], src, tmp, src_stride, src_x, src_y, b_w, b_h, rt, plane_index, w, h);
2694     }
2695
2696     if(same_block(lt, lb)){
2697         block[2]= block[0];
2698     }else if(same_block(rt, lb)){
2699         block[2]= block[1];
2700     }else{
2701         block[2]= ptmp;
2702         ptmp+=tmp_step;
2703         pred_block(s, block[2], src, tmp, src_stride, src_x, src_y, b_w, b_h, lb, plane_index, w, h);
2704     }
2705
2706     if(same_block(lt, rb) ){
2707         block[3]= block[0];
2708     }else if(same_block(rt, rb)){
2709         block[3]= block[1];
2710     }else if(same_block(lb, rb)){
2711         block[3]= block[2];
2712     }else{
2713         block[3]= ptmp;
2714         pred_block(s, block[3], src, tmp, src_stride, src_x, src_y, b_w, b_h, rb, plane_index, w, h);
2715     }
2716 #if 0
2717     for(y=0; y<b_h; y++){
2718         for(x=0; x<b_w; x++){
2719             int v=   obmc [x + y*obmc_stride] * block[3][x + y*src_stride] * (256/OBMC_MAX);
2720             if(add) dst[x + y*dst_stride] += v;
2721             else    dst[x + y*dst_stride] -= v;
2722         }
2723     }
2724     for(y=0; y<b_h; y++){
2725         uint8_t *obmc2= obmc + (obmc_stride>>1);
2726         for(x=0; x<b_w; x++){
2727             int v=   obmc2[x + y*obmc_stride] * block[2][x + y*src_stride] * (256/OBMC_MAX);
2728             if(add) dst[x + y*dst_stride] += v;
2729             else    dst[x + y*dst_stride] -= v;
2730         }
2731     }
2732     for(y=0; y<b_h; y++){
2733         uint8_t *obmc3= obmc + obmc_stride*(obmc_stride>>1);
2734         for(x=0; x<b_w; x++){
2735             int v=   obmc3[x + y*obmc_stride] * block[1][x + y*src_stride] * (256/OBMC_MAX);
2736             if(add) dst[x + y*dst_stride] += v;
2737             else    dst[x + y*dst_stride] -= v;
2738         }
2739     }
2740     for(y=0; y<b_h; y++){
2741         uint8_t *obmc3= obmc + obmc_stride*(obmc_stride>>1);
2742         uint8_t *obmc4= obmc3+ (obmc_stride>>1);
2743         for(x=0; x<b_w; x++){
2744             int v=   obmc4[x + y*obmc_stride] * block[0][x + y*src_stride] * (256/OBMC_MAX);
2745             if(add) dst[x + y*dst_stride] += v;
2746             else    dst[x + y*dst_stride] -= v;
2747         }
2748     }
2749 #else
2750     for(y=0; y<b_h; y++){
2751         //FIXME ugly missue of obmc_stride
2752         uint8_t *obmc1= obmc + y*obmc_stride;
2753         uint8_t *obmc2= obmc1+ (obmc_stride>>1);
2754         uint8_t *obmc3= obmc1+ obmc_stride*(obmc_stride>>1);
2755         uint8_t *obmc4= obmc3+ (obmc_stride>>1);
2756         for(x=0; x<b_w; x++){
2757             int v=   obmc1[x] * block[3][x + y*src_stride]
2758                     +obmc2[x] * block[2][x + y*src_stride]
2759                     +obmc3[x] * block[1][x + y*src_stride]
2760                     +obmc4[x] * block[0][x + y*src_stride];
2761
2762             v <<= 8 - LOG2_OBMC_MAX;
2763             if(FRAC_BITS != 8){
2764                 v += 1<<(7 - FRAC_BITS);
2765                 v >>= 8 - FRAC_BITS;
2766             }
2767             if(add){
2768                 v += dst[x + y*dst_stride];
2769                 v = (v + (1<<(FRAC_BITS-1))) >> FRAC_BITS;
2770                 if(v&(~255)) v= ~(v>>31);
2771                 dst8[x + y*src_stride] = v;
2772             }else{
2773                 dst[x + y*dst_stride] -= v;
2774             }
2775         }
2776     }
2777 #endif
2778 }
2779
2780 static always_inline void predict_slice_buffered(SnowContext *s, slice_buffer * sb, DWTELEM * old_buffer, int plane_index, int add, int mb_y){
2781     Plane *p= &s->plane[plane_index];
2782     const int mb_w= s->b_width  << s->block_max_depth;
2783     const int mb_h= s->b_height << s->block_max_depth;
2784     int x, y, mb_x;
2785     int block_size = MB_SIZE >> s->block_max_depth;
2786     int block_w    = plane_index ? block_size/2 : block_size;
2787     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];
2788     int obmc_stride= plane_index ? block_size : 2*block_size;
2789     int ref_stride= s->current_picture.linesize[plane_index];
2790     uint8_t *ref  = s->last_picture.data[plane_index];
2791     uint8_t *dst8= s->current_picture.data[plane_index];
2792     int w= p->width;
2793     int h= p->height;
2794     START_TIMER
2795
2796     if(s->keyframe || (s->avctx->debug&512)){
2797         if(mb_y==mb_h)
2798             return;
2799
2800         if(add){
2801             for(y=block_w*mb_y; y<FFMIN(h,block_w*(mb_y+1)); y++)
2802             {
2803 //                DWTELEM * line = slice_buffer_get_line(sb, y);
2804                 DWTELEM * line = sb->line[y];
2805                 for(x=0; x<w; x++)
2806                 {
2807 //                    int v= buf[x + y*w] + (128<<FRAC_BITS) + (1<<(FRAC_BITS-1));
2808                     int v= line[x] + (128<<FRAC_BITS) + (1<<(FRAC_BITS-1));
2809                     v >>= FRAC_BITS;
2810                     if(v&(~255)) v= ~(v>>31);
2811                     dst8[x + y*ref_stride]= v;
2812                 }
2813             }
2814         }else{
2815             for(y=block_w*mb_y; y<FFMIN(h,block_w*(mb_y+1)); y++)
2816             {
2817 //                DWTELEM * line = slice_buffer_get_line(sb, y);
2818                 DWTELEM * line = sb->line[y];
2819                 for(x=0; x<w; x++)
2820                 {
2821                     line[x] -= 128 << FRAC_BITS;
2822 //                    buf[x + y*w]-= 128<<FRAC_BITS;
2823                 }
2824             }
2825         }
2826
2827         return;
2828     }
2829
2830         for(mb_x=0; mb_x<=mb_w; mb_x++){
2831             START_TIMER
2832
2833             add_yblock_buffered(s, sb, old_buffer, dst8, ref, obmc,
2834                        block_w*mb_x - block_w/2,
2835                        block_w*mb_y - block_w/2,
2836                        block_w, block_w,
2837                        w, h,
2838                        w, ref_stride, obmc_stride,
2839                        mb_x - 1, mb_y - 1,
2840                        add, plane_index);
2841
2842             STOP_TIMER("add_yblock")
2843         }
2844
2845     STOP_TIMER("predict_slice")
2846 }
2847
2848 static always_inline void predict_slice(SnowContext *s, DWTELEM *buf, int plane_index, int add, int mb_y){
2849     Plane *p= &s->plane[plane_index];
2850     const int mb_w= s->b_width  << s->block_max_depth;
2851     const int mb_h= s->b_height << s->block_max_depth;
2852     int x, y, mb_x;
2853     int block_size = MB_SIZE >> s->block_max_depth;
2854     int block_w    = plane_index ? block_size/2 : block_size;
2855     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];
2856     const int obmc_stride= plane_index ? block_size : 2*block_size;
2857     int ref_stride= s->current_picture.linesize[plane_index];
2858     uint8_t *ref  = s->last_picture.data[plane_index];
2859     uint8_t *dst8= s->current_picture.data[plane_index];
2860     int w= p->width;
2861     int h= p->height;
2862     START_TIMER
2863
2864     if(s->keyframe || (s->avctx->debug&512)){
2865         if(mb_y==mb_h)
2866             return;
2867
2868         if(add){
2869             for(y=block_w*mb_y; y<FFMIN(h,block_w*(mb_y+1)); y++){
2870                 for(x=0; x<w; x++){
2871                     int v= buf[x + y*w] + (128<<FRAC_BITS) + (1<<(FRAC_BITS-1));
2872                     v >>= FRAC_BITS;
2873                     if(v&(~255)) v= ~(v>>31);
2874                     dst8[x + y*ref_stride]= v;
2875                 }
2876             }
2877         }else{
2878             for(y=block_w*mb_y; y<FFMIN(h,block_w*(mb_y+1)); y++){
2879                 for(x=0; x<w; x++){
2880                     buf[x + y*w]-= 128<<FRAC_BITS;
2881                 }
2882             }
2883         }
2884
2885         return;
2886     }
2887
2888         for(mb_x=0; mb_x<=mb_w; mb_x++){
2889             START_TIMER
2890
2891             add_yblock(s, buf, dst8, ref, obmc,
2892                        block_w*mb_x - block_w/2,
2893                        block_w*mb_y - block_w/2,
2894                        block_w, block_w,
2895                        w, h,
2896                        w, ref_stride, obmc_stride,
2897                        mb_x - 1, mb_y - 1,
2898                        add, 1, plane_index);
2899
2900             STOP_TIMER("add_yblock")
2901         }
2902
2903     STOP_TIMER("predict_slice")
2904 }
2905
2906 static always_inline void predict_plane(SnowContext *s, DWTELEM *buf, int plane_index, int add){
2907     const int mb_h= s->b_height << s->block_max_depth;
2908     int mb_y;
2909     for(mb_y=0; mb_y<=mb_h; mb_y++)
2910         predict_slice(s, buf, plane_index, add, mb_y);
2911 }
2912
2913 static int get_dc(SnowContext *s, int mb_x, int mb_y, int plane_index){
2914     int i, x2, y2;
2915     Plane *p= &s->plane[plane_index];
2916     const int block_size = MB_SIZE >> s->block_max_depth;
2917     const int block_w    = plane_index ? block_size/2 : block_size;
2918     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];
2919     const int obmc_stride= plane_index ? block_size : 2*block_size;
2920     const int ref_stride= s->current_picture.linesize[plane_index];
2921     uint8_t *ref= s->   last_picture.data[plane_index];
2922     uint8_t *src= s-> input_picture.data[plane_index];
2923     DWTELEM *dst= (DWTELEM*)s->m.obmc_scratchpad + plane_index*block_size*block_size*4;
2924     const int b_stride = s->b_width << s->block_max_depth;
2925     const int w= p->width;
2926     const int h= p->height;
2927     int index= mb_x + mb_y*b_stride;
2928     BlockNode *b= &s->block[index];
2929     BlockNode backup= *b;
2930     int ab=0;
2931     int aa=0;
2932
2933     b->type|= BLOCK_INTRA;
2934     b->color[plane_index]= 0;
2935     memset(dst, 0, obmc_stride*obmc_stride*sizeof(DWTELEM));
2936
2937     for(i=0; i<4; i++){
2938         int mb_x2= mb_x + (i &1) - 1;
2939         int mb_y2= mb_y + (i>>1) - 1;
2940         int x= block_w*mb_x2 + block_w/2;
2941         int y= block_w*mb_y2 + block_w/2;
2942
2943         add_yblock(s, dst + ((i&1)+(i>>1)*obmc_stride)*block_w, NULL, ref, obmc,
2944                     x, y, block_w, block_w, w, h, obmc_stride, ref_stride, obmc_stride, mb_x2, mb_y2, 0, 0, plane_index);
2945
2946         for(y2= FFMAX(y, 0); y2<FFMIN(h, y+block_w); y2++){
2947             for(x2= FFMAX(x, 0); x2<FFMIN(w, x+block_w); x2++){
2948                 int index= x2-(block_w*mb_x - block_w/2) + (y2-(block_w*mb_y - block_w/2))*obmc_stride;
2949                 int obmc_v= obmc[index];
2950                 int d;
2951                 if(y<0) obmc_v += obmc[index + block_w*obmc_stride];
2952                 if(x<0) obmc_v += obmc[index + block_w];
2953                 if(y+block_w>h) obmc_v += obmc[index - block_w*obmc_stride];
2954                 if(x+block_w>w) obmc_v += obmc[index - block_w];
2955                 //FIXME precalc this or simplify it somehow else
2956
2957                 d = -dst[index] + (1<<(FRAC_BITS-1));
2958                 dst[index] = d;
2959                 ab += (src[x2 + y2*ref_stride] - (d>>FRAC_BITS)) * obmc_v;
2960                 aa += obmc_v * obmc_v; //FIXME precalclate this
2961             }
2962         }
2963     }
2964     *b= backup;
2965
2966     return clip(((ab<<LOG2_OBMC_MAX) + aa/2)/aa, 0, 255); //FIXME we shouldnt need cliping
2967 }
2968
2969 static inline int get_block_bits(SnowContext *s, int x, int y, int w){
2970     const int b_stride = s->b_width << s->block_max_depth;
2971     const int b_height = s->b_height<< s->block_max_depth;
2972     int index= x + y*b_stride;
2973     BlockNode *b     = &s->block[index];
2974     BlockNode *left  = x ? &s->block[index-1] : &null_block;
2975     BlockNode *top   = y ? &s->block[index-b_stride] : &null_block;
2976     BlockNode *tl    = y && x ? &s->block[index-b_stride-1] : left;
2977     BlockNode *tr    = y && x+w<b_stride ? &s->block[index-b_stride+w] : tl;
2978     int dmx, dmy;
2979 //  int mx_context= av_log2(2*ABS(left->mx - top->mx));
2980 //  int my_context= av_log2(2*ABS(left->my - top->my));
2981
2982     if(x<0 || x>=b_stride || y>=b_height)
2983         return 0;
2984     dmx= b->mx - mid_pred(left->mx, top->mx, tr->mx);
2985     dmy= b->my - mid_pred(left->my, top->my, tr->my);
2986 /*
2987 1            0      0
2988 01X          1-2    1
2989 001XX        3-6    2-3
2990 0001XXX      7-14   4-7
2991 00001XXXX   15-30   8-15
2992 */
2993 //FIXME try accurate rate
2994 //FIXME intra and inter predictors if surrounding blocks arent the same type
2995     if(b->type & BLOCK_INTRA){
2996         return 3+2*( av_log2(2*ABS(left->color[0] - b->color[0]))
2997                    + av_log2(2*ABS(left->color[1] - b->color[1]))
2998                    + av_log2(2*ABS(left->color[2] - b->color[2])));
2999     }else
3000         return 2*(1 + av_log2(2*ABS(dmx))
3001                     + av_log2(2*ABS(dmy))); //FIXME kill the 2* can be merged in lambda
3002 }
3003
3004 static int get_block_rd(SnowContext *s, int mb_x, int mb_y, int plane_index, const uint8_t *obmc_edged){
3005     Plane *p= &s->plane[plane_index];
3006     const int block_size = MB_SIZE >> s->block_max_depth;
3007     const int block_w    = plane_index ? block_size/2 : block_size;
3008     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];
3009     const int obmc_stride= plane_index ? block_size : 2*block_size;
3010     const int ref_stride= s->current_picture.linesize[plane_index];
3011     uint8_t *ref= s->   last_picture.data[plane_index];
3012     uint8_t *dst= s->current_picture.data[plane_index];
3013     uint8_t *src= s->  input_picture.data[plane_index];
3014     DWTELEM *pred= (DWTELEM*)s->m.obmc_scratchpad + plane_index*block_size*block_size*4;
3015     uint8_t cur[ref_stride*2*MB_SIZE]; //FIXME alignment
3016     uint8_t tmp[ref_stride*(2*MB_SIZE+5)];
3017     const int b_stride = s->b_width << s->block_max_depth;
3018     const int b_height = s->b_height<< s->block_max_depth;
3019     const int w= p->width;
3020     const int h= p->height;
3021     int distortion;
3022     int rate= 0;
3023     const int penalty_factor= get_penalty_factor(s->lambda, s->lambda2, s->avctx->me_cmp);
3024     int sx= block_w*mb_x - block_w/2;
3025     int sy= block_w*mb_y - block_w/2;
3026     int x0= FFMAX(0,-sx);
3027     int y0= FFMAX(0,-sy);
3028     int x1= FFMIN(block_w*2, w-sx);
3029     int y1= FFMIN(block_w*2, h-sy);
3030     int i,x,y;
3031
3032     pred_block(s, cur, ref, tmp, ref_stride, sx, sy, block_w*2, block_w*2, &s->block[mb_x + mb_y*b_stride], plane_index, w, h);
3033
3034     for(y=y0; y<y1; y++){
3035         const uint8_t *obmc1= obmc_edged + y*obmc_stride;
3036         const DWTELEM *pred1 = pred + y*obmc_stride;
3037         uint8_t *cur1 = cur + y*ref_stride;
3038         uint8_t *dst1 = dst + sx + (sy+y)*ref_stride;
3039         for(x=x0; x<x1; x++){
3040             int v = (cur1[x] * obmc1[x]) << (FRAC_BITS - LOG2_OBMC_MAX);
3041             v = (v + pred1[x]) >> FRAC_BITS;
3042             if(v&(~255)) v= ~(v>>31);
3043             dst1[x] = v;
3044         }
3045     }
3046
3047     /* copy the regions where obmc[] = (uint8_t)256 */
3048     if(LOG2_OBMC_MAX == 8
3049         && (mb_x == 0 || mb_x == b_stride-1)
3050         && (mb_y == 0 || mb_y == b_height-1)){
3051         if(mb_x == 0)
3052             x1 = block_w;
3053         else
3054             x0 = block_w;
3055         if(mb_y == 0)
3056             y1 = block_w;
3057         else
3058             y0 = block_w;
3059         for(y=y0; y<y1; y++)
3060             memcpy(dst + sx+x0 + (sy+y)*ref_stride, cur + x0 + y*ref_stride, x1-x0);
3061     }
3062
3063     //FIXME sad/ssd can be broken up, but wavelet cmp should be one 32x32 block
3064     if(block_w==16){
3065         distortion = 0;
3066         for(i=0; i<4; i++){
3067             int off = sx+16*(i&1) + (sy+16*(i>>1))*ref_stride;
3068             distortion += s->dsp.me_cmp[0](&s->m, src + off, dst + off, ref_stride, 16);
3069         }
3070     }else{
3071         assert(block_w==8);
3072         distortion = s->dsp.me_cmp[0](&s->m, src + sx + sy*ref_stride, dst + sx + sy*ref_stride, ref_stride, block_w*2);
3073     }
3074
3075     if(plane_index==0){
3076         for(i=0; i<4; i++){
3077 /* ..RRr
3078  * .RXx.
3079  * rxx..
3080  */
3081             rate += get_block_bits(s, mb_x + (i&1) - (i>>1), mb_y + (i>>1), 1);
3082         }
3083         if(mb_x == b_stride-2)
3084             rate += get_block_bits(s, mb_x + 1, mb_y + 1, 1);
3085     }
3086     return distortion + rate*penalty_factor;
3087 }
3088
3089 static int get_4block_rd(SnowContext *s, int mb_x, int mb_y, int plane_index){
3090     int i, y2;
3091     Plane *p= &s->plane[plane_index];
3092     const int block_size = MB_SIZE >> s->block_max_depth;
3093     const int block_w    = plane_index ? block_size/2 : block_size;
3094     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];
3095     const int obmc_stride= plane_index ? block_size : 2*block_size;
3096     const int ref_stride= s->current_picture.linesize[plane_index];
3097     uint8_t *ref= s->   last_picture.data[plane_index];
3098     uint8_t *dst= s->current_picture.data[plane_index];
3099     uint8_t *src= s-> input_picture.data[plane_index];
3100     const static DWTELEM zero_dst[4096]; //FIXME
3101     const int b_stride = s->b_width << s->block_max_depth;
3102     const int b_height = s->b_height<< s->block_max_depth;
3103     const int w= p->width;
3104     const int h= p->height;
3105     int distortion= 0;
3106     int rate= 0;
3107     const int penalty_factor= get_penalty_factor(s->lambda, s->lambda2, s->avctx->me_cmp);
3108
3109     for(i=0; i<9; i++){
3110         int mb_x2= mb_x + (i%3) - 1;
3111         int mb_y2= mb_y + (i/3) - 1;
3112         int x= block_w*mb_x2 + block_w/2;
3113         int y= block_w*mb_y2 + block_w/2;
3114
3115         add_yblock(s, zero_dst, dst, ref, obmc,
3116                    x, y, block_w, block_w, w, h, /*dst_stride*/0, ref_stride, obmc_stride, mb_x2, mb_y2, 1, 1, plane_index);
3117
3118         //FIXME find a cleaner/simpler way to skip the outside stuff
3119         for(y2= y; y2<0; y2++)
3120             memcpy(dst + x + y2*ref_stride, src + x + y2*ref_stride, block_w);
3121         for(y2= h; y2<y+block_w; y2++)
3122             memcpy(dst + x + y2*ref_stride, src + x + y2*ref_stride, block_w);
3123         if(x<0){
3124             for(y2= y; y2<y+block_w; y2++)
3125                 memcpy(dst + x + y2*ref_stride, src + x + y2*ref_stride, -x);
3126         }
3127         if(x+block_w > w){
3128             for(y2= y; y2<y+block_w; y2++)
3129                 memcpy(dst + w + y2*ref_stride, src + w + y2*ref_stride, x+block_w - w);
3130         }
3131
3132         assert(block_w== 8 || block_w==16);
3133         distortion += s->dsp.me_cmp[block_w==8](&s->m, src + x + y*ref_stride, dst + x + y*ref_stride, ref_stride, block_w);
3134     }
3135
3136     if(plane_index==0){
3137         BlockNode *b= &s->block[mb_x+mb_y*b_stride];
3138         int merged= same_block(b,b+1) && same_block(b,b+b_stride) && same_block(b,b+b_stride+1);
3139
3140 /* ..RRRr
3141  * .RXXx.
3142  * .RXXx.
3143  * rxxx.
3144  */
3145         if(merged)
3146             rate = get_block_bits(s, mb_x, mb_y, 2);
3147         for(i=merged?4:0; i<9; i++){
3148             static const int dxy[9][2] = {{0,0},{1,0},{0,1},{1,1},{2,0},{2,1},{-1,2},{0,2},{1,2}};
3149             rate += get_block_bits(s, mb_x + dxy[i][0], mb_y + dxy[i][1], 1);
3150         }
3151     }
3152     return distortion + rate*penalty_factor;
3153 }
3154
3155 static 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){
3156     const int b_stride= s->b_width << s->block_max_depth;
3157     BlockNode *block= &s->block[mb_x + mb_y * b_stride];
3158     BlockNode backup= *block;
3159     int rd, index, value;
3160
3161     assert(mb_x>=0 && mb_y>=0);
3162     assert(mb_x<b_stride);
3163
3164     if(intra){
3165         block->color[0] = p[0];
3166         block->color[1] = p[1];
3167         block->color[2] = p[2];
3168         block->type |= BLOCK_INTRA;
3169     }else{
3170         index= (p[0] + 31*p[1]) & (ME_CACHE_SIZE-1);
3171         value= s->me_cache_generation + (p[0]>>10) + (p[1]<<6);
3172         if(s->me_cache[index] == value)
3173             return 0;
3174         s->me_cache[index]= value;
3175
3176         block->mx= p[0];
3177         block->my= p[1];
3178         block->type &= ~BLOCK_INTRA;
3179     }
3180
3181     rd= get_block_rd(s, mb_x, mb_y, 0, obmc_edged);
3182
3183 //FIXME chroma
3184     if(rd < *best_rd){
3185         *best_rd= rd;
3186         return 1;
3187     }else{
3188         *block= backup;
3189         return 0;
3190     }
3191 }
3192
3193 /* special case for int[2] args we discard afterward, fixes compilation prob with gcc 2.95 */
3194 static 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){
3195     int p[2] = {p0, p1};
3196     return check_block(s, mb_x, mb_y, p, 0, obmc_edged, best_rd);
3197 }
3198
3199 static always_inline int check_4block_inter(SnowContext *s, int mb_x, int mb_y, int p0, int p1, int *best_rd){
3200     const int b_stride= s->b_width << s->block_max_depth;
3201     BlockNode *block= &s->block[mb_x + mb_y * b_stride];
3202     BlockNode backup[4]= {block[0], block[1], block[b_stride], block[b_stride+1]};
3203     int rd, index, value;
3204
3205     assert(mb_x>=0 && mb_y>=0);
3206     assert(mb_x<b_stride);
3207     assert(((mb_x|mb_y)&1) == 0);
3208
3209     index= (p0 + 31*p1) & (ME_CACHE_SIZE-1);
3210     value= s->me_cache_generation + (p0>>10) + (p1<<6);
3211     if(s->me_cache[index] == value)
3212         return 0;
3213     s->me_cache[index]= value;
3214
3215     block->mx= p0;
3216     block->my= p1;
3217     block->type &= ~BLOCK_INTRA;
3218     block[1]= block[b_stride]= block[b_stride+1]= *block;
3219
3220     rd= get_4block_rd(s, mb_x, mb_y, 0);
3221
3222 //FIXME chroma
3223     if(rd < *best_rd){
3224         *best_rd= rd;
3225         return 1;
3226     }else{
3227         block[0]= backup[0];
3228         block[1]= backup[1];
3229         block[b_stride]= backup[2];
3230         block[b_stride+1]= backup[3];
3231         return 0;
3232     }
3233 }
3234
3235 static void iterative_me(SnowContext *s){
3236     int pass, mb_x, mb_y;
3237     const int b_width = s->b_width  << s->block_max_depth;
3238     const int b_height= s->b_height << s->block_max_depth;
3239     const int b_stride= b_width;
3240     int color[3];
3241
3242     {
3243         RangeCoder r = s->c;
3244         uint8_t state[sizeof(s->block_state)];
3245         memcpy(state, s->block_state, sizeof(s->block_state));
3246         for(mb_y= 0; mb_y<s->b_height; mb_y++)
3247             for(mb_x= 0; mb_x<s->b_width; mb_x++)
3248                 encode_q_branch(s, 0, mb_x, mb_y);
3249         s->c = r;
3250         memcpy(s->block_state, state, sizeof(s->block_state));
3251     }
3252
3253     for(pass=0; pass<50; pass++){
3254         int change= 0;
3255
3256         for(mb_y= 0; mb_y<b_height; mb_y++){
3257             for(mb_x= 0; mb_x<b_width; mb_x++){
3258                 int dia_change, i, j;
3259                 int best_rd= INT_MAX;
3260                 BlockNode backup;
3261                 const int index= mb_x + mb_y * b_stride;
3262                 BlockNode *block= &s->block[index];
3263                 BlockNode *tb =                   mb_y            ? &s->block[index-b_stride  ] : &null_block;
3264                 BlockNode *lb = mb_x                              ? &s->block[index         -1] : &null_block;
3265                 BlockNode *rb = mb_x+1<b_width                    ? &s->block[index         +1] : &null_block;
3266                 BlockNode *bb =                   mb_y+1<b_height ? &s->block[index+b_stride  ] : &null_block;
3267                 BlockNode *tlb= mb_x           && mb_y            ? &s->block[index-b_stride-1] : &null_block;
3268                 BlockNode *trb= mb_x+1<b_width && mb_y            ? &s->block[index-b_stride+1] : &null_block;
3269                 BlockNode *blb= mb_x           && mb_y+1<b_height ? &s->block[index+b_stride-1] : &null_block;
3270                 BlockNode *brb= mb_x+1<b_width && mb_y+1<b_height ? &s->block[index+b_stride+1] : &null_block;
3271                 const int b_w= (MB_SIZE >> s->block_max_depth);
3272                 uint8_t obmc_edged[b_w*2][b_w*2];
3273
3274                 if(pass && (block->type & BLOCK_OPT))
3275                     continue;
3276                 block->type |= BLOCK_OPT;
3277
3278                 backup= *block;
3279
3280                 if(!s->me_cache_generation)
3281                     memset(s->me_cache, 0, sizeof(s->me_cache));
3282                 s->me_cache_generation += 1<<22;
3283
3284                 //FIXME precalc
3285                 {
3286                     int x, y;
3287                     memcpy(obmc_edged, obmc_tab[s->block_max_depth], b_w*b_w*4);
3288                     if(mb_x==0)
3289                         for(y=0; y<b_w*2; y++)
3290                             memset(obmc_edged[y], obmc_edged[y][0] + obmc_edged[y][b_w-1], b_w);
3291                     if(mb_x==b_stride-1)
3292                         for(y=0; y<b_w*2; y++)
3293                             memset(obmc_edged[y]+b_w, obmc_edged[y][b_w] + obmc_edged[y][b_w*2-1], b_w);
3294                     if(mb_y==0){
3295                         for(x=0; x<b_w*2; x++)
3296                             obmc_edged[0][x] += obmc_edged[b_w-1][x];
3297                         for(y=1; y<b_w; y++)
3298                             memcpy(obmc_edged[y], obmc_edged[0], b_w*2);
3299                     }
3300                     if(mb_y==b_height-1){
3301                         for(x=0; x<b_w*2; x++)
3302                             obmc_edged[b_w*2-1][x] += obmc_edged[b_w][x];
3303                         for(y=b_w; y<b_w*2-1; y++)
3304                             memcpy(obmc_edged[y], obmc_edged[b_w*2-1], b_w*2);
3305                     }
3306                 }
3307
3308                 //skip stuff outside the picture
3309                 if(mb_x==0 || mb_y==0 || mb_x==b_width-1 || mb_y==b_height-1)
3310                 {
3311                     uint8_t *src= s->  input_picture.data[0];
3312                     uint8_t *dst= s->current_picture.data[0];
3313                     const int stride= s->current_picture.linesize[0];
3314                     const int block_w= MB_SIZE >> s->block_max_depth;
3315                     const int sx= block_w*mb_x - block_w/2;
3316                     const int sy= block_w*mb_y - block_w/2;
3317                     const int w= s->plane[0].width;
3318                     const int h= s->plane[0].height;
3319                     int y;
3320
3321                     for(y=sy; y<0; y++)
3322                         memcpy(dst + sx + y*stride, src + sx + y*stride, block_w*2);
3323                     for(y=h; y<sy+block_w*2; y++)
3324                         memcpy(dst + sx + y*stride, src + sx + y*stride, block_w*2);
3325                     if(sx<0){
3326                         for(y=sy; y<sy+block_w*2; y++)
3327                             memcpy(dst + sx + y*stride, src + sx + y*stride, -sx);
3328                     }
3329                     if(sx+block_w*2 > w){
3330                         for(y=sy; y<sy+block_w*2; y++)
3331                             memcpy(dst + w + y*stride, src + w + y*stride, sx+block_w*2 - w);
3332                     }
3333                 }
3334
3335                 // intra(black) = neighbors' contribution to the current block
3336                 for(i=0; i<3; i++)
3337                     color[i]= get_dc(s, mb_x, mb_y, i);
3338
3339                 // get previous score (cant be cached due to OBMC)
3340                 if(pass > 0 && (block->type&BLOCK_INTRA)){
3341                     int color0[3]= {block->color[0], block->color[1], block->color[2]};
3342                     check_block(s, mb_x, mb_y, color0, 1, *obmc_edged, &best_rd);
3343                 }else
3344                     check_block_inter(s, mb_x, mb_y, block->mx, block->my, *obmc_edged, &best_rd);
3345
3346                 check_block_inter(s, mb_x, mb_y, 0, 0, *obmc_edged, &best_rd);
3347                 check_block_inter(s, mb_x, mb_y, tb->mx, tb->my, *obmc_edged, &best_rd);
3348                 check_block_inter(s, mb_x, mb_y, lb->mx, lb->my, *obmc_edged, &best_rd);
3349                 check_block_inter(s, mb_x, mb_y, rb->mx, rb->my, *obmc_edged, &best_rd);
3350                 check_block_inter(s, mb_x, mb_y, bb->mx, bb->my, *obmc_edged, &best_rd);
3351
3352                 /* fullpel ME */
3353                 //FIXME avoid subpel interpol / round to nearest integer
3354                 do{
3355                     dia_change=0;
3356                     for(i=0; i<FFMAX(s->avctx->dia_size, 1); i++){
3357                         for(j=0; j<i; j++){
3358                             dia_change |= check_block_inter(s, mb_x, mb_y, block->mx+4*(i-j), block->my+(4*j), *obmc_edged, &best_rd);
3359                             dia_change |= check_block_inter(s, mb_x, mb_y, block->mx-4*(i-j), block->my-(4*j), *obmc_edged, &best_rd);
3360                             dia_change |= check_block_inter(s, mb_x, mb_y, block->mx+4*(i-j), block->my-(4*j), *obmc_edged, &best_rd);
3361                             dia_change |= check_block_inter(s, mb_x, mb_y, block->mx-4*(i-j), block->my+(4*j), *obmc_edged, &best_rd);
3362                         }
3363                     }
3364                 }while(dia_change);
3365                 /* subpel ME */
3366                 do{
3367                     static const int square[8][2]= {{+1, 0},{-1, 0},{ 0,+1},{ 0,-1},{+1,+1},{-1,-1},{+1,-1},{-1,+1},};
3368                     dia_change=0;
3369                     for(i=0; i<8; i++)
3370                         dia_change |= check_block_inter(s, mb_x, mb_y, block->mx+square[i][0], block->my+square[i][1], *obmc_edged, &best_rd);
3371                 }while(dia_change);
3372                 //FIXME or try the standard 2 pass qpel or similar
3373 #if 1
3374                 check_block(s, mb_x, mb_y, color, 1, *obmc_edged, &best_rd);
3375                 //FIXME RD style color selection
3376 #endif
3377                 if(!same_block(block, &backup)){
3378                     if(tb != &null_block) tb ->type &= ~BLOCK_OPT;
3379                     if(lb != &null_block) lb ->type &= ~BLOCK_OPT;
3380                     if(rb != &null_block) rb ->type &= ~BLOCK_OPT;
3381                     if(bb != &null_block) bb ->type &= ~BLOCK_OPT;
3382                     if(tlb!= &null_block) tlb->type &= ~BLOCK_OPT;
3383                     if(trb!= &null_block) trb->type &= ~BLOCK_OPT;
3384                     if(blb!= &null_block) blb->type &= ~BLOCK_OPT;
3385                     if(brb!= &null_block) brb->type &= ~BLOCK_OPT;
3386                     change ++;
3387                 }
3388             }
3389         }
3390         av_log(NULL, AV_LOG_ERROR, "pass:%d changed:%d\n", pass, change);
3391         if(!change)
3392             break;
3393     }
3394
3395     if(s->block_max_depth == 1){
3396         int change= 0;
3397         for(mb_y= 0; mb_y<b_height; mb_y+=2){
3398             for(mb_x= 0; mb_x<b_width; mb_x+=2){
3399                 int dia_change, i, j;
3400                 int best_rd, init_rd;
3401                 const int index= mb_x + mb_y * b_stride;
3402                 BlockNode *b[4];
3403
3404                 b[0]= &s->block[index];
3405                 b[1]= b[0]+1;
3406                 b[2]= b[0]+b_stride;
3407                 b[3]= b[2]+1;
3408                 if(same_block(b[0], b[1]) &&
3409                    same_block(b[0], b[2]) &&
3410                    same_block(b[0], b[3]))
3411                     continue;
3412
3413                 if(!s->me_cache_generation)
3414                     memset(s->me_cache, 0, sizeof(s->me_cache));
3415                 s->me_cache_generation += 1<<22;
3416
3417                 init_rd= best_rd= get_4block_rd(s, mb_x, mb_y, 0);
3418
3419                 check_4block_inter(s, mb_x, mb_y,
3420                                    (b[0]->mx + b[1]->mx + b[2]->mx + b[3]->mx + 2) >> 2,
3421                                    (b[0]->my + b[1]->my + b[2]->my + b[3]->my + 2) >> 2, &best_rd);
3422
3423                 for(i=0; i<4; i++)
3424                     if(!(b[i]->type&BLOCK_INTRA))
3425                         check_4block_inter(s, mb_x, mb_y, b[i]->mx, b[i]->my, &best_rd);
3426
3427                 if(init_rd != best_rd)
3428                     change++;
3429             }
3430         }
3431         av_log(NULL, AV_LOG_ERROR, "pass:4mv changed:%d\n", change*4);
3432     }
3433 }
3434
3435 static void quantize(SnowContext *s, SubBand *b, DWTELEM *src, int stride, int bias){
3436     const int level= b->level;
3437     const int w= b->width;
3438     const int h= b->height;
3439     const int qlog= clip(s->qlog + b->qlog, 0, QROOT*16);
3440     const int qmul= qexp[qlog&(QROOT-1)]<<(qlog>>QSHIFT);
3441     int x,y, thres1, thres2;
3442 //    START_TIMER
3443
3444     if(s->qlog == LOSSLESS_QLOG) return;
3445
3446     bias= bias ? 0 : (3*qmul)>>3;
3447     thres1= ((qmul - bias)>>QEXPSHIFT) - 1;
3448     thres2= 2*thres1;
3449
3450     if(!bias){
3451         for(y=0; y<h; y++){
3452             for(x=0; x<w; x++){
3453                 int i= src[x + y*stride];
3454
3455                 if((unsigned)(i+thres1) > thres2){
3456                     if(i>=0){
3457                         i<<= QEXPSHIFT;
3458                         i/= qmul; //FIXME optimize
3459                         src[x + y*stride]=  i;
3460                     }else{
3461                         i= -i;
3462                         i<<= QEXPSHIFT;
3463                         i/= qmul; //FIXME optimize
3464                         src[x + y*stride]= -i;
3465                     }
3466                 }else
3467                     src[x + y*stride]= 0;
3468             }
3469         }
3470     }else{
3471         for(y=0; y<h; y++){
3472             for(x=0; x<w; x++){
3473                 int i= src[x + y*stride];
3474
3475                 if((unsigned)(i+thres1) > thres2){
3476                     if(i>=0){
3477                         i<<= QEXPSHIFT;
3478                         i= (i + bias) / qmul; //FIXME optimize
3479                         src[x + y*stride]=  i;
3480                     }else{
3481                         i= -i;
3482                         i<<= QEXPSHIFT;
3483                         i= (i + bias) / qmul; //FIXME optimize
3484                         src[x + y*stride]= -i;
3485                     }
3486                 }else
3487                     src[x + y*stride]= 0;
3488             }
3489         }
3490     }
3491     if(level+1 == s->spatial_decomposition_count){
3492 //        STOP_TIMER("quantize")
3493     }
3494 }
3495
3496 static void dequantize_slice_buffered(SnowContext *s, slice_buffer * sb, SubBand *b, DWTELEM *src, int stride, int start_y, int end_y){
3497     const int w= b->width;
3498     const int qlog= clip(s->qlog + b->qlog, 0, QROOT*16);
3499     const int qmul= qexp[qlog&(QROOT-1)]<<(qlog>>QSHIFT);
3500     const int qadd= (s->qbias*qmul)>>QBIAS_SHIFT;
3501     int x,y;
3502     START_TIMER
3503
3504     if(s->qlog == LOSSLESS_QLOG) return;
3505
3506     for(y=start_y; y<end_y; y++){
3507 //        DWTELEM * line = slice_buffer_get_line_from_address(sb, src + (y * stride));
3508         DWTELEM * line = slice_buffer_get_line(sb, (y * b->stride_line) + b->buf_y_offset) + b->buf_x_offset;
3509         for(x=0; x<w; x++){
3510             int i= line[x];
3511             if(i<0){
3512                 line[x]= -((-i*qmul + qadd)>>(QEXPSHIFT)); //FIXME try different bias
3513             }else if(i>0){
3514                 line[x]=  (( i*qmul + qadd)>>(QEXPSHIFT));
3515             }
3516         }
3517     }
3518     if(w > 200 /*level+1 == s->spatial_decomposition_count*/){
3519         STOP_TIMER("dquant")
3520     }
3521 }
3522
3523 static void dequantize(SnowContext *s, SubBand *b, DWTELEM *src, int stride){
3524     const int w= b->width;
3525     const int h= b->height;
3526     const int qlog= clip(s->qlog + b->qlog, 0, QROOT*16);
3527     const int qmul= qexp[qlog&(QROOT-1)]<<(qlog>>QSHIFT);
3528     const int qadd= (s->qbias*qmul)>>QBIAS_SHIFT;
3529     int x,y;
3530     START_TIMER
3531
3532     if(s->qlog == LOSSLESS_QLOG) return;
3533
3534     for(y=0; y<h; y++){
3535         for(x=0; x<w; x++){
3536             int i= src[x + y*stride];
3537             if(i<0){
3538                 src[x + y*stride]= -((-i*qmul + qadd)>>(QEXPSHIFT)); //FIXME try different bias
3539             }else if(i>0){
3540                 src[x + y*stride]=  (( i*qmul + qadd)>>(QEXPSHIFT));
3541             }
3542         }
3543     }
3544     if(w > 200 /*level+1 == s->spatial_decomposition_count*/){
3545         STOP_TIMER("dquant")
3546     }
3547 }
3548
3549 static void decorrelate(SnowContext *s, SubBand *b, DWTELEM *src, int stride, int inverse, int use_median){
3550     const int w= b->width;
3551     const int h= b->height;
3552     int x,y;
3553
3554     for(y=h-1; y>=0; y--){
3555         for(x=w-1; x>=0; x--){
3556             int i= x + y*stride;
3557
3558             if(x){
3559                 if(use_median){
3560                     if(y && x+1<w) src[i] -= mid_pred(src[i - 1], src[i - stride], src[i - stride + 1]);
3561                     else  src[i] -= src[i - 1];
3562                 }else{
3563                     if(y) src[i] -= mid_pred(src[i - 1], src[i - stride], src[i - 1] + src[i - stride] - src[i - 1 - stride]);
3564                     else  src[i] -= src[i - 1];
3565                 }
3566             }else{
3567                 if(y) src[i] -= src[i - stride];
3568             }
3569         }
3570     }
3571 }
3572
3573 static void correlate_slice_buffered(SnowContext *s, slice_buffer * sb, SubBand *b, DWTELEM *src, int stride, int inverse, int use_median, int start_y, int end_y){
3574     const int w= b->width;
3575     int x,y;
3576
3577 //    START_TIMER
3578
3579     DWTELEM * line;
3580     DWTELEM * prev;
3581
3582     if (start_y != 0)
3583         line = slice_buffer_get_line(sb, ((start_y - 1) * b->stride_line) + b->buf_y_offset) + b->buf_x_offset;
3584
3585     for(y=start_y; y<end_y; y++){
3586         prev = line;
3587 //        line = slice_buffer_get_line_from_address(sb, src + (y * stride));
3588         line = slice_buffer_get_line(sb, (y * b->stride_line) + b->buf_y_offset) + b->buf_x_offset;
3589         for(x=0; x<w; x++){
3590             if(x){
3591                 if(use_median){
3592                     if(y && x+1<w) line[x] += mid_pred(line[x - 1], prev[x], prev[x + 1]);
3593                     else  line[x] += line[x - 1];
3594                 }else{
3595                     if(y) line[x] += mid_pred(line[x - 1], prev[x], line[x - 1] + prev[x] - prev[x - 1]);
3596                     else  line[x] += line[x - 1];
3597                 }
3598             }else{
3599                 if(y) line[x] += prev[x];
3600             }
3601         }
3602     }
3603
3604 //    STOP_TIMER("correlate")
3605 }
3606
3607 static void correlate(SnowContext *s, SubBand *b, DWTELEM *src, int stride, int inverse, int use_median){
3608     const int w= b->width;
3609     const int h= b->height;
3610     int x,y;
3611
3612     for(y=0; y<h; y++){
3613         for(x=0; x<w; x++){
3614             int i= x + y*stride;
3615
3616             if(x){
3617                 if(use_median){
3618                     if(y && x+1<w) src[i] += mid_pred(src[i - 1], src[i - stride], src[i - stride + 1]);
3619                     else  src[i] += src[i - 1];
3620                 }else{
3621                     if(y) src[i] += mid_pred(src[i - 1], src[i - stride], src[i - 1] + src[i - stride] - src[i - 1 - stride]);
3622                     else  src[i] += src[i - 1];
3623                 }
3624             }else{
3625                 if(y) src[i] += src[i - stride];
3626             }
3627         }
3628     }
3629 }
3630
3631 static void encode_header(SnowContext *s){
3632     int plane_index, level, orientation;
3633     uint8_t kstate[32];
3634
3635     memset(kstate, MID_STATE, sizeof(kstate));
3636
3637     put_rac(&s->c, kstate, s->keyframe);
3638     if(s->keyframe || s->always_reset)
3639         reset_contexts(s);
3640     if(s->keyframe){
3641         put_symbol(&s->c, s->header_state, s->version, 0);
3642         put_rac(&s->c, s->header_state, s->always_reset);
3643         put_symbol(&s->c, s->header_state, s->temporal_decomposition_type, 0);
3644         put_symbol(&s->c, s->header_state, s->temporal_decomposition_count, 0);
3645         put_symbol(&s->c, s->header_state, s->spatial_decomposition_count, 0);
3646         put_symbol(&s->c, s->header_state, s->colorspace_type, 0);
3647         put_symbol(&s->c, s->header_state, s->chroma_h_shift, 0);
3648         put_symbol(&s->c, s->header_state, s->chroma_v_shift, 0);
3649         put_rac(&s->c, s->header_state, s->spatial_scalability);
3650 //        put_rac(&s->c, s->header_state, s->rate_scalability);
3651
3652         for(plane_index=0; plane_index<2; plane_index++){
3653             for(level=0; level<s->spatial_decomposition_count; level++){
3654                 for(orientation=level ? 1:0; orientation<4; orientation++){
3655                     if(orientation==2) continue;
3656                     put_symbol(&s->c, s->header_state, s->plane[plane_index].band[level][orientation].qlog, 1);
3657                 }
3658             }
3659         }
3660     }
3661     put_symbol(&s->c, s->header_state, s->spatial_decomposition_type, 0);
3662     put_symbol(&s->c, s->header_state, s->qlog, 1);
3663     put_symbol(&s->c, s->header_state, s->mv_scale, 0);
3664     put_symbol(&s->c, s->header_state, s->qbias, 1);
3665     put_symbol(&s->c, s->header_state, s->block_max_depth, 0);
3666 }
3667
3668 static int decode_header(SnowContext *s){
3669     int plane_index, level, orientation;
3670     uint8_t kstate[32];
3671
3672     memset(kstate, MID_STATE, sizeof(kstate));
3673
3674     s->keyframe= get_rac(&s->c, kstate);
3675     if(s->keyframe || s->always_reset)
3676         reset_contexts(s);
3677     if(s->keyframe){
3678         s->version= get_symbol(&s->c, s->header_state, 0);
3679         if(s->version>0){
3680             av_log(s->avctx, AV_LOG_ERROR, "version %d not supported", s->version);
3681             return -1;
3682         }
3683         s->always_reset= get_rac(&s->c, s->header_state);
3684         s->temporal_decomposition_type= get_symbol(&s->c, s->header_state, 0);
3685         s->temporal_decomposition_count= get_symbol(&s->c, s->header_state, 0);
3686         s->spatial_decomposition_count= get_symbol(&s->c, s->header_state, 0);
3687         s->colorspace_type= get_symbol(&s->c, s->header_state, 0);
3688         s->chroma_h_shift= get_symbol(&s->c, s->header_state, 0);
3689         s->chroma_v_shift= get_symbol(&s->c, s->header_state, 0);
3690         s->spatial_scalability= get_rac(&s->c, s->header_state);
3691 //        s->rate_scalability= get_rac(&s->c, s->header_state);
3692
3693         for(plane_index=0; plane_index<3; plane_index++){
3694             for(level=0; level<s->spatial_decomposition_count; level++){
3695                 for(orientation=level ? 1:0; orientation<4; orientation++){
3696                     int q;
3697                     if     (plane_index==2) q= s->plane[1].band[level][orientation].qlog;
3698                     else if(orientation==2) q= s->plane[plane_index].band[level][1].qlog;
3699                     else                    q= get_symbol(&s->c, s->header_state, 1);
3700                     s->plane[plane_index].band[level][orientation].qlog= q;
3701                 }
3702             }
3703         }
3704     }
3705
3706     s->spatial_decomposition_type= get_symbol(&s->c, s->header_state, 0);
3707     if(s->spatial_decomposition_type > 2){
3708         av_log(s->avctx, AV_LOG_ERROR, "spatial_decomposition_type %d not supported", s->spatial_decomposition_type);
3709         return -1;
3710     }
3711
3712     s->qlog= get_symbol(&s->c, s->header_state, 1);
3713     s->mv_scale= get_symbol(&s->c, s->header_state, 0);
3714     s->qbias= get_symbol(&s->c, s->header_state, 1);
3715     s->block_max_depth= get_symbol(&s->c, s->header_state, 0);
3716     if(s->block_max_depth > 1 || s->block_max_depth < 0){
3717         av_log(s->avctx, AV_LOG_ERROR, "block_max_depth= %d is too large", s->block_max_depth);
3718         s->block_max_depth= 0;
3719         return -1;
3720     }
3721
3722     return 0;
3723 }
3724
3725 static void init_qexp(void){
3726     int i;
3727     double v=128;
3728
3729     for(i=0; i<QROOT; i++){
3730         qexp[i]= lrintf(v);
3731         v *= pow(2, 1.0 / QROOT);
3732     }
3733 }
3734
3735 static int common_init(AVCodecContext *avctx){
3736     SnowContext *s = avctx->priv_data;
3737     int width, height;
3738     int level, orientation, plane_index, dec;
3739
3740     s->avctx= avctx;
3741
3742     dsputil_init(&s->dsp, avctx);
3743
3744 #define mcf(dx,dy)\
3745     s->dsp.put_qpel_pixels_tab       [0][dy+dx/4]=\
3746     s->dsp.put_no_rnd_qpel_pixels_tab[0][dy+dx/4]=\
3747         s->dsp.put_h264_qpel_pixels_tab[0][dy+dx/4];\
3748     s->dsp.put_qpel_pixels_tab       [1][dy+dx/4]=\
3749     s->dsp.put_no_rnd_qpel_pixels_tab[1][dy+dx/4]=\
3750         s->dsp.put_h264_qpel_pixels_tab[1][dy+dx/4];
3751
3752     mcf( 0, 0)
3753     mcf( 4, 0)
3754     mcf( 8, 0)
3755     mcf(12, 0)
3756     mcf( 0, 4)
3757     mcf( 4, 4)
3758     mcf( 8, 4)
3759     mcf(12, 4)
3760     mcf( 0, 8)
3761     mcf( 4, 8)
3762     mcf( 8, 8)
3763     mcf(12, 8)
3764     mcf( 0,12)
3765     mcf( 4,12)
3766     mcf( 8,12)
3767     mcf(12,12)
3768
3769 #define mcfh(dx,dy)\
3770     s->dsp.put_pixels_tab       [0][dy/4+dx/8]=\
3771     s->dsp.put_no_rnd_pixels_tab[0][dy/4+dx/8]=\
3772         mc_block_hpel ## dx ## dy ## 16;\
3773     s->dsp.put_pixels_tab       [1][dy/4+dx/8]=\
3774     s->dsp.put_no_rnd_pixels_tab[1][dy/4+dx/8]=\
3775         mc_block_hpel ## dx ## dy ## 8;
3776
3777     mcfh(0, 0)
3778     mcfh(8, 0)
3779     mcfh(0, 8)
3780     mcfh(8, 8)
3781
3782     if(!qexp[0])
3783         init_qexp();
3784
3785     dec= s->spatial_decomposition_count= 5;
3786     s->spatial_decomposition_type= avctx->prediction_method; //FIXME add decorrelator type r transform_type
3787
3788     s->chroma_h_shift= 1; //FIXME XXX
3789     s->chroma_v_shift= 1;
3790
3791 //    dec += FFMAX(s->chroma_h_shift, s->chroma_v_shift);
3792
3793     width= s->avctx->width;
3794     height= s->avctx->height;
3795
3796     s->spatial_dwt_buffer= av_mallocz(width*height*sizeof(DWTELEM));
3797
3798     s->mv_scale= (s->avctx->flags & CODEC_FLAG_QPEL) ? 2 : 4;
3799     s->block_max_depth= (s->avctx->flags & CODEC_FLAG_4MV) ? 1 : 0;
3800
3801     for(plane_index=0; plane_index<3; plane_index++){
3802         int w= s->avctx->width;
3803         int h= s->avctx->height;
3804
3805         if(plane_index){
3806             w>>= s->chroma_h_shift;
3807             h>>= s->chroma_v_shift;
3808         }
3809         s->plane[plane_index].width = w;
3810         s->plane[plane_index].height= h;
3811 //av_log(NULL, AV_LOG_DEBUG, "%d %d\n", w, h);
3812         for(level=s->spatial_decomposition_count-1; level>=0; level--){
3813             for(orientation=level ? 1 : 0; orientation<4; orientation++){
3814                 SubBand *b= &s->plane[plane_index].band[level][orientation];
3815
3816                 b->buf= s->spatial_dwt_buffer;
3817                 b->level= level;
3818                 b->stride= s->plane[plane_index].width << (s->spatial_decomposition_count - level);
3819                 b->width = (w + !(orientation&1))>>1;
3820                 b->height= (h + !(orientation>1))>>1;
3821
3822                 b->stride_line = 1 << (s->spatial_decomposition_count - level);
3823                 b->buf_x_offset = 0;
3824                 b->buf_y_offset = 0;
3825
3826                 if(orientation&1){
3827                     b->buf += (w+1)>>1;
3828                     b->buf_x_offset = (w+1)>>1;
3829                 }
3830                 if(orientation>1){
3831                     b->buf += b->stride>>1;
3832                     b->buf_y_offset = b->stride_line >> 1;
3833                 }
3834
3835                 if(level)
3836                     b->parent= &s->plane[plane_index].band[level-1][orientation];
3837                 b->x_coeff=av_mallocz(((b->width+1) * b->height+1)*sizeof(x_and_coeff));
3838             }
3839             w= (w+1)>>1;
3840             h= (h+1)>>1;
3841         }
3842     }
3843
3844     reset_contexts(s);
3845 /*
3846     width= s->width= avctx->width;
3847     height= s->height= avctx->height;
3848
3849     assert(width && height);
3850 */
3851     s->avctx->get_buffer(s->avctx, &s->mconly_picture);
3852
3853     return 0;
3854 }
3855
3856 static void ratecontrol_1pass(SnowContext *s, AVFrame *pict)
3857 {
3858     /* estimate the frame's complexity as a sum of weighted dwt coefs.
3859      * FIXME we know exact mv bits at this point,
3860      * but ratecontrol isn't set up to include them. */
3861     uint32_t coef_sum= 0;
3862     int level, orientation;
3863
3864     for(level=0; level<s->spatial_decomposition_count; level++){
3865         for(orientation=level ? 1 : 0; orientation<4; orientation++){
3866             SubBand *b= &s->plane[0].band[level][orientation];
3867             DWTELEM *buf= b->buf;
3868             const int w= b->width;
3869             const int h= b->height;
3870             const int stride= b->stride;
3871             const int qlog= clip(2*QROOT + b->qlog, 0, QROOT*16);
3872             const int qmul= qexp[qlog&(QROOT-1)]<<(qlog>>QSHIFT);
3873             const int qdiv= (1<<16)/qmul;
3874             int x, y;
3875             if(orientation==0)
3876                 decorrelate(s, b, buf, stride, 1, 0);
3877             for(y=0; y<h; y++)
3878                 for(x=0; x<w; x++)
3879                     coef_sum+= abs(buf[x+y*stride]) * qdiv >> 16;
3880             if(orientation==0)
3881                 correlate(s, b, buf, stride, 1, 0);
3882         }
3883     }
3884
3885     /* ugly, ratecontrol just takes a sqrt again */
3886     coef_sum = (uint64_t)coef_sum * coef_sum >> 16;
3887     assert(coef_sum < INT_MAX);
3888
3889     if(pict->pict_type == I_TYPE){
3890         s->m.current_picture.mb_var_sum= coef_sum;
3891         s->m.current_picture.mc_mb_var_sum= 0;
3892     }else{
3893         s->m.current_picture.mc_mb_var_sum= coef_sum;
3894         s->m.current_picture.mb_var_sum= 0;
3895     }
3896
3897     pict->quality= ff_rate_estimate_qscale(&s->m, 1);
3898     s->lambda= pict->quality * 3/2;
3899 }
3900
3901 static void calculate_vissual_weight(SnowContext *s, Plane *p){
3902     int width = p->width;
3903     int height= p->height;
3904     int level, orientation, x, y;
3905
3906     for(level=0; level<s->spatial_decomposition_count; level++){
3907         for(orientation=level ? 1 : 0; orientation<4; orientation++){
3908             SubBand *b= &p->band[level][orientation];
3909             DWTELEM *buf= b->buf;
3910             int64_t error=0;
3911
3912             memset(s->spatial_dwt_buffer, 0, sizeof(int)*width*height);
3913             buf[b->width/2 + b->height/2*b->stride]= 256*256;
3914             ff_spatial_idwt(s->spatial_dwt_buffer, width, height, width, s->spatial_decomposition_type, s->spatial_decomposition_count);
3915             for(y=0; y<height; y++){
3916                 for(x=0; x<width; x++){
3917                     int64_t d= s->spatial_dwt_buffer[x + y*width];
3918                     error += d*d;
3919                 }
3920             }
3921
3922             b->qlog= (int)(log(352256.0/sqrt(error)) / log(pow(2.0, 1.0/QROOT))+0.5);
3923 //            av_log(NULL, AV_LOG_DEBUG, "%d %d %d\n", level, orientation, b->qlog/*, sqrt(error)*/);
3924         }
3925     }
3926 }
3927
3928 static int encode_init(AVCodecContext *avctx)
3929 {
3930     SnowContext *s = avctx->priv_data;
3931     int plane_index;
3932
3933     if(avctx->strict_std_compliance > FF_COMPLIANCE_EXPERIMENTAL){
3934         av_log(avctx, AV_LOG_ERROR, "this codec is under development, files encoded with it may not be decodable with future versions!!!\n"
3935                "use vstrict=-2 / -strict -2 to use it anyway\n");
3936         return -1;
3937     }
3938
3939     common_init(avctx);
3940     alloc_blocks(s);
3941
3942     s->version=0;
3943
3944     s->m.avctx   = avctx;
3945     s->m.flags   = avctx->flags;
3946     s->m.bit_rate= avctx->bit_rate;
3947
3948     s->m.me.scratchpad= av_mallocz((avctx->width+64)*2*16*2*sizeof(uint8_t));
3949     s->m.me.map       = av_mallocz(ME_MAP_SIZE*sizeof(uint32_t));
3950     s->m.me.score_map = av_mallocz(ME_MAP_SIZE*sizeof(uint32_t));
3951     s->m.obmc_scratchpad= av_mallocz(MB_SIZE*MB_SIZE*12*sizeof(uint32_t));
3952     h263_encode_init(&s->m); //mv_penalty
3953
3954     if(avctx->flags&CODEC_FLAG_PASS1){
3955         if(!avctx->stats_out)
3956             avctx->stats_out = av_mallocz(256);
3957     }
3958     if(!(avctx->flags&CODEC_FLAG_QSCALE)){
3959         if(ff_rate_control_init(&s->m) < 0)
3960             return -1;
3961     }
3962     s->pass1_rc= !(avctx->flags & (CODEC_FLAG_QSCALE|CODEC_FLAG_PASS2));
3963
3964     for(plane_index=0; plane_index<3; plane_index++){
3965         calculate_vissual_weight(s, &s->plane[plane_index]);
3966     }
3967
3968
3969     avctx->coded_frame= &s->current_picture;
3970     switch(avctx->pix_fmt){
3971 //    case PIX_FMT_YUV444P:
3972 //    case PIX_FMT_YUV422P:
3973     case PIX_FMT_YUV420P:
3974     case PIX_FMT_GRAY8:
3975 //    case PIX_FMT_YUV411P:
3976 //    case PIX_FMT_YUV410P:
3977         s->colorspace_type= 0;
3978         break;
3979 /*    case PIX_FMT_RGBA32:
3980         s->colorspace= 1;
3981         break;*/
3982     default:
3983         av_log(avctx, AV_LOG_ERROR, "format not supported\n");
3984         return -1;
3985     }
3986 //    avcodec_get_chroma_sub_sample(avctx->pix_fmt, &s->chroma_h_shift, &s->chroma_v_shift);
3987     s->chroma_h_shift= 1;
3988     s->chroma_v_shift= 1;
3989
3990     ff_set_cmp(&s->dsp, s->dsp.me_cmp, s->avctx->me_cmp);
3991     ff_set_cmp(&s->dsp, s->dsp.me_sub_cmp, s->avctx->me_sub_cmp);
3992
3993     s->avctx->get_buffer(s->avctx, &s->input_picture);
3994
3995     return 0;
3996 }
3997
3998 static int frame_start(SnowContext *s){
3999    AVFrame tmp;
4000    int w= s->avctx->width; //FIXME round up to x16 ?
4001    int h= s->avctx->height;
4002
4003     if(s->current_picture.data[0]){
4004         draw_edges(s->current_picture.data[0], s->current_picture.linesize[0], w   , h   , EDGE_WIDTH  );
4005         draw_edges(s->current_picture.data[1], s->current_picture.linesize[1], w>>1, h>>1, EDGE_WIDTH/2);
4006         draw_edges(s->current_picture.data[2], s->current_picture.linesize[2], w>>1, h>>1, EDGE_WIDTH/2);
4007     }
4008
4009     tmp= s->last_picture;
4010     s->last_picture= s->current_picture;
4011     s->current_picture= tmp;
4012
4013     s->current_picture.reference= 1;
4014     if(s->avctx->get_buffer(s->avctx, &s->current_picture) < 0){
4015         av_log(s->avctx, AV_LOG_ERROR, "get_buffer() failed\n");
4016         return -1;
4017     }
4018
4019     return 0;
4020 }
4021
4022 static int encode_frame(AVCodecContext *avctx, unsigned char *buf, int buf_size, void *data){
4023     SnowContext *s = avctx->priv_data;
4024     RangeCoder * const c= &s->c;
4025     AVFrame *pict = data;
4026     const int width= s->avctx->width;
4027     const int height= s->avctx->height;
4028     int level, orientation, plane_index, i, y;
4029
4030     ff_init_range_encoder(c, buf, buf_size);
4031     ff_build_rac_states(c, 0.05*(1LL<<32), 256-8);
4032
4033     for(i=0; i<3; i++){
4034         int shift= !!i;
4035         for(y=0; y<(height>>shift); y++)
4036             memcpy(&s->input_picture.data[i][y * s->input_picture.linesize[i]],
4037                    &pict->data[i][y * pict->linesize[i]],
4038                    width>>shift);
4039     }
4040     s->new_picture = *pict;
4041
4042     s->m.picture_number= avctx->frame_number;
4043     if(avctx->flags&CODEC_FLAG_PASS2){
4044         s->m.pict_type =
4045         pict->pict_type= s->m.rc_context.entry[avctx->frame_number].new_pict_type;
4046         s->keyframe= pict->pict_type==FF_I_TYPE;
4047         if(!(avctx->flags&CODEC_FLAG_QSCALE))
4048             pict->quality= ff_rate_estimate_qscale(&s->m, 0);
4049     }else{
4050         s->keyframe= avctx->gop_size==0 || avctx->frame_number % avctx->gop_size == 0;
4051         s->m.pict_type=
4052         pict->pict_type= s->keyframe ? FF_I_TYPE : FF_P_TYPE;
4053     }
4054
4055     if(s->pass1_rc && avctx->frame_number == 0)
4056         pict->quality= 2*FF_QP2LAMBDA;
4057     if(pict->quality){
4058         s->qlog= rint(QROOT*log(pict->quality / (float)FF_QP2LAMBDA)/log(2));
4059         //<64 >60
4060         s->qlog += 61*QROOT/8;
4061         s->lambda = pict->quality * 3/2;
4062     }
4063     if(s->qlog < 0 || (!pict->quality && (avctx->flags & CODEC_FLAG_QSCALE))){
4064         s->qlog= LOSSLESS_QLOG;
4065         s->lambda = 0;
4066     }//else keep previous frame's qlog until after motion est
4067
4068     frame_start(s);
4069     s->current_picture.key_frame= s->keyframe;
4070
4071     s->m.current_picture_ptr= &s->m.current_picture;
4072     if(pict->pict_type == P_TYPE){
4073         int block_width = (width +15)>>4;
4074         int block_height= (height+15)>>4;
4075         int stride= s->current_picture.linesize[0];
4076
4077         assert(s->current_picture.data[0]);
4078         assert(s->last_picture.data[0]);
4079
4080         s->m.avctx= s->avctx;
4081         s->m.current_picture.data[0]= s->current_picture.data[0];
4082         s->m.   last_picture.data[0]= s->   last_picture.data[0];
4083         s->m.    new_picture.data[0]= s->  input_picture.data[0];
4084         s->m.   last_picture_ptr= &s->m.   last_picture;
4085         s->m.linesize=
4086         s->m.   last_picture.linesize[0]=
4087         s->m.    new_picture.linesize[0]=
4088         s->m.current_picture.linesize[0]= stride;
4089         s->m.uvlinesize= s->current_picture.linesize[1];
4090         s->m.width = width;
4091         s->m.height= height;
4092         s->m.mb_width = block_width;
4093         s->m.mb_height= block_height;
4094         s->m.mb_stride=   s->m.mb_width+1;
4095         s->m.b8_stride= 2*s->m.mb_width+1;
4096         s->m.f_code=1;
4097         s->m.pict_type= pict->pict_type;
4098         s->m.me_method= s->avctx->me_method;
4099         s->m.me.scene_change_score=0;
4100         s->m.flags= s->avctx->flags;
4101         s->m.quarter_sample= (s->avctx->flags & CODEC_FLAG_QPEL)!=0;
4102         s->m.out_format= FMT_H263;
4103         s->m.unrestricted_mv= 1;
4104
4105         s->m.lambda = s->lambda;
4106         s->m.qscale= (s->m.lambda*139 + FF_LAMBDA_SCALE*64) >> (FF_LAMBDA_SHIFT + 7);
4107         s->lambda2= s->m.lambda2= (s->m.lambda*s->m.lambda + FF_LAMBDA_SCALE/2) >> FF_LAMBDA_SHIFT;
4108
4109         s->m.dsp= s->dsp; //move
4110         ff_init_me(&s->m);
4111         s->dsp= s->m.dsp;
4112     }
4113
4114 redo_frame:
4115
4116     s->m.pict_type = pict->pict_type;
4117     s->qbias= pict->pict_type == P_TYPE ? 2 : 0;
4118
4119     encode_header(s);
4120     s->m.misc_bits = 8*(s->c.bytestream - s->c.bytestream_start);
4121     encode_blocks(s);
4122     s->m.mv_bits = 8*(s->c.bytestream - s->c.bytestream_start) - s->m.misc_bits;
4123
4124     for(plane_index=0; plane_index<3; plane_index++){
4125         Plane *p= &s->plane[plane_index];
4126         int w= p->width;
4127         int h= p->height;
4128         int x, y;
4129 //        int bits= put_bits_count(&s->c.pb);
4130
4131         //FIXME optimize
4132      if(pict->data[plane_index]) //FIXME gray hack
4133         for(y=0; y<h; y++){
4134             for(x=0; x<w; x++){
4135                 s->spatial_dwt_buffer[y*w + x]= pict->data[plane_index][y*pict->linesize[plane_index] + x]<<FRAC_BITS;
4136             }
4137         }
4138         predict_plane(s, s->spatial_dwt_buffer, plane_index, 0);
4139
4140         if(   plane_index==0
4141            && pict->pict_type == P_TYPE
4142            && !(avctx->flags&CODEC_FLAG_PASS2)
4143            && s->m.me.scene_change_score > s->avctx->scenechange_threshold){
4144             ff_init_range_encoder(c, buf, buf_size);
4145             ff_build_rac_states(c, 0.05*(1LL<<32), 256-8);
4146             pict->pict_type= FF_I_TYPE;
4147             s->keyframe=1;
4148             reset_contexts(s);
4149             goto redo_frame;
4150         }
4151
4152         if(s->qlog == LOSSLESS_QLOG){
4153             for(y=0; y<h; y++){
4154                 for(x=0; x<w; x++){
4155                     s->spatial_dwt_buffer[y*w + x]= (s->spatial_dwt_buffer[y*w + x] + (1<<(FRAC_BITS-1))-1)>>FRAC_BITS;
4156                 }
4157             }
4158         }
4159
4160         ff_spatial_dwt(s->spatial_dwt_buffer, w, h, w, s->spatial_decomposition_type, s->spatial_decomposition_count);
4161
4162         if(s->pass1_rc && plane_index==0)
4163             ratecontrol_1pass(s, pict);
4164
4165         for(level=0; level<s->spatial_decomposition_count; level++){
4166             for(orientation=level ? 1 : 0; orientation<4; orientation++){
4167                 SubBand *b= &p->band[level][orientation];
4168
4169                 quantize(s, b, b->buf, b->stride, s->qbias);
4170                 if(orientation==0)
4171                     decorrelate(s, b, b->buf, b->stride, pict->pict_type == P_TYPE, 0);
4172                 encode_subband(s, b, b->buf, b->parent ? b->parent->buf : NULL, b->stride, orientation);
4173                 assert(b->parent==NULL || b->parent->stride == b->stride*2);
4174                 if(orientation==0)
4175                     correlate(s, b, b->buf, b->stride, 1, 0);
4176             }
4177         }
4178 //        av_log(NULL, AV_LOG_DEBUG, "plane:%d bits:%d\n", plane_index, put_bits_count(&s->c.pb) - bits);
4179
4180         for(level=0; level<s->spatial_decomposition_count; level++){
4181             for(orientation=level ? 1 : 0; orientation<4; orientation++){
4182                 SubBand *b= &p->band[level][orientation];
4183
4184                 dequantize(s, b, b->buf, b->stride);
4185             }
4186         }
4187
4188         ff_spatial_idwt(s->spatial_dwt_buffer, w, h, w, s->spatial_decomposition_type, s->spatial_decomposition_count);
4189         if(s->qlog == LOSSLESS_QLOG){
4190             for(y=0; y<h; y++){
4191                 for(x=0; x<w; x++){
4192                     s->spatial_dwt_buffer[y*w + x]<<=FRAC_BITS;
4193                 }
4194             }
4195         }
4196 {START_TIMER
4197         predict_plane(s, s->spatial_dwt_buffer, plane_index, 1);
4198 STOP_TIMER("pred-conv")}
4199         if(s->avctx->flags&CODEC_FLAG_PSNR){
4200             int64_t error= 0;
4201
4202     if(pict->data[plane_index]) //FIXME gray hack
4203             for(y=0; y<h; y++){
4204                 for(x=0; x<w; x++){
4205                     int d= s->current_picture.data[plane_index][y*s->current_picture.linesize[plane_index] + x] - pict->data[plane_index][y*pict->linesize[plane_index] + x];
4206                     error += d*d;
4207                 }
4208             }
4209             s->avctx->error[plane_index] += error;
4210             s->current_picture.error[plane_index] = error;
4211         }
4212     }
4213
4214     if(s->last_picture.data[0])
4215         avctx->release_buffer(avctx, &s->last_picture);
4216
4217     s->current_picture.coded_picture_number = avctx->frame_number;
4218     s->current_picture.pict_type = pict->pict_type;
4219     s->current_picture.quality = pict->quality;
4220     s->m.frame_bits = 8*(s->c.bytestream - s->c.bytestream_start);
4221     s->m.p_tex_bits = s->m.frame_bits - s->m.misc_bits - s->m.mv_bits;
4222     s->m.current_picture.display_picture_number =
4223     s->m.current_picture.coded_picture_number = avctx->frame_number;
4224     s->m.current_picture.quality = pict->quality;
4225     s->m.total_bits += 8*(s->c.bytestream - s->c.bytestream_start);
4226     if(s->pass1_rc)
4227         ff_rate_estimate_qscale(&s->m, 0);
4228     if(avctx->flags&CODEC_FLAG_PASS1)
4229         ff_write_pass1_stats(&s->m);
4230     s->m.last_pict_type = s->m.pict_type;
4231
4232     emms_c();
4233
4234     return ff_rac_terminate(c);
4235 }
4236
4237 static void common_end(SnowContext *s){
4238     int plane_index, level, orientation;
4239
4240     av_freep(&s->spatial_dwt_buffer);
4241
4242     av_freep(&s->m.me.scratchpad);
4243     av_freep(&s->m.me.map);
4244     av_freep(&s->m.me.score_map);
4245     av_freep(&s->m.obmc_scratchpad);
4246
4247     av_freep(&s->block);
4248
4249     for(plane_index=0; plane_index<3; plane_index++){
4250         for(level=s->spatial_decomposition_count-1; level>=0; level--){
4251             for(orientation=level ? 1 : 0; orientation<4; orientation++){
4252                 SubBand *b= &s->plane[plane_index].band[level][orientation];
4253
4254                 av_freep(&b->x_coeff);
4255             }
4256         }
4257     }
4258 }
4259
4260 static int encode_end(AVCodecContext *avctx)
4261 {
4262     SnowContext *s = avctx->priv_data;
4263
4264     common_end(s);
4265     av_free(avctx->stats_out);
4266
4267     return 0;
4268 }
4269
4270 static int decode_init(AVCodecContext *avctx)
4271 {
4272     SnowContext *s = avctx->priv_data;
4273     int block_size;
4274
4275     avctx->pix_fmt= PIX_FMT_YUV420P;
4276