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