optimize quantizaton (about 3x faster)
[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 #include "cabac.h"
23
24 #include "mpegvideo.h"
25
26 #undef NDEBUG
27 #include <assert.h>
28
29 #define MAX_DECOMPOSITIONS 8
30 #define MAX_PLANES 4
31 #define DWTELEM int
32 #define QROOT 8 
33
34 static const int8_t quant3[256]={
35  0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
36  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
37  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
38  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
39  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
40  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
41  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
42  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
43 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
44 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
45 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
46 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
47 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
48 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-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, 0,
51 };
52 static const int8_t quant3b[256]={
53  0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
54  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
55  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
56  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
57  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
58  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
59  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
60  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
61 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
62 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
63 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
64 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
65 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
66 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
67 -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 };
70 static const int8_t quant5[256]={
71  0, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
72  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
73  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
74  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
75  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
76  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
77  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
78  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
79 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
80 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
81 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
82 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
83 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
84 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
85 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,
86 -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-1,-1,-1,
87 };
88 static const int8_t quant7[256]={
89  0, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
90  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
91  2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
92  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
93  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
94  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
95  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
96  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
97 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
98 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
99 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
100 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
101 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,
102 -3,-3,-3,-3,-3,-3,-3,-3,-3,-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,-1,-1,
105 };
106 static const int8_t quant9[256]={
107  0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3,
108  3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
109  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
110  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
111  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
112  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
113  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
114  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
115 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
116 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
117 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
118 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
119 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
120 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,
121 -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-3,-3,-3,-3,
122 -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-2,-2,-2,-2,-1,-1,
123 };
124 static const int8_t quant11[256]={
125  0, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4,
126  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
127  4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
128  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
129  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
130  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
131  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
132  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
133 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
134 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
135 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
136 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
137 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,
138 -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-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,-3,-3,-3,-3,-3,-3,-3,-2,-2,-2,-1,
141 };
142 static const int8_t quant13[256]={
143  0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4,
144  4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
145  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
146  5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
147  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
148  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
149  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
150  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
151 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,
152 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,
153 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,
154 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,
155 -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-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 -4,-4,-4,-4,-4,-4,-4,-4,-4,-3,-3,-3,-3,-2,-2,-1,
159 };
160
161 #define OBMC_MAX 64
162 #if 0 //64*cubic
163 static const uint8_t obmc32[1024]={
164  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,
165  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,
166  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,
167  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,
168  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,
169  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,
170  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,
171  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,
172  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,
173  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,
174  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,
175  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,
176  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,
177  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,
178  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,
179  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,
180  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,
181  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,
182  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,
183  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,
184  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,
185  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,
186  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,
187  0, 1, 2, 4, 7,10,13,16,19,22,25,28,31,33,34,35,35,34,33,31,28,25,22,19,16,13,10, 7, 4, 2, 1, 0,
188  0, 1, 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,
189  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,
190  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,
191  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,
192  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,
193  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,
194  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,
195  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,
196 //error:0.000022
197 };
198 static const uint8_t obmc16[256]={
199  0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
200  0, 1, 1, 2, 4, 5, 5, 6, 6, 5, 5, 4, 2, 1, 1, 0,
201  0, 1, 4, 6, 9,11,13,15,15,13,11, 9, 6, 4, 1, 0,
202  0, 2, 6,11,15,20,24,26,26,24,20,15,11, 6, 2, 0,
203  0, 4, 9,15,23,29,34,38,38,34,29,23,15, 9, 4, 0,
204  0, 5,11,20,29,38,45,49,49,45,38,29,20,11, 5, 0,
205  1, 5,13,24,34,45,53,57,57,53,45,34,24,13, 5, 1,
206  1, 6,15,26,38,49,57,62,62,57,49,38,26,15, 6, 1,
207  1, 6,15,26,38,49,57,62,62,57,49,38,26,15, 6, 1,
208  1, 5,13,24,34,45,53,57,57,53,45,34,24,13, 5, 1,
209  0, 5,11,20,29,38,45,49,49,45,38,29,20,11, 5, 0,
210  0, 4, 9,15,23,29,34,38,38,34,29,23,15, 9, 4, 0,
211  0, 2, 6,11,15,20,24,26,26,24,20,15,11, 6, 2, 0,
212  0, 1, 4, 6, 9,11,13,15,15,13,11, 9, 6, 4, 1, 0,
213  0, 1, 1, 2, 4, 5, 5, 6, 6, 5, 5, 4, 2, 1, 1, 0,
214  0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
215 //error:0.000033
216 };
217 #elif 1 // 64*linear
218 static const uint8_t obmc32[1024]={
219  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,
220  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,
221  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,
222  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,
223  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,
224  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,
225  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,
226  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,
227  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,
228  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,
229  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,
230  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,
231  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,
232  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,
233  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,
234  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,
235  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,
236  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,
237  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,
238  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,
239  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,
240  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,
241  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,
242  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,
243  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,
244  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,
245  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,
246  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,
247  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,
248  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,
249  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,
250  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,
251  //error:0.000020
252 };
253 static const uint8_t obmc16[256]={
254  0, 1, 1, 2, 2, 3, 3, 4, 4, 3, 3, 2, 2, 1, 1, 0,
255  1, 2, 4, 5, 7, 8,10,11,11,10, 8, 7, 5, 4, 2, 1,
256  1, 4, 6, 9,11,14,16,19,19,16,14,11, 9, 6, 4, 1,
257  2, 5, 9,12,16,19,23,26,26,23,19,16,12, 9, 5, 2,
258  2, 7,11,16,20,25,29,34,34,29,25,20,16,11, 7, 2,
259  3, 8,14,19,25,30,36,41,41,36,30,25,19,14, 8, 3,
260  3,10,16,23,29,36,42,49,49,42,36,29,23,16,10, 3,
261  4,11,19,26,34,41,49,56,56,49,41,34,26,19,11, 4,
262  4,11,19,26,34,41,49,56,56,49,41,34,26,19,11, 4,
263  3,10,16,23,29,36,42,49,49,42,36,29,23,16,10, 3,
264  3, 8,14,19,25,30,36,41,41,36,30,25,19,14, 8, 3,
265  2, 7,11,16,20,25,29,34,34,29,25,20,16,11, 7, 2,
266  2, 5, 9,12,16,19,23,26,26,23,19,16,12, 9, 5, 2,
267  1, 4, 6, 9,11,14,16,19,19,16,14,11, 9, 6, 4, 1,
268  1, 2, 4, 5, 7, 8,10,11,11,10, 8, 7, 5, 4, 2, 1,
269  0, 1, 1, 2, 2, 3, 3, 4, 4, 3, 3, 2, 2, 1, 1, 0,
270 //error:0.000015
271 };
272 #else //64*cos
273 static const uint8_t obmc32[1024]={
274  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,
275  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,
276  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,
277  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,
278  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,
279  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,
280  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,
281  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,
282  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,
283  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,
284  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,
285  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,
286  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,
287  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,
288  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,
289  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,
290  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,
291  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,
292  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,
293  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,
294  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,
295  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,
296  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,
297  0, 1, 2, 4, 6, 9,12,16,19,23,26,29,31,33,34,35,35,34,33,31,29,26,23,19,16,12, 9, 6, 4, 2, 1, 0,
298  0, 1, 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,
299  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,
300  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,
301  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,
302  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,
303  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,
304  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,
305  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,
306 //error:0.000022
307 };
308 static const uint8_t obmc16[256]={
309  0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0,
310  0, 0, 1, 2, 3, 4, 5, 5, 5, 5, 4, 3, 2, 1, 0, 0,
311  0, 1, 3, 6, 8,11,13,14,14,13,11, 8, 6, 3, 1, 0,
312  0, 2, 6,10,15,20,24,26,26,24,20,15,10, 6, 2, 0,
313  0, 3, 8,16,23,30,35,38,38,35,30,23,16, 8, 3, 0,
314  1, 4,11,20,30,39,46,49,49,46,39,30,20,11, 4, 1,
315  1, 5,13,24,35,46,54,58,58,54,46,35,24,13, 5, 1,
316  0, 5,14,26,38,49,58,63,63,58,49,38,26,14, 5, 0,
317  0, 5,14,26,38,49,58,63,63,58,49,38,26,14, 5, 0,
318  1, 5,13,24,35,46,54,58,58,54,46,35,24,13, 5, 1,
319  1, 4,11,20,30,39,46,49,49,46,39,30,20,11, 4, 1,
320  0, 3, 8,16,23,30,35,38,38,35,30,23,16, 8, 3, 0,
321  0, 2, 6,10,15,20,24,26,26,24,20,15,10, 6, 2, 0,
322  0, 1, 3, 6, 8,11,13,14,14,13,11, 8, 6, 3, 1, 0,
323  0, 0, 1, 2, 3, 4, 5, 5, 5, 5, 4, 3, 2, 1, 0, 0,
324  0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0,
325 //error:0.000022
326 };
327 #endif
328
329
330 typedef struct SubBand{
331     int level;
332     int stride;
333     int width;
334     int height;
335     int qlog;                                   ///< log(qscale)/log[2^(1/6)]
336     DWTELEM *buf;
337     struct SubBand *parent;
338     uint8_t state[/*7*2*/ 7 + 512][32];
339 }SubBand;
340
341 typedef struct Plane{
342     int width;
343     int height;
344     SubBand band[MAX_DECOMPOSITIONS][4];
345 }Plane;
346
347 typedef struct SnowContext{
348 //    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)
349
350     AVCodecContext *avctx;
351     CABACContext c;
352     DSPContext dsp;
353     AVFrame input_picture;
354     AVFrame current_picture;
355     AVFrame last_picture;
356     AVFrame mconly_picture;
357 //     uint8_t q_context[16];
358     uint8_t header_state[32];
359     int keyframe;
360     int version;
361     int spatial_decomposition_type;
362     int temporal_decomposition_type;
363     int spatial_decomposition_count;
364     int temporal_decomposition_count;
365     DWTELEM *spatial_dwt_buffer;
366     DWTELEM *pred_buffer;
367     int colorspace_type;
368     int chroma_h_shift;
369     int chroma_v_shift;
370     int spatial_scalability;
371     int qlog;
372     int mv_scale;
373     int qbias;
374 #define QBIAS_SHIFT 3
375     int b_width; //FIXME remove?
376     int b_height; //FIXME remove?
377     Plane plane[MAX_PLANES];
378     SubBand mb_band;
379     SubBand mv_band[2];
380     
381     uint16_t *mb_type;
382     uint8_t *mb_mean;
383     uint32_t *dummy;
384     int16_t (*motion_val8)[2];
385     int16_t (*motion_val16)[2];
386     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)
387 }SnowContext;
388
389 #define QEXPSHIFT 7 //FIXME try to change this to 0
390 static const uint8_t qexp[8]={
391     128, 140, 152, 166, 181, 197, 215, 235
392 //   64,  70,  76,  83,  91,  99, 108, 117
393 //   32,  35,  38,  41,  45,  49,  54,  59
394 //   16,  17,  19,  21,  23,  25,  27,  29
395 //    8,   9,  10,  10,  11,  12,  13,  15
396 };
397
398 static inline int mirror(int v, int m){
399     if     (v<0) return -v;
400     else if(v>m) return 2*m-v;
401     else         return v;
402 }
403
404 static inline void put_symbol(CABACContext *c, uint8_t *state, int v, int is_signed){
405     int i;
406
407     if(v){
408         const int a= ABS(v);
409         const int e= av_log2(a);
410 #if 1
411         const int el= FFMIN(e, 10);   
412         put_cabac(c, state+0, 0);
413
414         for(i=0; i<el; i++){
415             put_cabac(c, state+1+i, 1);  //1..10
416         }
417         for(; i<e; i++){
418             put_cabac(c, state+1+9, 1);  //1..10
419         }
420         put_cabac(c, state+1+FFMIN(i,9), 0);
421
422         for(i=e-1; i>=el; i--){
423             put_cabac(c, state+22+9, (a>>i)&1); //22..31
424         }
425         for(; i>=0; i--){
426             put_cabac(c, state+22+i, (a>>i)&1); //22..31
427         }
428
429         if(is_signed)
430             put_cabac(c, state+11 + el, v < 0); //11..21
431 #else
432         
433         put_cabac(c, state+0, 0);
434         if(e<=9){
435             for(i=0; i<e; i++){
436                 put_cabac(c, state+1+i, 1);  //1..10
437             }
438             put_cabac(c, state+1+i, 0);
439
440             for(i=e-1; i>=0; i--){
441                 put_cabac(c, state+22+i, (a>>i)&1); //22..31
442             }
443
444             if(is_signed)
445                 put_cabac(c, state+11 + e, v < 0); //11..21
446         }else{
447             for(i=0; i<e; i++){
448                 put_cabac(c, state+1+FFMIN(i,9), 1);  //1..10
449             }
450             put_cabac(c, state+1+FFMIN(i,9), 0);
451
452             for(i=e-1; i>=0; i--){
453                 put_cabac(c, state+22+FFMIN(i,9), (a>>i)&1); //22..31
454             }
455
456             if(is_signed)
457                 put_cabac(c, state+11 + FFMIN(e,10), v < 0); //11..21
458         }
459 #endif
460     }else{
461         put_cabac(c, state+0, 1);
462     }
463 }
464
465 static inline int get_symbol(CABACContext *c, uint8_t *state, int is_signed){
466     if(get_cabac(c, state+0))
467         return 0;
468     else{
469         int i, e, a, el;
470  //FIXME try to merge loops with FFMIN() maybe they are equally fast and they are surly cuter
471         for(e=0; e<10; e++){ 
472             if(get_cabac(c, state + 1 + e)==0) // 1..10
473                 break;
474         }
475         el= e;
476  
477         if(e==10){
478             while(get_cabac(c, state + 1 + 9)) //10
479                 e++;
480         }
481         a= 1;
482         for(i=e-1; i>=el; i--){
483             a += a + get_cabac(c, state+22+9); //31
484         }
485         for(; i>=0; i--){
486             a += a + get_cabac(c, state+22+i); //22..31
487         }
488
489         if(is_signed && get_cabac(c, state+11 + el)) //11..21
490             return -a;
491         else
492             return a;
493     }
494 }
495
496 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){
497     const int mirror_left= !highpass;
498     const int mirror_right= (width&1) ^ highpass;
499     const int w= (width>>1) - 1 + (highpass & width);
500     int i;
501
502 #define LIFT(src, ref, inv) ((src) + ((inv) ? - (ref) : + (ref)))
503     if(mirror_left){
504         dst[0] = LIFT(src[0], ((mul*2*ref[0]+add)>>shift), inverse);
505         dst += dst_step;
506         src += src_step;
507     }
508     
509     for(i=0; i<w; i++){
510         dst[i*dst_step] = LIFT(src[i*src_step], ((mul*(ref[i*ref_step] + ref[(i+1)*ref_step])+add)>>shift), inverse);
511     }
512     
513     if(mirror_right){
514         dst[w*dst_step] = LIFT(src[w*src_step], ((mul*2*ref[w*ref_step]+add)>>shift), inverse);
515     }
516 }
517
518 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){
519     const int mirror_left= !highpass;
520     const int mirror_right= (width&1) ^ highpass;
521     const int w= (width>>1) - 1 + (highpass & width);
522     int i;
523
524     if(mirror_left){
525         int r= 3*2*ref[0];
526         r += r>>4;
527         r += r>>8;
528         dst[0] = LIFT(src[0], ((r+add)>>shift), inverse);
529         dst += dst_step;
530         src += src_step;
531     }
532     
533     for(i=0; i<w; i++){
534         int r= 3*(ref[i*ref_step] + ref[(i+1)*ref_step]);
535         r += r>>4;
536         r += r>>8;
537         dst[i*dst_step] = LIFT(src[i*src_step], ((r+add)>>shift), inverse);
538     }
539     
540     if(mirror_right){
541         int r= 3*2*ref[w*ref_step];
542         r += r>>4;
543         r += r>>8;
544         dst[w*dst_step] = LIFT(src[w*src_step], ((r+add)>>shift), inverse);
545     }
546 }
547
548
549 static void inplace_lift(int *dst, int width, int *coeffs, int n, int shift, int start, int inverse){
550     int x, i;
551     
552     for(x=start; x<width; x+=2){
553         int64_t sum=0;
554
555         for(i=0; i<n; i++){
556             int x2= x + 2*i - n + 1;
557             if     (x2<     0) x2= -x2;
558             else if(x2>=width) x2= 2*width-x2-2;
559             sum += coeffs[i]*(int64_t)dst[x2];
560         }
561         if(inverse) dst[x] -= (sum + (1<<shift)/2)>>shift;
562         else        dst[x] += (sum + (1<<shift)/2)>>shift;
563     }
564 }
565
566 static void inplace_liftV(int *dst, int width, int height, int stride, int *coeffs, int n, int shift, int start, int inverse){
567     int x, y, i;
568     for(y=start; y<height; y+=2){
569         for(x=0; x<width; x++){
570             int64_t sum=0;
571     
572             for(i=0; i<n; i++){
573                 int y2= y + 2*i - n + 1;
574                 if     (y2<      0) y2= -y2;
575                 else if(y2>=height) y2= 2*height-y2-2;
576                 sum += coeffs[i]*(int64_t)dst[x + y2*stride];
577             }
578             if(inverse) dst[x + y*stride] -= (sum + (1<<shift)/2)>>shift;
579             else        dst[x + y*stride] += (sum + (1<<shift)/2)>>shift;
580         }
581     }
582 }
583
584 #define SCALEX 1
585 #define LX0 0
586 #define LX1 1
587
588 #if 0 // more accurate 9/7
589 #define N1 2
590 #define SHIFT1 14
591 #define COEFFS1 (int[]){-25987,-25987}
592 #define N2 2
593 #define SHIFT2 19
594 #define COEFFS2 (int[]){-27777,-27777}
595 #define N3 2
596 #define SHIFT3 15
597 #define COEFFS3 (int[]){28931,28931}
598 #define N4 2
599 #define SHIFT4 15
600 #define COEFFS4 (int[]){14533,14533}
601 #elif 1 // 13/7 CRF
602 #define N1 4
603 #define SHIFT1 4
604 #define COEFFS1 (int[]){1,-9,-9,1}
605 #define N2 4
606 #define SHIFT2 4
607 #define COEFFS2 (int[]){-1,5,5,-1}
608 #define N3 0
609 #define SHIFT3 1
610 #define COEFFS3 NULL
611 #define N4 0
612 #define SHIFT4 1
613 #define COEFFS4 NULL
614 #elif 1 // 3/5
615 #define LX0 1
616 #define LX1 0
617 #define SCALEX 0.5
618 #define N1 2
619 #define SHIFT1 1
620 #define COEFFS1 (int[]){1,1}
621 #define N2 2
622 #define SHIFT2 2
623 #define COEFFS2 (int[]){-1,-1}
624 #define N3 0
625 #define SHIFT3 0
626 #define COEFFS3 NULL
627 #define N4 0
628 #define SHIFT4 0
629 #define COEFFS4 NULL
630 #elif 1 // 11/5 
631 #define N1 0
632 #define SHIFT1 1
633 #define COEFFS1 NULL
634 #define N2 2
635 #define SHIFT2 2
636 #define COEFFS2 (int[]){-1,-1}
637 #define N3 2
638 #define SHIFT3 0
639 #define COEFFS3 (int[]){-1,-1}
640 #define N4 4
641 #define SHIFT4 7
642 #define COEFFS4 (int[]){-5,29,29,-5}
643 #define SCALEX 4
644 #elif 1 // 9/7 CDF
645 #define N1 2
646 #define SHIFT1 7
647 #define COEFFS1 (int[]){-203,-203}
648 #define N2 2
649 #define SHIFT2 12
650 #define COEFFS2 (int[]){-217,-217}
651 #define N3 2
652 #define SHIFT3 7
653 #define COEFFS3 (int[]){113,113}
654 #define N4 2
655 #define SHIFT4 9
656 #define COEFFS4 (int[]){227,227}
657 #define SCALEX 1
658 #elif 1 // 7/5 CDF
659 #define N1 0
660 #define SHIFT1 1
661 #define COEFFS1 NULL
662 #define N2 2
663 #define SHIFT2 2
664 #define COEFFS2 (int[]){-1,-1}
665 #define N3 2
666 #define SHIFT3 0
667 #define COEFFS3 (int[]){-1,-1}
668 #define N4 2
669 #define SHIFT4 4
670 #define COEFFS4 (int[]){3,3}
671 #elif 1 // 9/7 MN
672 #define N1 4
673 #define SHIFT1 4
674 #define COEFFS1 (int[]){1,-9,-9,1}
675 #define N2 2
676 #define SHIFT2 2
677 #define COEFFS2 (int[]){1,1}
678 #define N3 0
679 #define SHIFT3 1
680 #define COEFFS3 NULL
681 #define N4 0
682 #define SHIFT4 1
683 #define COEFFS4 NULL
684 #else // 13/7 CRF
685 #define N1 4
686 #define SHIFT1 4
687 #define COEFFS1 (int[]){1,-9,-9,1}
688 #define N2 4
689 #define SHIFT2 4
690 #define COEFFS2 (int[]){-1,5,5,-1}
691 #define N3 0
692 #define SHIFT3 1
693 #define COEFFS3 NULL
694 #define N4 0
695 #define SHIFT4 1
696 #define COEFFS4 NULL
697 #endif
698 static void horizontal_decomposeX(int *b, int width){
699     int temp[width];
700     const int width2= width>>1;
701     const int w2= (width+1)>>1;
702     int A1,A2,A3,A4, x;
703
704     inplace_lift(b, width, COEFFS1, N1, SHIFT1, LX1, 0);
705     inplace_lift(b, width, COEFFS2, N2, SHIFT2, LX0, 0);
706     inplace_lift(b, width, COEFFS3, N3, SHIFT3, LX1, 0);
707     inplace_lift(b, width, COEFFS4, N4, SHIFT4, LX0, 0);
708     
709     for(x=0; x<width2; x++){
710         temp[x   ]= b[2*x    ];
711         temp[x+w2]= b[2*x + 1];
712     }
713     if(width&1)
714         temp[x   ]= b[2*x    ];
715     memcpy(b, temp, width*sizeof(int));
716 }
717
718 static void horizontal_composeX(int *b, int width){
719     int temp[width];
720     const int width2= width>>1;
721     int A1,A2,A3,A4, x;
722     const int w2= (width+1)>>1;
723
724     memcpy(temp, b, width*sizeof(int));
725     for(x=0; x<width2; x++){
726         b[2*x    ]= temp[x   ];
727         b[2*x + 1]= temp[x+w2];
728     }
729     if(width&1)
730         b[2*x    ]= temp[x   ];
731
732     inplace_lift(b, width, COEFFS4, N4, SHIFT4, LX0, 1);
733     inplace_lift(b, width, COEFFS3, N3, SHIFT3, LX1, 1);
734     inplace_lift(b, width, COEFFS2, N2, SHIFT2, LX0, 1);
735     inplace_lift(b, width, COEFFS1, N1, SHIFT1, LX1, 1);
736 }
737
738 static void spatial_decomposeX(int *buffer, int width, int height, int stride){
739     int x, y;
740   
741     for(y=0; y<height; y++){
742         for(x=0; x<width; x++){
743             buffer[y*stride + x] *= SCALEX;
744         }
745     }
746
747     for(y=0; y<height; y++){
748         horizontal_decomposeX(buffer + y*stride, width);
749     }
750     
751     inplace_liftV(buffer, width, height, stride, COEFFS1, N1, SHIFT1, LX1, 0);
752     inplace_liftV(buffer, width, height, stride, COEFFS2, N2, SHIFT2, LX0, 0);
753     inplace_liftV(buffer, width, height, stride, COEFFS3, N3, SHIFT3, LX1, 0);
754     inplace_liftV(buffer, width, height, stride, COEFFS4, N4, SHIFT4, LX0, 0);    
755 }
756
757 static void spatial_composeX(int *buffer, int width, int height, int stride){
758     int x, y;
759   
760     inplace_liftV(buffer, width, height, stride, COEFFS4, N4, SHIFT4, LX0, 1);
761     inplace_liftV(buffer, width, height, stride, COEFFS3, N3, SHIFT3, LX1, 1);
762     inplace_liftV(buffer, width, height, stride, COEFFS2, N2, SHIFT2, LX0, 1);
763     inplace_liftV(buffer, width, height, stride, COEFFS1, N1, SHIFT1, LX1, 1);
764
765     for(y=0; y<height; y++){
766         horizontal_composeX(buffer + y*stride, width);
767     }
768
769     for(y=0; y<height; y++){
770         for(x=0; x<width; x++){
771             buffer[y*stride + x] /= SCALEX;
772         }
773     }
774 }
775
776 static void horizontal_decompose53i(int *b, int width){
777     int temp[width];
778     const int width2= width>>1;
779     int A1,A2,A3,A4, x;
780     const int w2= (width+1)>>1;
781
782     for(x=0; x<width2; x++){
783         temp[x   ]= b[2*x    ];
784         temp[x+w2]= b[2*x + 1];
785     }
786     if(width&1)
787         temp[x   ]= b[2*x    ];
788 #if 0
789     A2= temp[1       ];
790     A4= temp[0       ];
791     A1= temp[0+width2];
792     A1 -= (A2 + A4)>>1;
793     A4 += (A1 + 1)>>1;
794     b[0+width2] = A1;
795     b[0       ] = A4;
796     for(x=1; x+1<width2; x+=2){
797         A3= temp[x+width2];
798         A4= temp[x+1     ];
799         A3 -= (A2 + A4)>>1;
800         A2 += (A1 + A3 + 2)>>2;
801         b[x+width2] = A3;
802         b[x       ] = A2;
803
804         A1= temp[x+1+width2];
805         A2= temp[x+2       ];
806         A1 -= (A2 + A4)>>1;
807         A4 += (A1 + A3 + 2)>>2;
808         b[x+1+width2] = A1;
809         b[x+1       ] = A4;
810     }
811     A3= temp[width-1];
812     A3 -= A2;
813     A2 += (A1 + A3 + 2)>>2;
814     b[width -1] = A3;
815     b[width2-1] = A2;
816 #else        
817     lift(b+w2, temp+w2, temp, 1, 1, 1, width, -1, 0, 1, 1, 0);
818     lift(b   , temp   , b+w2, 1, 1, 1, width,  1, 2, 2, 0, 0);
819 #endif
820 }
821
822 static void vertical_decompose53iH0(int *b0, int *b1, int *b2, int width){
823     int i;
824     
825     for(i=0; i<width; i++){
826         b1[i] -= (b0[i] + b2[i])>>1;
827     }
828 }
829
830 static void vertical_decompose53iL0(int *b0, int *b1, int *b2, int width){
831     int i;
832     
833     for(i=0; i<width; i++){
834         b1[i] += (b0[i] + b2[i] + 2)>>2;
835     }
836 }
837
838 static void spatial_decompose53i(int *buffer, int width, int height, int stride){
839     int x, y;
840     DWTELEM *b0= buffer + mirror(-2-1, height-1)*stride;
841     DWTELEM *b1= buffer + mirror(-2  , height-1)*stride;
842   
843     for(y=-2; y<height; y+=2){
844         DWTELEM *b2= buffer + mirror(y+1, height-1)*stride;
845         DWTELEM *b3= buffer + mirror(y+2, height-1)*stride;
846
847 {START_TIMER
848         if(b1 <= b3)     horizontal_decompose53i(b2, width);
849         if(y+2 < height) horizontal_decompose53i(b3, width);
850 STOP_TIMER("horizontal_decompose53i")}
851         
852 {START_TIMER
853         if(b1 <= b3) vertical_decompose53iH0(b1, b2, b3, width);
854         if(b0 <= b2) vertical_decompose53iL0(b0, b1, b2, width);
855 STOP_TIMER("vertical_decompose53i*")}
856         
857         b0=b2;
858         b1=b3;
859     }
860 }
861
862 #define lift5 lift
863 #if 1
864 #define W_AM 3
865 #define W_AO 0
866 #define W_AS 1
867
868 #define W_BM 1
869 #define W_BO 8
870 #define W_BS 4
871
872 #undef lift5
873 #define W_CM 9999
874 #define W_CO 2
875 #define W_CS 2
876
877 #define W_DM 15
878 #define W_DO 16
879 #define W_DS 5
880 #elif 0
881 #define W_AM 55
882 #define W_AO 16
883 #define W_AS 5
884
885 #define W_BM 3
886 #define W_BO 32
887 #define W_BS 6
888
889 #define W_CM 127
890 #define W_CO 64
891 #define W_CS 7
892
893 #define W_DM 7
894 #define W_DO 8
895 #define W_DS 4
896 #elif 0
897 #define W_AM 97
898 #define W_AO 32
899 #define W_AS 6
900
901 #define W_BM 63
902 #define W_BO 512
903 #define W_BS 10
904
905 #define W_CM 13
906 #define W_CO 8
907 #define W_CS 4
908
909 #define W_DM 15
910 #define W_DO 16
911 #define W_DS 5
912
913 #else
914
915 #define W_AM 203
916 #define W_AO 64
917 #define W_AS 7
918
919 #define W_BM 217
920 #define W_BO 2048
921 #define W_BS 12
922
923 #define W_CM 113
924 #define W_CO 64
925 #define W_CS 7
926
927 #define W_DM 227
928 #define W_DO 128
929 #define W_DS 9
930 #endif
931 static void horizontal_decompose97i(int *b, int width){
932     int temp[width];
933     const int w2= (width+1)>>1;
934
935     lift (temp+w2, b    +1, b      , 1, 2, 2, width, -W_AM, W_AO, W_AS, 1, 0);
936     lift (temp   , b      , temp+w2, 1, 2, 1, width, -W_BM, W_BO, W_BS, 0, 0);
937     lift5(b   +w2, temp+w2, temp   , 1, 1, 1, width,  W_CM, W_CO, W_CS, 1, 0);
938     lift (b      , temp   , b   +w2, 1, 1, 1, width,  W_DM, W_DO, W_DS, 0, 0);
939 }
940
941
942 static void vertical_decompose97iH0(int *b0, int *b1, int *b2, int width){
943     int i;
944     
945     for(i=0; i<width; i++){
946         b1[i] -= (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
947     }
948 }
949
950 static void vertical_decompose97iH1(int *b0, int *b1, int *b2, int width){
951     int i;
952     
953     for(i=0; i<width; i++){
954 #ifdef lift5
955         b1[i] += (W_CM*(b0[i] + b2[i])+W_CO)>>W_CS;
956 #else
957         int r= 3*(b0[i] + b2[i]);
958         r+= r>>4;
959         r+= r>>8;
960         b1[i] += (r+W_CO)>>W_CS;
961 #endif
962     }
963 }
964
965 static void vertical_decompose97iL0(int *b0, int *b1, int *b2, int width){
966     int i;
967     
968     for(i=0; i<width; i++){
969         b1[i] -= (W_BM*(b0[i] + b2[i])+W_BO)>>W_BS;
970     }
971 }
972
973 static void vertical_decompose97iL1(int *b0, int *b1, int *b2, int width){
974     int i;
975     
976     for(i=0; i<width; i++){
977         b1[i] += (W_DM*(b0[i] + b2[i])+W_DO)>>W_DS;
978     }
979 }
980
981 static void spatial_decompose97i(int *buffer, int width, int height, int stride){
982     int x, y;
983     DWTELEM *b0= buffer + mirror(-4-1, height-1)*stride;
984     DWTELEM *b1= buffer + mirror(-4  , height-1)*stride;
985     DWTELEM *b2= buffer + mirror(-4+1, height-1)*stride;
986     DWTELEM *b3= buffer + mirror(-4+2, height-1)*stride;
987   
988     for(y=-4; y<height; y+=2){
989         DWTELEM *b4= buffer + mirror(y+3, height-1)*stride;
990         DWTELEM *b5= buffer + mirror(y+4, height-1)*stride;
991
992 {START_TIMER
993         if(b3 <= b5)     horizontal_decompose97i(b4, width);
994         if(y+4 < height) horizontal_decompose97i(b5, width);
995 if(width>400){
996 STOP_TIMER("horizontal_decompose97i")
997 }}
998         
999 {START_TIMER
1000         if(b3 <= b5) vertical_decompose97iH0(b3, b4, b5, width);
1001         if(b2 <= b4) vertical_decompose97iL0(b2, b3, b4, width);
1002         if(b1 <= b3) vertical_decompose97iH1(b1, b2, b3, width);
1003         if(b0 <= b2) vertical_decompose97iL1(b0, b1, b2, width);
1004
1005 if(width>400){
1006 STOP_TIMER("vertical_decompose97i")
1007 }}
1008         
1009         b0=b2;
1010         b1=b3;
1011         b2=b4;
1012         b3=b5;
1013     }
1014 }
1015
1016 static void spatial_dwt(SnowContext *s, int *buffer, int width, int height, int stride){
1017     int level;
1018     
1019     for(level=0; level<s->spatial_decomposition_count; level++){
1020         switch(s->spatial_decomposition_type){
1021         case 0: spatial_decompose97i(buffer, width>>level, height>>level, stride<<level); break;
1022         case 1: spatial_decompose53i(buffer, width>>level, height>>level, stride<<level); break;
1023         case 2: spatial_decomposeX  (buffer, width>>level, height>>level, stride<<level); break;
1024         }
1025     }
1026 }
1027
1028 static void horizontal_compose53i(int *b, int width){
1029     int temp[width];
1030     const int width2= width>>1;
1031     const int w2= (width+1)>>1;
1032     int A1,A2,A3,A4, x;
1033
1034 #if 0
1035     A2= temp[1       ];
1036     A4= temp[0       ];
1037     A1= temp[0+width2];
1038     A1 -= (A2 + A4)>>1;
1039     A4 += (A1 + 1)>>1;
1040     b[0+width2] = A1;
1041     b[0       ] = A4;
1042     for(x=1; x+1<width2; x+=2){
1043         A3= temp[x+width2];
1044         A4= temp[x+1     ];
1045         A3 -= (A2 + A4)>>1;
1046         A2 += (A1 + A3 + 2)>>2;
1047         b[x+width2] = A3;
1048         b[x       ] = A2;
1049
1050         A1= temp[x+1+width2];
1051         A2= temp[x+2       ];
1052         A1 -= (A2 + A4)>>1;
1053         A4 += (A1 + A3 + 2)>>2;
1054         b[x+1+width2] = A1;
1055         b[x+1       ] = A4;
1056     }
1057     A3= temp[width-1];
1058     A3 -= A2;
1059     A2 += (A1 + A3 + 2)>>2;
1060     b[width -1] = A3;
1061     b[width2-1] = A2;
1062 #else   
1063     lift(temp   , b   , b+w2, 1, 1, 1, width,  1, 2, 2, 0, 1);
1064     lift(temp+w2, b+w2, temp, 1, 1, 1, width, -1, 0, 1, 1, 1);
1065 #endif
1066     for(x=0; x<width2; x++){
1067         b[2*x    ]= temp[x   ];
1068         b[2*x + 1]= temp[x+w2];
1069     }
1070     if(width&1)
1071         b[2*x    ]= temp[x   ];
1072 }
1073
1074 static void vertical_compose53iH0(int *b0, int *b1, int *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_compose53iL0(int *b0, int *b1, int *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_compose53i(int *buffer, int width, int height, int stride){
1091     int x, y;
1092     DWTELEM *b0= buffer + mirror(-1-1, height-1)*stride;
1093     DWTELEM *b1= buffer + mirror(-1  , height-1)*stride;
1094   
1095     for(y=-1; 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(b1 <= b3) vertical_compose53iL0(b1, b2, b3, width);
1101         if(b0 <= b2) vertical_compose53iH0(b0, b1, b2, width);
1102 STOP_TIMER("vertical_compose53i*")}
1103
1104 {START_TIMER
1105         if(y-1 >= 0) horizontal_compose53i(b0, width);
1106         if(b0 <= b2) horizontal_compose53i(b1, width);
1107 STOP_TIMER("horizontal_compose53i")}
1108
1109         b0=b2;
1110         b1=b3;
1111     }
1112 }   
1113
1114  
1115 static void horizontal_compose97i(int *b, int width){
1116     int temp[width];
1117     const int w2= (width+1)>>1;
1118
1119     lift (temp   , b      , b   +w2, 1, 1, 1, width,  W_DM, W_DO, W_DS, 0, 1);
1120     lift5(temp+w2, b   +w2, temp   , 1, 1, 1, width,  W_CM, W_CO, W_CS, 1, 1);
1121     lift (b      , temp   , temp+w2, 2, 1, 1, width, -W_BM, W_BO, W_BS, 0, 1);
1122     lift (b+1    , temp+w2, b      , 2, 1, 2, width, -W_AM, W_AO, W_AS, 1, 1);
1123 }
1124
1125 static void vertical_compose97iH0(int *b0, int *b1, int *b2, int width){
1126     int i;
1127     
1128     for(i=0; i<width; i++){
1129         b1[i] += (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
1130     }
1131 }
1132
1133 static void vertical_compose97iH1(int *b0, int *b1, int *b2, int width){
1134     int i;
1135     
1136     for(i=0; i<width; i++){
1137 #ifdef lift5
1138         b1[i] -= (W_CM*(b0[i] + b2[i])+W_CO)>>W_CS;
1139 #else
1140         int r= 3*(b0[i] + b2[i]);
1141         r+= r>>4;
1142         r+= r>>8;
1143         b1[i] -= (r+W_CO)>>W_CS;
1144 #endif
1145     }
1146 }
1147
1148 static void vertical_compose97iL0(int *b0, int *b1, int *b2, int width){
1149     int i;
1150     
1151     for(i=0; i<width; i++){
1152         b1[i] += (W_BM*(b0[i] + b2[i])+W_BO)>>W_BS;
1153     }
1154 }
1155
1156 static void vertical_compose97iL1(int *b0, int *b1, int *b2, int width){
1157     int i;
1158     
1159     for(i=0; i<width; i++){
1160         b1[i] -= (W_DM*(b0[i] + b2[i])+W_DO)>>W_DS;
1161     }
1162 }
1163
1164 static void spatial_compose97i(int *buffer, int width, int height, int stride){
1165     int x, y;
1166     DWTELEM *b0= buffer + mirror(-3-1, height-1)*stride;
1167     DWTELEM *b1= buffer + mirror(-3  , height-1)*stride;
1168     DWTELEM *b2= buffer + mirror(-3+1, height-1)*stride;
1169     DWTELEM *b3= buffer + mirror(-3+2, height-1)*stride;
1170
1171     for(y=-3; y<=height; y+=2){
1172         DWTELEM *b4= buffer + mirror(y+3, height-1)*stride;
1173         DWTELEM *b5= buffer + mirror(y+4, height-1)*stride;
1174
1175         if(stride == width && y+4 < height && 0){ 
1176             int x;
1177             for(x=0; x<width/2; x++)
1178                 b5[x] += 64*2;
1179             for(; x<width; x++)
1180                 b5[x] += 169*2;
1181         }
1182         
1183 {START_TIMER
1184         if(b3 <= b5) vertical_compose97iL1(b3, b4, b5, width);
1185         if(b2 <= b4) vertical_compose97iH1(b2, b3, b4, width);
1186         if(b1 <= b3) vertical_compose97iL0(b1, b2, b3, width);
1187         if(b0 <= b2) vertical_compose97iH0(b0, b1, b2, width);
1188 if(width>400){
1189 STOP_TIMER("vertical_compose97i")}}
1190
1191 {START_TIMER
1192         if(y-1>=  0) horizontal_compose97i(b0, width);
1193         if(b0 <= b2) horizontal_compose97i(b1, width);
1194 if(width>400 && b0 <= b2){
1195 STOP_TIMER("horizontal_compose97i")}}
1196         
1197         b0=b2;
1198         b1=b3;
1199         b2=b4;
1200         b3=b5;
1201     }
1202 }
1203
1204 static void spatial_idwt(SnowContext *s, int *buffer, int width, int height, int stride){
1205     int level;
1206
1207     for(level=s->spatial_decomposition_count-1; level>=0; level--){
1208         switch(s->spatial_decomposition_type){
1209         case 0: spatial_compose97i(buffer, width>>level, height>>level, stride<<level); break;
1210         case 1: spatial_compose53i(buffer, width>>level, height>>level, stride<<level); break;
1211         case 2: spatial_composeX  (buffer, width>>level, height>>level, stride<<level); break;
1212         }
1213     }
1214 }
1215
1216 static const int hilbert[16][2]={
1217     {0,0}, {1,0}, {1,1}, {0,1},
1218     {0,2}, {0,3}, {1,3}, {1,2},
1219     {2,2}, {2,3}, {3,3}, {3,2},
1220     {3,1}, {2,1}, {2,0}, {3,0},
1221 };
1222 #if 0
1223 -o o-
1224  | |
1225  o-o
1226  
1227 -o-o o-o-
1228    | | 
1229  o-o o-o
1230  |     |
1231  o o-o o
1232  | | | |
1233  o-o o-o
1234  
1235  0112122312232334122323342334
1236  0123456789ABCDEF0123456789AB
1237  RLLRMRRLLRRMRLLMLRRLMLLRRLLM
1238  
1239  4  B  F 14 1B
1240  4 11 15 20 27
1241  
1242 -o o-o-o o-o-o o-
1243  | |   | |   | |
1244  o-o o-o o-o o-o
1245      |     |
1246  o-o o-o o-o o-o
1247  | |   | |   | |
1248  o o-o-o o-o-o o
1249  |             |
1250  o-o o-o-o-o o-o
1251    | |     | | 
1252  o-o o-o o-o o-o
1253  |     | |     |
1254  o o-o o o o-o o
1255  | | | | | | | |
1256  o-o o-o o-o o-o
1257
1258 #endif
1259
1260 #define SVI(a, i, x, y) \
1261 {\
1262     a[i][0]= x;\
1263     a[i][1]= y;\
1264     i++;\
1265 }
1266
1267 static int sig_cmp(const void *a, const void *b){
1268     const int16_t* da = (const int16_t *) a;
1269     const int16_t* db = (const int16_t *) b;
1270     
1271     if(da[1] != db[1]) return da[1] - db[1];
1272     else               return da[0] - db[0];
1273 }
1274
1275
1276 static void encode_subband(SnowContext *s, SubBand *b, DWTELEM *src, DWTELEM *parent, int stride, int orientation){
1277     const int level= b->level;
1278     const int w= b->width;
1279     const int h= b->height;
1280     int x, y;
1281     
1282 #if 0
1283     if(orientation==3 && parent && 0){
1284         int16_t candidate[w*h][2];
1285         uint8_t state[w*h];
1286         int16_t boarder[3][w*h*4][2];
1287         int16_t significant[w*h][2];
1288         int candidate_count=0;
1289         int boarder_count[3]={0,0,0};
1290         int significant_count=0;
1291         int rle_pos=0;
1292         int v, last_v;
1293         int primary= orientation==1;
1294         
1295         memset(candidate, 0, sizeof(candidate));
1296         memset(state, 0, sizeof(state));
1297         memset(boarder, 0, sizeof(boarder));
1298         
1299         for(y=0; y<h; y++){
1300             for(x=0; x<w; x++){
1301                 if(parent[(x>>1) + (y>>1)*2*stride])
1302                     SVI(candidate, candidate_count, x, y)
1303             }
1304         }
1305
1306         for(;;){
1307             while(candidate_count && !boarder_count[0] && !boarder_count[1] && !boarder_count[2]){
1308                 candidate_count--;
1309                 x= candidate[ candidate_count][0];
1310                 y= candidate[ candidate_count][1];
1311                 if(state[x + y*w])
1312                     continue;
1313                 state[x + y*w]= 1;
1314                 v= !!src[x + y*stride];
1315                 put_cabac(&s->c, &b->state[0][0], v);
1316                 if(v){
1317                     SVI(significant, significant_count, x,y)
1318                     if(x     && !state[x - 1 +  y   *w]) SVI(boarder[0],boarder_count[0],x-1,y  )
1319                     if(y     && !state[x     + (y-1)*w]) SVI(boarder[1],boarder_count[1],x  ,y-1)
1320                     if(x+1<w && !state[x + 1 +  y   *w]) SVI(boarder[0],boarder_count[0],x+1,y  )
1321                     if(y+1<h && !state[x     + (y+1)*w]) SVI(boarder[1],boarder_count[1],x  ,y+1)
1322                     if(x     && y     && !state[x - 1 + (y-1)*w]) SVI(boarder[2],boarder_count[2],x-1,y-1)
1323                     if(x     && y+1<h && !state[x - 1 + (y+1)*w]) SVI(boarder[2],boarder_count[2],x-1,y+1)
1324                     if(x+1<w && y+1<h && !state[x + 1 + (y+1)*w]) SVI(boarder[2],boarder_count[2],x+1,y+1)
1325                     if(x+1<w && y     && !state[x + 1 + (y-1)*w]) SVI(boarder[2],boarder_count[2],x+1,y-1)
1326                 }
1327             }
1328             while(!boarder_count[0] && !boarder_count[1] && !boarder_count[2] && rle_pos < w*h){
1329                 int run=0;
1330                 for(; rle_pos < w*h;){
1331                     x= rle_pos % w; //FIXME speed
1332                     y= rle_pos / w;
1333                     rle_pos++;
1334                     if(state[x + y*w])
1335                         continue;
1336                     state[x + y*w]= 1;
1337                     v= !!src[x + y*stride];
1338                     if(v){
1339                         put_symbol(&s->c, b->state[1], run, 0);
1340                         SVI(significant, significant_count, x,y)
1341                         if(x     && !state[x - 1 +  y   *w]) SVI(boarder[0],boarder_count[0],x-1,y  )
1342                         if(y     && !state[x     + (y-1)*w]) SVI(boarder[1],boarder_count[1],x  ,y-1)
1343                         if(x+1<w && !state[x + 1 +  y   *w]) SVI(boarder[0],boarder_count[0],x+1,y  )
1344                         if(y+1<h && !state[x     + (y+1)*w]) SVI(boarder[1],boarder_count[1],x  ,y+1)
1345                         if(x     && y     && !state[x - 1 + (y-1)*w]) SVI(boarder[2],boarder_count[2],x-1,y-1)
1346                         if(x     && y+1<h && !state[x - 1 + (y+1)*w]) SVI(boarder[2],boarder_count[2],x-1,y+1)
1347                         if(x+1<w && y+1<h && !state[x + 1 + (y+1)*w]) SVI(boarder[2],boarder_count[2],x+1,y+1)
1348                         if(x+1<w && y     && !state[x + 1 + (y-1)*w]) SVI(boarder[2],boarder_count[2],x+1,y-1)
1349                         break;
1350 //FIXME                note only right & down can be boarders
1351                     }
1352                     run++;
1353                 }
1354             }
1355             if(!boarder_count[0] && !boarder_count[1] && !boarder_count[2])
1356                 break;
1357             
1358             while(boarder_count[0] || boarder_count[1] || boarder_count[2]){
1359                 int index;
1360                 
1361                 if     (boarder_count[  primary]) index=  primary;
1362                 else if(boarder_count[1-primary]) index=1-primary;
1363                 else                              index=2;
1364                 
1365                 boarder_count[index]--;
1366                 x= boarder[index][ boarder_count[index] ][0];
1367                 y= boarder[index][ boarder_count[index] ][1];
1368                 if(state[x + y*w]) //FIXME maybe check earlier
1369                     continue;
1370                 state[x + y*w]= 1;
1371                 v= !!src[x + y*stride];
1372                 put_cabac(&s->c, &b->state[0][index+1], v);
1373                 if(v){
1374                     SVI(significant, significant_count, x,y)
1375                     if(x     && !state[x - 1 +  y   *w]) SVI(boarder[0],boarder_count[0],x-1,y  )
1376                     if(y     && !state[x     + (y-1)*w]) SVI(boarder[1],boarder_count[1],x  ,y-1)
1377                     if(x+1<w && !state[x + 1 +  y   *w]) SVI(boarder[0],boarder_count[0],x+1,y  )
1378                     if(y+1<h && !state[x     + (y+1)*w]) SVI(boarder[1],boarder_count[1],x  ,y+1)
1379                     if(x     && y     && !state[x - 1 + (y-1)*w]) SVI(boarder[2],boarder_count[2],x-1,y-1)
1380                     if(x     && y+1<h && !state[x - 1 + (y+1)*w]) SVI(boarder[2],boarder_count[2],x-1,y+1)
1381                     if(x+1<w && y+1<h && !state[x + 1 + (y+1)*w]) SVI(boarder[2],boarder_count[2],x+1,y+1)
1382                     if(x+1<w && y     && !state[x + 1 + (y-1)*w]) SVI(boarder[2],boarder_count[2],x+1,y-1)
1383                 }
1384             }
1385         }
1386         //FIXME sort significant coeffs maybe
1387         if(1){
1388             qsort(significant, significant_count, sizeof(int16_t[2]), sig_cmp);
1389         }
1390         
1391         last_v=1;
1392         while(significant_count){
1393             int context= 3 + quant7[last_v&0xFF]; //use significance of suroundings
1394             significant_count--;
1395             x= significant[significant_count][0];//FIXME try opposit direction
1396             y= significant[significant_count][1];
1397             v= src[x + y*stride];
1398             put_symbol(&s->c, b->state[context + 2], v, 1); //FIXME try to avoid first bit, try this with the old code too!!
1399             last_v= v;
1400         }
1401     }
1402 #endif
1403     if(1){
1404         int run=0;
1405         int runs[w*h];
1406         int run_index=0;
1407                 
1408         for(y=0; y<h; y++){
1409             for(x=0; x<w; x++){
1410                 int v, p=0;
1411                 int /*ll=0, */l=0, lt=0, t=0, rt=0;
1412                 v= src[x + y*stride];
1413
1414                 if(y){
1415                     t= src[x + (y-1)*stride];
1416                     if(x){
1417                         lt= src[x - 1 + (y-1)*stride];
1418                     }
1419                     if(x + 1 < w){
1420                         rt= src[x + 1 + (y-1)*stride];
1421                     }
1422                 }
1423                 if(x){
1424                     l= src[x - 1 + y*stride];
1425                     /*if(x > 1){
1426                         if(orientation==1) ll= src[y + (x-2)*stride];
1427                         else               ll= src[x - 2 + y*stride];
1428                     }*/
1429                 }
1430                 if(parent){
1431                     int px= x>>1;
1432                     int py= y>>1;
1433                     if(px<b->parent->width && py<b->parent->height) 
1434                         p= parent[px + py*2*stride];
1435                 }
1436                 if(!(/*ll|*/l|lt|t|rt|p)){
1437                     if(v){
1438                         runs[run_index++]= run;
1439                         run=0;
1440                     }else{
1441                         run++;
1442                     }
1443                 }
1444             }
1445         }
1446         runs[run_index++]= run;
1447         run_index=0;
1448         run= runs[run_index++];
1449
1450         put_symbol(&s->c, b->state[1], run, 0);
1451         
1452         for(y=0; y<h; y++){
1453             for(x=0; x<w; x++){
1454                 int v, p=0;
1455                 int /*ll=0, */l=0, lt=0, t=0, rt=0;
1456                 v= src[x + y*stride];
1457
1458                 if(y){
1459                     t= src[x + (y-1)*stride];
1460                     if(x){
1461                         lt= src[x - 1 + (y-1)*stride];
1462                     }
1463                     if(x + 1 < w){
1464                         rt= src[x + 1 + (y-1)*stride];
1465                     }
1466                 }
1467                 if(x){
1468                     l= src[x - 1 + y*stride];
1469                     /*if(x > 1){
1470                         if(orientation==1) ll= src[y + (x-2)*stride];
1471                         else               ll= src[x - 2 + y*stride];
1472                     }*/
1473                 }
1474                 if(parent){
1475                     int px= x>>1;
1476                     int py= y>>1;
1477                     if(px<b->parent->width && py<b->parent->height) 
1478                         p= parent[px + py*2*stride];
1479                 }
1480                 if(/*ll|*/l|lt|t|rt|p){
1481                     int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));
1482
1483                     put_cabac(&s->c, &b->state[0][context], !!v);
1484                 }else{
1485                     if(!run){
1486                         run= runs[run_index++];
1487                         put_symbol(&s->c, b->state[1], run, 0);
1488                         assert(v);
1489                     }else{
1490                         run--;
1491                         assert(!v);
1492                     }
1493                 }
1494                 if(v){
1495                     int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));
1496
1497                     put_symbol(&s->c, b->state[context + 2], ABS(v)-1, 0);
1498                     put_cabac(&s->c, &b->state[0][16 + 1 + 3 + quant3b[l&0xFF] + 3*quant3b[t&0xFF]], v<0);
1499                 }
1500             }
1501         }
1502         return;
1503     }
1504 }
1505
1506 static inline void decode_subband(SnowContext *s, SubBand *b, DWTELEM *src, DWTELEM *parent, int stride, int orientation){
1507     const int level= b->level;
1508     const int w= b->width;
1509     const int h= b->height;
1510     int x,y;
1511
1512     START_TIMER
1513     
1514     if(1){
1515         int run;
1516                 
1517         for(y=0; y<b->height; y++)
1518             memset(&src[y*stride], 0, b->width*sizeof(DWTELEM));
1519
1520         run= get_symbol(&s->c, b->state[1], 0);
1521         for(y=0; y<h; y++){
1522             for(x=0; x<w; x++){
1523                 int v, p=0;
1524                 int /*ll=0, */l=0, lt=0, t=0, rt=0;
1525
1526                 if(y){
1527                     t= src[x + (y-1)*stride];
1528                     if(x){
1529                         lt= src[x - 1 + (y-1)*stride];
1530                     }
1531                     if(x + 1 < w){
1532                         rt= src[x + 1 + (y-1)*stride];
1533                     }
1534                 }
1535                 if(x){
1536                     l= src[x - 1 + y*stride];
1537                     /*if(x > 1){
1538                         if(orientation==1) ll= src[y + (x-2)*stride];
1539                         else               ll= src[x - 2 + y*stride];
1540                     }*/
1541                 }
1542                 if(parent){
1543                     int px= x>>1;
1544                     int py= y>>1;
1545                     if(px<b->parent->width && py<b->parent->height) 
1546                         p= parent[px + py*2*stride];
1547                 }
1548                 if(/*ll|*/l|lt|t|rt|p){
1549                     int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));
1550
1551                     v=get_cabac(&s->c, &b->state[0][context]);
1552                 }else{
1553                     if(!run){
1554                         run= get_symbol(&s->c, b->state[1], 0);
1555                         //FIXME optimize this here
1556                         //FIXME try to store a more naive run
1557                         v=1;
1558                     }else{
1559                         run--;
1560                         v=0;
1561                     }
1562                 }
1563                 if(v){
1564                     int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));
1565                     v= get_symbol(&s->c, b->state[context + 2], 0) + 1;
1566                     if(get_cabac(&s->c, &b->state[0][16 + 1 + 3 + quant3b[l&0xFF] + 3*quant3b[t&0xFF]]))
1567                         v= -v;
1568                     src[x + y*stride]= v;
1569                 }
1570             }
1571         }
1572         if(level+1 == s->spatial_decomposition_count){
1573             STOP_TIMER("decode_subband")
1574         }
1575         
1576         return;
1577     }
1578 }
1579
1580 static void reset_contexts(SnowContext *s){
1581     int plane_index, level, orientation;
1582
1583     for(plane_index=0; plane_index<2; plane_index++){
1584         for(level=0; level<s->spatial_decomposition_count; level++){
1585             for(orientation=level ? 1:0; orientation<4; orientation++){
1586                 memset(s->plane[plane_index].band[level][orientation].state, 0, sizeof(s->plane[plane_index].band[level][orientation].state));
1587             }
1588         }
1589     }
1590     memset(s->mb_band.state, 0, sizeof(s->mb_band.state));
1591     memset(s->mv_band[0].state, 0, sizeof(s->mv_band[0].state));
1592     memset(s->mv_band[1].state, 0, sizeof(s->mv_band[1].state));
1593     memset(s->header_state, 0, sizeof(s->header_state));
1594 }
1595
1596 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){
1597     int x, y;
1598
1599     for(y=0; y < b_h+5; y++){
1600         for(x=0; x < b_w; x++){
1601             int a0= src[x     + y*stride];
1602             int a1= src[x + 1 + y*stride];
1603             int a2= src[x + 2 + y*stride];
1604             int a3= src[x + 3 + y*stride];
1605             int a4= src[x + 4 + y*stride];
1606             int a5= src[x + 5 + y*stride];
1607 //            int am= 9*(a1+a2) - (a0+a3);
1608             int am= 20*(a2+a3) - 5*(a1+a4) + (a0+a5);
1609 //            int am= 18*(a2+a3) - 2*(a1+a4);
1610 //             int aL= (-7*a0 + 105*a1 + 35*a2 - 5*a3)>>3;
1611 //             int aR= (-7*a3 + 105*a2 + 35*a1 - 5*a0)>>3;
1612
1613 //            if(b_w==16) am= 8*(a1+a2);
1614
1615             if(dx<8) tmp[x + y*stride]= (32*a2*( 8-dx) +    am* dx    + 128)>>8;
1616             else     tmp[x + y*stride]= (   am*(16-dx) + 32*a3*(dx-8) + 128)>>8;
1617
1618 /*            if     (dx< 4) tmp[x + y*stride]= (16*a1*( 4-dx) +    aL* dx     + 32)>>6;
1619             else if(dx< 8) tmp[x + y*stride]= (   aL*( 8-dx) +    am*(dx- 4) + 32)>>6;
1620             else if(dx<12) tmp[x + y*stride]= (   am*(12-dx) +    aR*(dx- 8) + 32)>>6;
1621             else           tmp[x + y*stride]= (   aR*(16-dx) + 16*a2*(dx-12) + 32)>>6;*/
1622         }
1623     }
1624     for(y=0; y < b_h; y++){
1625         for(x=0; x < b_w; x++){
1626             int a0= tmp[x +  y     *stride];
1627             int a1= tmp[x + (y + 1)*stride];
1628             int a2= tmp[x + (y + 2)*stride];
1629             int a3= tmp[x + (y + 3)*stride];
1630             int a4= tmp[x + (y + 4)*stride];
1631             int a5= tmp[x + (y + 5)*stride];
1632             int am= 20*(a2+a3) - 5*(a1+a4) + (a0+a5);
1633 //            int am= 18*(a2+a3) - 2*(a1+a4);
1634 /*            int aL= (-7*a0 + 105*a1 + 35*a2 - 5*a3)>>3;
1635             int aR= (-7*a3 + 105*a2 + 35*a1 - 5*a0)>>3;*/
1636             
1637 //            if(b_w==16) am= 8*(a1+a2);
1638
1639             if(dy<8) dst[x + y*stride]= (32*a2*( 8-dy) +    am* dy    + 128)>>8;
1640             else     dst[x + y*stride]= (   am*(16-dy) + 32*a3*(dy-8) + 128)>>8;
1641
1642 /*            if     (dy< 4) tmp[x + y*stride]= (16*a1*( 4-dy) +    aL* dy     + 32)>>6;
1643             else if(dy< 8) tmp[x + y*stride]= (   aL*( 8-dy) +    am*(dy- 4) + 32)>>6;
1644             else if(dy<12) tmp[x + y*stride]= (   am*(12-dy) +    aR*(dy- 8) + 32)>>6;
1645             else           tmp[x + y*stride]= (   aR*(16-dy) + 16*a2*(dy-12) + 32)>>6;*/
1646         }
1647     }
1648 }
1649
1650 #define mcb(dx,dy,b_w)\
1651 static void mc_block ## dx ## dy(uint8_t *dst, uint8_t *src, int stride){\
1652     uint8_t tmp[stride*(b_w+5)];\
1653     mc_block(dst, src-2-2*stride, tmp, stride, b_w, b_w, dx, dy);\
1654 }
1655
1656 mcb( 0, 0,16)
1657 mcb( 4, 0,16)
1658 mcb( 8, 0,16)
1659 mcb(12, 0,16)
1660 mcb( 0, 4,16)
1661 mcb( 4, 4,16)
1662 mcb( 8, 4,16)
1663 mcb(12, 4,16)
1664 mcb( 0, 8,16)
1665 mcb( 4, 8,16)
1666 mcb( 8, 8,16)
1667 mcb(12, 8,16)
1668 mcb( 0,12,16)
1669 mcb( 4,12,16)
1670 mcb( 8,12,16)
1671 mcb(12,12,16)
1672
1673 #define mca(dx,dy,b_w)\
1674 static void mc_block_hpel ## dx ## dy(uint8_t *dst, uint8_t *src, int stride, int h){\
1675     uint8_t tmp[stride*(b_w+5)];\
1676     assert(h==b_w);\
1677     mc_block(dst, src-2-2*stride, tmp, stride, b_w, b_w, dx, dy);\
1678 }
1679
1680 mca( 0, 0,16)
1681 mca( 8, 0,16)
1682 mca( 0, 8,16)
1683 mca( 8, 8,16)
1684
1685 static void add_xblock(DWTELEM *dst, uint8_t *src, uint8_t *obmc, int s_x, int s_y, int b_w, int b_h, int mv_x, int mv_y, int w, int h, int dst_stride, int src_stride, int obmc_stride, int mb_type, int add){
1686     uint8_t tmp[src_stride*(b_h+5)]; //FIXME move to context to gurantee alignment
1687     int x,y;
1688
1689     if(s_x<0){
1690         obmc -= s_x;
1691         b_w += s_x;
1692         s_x=0;
1693     }else if(s_x + b_w > w){
1694         b_w = w - s_x;
1695     }
1696     if(s_y<0){
1697         obmc -= s_y*obmc_stride;
1698         b_h += s_y;
1699         s_y=0;
1700     }else if(s_y + b_h> h){
1701         b_h = h - s_y;
1702     }
1703
1704     dst += s_x + s_y*dst_stride;
1705     
1706     if(mb_type==1){
1707         src += s_x + s_y*src_stride;
1708         for(y=0; y < b_h; y++){
1709             for(x=0; x < b_w; x++){
1710                 if(add) dst[x + y*dst_stride] += obmc[x + y*obmc_stride] * 128 * (256/OBMC_MAX);
1711                 else    dst[x + y*dst_stride] -= obmc[x + y*obmc_stride] * 128 * (256/OBMC_MAX);
1712             }
1713         }
1714     }else{
1715         int dx= mv_x&15;
1716         int dy= mv_y&15;
1717 //        int dxy= (mv_x&1) + 2*(mv_y&1);
1718
1719         s_x += (mv_x>>4) - 2;
1720         s_y += (mv_y>>4) - 2;
1721         src += s_x + s_y*src_stride;
1722         //use dsputil
1723     
1724         if(   (unsigned)s_x >= w - b_w - 4
1725            || (unsigned)s_y >= h - b_h - 4){
1726             ff_emulated_edge_mc(tmp + 32, src, src_stride, b_w+5, b_h+5, s_x, s_y, w, h);
1727             src= tmp + 32;
1728         }
1729
1730         if(mb_type==0){
1731             mc_block(tmp, src, tmp + 64+8, src_stride, b_w, b_h, dx, dy);
1732         }else{
1733             int sum=0;
1734             for(y=0; y < b_h; y++){
1735                 for(x=0; x < b_w; x++){
1736                     sum += src[x+  y*src_stride];
1737                 }
1738             }
1739             sum= (sum + b_h*b_w/2) / (b_h*b_w);
1740             for(y=0; y < b_h; y++){
1741                 for(x=0; x < b_w; x++){
1742                     tmp[x + y*src_stride]= sum; 
1743                 }
1744             }
1745         }
1746
1747         for(y=0; y < b_h; y++){
1748             for(x=0; x < b_w; x++){
1749                 if(add) dst[x + y*dst_stride] += obmc[x + y*obmc_stride] * tmp[x + y*src_stride] * (256/OBMC_MAX);
1750                 else    dst[x + y*dst_stride] -= obmc[x + y*obmc_stride] * tmp[x + y*src_stride] * (256/OBMC_MAX);
1751             }
1752         }
1753     }
1754 }
1755
1756 static void predict_plane(SnowContext *s, DWTELEM *buf, int plane_index, int add){
1757     Plane *p= &s->plane[plane_index];
1758     const int mb_w= s->mb_band.width;
1759     const int mb_h= s->mb_band.height;
1760     const int mb_stride= s->mb_band.stride;
1761     int x, y, mb_x, mb_y;
1762     int scale      = plane_index ?  s->mv_scale : 2*s->mv_scale;
1763     int block_w    = plane_index ?  8 : 16;
1764     uint8_t *obmc  = plane_index ? obmc16 : obmc32;
1765     int obmc_stride= plane_index ? 16 : 32;
1766     int ref_stride= s->last_picture.linesize[plane_index];
1767     uint8_t *ref  = s->last_picture.data[plane_index];
1768     int w= p->width;
1769     int h= p->height;
1770     
1771 if(s->avctx->debug&512){
1772     for(y=0; y<h; y++){
1773         for(x=0; x<w; x++){
1774             if(add) buf[x + y*w]+= 128*256;
1775             else    buf[x + y*w]-= 128*256;
1776         }
1777     }
1778     
1779     return;
1780 }
1781     for(mb_y=-1; mb_y<=mb_h; mb_y++){
1782         for(mb_x=-1; mb_x<=mb_w; mb_x++){
1783             int index= clip(mb_x, 0, mb_w-1) + clip(mb_y, 0, mb_h-1)*mb_stride;
1784
1785             add_xblock(buf, ref, obmc, 
1786                        block_w*mb_x - block_w/2, 
1787                        block_w*mb_y - block_w/2,
1788                        2*block_w, 2*block_w,
1789                        s->mv_band[0].buf[index]*scale, s->mv_band[1].buf[index]*scale,
1790                        w, h,
1791                        w, ref_stride, obmc_stride, 
1792                        s->mb_band.buf[index], add);
1793
1794         }
1795     }
1796 }
1797
1798 static void quantize(SnowContext *s, SubBand *b, DWTELEM *src, int stride, int bias){
1799     const int level= b->level;
1800     const int w= b->width;
1801     const int h= b->height;
1802     const int qlog= clip(s->qlog + b->qlog, 0, 128);
1803     const int qmul= qexp[qlog&7]<<(qlog>>3);
1804     int x,y, thres1, thres2;
1805     START_TIMER
1806
1807     assert(QROOT==8);
1808
1809     bias= bias ? 0 : (3*qmul)>>3;
1810     thres1= ((qmul - bias)>>QEXPSHIFT) - 1;
1811     thres2= 2*thres1;
1812     
1813     if(!bias){
1814         for(y=0; y<h; y++){
1815             for(x=0; x<w; x++){
1816                 int i= src[x + y*stride];
1817                 
1818                 if((unsigned)(i+thres1) > thres2){
1819                     if(i>=0){
1820                         i<<= QEXPSHIFT;
1821                         i/= qmul; //FIXME optimize
1822                         src[x + y*stride]=  i;
1823                     }else{
1824                         i= -i;
1825                         i<<= QEXPSHIFT;
1826                         i/= qmul; //FIXME optimize
1827                         src[x + y*stride]= -i;
1828                     }
1829                 }else
1830                     src[x + y*stride]= 0;
1831             }
1832         }
1833     }else{
1834         for(y=0; y<h; y++){
1835             for(x=0; x<w; x++){
1836                 int i= src[x + y*stride]; 
1837                 
1838                 if((unsigned)(i+thres1) > thres2){
1839                     if(i>=0){
1840                         i<<= QEXPSHIFT;
1841                         i= (i + bias) / qmul; //FIXME optimize
1842                         src[x + y*stride]=  i;
1843                     }else{
1844                         i= -i;
1845                         i<<= QEXPSHIFT;
1846                         i= (i + bias) / qmul; //FIXME optimize
1847                         src[x + y*stride]= -i;
1848                     }
1849                 }else
1850                     src[x + y*stride]= 0;
1851             }
1852         }
1853     }
1854     if(level+1 == s->spatial_decomposition_count){
1855 //        STOP_TIMER("quantize")
1856     }
1857 }
1858
1859 static void dequantize(SnowContext *s, SubBand *b, DWTELEM *src, int stride){
1860     const int level= b->level;
1861     const int w= b->width;
1862     const int h= b->height;
1863     const int qlog= clip(s->qlog + b->qlog, 0, 128);
1864     const int qmul= qexp[qlog&7]<<(qlog>>3);
1865     const int qadd= (s->qbias*qmul)>>QBIAS_SHIFT;
1866     int x,y;
1867     
1868     assert(QROOT==8);
1869
1870     for(y=0; y<h; y++){
1871         for(x=0; x<w; x++){
1872             int i= src[x + y*stride];
1873             if(i<0){
1874                 src[x + y*stride]= -((-i*qmul + qadd)>>(QEXPSHIFT)); //FIXME try different bias
1875             }else if(i>0){
1876                 src[x + y*stride]=  (( i*qmul + qadd)>>(QEXPSHIFT));
1877             }
1878         }
1879     }
1880 }
1881
1882 static void decorrelate(SnowContext *s, SubBand *b, DWTELEM *src, int stride, int inverse, int use_median){
1883     const int w= b->width;
1884     const int h= b->height;
1885     int x,y;
1886     
1887     for(y=h-1; y>=0; y--){
1888         for(x=w-1; x>=0; x--){
1889             int i= x + y*stride;
1890             
1891             if(x){
1892                 if(use_median){
1893                     if(y && x+1<w) src[i] -= mid_pred(src[i - 1], src[i - stride], src[i - stride + 1]);
1894                     else  src[i] -= src[i - 1];
1895                 }else{
1896                     if(y) src[i] -= mid_pred(src[i - 1], src[i - stride], src[i - 1] + src[i - stride] - src[i - 1 - stride]);
1897                     else  src[i] -= src[i - 1];
1898                 }
1899             }else{
1900                 if(y) src[i] -= src[i - stride];
1901             }
1902         }
1903     }
1904 }
1905
1906 static void correlate(SnowContext *s, SubBand *b, DWTELEM *src, int stride, int inverse, int use_median){
1907     const int w= b->width;
1908     const int h= b->height;
1909     int x,y;
1910     
1911     for(y=0; y<h; y++){
1912         for(x=0; x<w; x++){
1913             int i= x + y*stride;
1914             
1915             if(x){
1916                 if(use_median){
1917                     if(y && x+1<w) src[i] += mid_pred(src[i - 1], src[i - stride], src[i - stride + 1]);
1918                     else  src[i] += src[i - 1];
1919                 }else{
1920                     if(y) src[i] += mid_pred(src[i - 1], src[i - stride], src[i - 1] + src[i - stride] - src[i - 1 - stride]);
1921                     else  src[i] += src[i - 1];
1922                 }
1923             }else{
1924                 if(y) src[i] += src[i - stride];
1925             }
1926         }
1927     }
1928 }
1929
1930 static void encode_header(SnowContext *s){
1931     int plane_index, level, orientation;
1932
1933     put_cabac(&s->c, s->header_state, s->keyframe); // state clearing stuff?
1934     if(s->keyframe){
1935         put_symbol(&s->c, s->header_state, s->version, 0);
1936         put_symbol(&s->c, s->header_state, s->temporal_decomposition_type, 0);
1937         put_symbol(&s->c, s->header_state, s->temporal_decomposition_count, 0);
1938         put_symbol(&s->c, s->header_state, s->spatial_decomposition_count, 0);
1939         put_symbol(&s->c, s->header_state, s->colorspace_type, 0);
1940         put_symbol(&s->c, s->header_state, s->b_width, 0);
1941         put_symbol(&s->c, s->header_state, s->b_height, 0);
1942         put_symbol(&s->c, s->header_state, s->chroma_h_shift, 0);
1943         put_symbol(&s->c, s->header_state, s->chroma_v_shift, 0);
1944         put_cabac(&s->c, s->header_state, s->spatial_scalability);
1945 //        put_cabac(&s->c, s->header_state, s->rate_scalability);
1946
1947         for(plane_index=0; plane_index<2; plane_index++){
1948             for(level=0; level<s->spatial_decomposition_count; level++){
1949                 for(orientation=level ? 1:0; orientation<4; orientation++){
1950                     if(orientation==2) continue;
1951                     put_symbol(&s->c, s->header_state, s->plane[plane_index].band[level][orientation].qlog, 1);
1952                 }
1953             }
1954         }
1955     }
1956     put_symbol(&s->c, s->header_state, s->spatial_decomposition_type, 0);
1957     put_symbol(&s->c, s->header_state, s->qlog, 1); 
1958     put_symbol(&s->c, s->header_state, s->mv_scale, 0); 
1959     put_symbol(&s->c, s->header_state, s->qbias, 1);
1960 }
1961
1962 static int decode_header(SnowContext *s){
1963     int plane_index, level, orientation;
1964
1965     s->keyframe= get_cabac(&s->c, s->header_state);
1966     if(s->keyframe){
1967         s->version= get_symbol(&s->c, s->header_state, 0);
1968         if(s->version>0){
1969             av_log(s->avctx, AV_LOG_ERROR, "version %d not supported", s->version);
1970             return -1;
1971         }
1972         s->temporal_decomposition_type= get_symbol(&s->c, s->header_state, 0);
1973         s->temporal_decomposition_count= get_symbol(&s->c, s->header_state, 0);
1974         s->spatial_decomposition_count= get_symbol(&s->c, s->header_state, 0);
1975         s->colorspace_type= get_symbol(&s->c, s->header_state, 0);
1976         s->b_width= get_symbol(&s->c, s->header_state, 0);
1977         s->b_height= get_symbol(&s->c, s->header_state, 0);
1978         s->chroma_h_shift= get_symbol(&s->c, s->header_state, 0);
1979         s->chroma_v_shift= get_symbol(&s->c, s->header_state, 0);
1980         s->spatial_scalability= get_cabac(&s->c, s->header_state);
1981 //        s->rate_scalability= get_cabac(&s->c, s->header_state);
1982
1983         for(plane_index=0; plane_index<3; plane_index++){
1984             for(level=0; level<s->spatial_decomposition_count; level++){
1985                 for(orientation=level ? 1:0; orientation<4; orientation++){
1986                     int q;
1987                     if     (plane_index==2) q= s->plane[1].band[level][orientation].qlog;
1988                     else if(orientation==2) q= s->plane[plane_index].band[level][1].qlog;
1989                     else                    q= get_symbol(&s->c, s->header_state, 1);
1990                     s->plane[plane_index].band[level][orientation].qlog= q;
1991                 }
1992             }
1993         }
1994     }
1995     
1996     s->spatial_decomposition_type= get_symbol(&s->c, s->header_state, 0);
1997     if(s->spatial_decomposition_type > 2){
1998         av_log(s->avctx, AV_LOG_ERROR, "spatial_decomposition_type %d not supported", s->spatial_decomposition_type);
1999         return -1;
2000     }
2001     
2002     s->qlog= get_symbol(&s->c, s->header_state, 1);
2003     s->mv_scale= get_symbol(&s->c, s->header_state, 0);
2004     s->qbias= get_symbol(&s->c, s->header_state, 1);
2005
2006     return 0;
2007 }
2008
2009 static int common_init(AVCodecContext *avctx){
2010     SnowContext *s = avctx->priv_data;
2011     int width, height;
2012     int level, orientation, plane_index, dec;
2013
2014     s->avctx= avctx;
2015         
2016     dsputil_init(&s->dsp, avctx);
2017
2018 #define mcf(dx,dy)\
2019     s->dsp.put_qpel_pixels_tab       [0][dy+dx/4]=\
2020     s->dsp.put_no_rnd_qpel_pixels_tab[0][dy+dx/4]=\
2021         mc_block ## dx ## dy;
2022
2023     mcf( 0, 0)
2024     mcf( 4, 0)
2025     mcf( 8, 0)
2026     mcf(12, 0)
2027     mcf( 0, 4)
2028     mcf( 4, 4)
2029     mcf( 8, 4)
2030     mcf(12, 4)
2031     mcf( 0, 8)
2032     mcf( 4, 8)
2033     mcf( 8, 8)
2034     mcf(12, 8)
2035     mcf( 0,12)
2036     mcf( 4,12)
2037     mcf( 8,12)
2038     mcf(12,12)
2039
2040 #define mcfh(dx,dy)\
2041     s->dsp.put_pixels_tab       [0][dy/4+dx/8]=\
2042     s->dsp.put_no_rnd_pixels_tab[0][dy/4+dx/8]=\
2043         mc_block_hpel ## dx ## dy;
2044
2045     mcfh(0, 0)
2046     mcfh(8, 0)
2047     mcfh(0, 8)
2048     mcfh(8, 8)
2049         
2050     dec= s->spatial_decomposition_count= 5;
2051     s->spatial_decomposition_type= avctx->prediction_method; //FIXME add decorrelator type r transform_type
2052     
2053     s->chroma_h_shift= 1; //FIXME XXX
2054     s->chroma_v_shift= 1;
2055     
2056 //    dec += FFMAX(s->chroma_h_shift, s->chroma_v_shift);
2057     
2058     s->b_width = (s->avctx->width +(1<<dec)-1)>>dec;
2059     s->b_height= (s->avctx->height+(1<<dec)-1)>>dec;
2060     
2061     s->spatial_dwt_buffer= av_mallocz(s->b_width*s->b_height*sizeof(DWTELEM)<<(2*dec));
2062     s->pred_buffer= av_mallocz(s->b_width*s->b_height*sizeof(DWTELEM)<<(2*dec));
2063     
2064     s->mv_scale= (s->avctx->flags & CODEC_FLAG_QPEL) ? 2 : 4;
2065     
2066     for(plane_index=0; plane_index<3; plane_index++){    
2067         int w= s->avctx->width;
2068         int h= s->avctx->height;
2069
2070         if(plane_index){
2071             w>>= s->chroma_h_shift;
2072             h>>= s->chroma_v_shift;
2073         }
2074         s->plane[plane_index].width = w;
2075         s->plane[plane_index].height= h;
2076 av_log(NULL, AV_LOG_DEBUG, "%d %d\n", w, h);
2077         for(level=s->spatial_decomposition_count-1; level>=0; level--){
2078             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2079                 SubBand *b= &s->plane[plane_index].band[level][orientation];
2080                 
2081                 b->buf= s->spatial_dwt_buffer;
2082                 b->level= level;
2083                 b->stride= s->plane[plane_index].width << (s->spatial_decomposition_count - level);
2084                 b->width = (w + !(orientation&1))>>1;
2085                 b->height= (h + !(orientation>1))>>1;
2086                 
2087                 if(orientation&1) b->buf += (w+1)>>1;
2088                 if(orientation>1) b->buf += b->stride>>1;
2089                 
2090                 if(level)
2091                     b->parent= &s->plane[plane_index].band[level-1][orientation];
2092             }
2093             w= (w+1)>>1;
2094             h= (h+1)>>1;
2095         }
2096     }
2097     
2098     //FIXME init_subband() ?
2099     s->mb_band.stride= s->mv_band[0].stride= s->mv_band[1].stride=
2100     s->mb_band.width = s->mv_band[0].width = s->mv_band[1].width = (s->avctx->width + 15)>>4;
2101     s->mb_band.height= s->mv_band[0].height= s->mv_band[1].height= (s->avctx->height+ 15)>>4;
2102     s->mb_band   .buf= av_mallocz(s->mb_band   .stride * s->mb_band   .height*sizeof(DWTELEM));
2103     s->mv_band[0].buf= av_mallocz(s->mv_band[0].stride * s->mv_band[0].height*sizeof(DWTELEM));
2104     s->mv_band[1].buf= av_mallocz(s->mv_band[1].stride * s->mv_band[1].height*sizeof(DWTELEM));
2105
2106     reset_contexts(s);
2107 /*    
2108     width= s->width= avctx->width;
2109     height= s->height= avctx->height;
2110     
2111     assert(width && height);
2112 */
2113     s->avctx->get_buffer(s->avctx, &s->mconly_picture);
2114     
2115     return 0;
2116 }
2117
2118
2119 static void calculate_vissual_weight(SnowContext *s, Plane *p){
2120     int width = p->width;
2121     int height= p->height;
2122     int i, level, orientation, x, y;
2123
2124     for(level=0; level<s->spatial_decomposition_count; level++){
2125         for(orientation=level ? 1 : 0; orientation<4; orientation++){
2126             SubBand *b= &p->band[level][orientation];
2127             DWTELEM *buf= b->buf;
2128             int64_t error=0;
2129             
2130             memset(s->spatial_dwt_buffer, 0, sizeof(int)*width*height);
2131             buf[b->width/2 + b->height/2*b->stride]= 256*256;
2132             spatial_idwt(s, s->spatial_dwt_buffer, width, height, width);
2133             for(y=0; y<height; y++){
2134                 for(x=0; x<width; x++){
2135                     int64_t d= s->spatial_dwt_buffer[x + y*width];
2136                     error += d*d;
2137                 }
2138             }
2139
2140             b->qlog= (int)(log(352256.0/sqrt(error)) / log(pow(2.0, 1.0/QROOT))+0.5);
2141             av_log(NULL, AV_LOG_DEBUG, "%d %d %d\n", level, orientation, b->qlog/*, sqrt(error)*/);
2142         }
2143     }
2144 }
2145
2146 static int encode_init(AVCodecContext *avctx)
2147 {
2148     SnowContext *s = avctx->priv_data;
2149     int i;
2150     int level, orientation, plane_index;
2151
2152     common_init(avctx);
2153  
2154     s->version=0;
2155     
2156     s->m.me.scratchpad= av_mallocz((avctx->width+64)*2*16*2*sizeof(uint8_t));
2157     s->m.me.map       = av_mallocz(ME_MAP_SIZE*sizeof(uint32_t));
2158     s->m.me.score_map = av_mallocz(ME_MAP_SIZE*sizeof(uint32_t));
2159     s->mb_type        = av_mallocz((s->mb_band.width+1)*s->mb_band.height*sizeof(int16_t));
2160     s->mb_mean        = av_mallocz((s->mb_band.width+1)*s->mb_band.height*sizeof(int8_t ));
2161     s->dummy          = av_mallocz((s->mb_band.width+1)*s->mb_band.height*sizeof(int32_t));
2162     h263_encode_init(&s->m); //mv_penalty
2163
2164     for(plane_index=0; plane_index<3; plane_index++){
2165         calculate_vissual_weight(s, &s->plane[plane_index]);
2166     }
2167     
2168     
2169     avctx->coded_frame= &s->current_picture;
2170     switch(avctx->pix_fmt){
2171 //    case PIX_FMT_YUV444P:
2172 //    case PIX_FMT_YUV422P:
2173     case PIX_FMT_YUV420P:
2174     case PIX_FMT_GRAY8:
2175 //    case PIX_FMT_YUV411P:
2176 //    case PIX_FMT_YUV410P:
2177         s->colorspace_type= 0;
2178         break;
2179 /*    case PIX_FMT_RGBA32:
2180         s->colorspace= 1;
2181         break;*/
2182     default:
2183         av_log(avctx, AV_LOG_ERROR, "format not supported\n");
2184         return -1;
2185     }
2186 //    avcodec_get_chroma_sub_sample(avctx->pix_fmt, &s->chroma_h_shift, &s->chroma_v_shift);
2187     s->chroma_h_shift= 1;
2188     s->chroma_v_shift= 1;
2189     return 0;
2190 }
2191
2192 static int frame_start(SnowContext *s){
2193    AVFrame tmp;
2194
2195    if(s->keyframe)
2196         reset_contexts(s);
2197  
2198     tmp= s->last_picture;
2199     s->last_picture= s->current_picture;
2200     s->current_picture= tmp;
2201     
2202     s->current_picture.reference= 1;
2203     if(s->avctx->get_buffer(s->avctx, &s->current_picture) < 0){
2204         av_log(s->avctx, AV_LOG_ERROR, "get_buffer() failed\n");
2205         return -1;
2206     }
2207     
2208     return 0;
2209 }
2210
2211 static int encode_frame(AVCodecContext *avctx, unsigned char *buf, int buf_size, void *data){
2212     SnowContext *s = avctx->priv_data;
2213     CABACContext * const c= &s->c;
2214     AVFrame *pict = data;
2215     const int width= s->avctx->width;
2216     const int height= s->avctx->height;
2217     int used_count= 0;
2218     int log2_threshold, level, orientation, plane_index, i;
2219
2220     if(avctx->strict_std_compliance >= 0){
2221         av_log(avctx, AV_LOG_ERROR, "this codec is under development, files encoded with it wont be decodeable with future versions!!!\n"
2222                "use vstrict=-1 to use it anyway\n");
2223         return -1;
2224     }
2225         
2226     ff_init_cabac_encoder(c, buf, buf_size);
2227     ff_init_cabac_states(c, ff_h264_lps_range, ff_h264_mps_state, ff_h264_lps_state, 64);
2228     
2229     s->input_picture = *pict;
2230
2231     memset(s->header_state, 0, sizeof(s->header_state));
2232
2233     s->keyframe=avctx->gop_size==0 || avctx->frame_number % avctx->gop_size == 0;
2234     pict->pict_type= s->keyframe ? FF_I_TYPE : FF_P_TYPE;
2235     
2236     s->qlog= rint(QROOT*log(pict->quality / (float)FF_QP2LAMBDA)/log(2));
2237     //<64 >60
2238     s->qlog += 61;
2239
2240     for(i=0; i<s->mb_band.stride * s->mb_band.height; i++){
2241         s->mb_band.buf[i]= s->keyframe;
2242     }
2243     
2244     frame_start(s);
2245
2246     if(pict->pict_type == P_TYPE){
2247         int block_width = (width +15)>>4;
2248         int block_height= (height+15)>>4;
2249         int stride= s->current_picture.linesize[0];
2250         uint8_t *src_plane= s->input_picture.data[0];
2251         int src_stride= s->input_picture.linesize[0];
2252         int x,y;
2253         
2254         assert(s->current_picture.data[0]);
2255         assert(s->last_picture.data[0]);
2256      
2257         s->m.avctx= s->avctx;
2258         s->m.current_picture.data[0]= s->current_picture.data[0];
2259         s->m.   last_picture.data[0]= s->   last_picture.data[0];
2260         s->m.    new_picture.data[0]= s->  input_picture.data[0];
2261         s->m.current_picture_ptr= &s->m.current_picture;
2262         s->m.   last_picture_ptr= &s->m.   last_picture;
2263         s->m.linesize=
2264         s->m.   last_picture.linesize[0]=
2265         s->m.    new_picture.linesize[0]=
2266         s->m.current_picture.linesize[0]= stride;
2267         s->m.width = width;
2268         s->m.height= height;
2269         s->m.mb_width = block_width;
2270         s->m.mb_height= block_height;
2271         s->m.mb_stride=   s->m.mb_width+1;
2272         s->m.b8_stride= 2*s->m.mb_width+1;
2273         s->m.f_code=1;
2274         s->m.pict_type= pict->pict_type;
2275         s->m.me_method= s->avctx->me_method;
2276         s->m.me.scene_change_score=0;
2277         s->m.flags= s->avctx->flags;
2278         s->m.quarter_sample= (s->avctx->flags & CODEC_FLAG_QPEL)!=0;
2279         s->m.out_format= FMT_H263;
2280         s->m.unrestricted_mv= 1;
2281
2282         s->m.lambda= pict->quality * 3/2; //FIXME bug somewhere else
2283         s->m.qscale= (s->m.lambda*139 + FF_LAMBDA_SCALE*64) >> (FF_LAMBDA_SHIFT + 7);
2284         s->m.lambda2= (s->m.lambda*s->m.lambda + FF_LAMBDA_SCALE/2) >> FF_LAMBDA_SHIFT;
2285
2286         if(!s->motion_val8){
2287             s->motion_val8 = av_mallocz(s->m.b8_stride*block_height*2*2*sizeof(int16_t));
2288             s->motion_val16= av_mallocz(s->m.mb_stride*block_height*2*sizeof(int16_t));
2289         }
2290         
2291         s->m.mb_type= s->mb_type;
2292         
2293         //dummies, to avoid segfaults
2294         s->m.current_picture.mb_mean  = s->mb_mean;
2295         s->m.current_picture.mb_var   = (int16_t*)s->dummy;
2296         s->m.current_picture.mc_mb_var= (int16_t*)s->dummy;
2297         s->m.current_picture.mb_type  = s->dummy;
2298         
2299         s->m.current_picture.motion_val[0]= s->motion_val8;
2300         s->m.p_mv_table= s->motion_val16;
2301         s->m.dsp= s->dsp; //move
2302         ff_init_me(&s->m);
2303     
2304         
2305         s->m.me.pre_pass=1;
2306         s->m.me.dia_size= s->avctx->pre_dia_size;
2307         s->m.first_slice_line=1;
2308         for(y= block_height-1; y >= 0; y--) {
2309             uint8_t src[stride*16];
2310
2311             s->m.new_picture.data[0]= src - y*16*stride; //ugly
2312             s->m.mb_y= y;
2313             for(i=0; i<16 && i + 16*y<height; i++){
2314                 memcpy(&src[i*stride], &src_plane[(i+16*y)*src_stride], width);
2315                 for(x=width; x<16*block_width; x++)
2316                     src[i*stride+x]= src[i*stride+x-1];
2317             }
2318             for(; i<16 && i + 16*y<16*block_height; i++)
2319                 memcpy(&src[i*stride], &src[(i-1)*stride], 16*block_width);
2320
2321             for(x=block_width-1; x >=0 ;x--) {
2322                 s->m.mb_x= x;
2323                 ff_init_block_index(&s->m);
2324                 ff_update_block_index(&s->m);
2325                 ff_pre_estimate_p_frame_motion(&s->m, x, y);
2326             }
2327             s->m.first_slice_line=0;
2328         }        
2329         s->m.me.pre_pass=0;
2330         
2331         
2332         s->m.me.dia_size= s->avctx->dia_size;
2333         s->m.first_slice_line=1;
2334         for (y = 0; y < block_height; y++) {
2335             uint8_t src[stride*16];
2336
2337             s->m.new_picture.data[0]= src - y*16*stride; //ugly
2338             s->m.mb_y= y;
2339             
2340             assert(width <= stride);
2341             assert(width <= 16*block_width);
2342     
2343             for(i=0; i<16 && i + 16*y<height; i++){
2344                 memcpy(&src[i*stride], &src_plane[(i+16*y)*src_stride], width);
2345                 for(x=width; x<16*block_width; x++)
2346                     src[i*stride+x]= src[i*stride+x-1];
2347             }
2348             for(; i<16 && i + 16*y<16*block_height; i++)
2349                 memcpy(&src[i*stride], &src[(i-1)*stride], 16*block_width);
2350     
2351             for (x = 0; x < block_width; x++) {
2352                 int mb_xy= x + y*(s->mb_band.stride);
2353                 s->m.mb_x= x;
2354                 ff_init_block_index(&s->m);
2355                 ff_update_block_index(&s->m);
2356                 
2357                 ff_estimate_p_frame_motion(&s->m, x, y);
2358                 
2359                 s->mb_band   .buf[mb_xy]= (s->m.mb_type[x + y*s->m.mb_stride]&CANDIDATE_MB_TYPE_INTER)
2360                  ? 0 : 2;
2361                 s->mv_band[0].buf[mb_xy]= s->motion_val16[x + y*s->m.mb_stride][0];
2362                 s->mv_band[1].buf[mb_xy]= s->motion_val16[x + y*s->m.mb_stride][1];
2363                 
2364                 if(s->mb_band   .buf[x + y*(s->mb_band.stride)]==2 && 0){
2365                     int dc0=128, dc1=128, dc, dc2, dir;
2366                     int offset= (s->avctx->flags & CODEC_FLAG_QPEL) ? 64 : 32;
2367                     
2368                     dc       =s->mb_mean[x +  y   *s->m.mb_stride    ];
2369                     if(x) dc0=s->mb_mean[x +  y   *s->m.mb_stride - 1];
2370                     if(y) dc1=s->mb_mean[x + (y-1)*s->m.mb_stride    ];
2371                     dc2= (dc0+dc1)>>1;
2372 #if 0
2373                     if     (ABS(dc0 - dc) < ABS(dc1 - dc) && ABS(dc0 - dc) < ABS(dc2 - dc))
2374                         dir= 1;
2375                     else if(ABS(dc0 - dc) >=ABS(dc1 - dc) && ABS(dc1 - dc) < ABS(dc2 - dc))
2376                         dir=-1;
2377                     else
2378                         dir=0;
2379 #endif                    
2380                     if(ABS(dc0 - dc) < ABS(dc1 - dc) && x){
2381                         s->mv_band[0].buf[mb_xy]= s->mv_band[0].buf[x + y*(s->mb_band.stride)-1] - offset;
2382                         s->mv_band[1].buf[mb_xy]= s->mv_band[1].buf[x + y*(s->mb_band.stride)-1];
2383                         s->mb_mean[x +  y   *s->m.mb_stride    ]= dc0;
2384                     }else if(y){
2385                         s->mv_band[0].buf[mb_xy]= s->mv_band[0].buf[x + (y-1)*(s->mb_band.stride)];
2386                         s->mv_band[1].buf[mb_xy]= s->mv_band[1].buf[x + (y-1)*(s->mb_band.stride)] - offset;
2387                         s->mb_mean[x +  y   *s->m.mb_stride    ]= dc1;
2388                     }
2389                 }
2390 //                s->mb_band   .buf[x + y*(s->mb_band.stride)]=1; //FIXME intra only test
2391             }
2392             s->m.first_slice_line=0;
2393         }
2394         assert(s->m.pict_type == P_TYPE);
2395         if(s->m.me.scene_change_score > s->avctx->scenechange_threshold){
2396             s->m.pict_type= 
2397             pict->pict_type =I_TYPE;
2398             for(i=0; i<s->mb_band.stride * s->mb_band.height; i++){
2399                 s->mb_band.buf[i]= 1;
2400                 s->mv_band[0].buf[i]=
2401                 s->mv_band[1].buf[i]= 0;
2402             }
2403     //printf("Scene change detected, encoding as I Frame %d %d\n", s->current_picture.mb_var_sum, s->current_picture.mc_mb_var_sum);
2404         }        
2405     }
2406         
2407     s->m.first_slice_line=1;
2408     
2409     s->qbias= pict->pict_type == P_TYPE ? 2 : 0;
2410
2411     encode_header(s);
2412     
2413     decorrelate(s, &s->mb_band   , s->mb_band   .buf, s->mb_band   .stride, 0, 1);
2414     decorrelate(s, &s->mv_band[0], s->mv_band[0].buf, s->mv_band[0].stride, 0, 1);
2415     decorrelate(s, &s->mv_band[1], s->mv_band[1].buf, s->mv_band[1].stride, 0 ,1);
2416     encode_subband(s, &s->mb_band   , s->mb_band   .buf, NULL, s->mb_band   .stride, 0);
2417     encode_subband(s, &s->mv_band[0], s->mv_band[0].buf, NULL, s->mv_band[0].stride, 0);
2418     encode_subband(s, &s->mv_band[1], s->mv_band[1].buf, NULL, s->mv_band[1].stride, 0);
2419     
2420 //FIXME avoid this
2421     correlate(s, &s->mb_band   , s->mb_band   .buf, s->mb_band   .stride, 1, 1);
2422     correlate(s, &s->mv_band[0], s->mv_band[0].buf, s->mv_band[0].stride, 1, 1);
2423     correlate(s, &s->mv_band[1], s->mv_band[1].buf, s->mv_band[1].stride, 1, 1);
2424     
2425     for(plane_index=0; plane_index<3; plane_index++){
2426         Plane *p= &s->plane[plane_index];
2427         int w= p->width;
2428         int h= p->height;
2429         int x, y;
2430         int bits= put_bits_count(&s->c.pb);
2431
2432         //FIXME optimize
2433 #if QPRED
2434         memset(s->pred_buffer, 0, sizeof(DWTELEM)*w*h);
2435         predict_plane(s, s->pred_buffer, plane_index, 1);
2436         spatial_dwt(s, s->pred_buffer, w, h, w);
2437         for(level=0; level<s->spatial_decomposition_count; level++){
2438             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2439                 SubBand *b= &p->band[level][orientation];
2440                 int delta= ((int)s->pred_buffer - (int)s->spatial_dwt_buffer)/sizeof(DWTELEM);
2441                 
2442                 quantize  (s, b, b->buf + delta, b->stride, s->qbias);
2443                 dequantize(s, b, b->buf + delta, b->stride);
2444             }
2445         }
2446         for(y=0; y<h; y++){
2447             for(x=0; x<w; x++){
2448                 s->spatial_dwt_buffer[y*w + x]= pict->data[plane_index][y*pict->linesize[plane_index] + x]<<8;
2449             }
2450         }
2451         spatial_dwt(s, s->spatial_dwt_buffer, w, h, w);
2452         for(y=0; y<h; y++){
2453             for(x=0; x<w; x++){
2454                 s->spatial_dwt_buffer[y*w + x]-= s->pred_buffer[y*w + x];
2455             }
2456         }
2457 #else
2458      if(pict->data[plane_index]) //FIXME gray hack
2459         for(y=0; y<h; y++){
2460             for(x=0; x<w; x++){
2461                 s->spatial_dwt_buffer[y*w + x]= pict->data[plane_index][y*pict->linesize[plane_index] + x]<<8;
2462             }
2463         }
2464         predict_plane(s, s->spatial_dwt_buffer, plane_index, 0);
2465         spatial_dwt(s, s->spatial_dwt_buffer, w, h, w);
2466 #endif
2467  
2468         for(level=0; level<s->spatial_decomposition_count; level++){
2469             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2470                 SubBand *b= &p->band[level][orientation];
2471                 
2472                 quantize(s, b, b->buf, b->stride, s->qbias);
2473                 if(orientation==0)
2474                     decorrelate(s, b, b->buf, b->stride, pict->pict_type == P_TYPE, 0);
2475                 encode_subband(s, b, b->buf, b->parent ? b->parent->buf : NULL, b->stride, orientation);
2476                 assert(b->parent==NULL || b->parent->stride == b->stride*2);
2477                 if(orientation==0)
2478                     correlate(s, b, b->buf, b->stride, 1, 0);
2479             }
2480         }
2481 //        av_log(NULL, AV_LOG_DEBUG, "plane:%d bits:%d\n", plane_index, put_bits_count(&s->c.pb) - bits);
2482
2483         for(level=0; level<s->spatial_decomposition_count; level++){
2484             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2485                 SubBand *b= &p->band[level][orientation];
2486
2487                 dequantize(s, b, b->buf, b->stride);
2488             }
2489         }
2490         
2491 #if QPRED
2492         for(y=0; y<h; y++){
2493             for(x=0; x<w; x++){
2494                 s->spatial_dwt_buffer[y*w + x]+= s->pred_buffer[y*w + x];
2495             }
2496         }
2497         spatial_idwt(s, s->spatial_dwt_buffer, w, h, w);
2498 #else
2499         spatial_idwt(s, s->spatial_dwt_buffer, w, h, w);
2500         predict_plane(s, s->spatial_dwt_buffer, plane_index, 1);
2501 #endif
2502         //FIXME optimize
2503         for(y=0; y<h; y++){
2504             for(x=0; x<w; x++){
2505                 int v= (s->spatial_dwt_buffer[y*w + x]+128)>>8;
2506                 if(v&(~255)) v= ~(v>>31);
2507                 s->current_picture.data[plane_index][y*s->current_picture.linesize[plane_index] + x]= v;
2508             }
2509         }
2510         if(s->avctx->flags&CODEC_FLAG_PSNR){
2511             int64_t error= 0;
2512             
2513     if(pict->data[plane_index]) //FIXME gray hack
2514             for(y=0; y<h; y++){
2515                 for(x=0; x<w; x++){
2516                     int d= s->spatial_dwt_buffer[y*w + x] - pict->data[plane_index][y*pict->linesize[plane_index] + x]*256;
2517                     error += d*d;
2518                 }
2519             }
2520             error= (error + 128*256)>>16;
2521             s->avctx->error[plane_index] += error;
2522             s->avctx->error[3] += error;
2523         }
2524     }
2525
2526     if(s->last_picture.data[0])
2527         avctx->release_buffer(avctx, &s->last_picture);
2528
2529     emms_c();
2530     
2531     return put_cabac_terminate(c, 1);
2532 }
2533
2534 static void common_end(SnowContext *s){
2535     av_freep(&s->spatial_dwt_buffer);
2536     av_freep(&s->mb_band.buf);
2537     av_freep(&s->mv_band[0].buf);
2538     av_freep(&s->mv_band[1].buf);
2539
2540     av_freep(&s->m.me.scratchpad);    
2541     av_freep(&s->m.me.map);
2542     av_freep(&s->m.me.score_map);
2543     av_freep(&s->mb_type);
2544     av_freep(&s->mb_mean);
2545     av_freep(&s->dummy);
2546     av_freep(&s->motion_val8);
2547     av_freep(&s->motion_val16);
2548 }
2549
2550 static int encode_end(AVCodecContext *avctx)
2551 {
2552     SnowContext *s = avctx->priv_data;
2553
2554     common_end(s);
2555
2556     return 0;
2557 }
2558
2559 static int decode_init(AVCodecContext *avctx)
2560 {
2561 //    SnowContext *s = avctx->priv_data;
2562
2563     common_init(avctx);
2564     
2565     return 0;
2566 }
2567
2568 static int decode_frame(AVCodecContext *avctx, void *data, int *data_size, uint8_t *buf, int buf_size){
2569     SnowContext *s = avctx->priv_data;
2570     CABACContext * const c= &s->c;
2571     const int width= s->avctx->width;
2572     const int height= s->avctx->height;
2573     int bytes_read;
2574     AVFrame *picture = data;
2575     int log2_threshold, level, orientation, plane_index;
2576     
2577
2578     /* no supplementary picture */
2579     if (buf_size == 0)
2580         return 0;
2581
2582     ff_init_cabac_decoder(c, buf, buf_size);
2583     ff_init_cabac_states(c, ff_h264_lps_range, ff_h264_mps_state, ff_h264_lps_state, 64);
2584
2585     memset(s->header_state, 0, sizeof(s->header_state));
2586
2587     s->current_picture.pict_type= FF_I_TYPE; //FIXME I vs. P
2588     decode_header(s);
2589
2590     frame_start(s);
2591     //keyframe flag dupliaction mess FIXME
2592     if(avctx->debug&FF_DEBUG_PICT_INFO)
2593         av_log(avctx, AV_LOG_ERROR, "keyframe:%d qlog:%d\n", s->keyframe, s->qlog);
2594     
2595     decode_subband(s, &s->mb_band   , s->mb_band   .buf, NULL, s->mb_band   .stride, 0);
2596     decode_subband(s, &s->mv_band[0], s->mv_band[0].buf, NULL, s->mv_band[0].stride, 0);
2597     decode_subband(s, &s->mv_band[1], s->mv_band[1].buf, NULL, s->mv_band[1].stride, 0);
2598     correlate(s, &s->mb_band   , s->mb_band   .buf, s->mb_band   .stride, 1, 1);
2599     correlate(s, &s->mv_band[0], s->mv_band[0].buf, s->mv_band[0].stride, 1, 1);
2600     correlate(s, &s->mv_band[1], s->mv_band[1].buf, s->mv_band[1].stride, 1, 1);
2601
2602     for(plane_index=0; plane_index<3; plane_index++){
2603         Plane *p= &s->plane[plane_index];
2604         int w= p->width;
2605         int h= p->height;
2606         int x, y;
2607         
2608 if(s->avctx->debug&2048){
2609         memset(s->spatial_dwt_buffer, 0, sizeof(DWTELEM)*w*h);
2610         predict_plane(s, s->spatial_dwt_buffer, plane_index, 1);
2611
2612         for(y=0; y<h; y++){
2613             for(x=0; x<w; x++){
2614                 int v= (s->spatial_dwt_buffer[y*w + x]+128)>>8;
2615                 if(v&(~255)) v= ~(v>>31);
2616                 s->mconly_picture.data[plane_index][y*s->mconly_picture.linesize[plane_index] + x]= v;
2617             }
2618         }
2619 }
2620         for(level=0; level<s->spatial_decomposition_count; level++){
2621             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2622                 SubBand *b= &p->band[level][orientation];
2623
2624                 decode_subband(s, b, b->buf, b->parent ? b->parent->buf : NULL, b->stride, orientation);
2625                 if(orientation==0)
2626                     correlate(s, b, b->buf, b->stride, 1, 0);
2627             }
2628         }
2629 if(!(s->avctx->debug&1024))
2630         for(level=0; level<s->spatial_decomposition_count; level++){
2631             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2632                 SubBand *b= &p->band[level][orientation];
2633
2634                 dequantize(s, b, b->buf, b->stride);
2635             }
2636         }
2637
2638 #if QPRED
2639         memset(s->pred_buffer, 0, sizeof(DWTELEM)*w*h);
2640         predict_plane(s, s->pred_buffer, plane_index, 1);
2641         spatial_dwt(s, s->pred_buffer, w, h, w);
2642         for(level=0; level<s->spatial_decomposition_count; level++){
2643             for(orientation=level ? 1 : 0; orientation<4; orientation++){
2644                 SubBand *b= &p->band[level][orientation];
2645                 int delta= ((int)s->pred_buffer - (int)s->spatial_dwt_buffer)/sizeof(DWTELEM);
2646                 
2647                 quantize  (s, b, b->buf + delta, b->stride, s->qbias);
2648                 dequantize(s, b, b->buf + delta, b->stride);
2649             }
2650         }
2651         for(y=0; y<h; y++){
2652             for(x=0; x<w; x++){
2653                 s->spatial_dwt_buffer[y*w + x]+= s->pred_buffer[y*w + x];
2654             }
2655         }
2656         spatial_idwt(s, s->spatial_dwt_buffer, w, h, w);
2657 #else
2658         spatial_idwt(s, s->spatial_dwt_buffer, w, h, w);
2659         predict_plane(s, s->spatial_dwt_buffer, plane_index, 1);
2660 #endif
2661
2662         //FIXME optimize
2663         for(y=0; y<h; y++){
2664             for(x=0; x<w; x++){
2665                 int v= (s->spatial_dwt_buffer[y*w + x]+128)>>8;
2666                 if(v&(~255)) v= ~(v>>31);
2667                 s->current_picture.data[plane_index][y*s->current_picture.linesize[plane_index] + x]= v;
2668             }
2669         }
2670     }
2671             
2672     emms_c();
2673
2674     if(s->last_picture.data[0])
2675         avctx->release_buffer(avctx, &s->last_picture);
2676
2677 if(!(s->avctx->debug&2048))        
2678     *picture= s->current_picture;
2679 else
2680     *picture= s->mconly_picture;
2681     
2682     *data_size = sizeof(AVFrame);
2683     
2684     bytes_read= get_cabac_terminate(c);
2685     if(bytes_read ==0) av_log(s->avctx, AV_LOG_ERROR, "error at end of frame\n");
2686
2687     return bytes_read;
2688 }
2689
2690 static int decode_end(AVCodecContext *avctx)
2691 {
2692     SnowContext *s = avctx->priv_data;
2693
2694     common_end(s);
2695
2696     return 0;
2697 }
2698
2699 AVCodec snow_decoder = {
2700     "snow",
2701     CODEC_TYPE_VIDEO,
2702     CODEC_ID_SNOW,
2703     sizeof(SnowContext),
2704     decode_init,
2705     NULL,
2706     decode_end,
2707     decode_frame,
2708     0 /*CODEC_CAP_DR1*/ /*| CODEC_CAP_DRAW_HORIZ_BAND*/,
2709     NULL
2710 };
2711
2712 AVCodec snow_encoder = {
2713     "snow",
2714     CODEC_TYPE_VIDEO,
2715     CODEC_ID_SNOW,
2716     sizeof(SnowContext),
2717     encode_init,
2718     encode_frame,
2719     encode_end,
2720 };
2721
2722
2723 #if 0
2724 #undef malloc
2725 #undef free
2726 #undef printf
2727
2728 int main(){
2729     int width=256;
2730     int height=256;
2731     int buffer[2][width*height];
2732     SnowContext s;
2733     int i;
2734     s.spatial_decomposition_count=6;
2735     s.spatial_decomposition_type=1;
2736     
2737     printf("testing 5/3 DWT\n");
2738     for(i=0; i<width*height; i++)
2739         buffer[0][i]= buffer[1][i]= random()%54321 - 12345;
2740     
2741     spatial_dwt(&s, buffer[0], width, height, width);
2742     spatial_idwt(&s, buffer[0], width, height, width);
2743     
2744     for(i=0; i<width*height; i++)
2745         if(buffer[0][i]!= buffer[1][i]) printf("fsck: %d %d %d\n",i, buffer[0][i], buffer[1][i]);
2746
2747     printf("testing 9/7 DWT\n");
2748     s.spatial_decomposition_type=0;
2749     for(i=0; i<width*height; i++)
2750         buffer[0][i]= buffer[1][i]= random()%54321 - 12345;
2751     
2752     spatial_dwt(&s, buffer[0], width, height, width);
2753     spatial_idwt(&s, buffer[0], width, height, width);
2754     
2755     for(i=0; i<width*height; i++)
2756         if(buffer[0][i]!= buffer[1][i]) printf("fsck: %d %d %d\n",i, buffer[0][i], buffer[1][i]);
2757         
2758     printf("testing AC coder\n");
2759     memset(s.header_state, 0, sizeof(s.header_state));
2760     ff_init_cabac_encoder(&s.c, buffer[0], 256*256);
2761     ff_init_cabac_states(&s.c, ff_h264_lps_range, ff_h264_mps_state, ff_h264_lps_state, 64);
2762         
2763     for(i=-256; i<256; i++){
2764 START_TIMER
2765         put_symbol(&s.c, s.header_state, i*i*i/3*ABS(i), 1);
2766 STOP_TIMER("put_symbol")
2767     }
2768     put_cabac_terminate(&s.c, 1);
2769
2770     memset(s.header_state, 0, sizeof(s.header_state));
2771     ff_init_cabac_decoder(&s.c, buffer[0], 256*256);
2772     ff_init_cabac_states(&s.c, ff_h264_lps_range, ff_h264_mps_state, ff_h264_lps_state, 64);
2773     
2774     for(i=-256; i<256; i++){
2775         int j;
2776 START_TIMER
2777         j= get_symbol(&s.c, s.header_state, 1);
2778 STOP_TIMER("get_symbol")
2779         if(j!=i*i*i/3*ABS(i)) printf("fsck: %d != %d\n", i, j);
2780     }
2781 {
2782 int level, orientation, x, y;
2783 int64_t errors[8][4];
2784 int64_t g=0;
2785
2786     memset(errors, 0, sizeof(errors));
2787     s.spatial_decomposition_count=3;
2788     s.spatial_decomposition_type=0;
2789     for(level=0; level<s.spatial_decomposition_count; level++){
2790         for(orientation=level ? 1 : 0; orientation<4; orientation++){
2791             int w= width  >> (s.spatial_decomposition_count-level);
2792             int h= height >> (s.spatial_decomposition_count-level);
2793             int stride= width  << (s.spatial_decomposition_count-level);
2794             DWTELEM *buf= buffer[0];
2795             int64_t error=0;
2796
2797             if(orientation&1) buf+=w;
2798             if(orientation>1) buf+=stride>>1;
2799             
2800             memset(buffer[0], 0, sizeof(int)*width*height);
2801             buf[w/2 + h/2*stride]= 256*256;
2802             spatial_idwt(&s, buffer[0], width, height, width);
2803             for(y=0; y<height; y++){
2804                 for(x=0; x<width; x++){
2805                     int64_t d= buffer[0][x + y*width];
2806                     error += d*d;
2807                     if(ABS(width/2-x)<9 && ABS(height/2-y)<9 && level==2) printf("%8lld ", d);
2808                 }
2809                 if(ABS(height/2-y)<9 && level==2) printf("\n");
2810             }
2811             error= (int)(sqrt(error)+0.5);
2812             errors[level][orientation]= error;
2813             if(g) g=ff_gcd(g, error);
2814             else g= error;
2815         }
2816     }
2817     printf("static int const visual_weight[][4]={\n");
2818     for(level=0; level<s.spatial_decomposition_count; level++){
2819         printf("  {");
2820         for(orientation=0; orientation<4; orientation++){
2821             printf("%8lld,", errors[level][orientation]/g);
2822         }
2823         printf("},\n");
2824     }
2825     printf("};\n");
2826     {
2827             int level=2;
2828             int orientation=3;
2829             int w= width  >> (s.spatial_decomposition_count-level);
2830             int h= height >> (s.spatial_decomposition_count-level);
2831             int stride= width  << (s.spatial_decomposition_count-level);
2832             DWTELEM *buf= buffer[0];
2833             int64_t error=0;
2834
2835             buf+=w;
2836             buf+=stride>>1;
2837             
2838             memset(buffer[0], 0, sizeof(int)*width*height);
2839 #if 1
2840             for(y=0; y<height; y++){
2841                 for(x=0; x<width; x++){
2842                     int tab[4]={0,2,3,1};
2843                     buffer[0][x+width*y]= 256*256*tab[(x&1) + 2*(y&1)];
2844                 }
2845             }
2846             spatial_dwt(&s, buffer[0], width, height, width);
2847 #else
2848             for(y=0; y<h; y++){
2849                 for(x=0; x<w; x++){
2850                     buf[x + y*stride  ]=169;
2851                     buf[x + y*stride-w]=64;
2852                 }
2853             }
2854             spatial_idwt(&s, buffer[0], width, height, width);
2855 #endif
2856             for(y=0; y<height; y++){
2857                 for(x=0; x<width; x++){
2858                     int64_t d= buffer[0][x + y*width];
2859                     error += d*d;
2860                     if(ABS(width/2-x)<9 && ABS(height/2-y)<9) printf("%8lld ", d);
2861                 }
2862                 if(ABS(height/2-y)<9) printf("\n");
2863             }
2864     }
2865
2866 }
2867     return 0;
2868 }
2869 #endif
2870