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