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