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