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