git.sesse.net Git - ffmpeg/blob - libavcodec/snow.c

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