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

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