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

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