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