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

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