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