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