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