falkTX
/
FFmpeg
mirror of https://github.com/falkTX/FFmpeg.git
1
0
Fork 0
Code Issues 0 Releases 338 Wiki Activity
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
3954 Commits
32 Branches
2.1GB
Tree: b43ff7e2b0
Branches Tags
${ item.name }
Create branch ${ searchTerm }
from 'b43ff7e2b0'
${ noResults }
FFmpeg/libavcodec/snow.c

3163 lines
106KB
Raw Blame History 

  1. /*

  2.  * Copyright (C) 2004 Michael Niedermayer <michaelni@gmx.at>

  3.  *

  4.  * This library is free software; you can redistribute it and/or

  5.  * modify it under the terms of the GNU Lesser General Public

  6.  * License as published by the Free Software Foundation; either

  7.  * version 2 of the License, or (at your option) any later version.

  8.  *

  9.  * This library is distributed in the hope that it will be useful,

  10.  * but WITHOUT ANY WARRANTY; without even the implied warranty of

  11.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

  12.  * Lesser General Public License for more details.

  13.  *

  14.  * You should have received a copy of the GNU Lesser General Public

  15.  * License along with this library; if not, write to the Free Software

  16.  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

  17.  */

  18. 

  19. #include "avcodec.h"

  20. #include "common.h"

  21. #include "dsputil.h"

  22. 

  23. #include "rangecoder.h"

  24. #define MID_STATE 128

  25. 

  26. #include "mpegvideo.h"

  27. 

  28. #undef NDEBUG

  29. #include <assert.h>

  30. 

  31. #define MAX_DECOMPOSITIONS 8

  32. #define MAX_PLANES 4

  33. #define DWTELEM int

  34. #define QROOT 8 

  35. #define LOSSLESS_QLOG -128

  36. #define FRAC_BITS 8

  37. 

  38. static const int8_t quant3[256]={

  39.  0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  40.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  41.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  42.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  43.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  44.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  45.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  46.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  47. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  48. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  49. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  50. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  51. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  52. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  53. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  54. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 0,

  55. };

  56. static const int8_t quant3b[256]={

  57.  0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  58.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  59.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  60.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  61.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  62.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  63.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  64.  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

  65. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  66. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  67. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  68. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  69. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  70. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  71. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  72. -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

  73. };

  74. static const int8_t quant5[256]={

  75.  0, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,

  76.  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,

  77.  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,

  78.  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,

  79.  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,

  80.  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,

  81.  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,

  82.  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,

  83. -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,

  84. -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,

  85. -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,

  86. -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,

  87. -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,

  88. -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,

  89. -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,

  90. -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-1,-1,-1,

  91. };

  92. static const int8_t quant7[256]={

  93.  0, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,

  94.  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,

  95.  2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,

  96.  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,

  97.  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,

  98.  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,

  99.  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,

  100.  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,

  101. -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,

  102. -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,

  103. -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,

  104. -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,

  105. -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,

  106. -3,-3,-3,-3,-3,-3,-3,-3,-3,-2,-2,-2,-2,-2,-2,-2,

  107. -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,

  108. -2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-1,-1,

  109. };

  110. static const int8_t quant9[256]={

  111.  0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3,

  112.  3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,

  113.  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,

  114.  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,

  115.  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,

  116.  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,

  117.  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,

  118.  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,

  119. -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,

  120. -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,

  121. -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,

  122. -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,

  123. -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,

  124. -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,

  125. -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-3,-3,-3,-3,

  126. -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-2,-2,-2,-2,-1,-1,

  127. };

  128. static const int8_t quant11[256]={

  129.  0, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4,

  130.  4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,

  131.  4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,

  132.  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,

  133.  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,

  134.  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,

  135.  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,

  136.  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,

  137. -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,

  138. -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,

  139. -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,

  140. -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,

  141. -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,

  142. -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-4,-4,

  143. -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,

  144. -4,-4,-4,-4,-4,-3,-3,-3,-3,-3,-3,-3,-2,-2,-2,-1,

  145. };

  146. static const int8_t quant13[256]={

  147.  0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4,

  148.  4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,

  149.  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,

  150.  5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,

  151.  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,

  152.  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,

  153.  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,

  154.  6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,

  155. -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,

  156. -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,

  157. -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,

  158. -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,

  159. -6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-6,-5,

  160. -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,

  161. -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,

  162. -4,-4,-4,-4,-4,-4,-4,-4,-4,-3,-3,-3,-3,-2,-2,-1,

  163. };

  164. 

  165. #define LOG2_OBMC_MAX 6

  166. #define OBMC_MAX (1<<(LOG2_OBMC_MAX))

  167. #if 0 //64*cubic

  168. static const uint8_t obmc32[1024]={

  169.  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

  170.  0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0,

  171.  0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 2, 2, 1, 1, 1, 1, 0, 0, 0,

  172.  0, 0, 1, 1, 2, 2, 3, 4, 4, 5, 6, 6, 7, 7, 8, 8, 8, 8, 7, 7, 6, 6, 5, 4, 4, 3, 2, 2, 1, 1, 0, 0,

  173.  0, 0, 1, 2, 2, 3, 4, 6, 7, 8, 9,10,11,12,12,12,12,12,12,11,10, 9, 8, 7, 6, 4, 3, 2, 2, 1, 0, 0,

  174.  0, 1, 1, 2, 3, 5, 6, 8,10,11,13,14,15,16,17,18,18,17,16,15,14,13,11,10, 8, 6, 5, 3, 2, 1, 1, 0,

  175.  0, 1, 1, 3, 4, 6, 8,10,13,15,17,19,20,22,22,23,23,22,22,20,19,17,15,13,10, 8, 6, 4, 3, 1, 1, 0,

  176.  0, 1, 2, 4, 6, 8,10,13,16,19,21,23,25,27,28,29,29,28,27,25,23,21,19,16,13,10, 8, 6, 4, 2, 1, 0,

  177.  0, 1, 2, 4, 7,10,13,16,19,22,25,28,31,33,34,35,35,34,33,31,28,25,22,19,16,13,10, 7, 4, 2, 1, 0,

  178.  0, 1, 3, 5, 8,11,15,19,22,26,30,33,36,38,40,41,41,40,38,36,33,30,26,22,19,15,11, 8, 5, 3, 1, 0,

  179.  0, 1, 3, 6, 9,12,17,21,25,30,34,38,41,44,45,46,46,45,44,41,38,34,30,25,21,17,12, 9, 6, 3, 1, 0,

  180.  0, 1, 3, 6,10,14,19,23,28,33,38,42,45,48,51,52,52,51,48,45,42,38,33,28,23,19,14,10, 6, 3, 1, 0,

  181.  0, 1, 4, 7,11,15,20,25,31,36,41,45,49,52,55,56,56,55,52,49,45,41,36,31,25,20,15,11, 7, 4, 1, 0,

  182.  0, 2, 4, 7,12,16,22,27,33,38,44,48,52,56,58,60,60,58,56,52,48,44,38,33,27,22,16,12, 7, 4, 2, 0,

  183.  0, 1, 4, 8,12,17,22,28,34,40,45,51,55,58,61,62,62,61,58,55,51,45,40,34,28,22,17,12, 8, 4, 1, 0,

  184.  0, 2, 4, 8,12,18,23,29,35,41,46,52,56,60,62,64,64,62,60,56,52,46,41,35,29,23,18,12, 8, 4, 2, 0,

  185.  0, 2, 4, 8,12,18,23,29,35,41,46,52,56,60,62,64,64,62,60,56,52,46,41,35,29,23,18,12, 8, 4, 2, 0,

  186.  0, 1, 4, 8,12,17,22,28,34,40,45,51,55,58,61,62,62,61,58,55,51,45,40,34,28,22,17,12, 8, 4, 1, 0,

  187.  0, 2, 4, 7,12,16,22,27,33,38,44,48,52,56,58,60,60,58,56,52,48,44,38,33,27,22,16,12, 7, 4, 2, 0,

  188.  0, 1, 4, 7,11,15,20,25,31,36,41,45,49,52,55,56,56,55,52,49,45,41,36,31,25,20,15,11, 7, 4, 1, 0,

  189.  0, 1, 3, 6,10,14,19,23,28,33,38,42,45,48,51,52,52,51,48,45,42,38,33,28,23,19,14,10, 6, 3, 1, 0,

  190.  0, 1, 3, 6, 9,12,17,21,25,30,34,38,41,44,45,46,46,45,44,41,38,34,30,25,21,17,12, 9, 6, 3, 1, 0,

  191.  0, 1, 3, 5, 8,11,15,19,22,26,30,33,36,38,40,41,41,40,38,36,33,30,26,22,19,15,11, 8, 5, 3, 1, 0,

  192.  0, 1, 2, 4, 7,10,13,16,19,22,25,28,31,33,34,35,35,34,33,31,28,25,22,19,16,13,10, 7, 4, 2, 1, 0,

  193.  0, 1, 2, 4, 6, 8,10,13,16,19,21,23,25,27,28,29,29,28,27,25,23,21,19,16,13,10, 8, 6, 4, 2, 1, 0,

  194.  0, 1, 1, 3, 4, 6, 8,10,13,15,17,19,20,22,22,23,23,22,22,20,19,17,15,13,10, 8, 6, 4, 3, 1, 1, 0,

  195.  0, 1, 1, 2, 3, 5, 6, 8,10,11,13,14,15,16,17,18,18,17,16,15,14,13,11,10, 8, 6, 5, 3, 2, 1, 1, 0,

  196.  0, 0, 1, 2, 2, 3, 4, 6, 7, 8, 9,10,11,12,12,12,12,12,12,11,10, 9, 8, 7, 6, 4, 3, 2, 2, 1, 0, 0,

  197.  0, 0, 1, 1, 2, 2, 3, 4, 4, 5, 6, 6, 7, 7, 8, 8, 8, 8, 7, 7, 6, 6, 5, 4, 4, 3, 2, 2, 1, 1, 0, 0,

  198.  0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 2, 2, 1, 1, 1, 1, 0, 0, 0,

  199.  0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0,

  200.  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

  201. //error:0.000022

  202. };

  203. static const uint8_t obmc16[256]={

  204.  0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,

  205.  0, 1, 1, 2, 4, 5, 5, 6, 6, 5, 5, 4, 2, 1, 1, 0,

  206.  0, 1, 4, 6, 9,11,13,15,15,13,11, 9, 6, 4, 1, 0,

  207.  0, 2, 6,11,15,20,24,26,26,24,20,15,11, 6, 2, 0,

  208.  0, 4, 9,15,23,29,34,38,38,34,29,23,15, 9, 4, 0,

  209.  0, 5,11,20,29,38,45,49,49,45,38,29,20,11, 5, 0,

  210.  1, 5,13,24,34,45,53,57,57,53,45,34,24,13, 5, 1,

  211.  1, 6,15,26,38,49,57,62,62,57,49,38,26,15, 6, 1,

  212.  1, 6,15,26,38,49,57,62,62,57,49,38,26,15, 6, 1,

  213.  1, 5,13,24,34,45,53,57,57,53,45,34,24,13, 5, 1,

  214.  0, 5,11,20,29,38,45,49,49,45,38,29,20,11, 5, 0,

  215.  0, 4, 9,15,23,29,34,38,38,34,29,23,15, 9, 4, 0,

  216.  0, 2, 6,11,15,20,24,26,26,24,20,15,11, 6, 2, 0,

  217.  0, 1, 4, 6, 9,11,13,15,15,13,11, 9, 6, 4, 1, 0,

  218.  0, 1, 1, 2, 4, 5, 5, 6, 6, 5, 5, 4, 2, 1, 1, 0,

  219.  0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,

  220. //error:0.000033

  221. };

  222. #elif 1 // 64*linear

  223. static const uint8_t obmc32[1024]={

  224.  0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0,

  225.  0, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 2, 2, 2, 1, 1, 1, 0,

  226.  0, 1, 2, 2, 3, 3, 4, 5, 5, 6, 7, 7, 8, 8, 9,10,10, 9, 8, 8, 7, 7, 6, 5, 5, 4, 3, 3, 2, 2, 1, 0,

  227.  0, 1, 2, 3, 4, 5, 6, 7, 7, 8, 9,10,11,12,13,14,14,13,12,11,10, 9, 8, 7, 7, 6, 5, 4, 3, 2, 1, 0,

  228.  1, 2, 3, 4, 5, 6, 7, 8,10,11,12,13,14,15,16,17,17,16,15,14,13,12,11,10, 8, 7, 6, 5, 4, 3, 2, 1,

  229.  1, 2, 3, 5, 6, 8, 9,10,12,13,14,16,17,19,20,21,21,20,19,17,16,14,13,12,10, 9, 8, 6, 5, 3, 2, 1,

  230.  1, 2, 4, 6, 7, 9,11,12,14,15,17,19,20,22,24,25,25,24,22,20,19,17,15,14,12,11, 9, 7, 6, 4, 2, 1,

  231.  1, 3, 5, 7, 8,10,12,14,16,18,20,22,23,25,27,29,29,27,25,23,22,20,18,16,14,12,10, 8, 7, 5, 3, 1,

  232.  1, 3, 5, 7,10,12,14,16,18,20,22,24,27,29,31,33,33,31,29,27,24,22,20,18,16,14,12,10, 7, 5, 3, 1,

  233.  1, 4, 6, 8,11,13,15,18,20,23,25,27,30,32,34,37,37,34,32,30,27,25,23,20,18,15,13,11, 8, 6, 4, 1,

  234.  1, 4, 7, 9,12,14,17,20,22,25,28,30,33,35,38,41,41,38,35,33,30,28,25,22,20,17,14,12, 9, 7, 4, 1,

  235.  1, 4, 7,10,13,16,19,22,24,27,30,33,36,39,42,45,45,42,39,36,33,30,27,24,22,19,16,13,10, 7, 4, 1,

  236.  2, 5, 8,11,14,17,20,23,27,30,33,36,39,42,45,48,48,45,42,39,36,33,30,27,23,20,17,14,11, 8, 5, 2,

  237.  2, 5, 8,12,15,19,22,25,29,32,35,39,42,46,49,52,52,49,46,42,39,35,32,29,25,22,19,15,12, 8, 5, 2,

  238.  2, 5, 9,13,16,20,24,27,31,34,38,42,45,49,53,56,56,53,49,45,42,38,34,31,27,24,20,16,13, 9, 5, 2,

  239.  2, 6,10,14,17,21,25,29,33,37,41,45,48,52,56,60,60,56,52,48,45,41,37,33,29,25,21,17,14,10, 6, 2,

  240.  2, 6,10,14,17,21,25,29,33,37,41,45,48,52,56,60,60,56,52,48,45,41,37,33,29,25,21,17,14,10, 6, 2,

  241.  2, 5, 9,13,16,20,24,27,31,34,38,42,45,49,53,56,56,53,49,45,42,38,34,31,27,24,20,16,13, 9, 5, 2,

  242.  2, 5, 8,12,15,19,22,25,29,32,35,39,42,46,49,52,52,49,46,42,39,35,32,29,25,22,19,15,12, 8, 5, 2,

  243.  2, 5, 8,11,14,17,20,23,27,30,33,36,39,42,45,48,48,45,42,39,36,33,30,27,23,20,17,14,11, 8, 5, 2,

  244.  1, 4, 7,10,13,16,19,22,24,27,30,33,36,39,42,45,45,42,39,36,33,30,27,24,22,19,16,13,10, 7, 4, 1,

  245.  1, 4, 7, 9,12,14,17,20,22,25,28,30,33,35,38,41,41,38,35,33,30,28,25,22,20,17,14,12, 9, 7, 4, 1,

  246.  1, 4, 6, 8,11,13,15,18,20,23,25,27,30,32,34,37,37,34,32,30,27,25,23,20,18,15,13,11, 8, 6, 4, 1,

  247.  1, 3, 5, 7,10,12,14,16,18,20,22,24,27,29,31,33,33,31,29,27,24,22,20,18,16,14,12,10, 7, 5, 3, 1,

  248.  1, 3, 5, 7, 8,10,12,14,16,18,20,22,23,25,27,29,29,27,25,23,22,20,18,16,14,12,10, 8, 7, 5, 3, 1,

  249.  1, 2, 4, 6, 7, 9,11,12,14,15,17,19,20,22,24,25,25,24,22,20,19,17,15,14,12,11, 9, 7, 6, 4, 2, 1,

  250.  1, 2, 3, 5, 6, 8, 9,10,12,13,14,16,17,19,20,21,21,20,19,17,16,14,13,12,10, 9, 8, 6, 5, 3, 2, 1,

  251.  1, 2, 3, 4, 5, 6, 7, 8,10,11,12,13,14,15,16,17,17,16,15,14,13,12,11,10, 8, 7, 6, 5, 4, 3, 2, 1,

  252.  0, 1, 2, 3, 4, 5, 6, 7, 7, 8, 9,10,11,12,13,14,14,13,12,11,10, 9, 8, 7, 7, 6, 5, 4, 3, 2, 1, 0,

  253.  0, 1, 2, 2, 3, 3, 4, 5, 5, 6, 7, 7, 8, 8, 9,10,10, 9, 8, 8, 7, 7, 6, 5, 5, 4, 3, 3, 2, 2, 1, 0,

  254.  0, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 2, 2, 2, 1, 1, 1, 0,

  255.  0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0,

  256.  //error:0.000020

  257. };

  258. static const uint8_t obmc16[256]={

  259.  0, 1, 1, 2, 2, 3, 3, 4, 4, 3, 3, 2, 2, 1, 1, 0,

  260.  1, 2, 4, 5, 7, 8,10,11,11,10, 8, 7, 5, 4, 2, 1,

  261.  1, 4, 6, 9,11,14,16,19,19,16,14,11, 9, 6, 4, 1,

  262.  2, 5, 9,12,16,19,23,26,26,23,19,16,12, 9, 5, 2,

  263.  2, 7,11,16,20,25,29,34,34,29,25,20,16,11, 7, 2,

  264.  3, 8,14,19,25,30,36,41,41,36,30,25,19,14, 8, 3,

  265.  3,10,16,23,29,36,42,49,49,42,36,29,23,16,10, 3,

  266.  4,11,19,26,34,41,49,56,56,49,41,34,26,19,11, 4,

  267.  4,11,19,26,34,41,49,56,56,49,41,34,26,19,11, 4,

  268.  3,10,16,23,29,36,42,49,49,42,36,29,23,16,10, 3,

  269.  3, 8,14,19,25,30,36,41,41,36,30,25,19,14, 8, 3,

  270.  2, 7,11,16,20,25,29,34,34,29,25,20,16,11, 7, 2,

  271.  2, 5, 9,12,16,19,23,26,26,23,19,16,12, 9, 5, 2,

  272.  1, 4, 6, 9,11,14,16,19,19,16,14,11, 9, 6, 4, 1,

  273.  1, 2, 4, 5, 7, 8,10,11,11,10, 8, 7, 5, 4, 2, 1,

  274.  0, 1, 1, 2, 2, 3, 3, 4, 4, 3, 3, 2, 2, 1, 1, 0,

  275. //error:0.000015

  276. };

  277. #else //64*cos

  278. static const uint8_t obmc32[1024]={

  279.  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

  280.  0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,

  281.  0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 1, 1, 1, 1, 0, 0, 0,

  282.  0, 0, 1, 1, 1, 2, 2, 3, 4, 5, 5, 6, 7, 7, 7, 7, 7, 7, 7, 7, 6, 5, 5, 4, 3, 2, 2, 1, 1, 1, 0, 0,

  283.  0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 9,10,11,11,12,12,12,12,11,11,10, 9, 7, 6, 5, 4, 3, 2, 1, 1, 0, 0,

  284.  0, 0, 1, 2, 3, 5, 6, 8, 9,11,12,14,15,16,17,17,17,17,16,15,14,12,11, 9, 8, 6, 5, 3, 2, 1, 0, 0,

  285.  0, 1, 1, 2, 4, 6, 8,10,12,15,17,19,20,21,22,23,23,22,21,20,19,17,15,12,10, 8, 6, 4, 2, 1, 1, 0,

  286.  0, 1, 2, 3, 5, 8,10,13,16,19,21,24,26,27,28,29,29,28,27,26,24,21,19,16,13,10, 8, 5, 3, 2, 1, 0,

  287.  0, 1, 2, 4, 6, 9,12,16,19,23,26,29,31,33,34,35,35,34,33,31,29,26,23,19,16,12, 9, 6, 4, 2, 1, 0,

  288.  0, 1, 3, 5, 7,11,15,19,23,26,30,34,37,39,40,41,41,40,39,37,34,30,26,23,19,15,11, 7, 5, 3, 1, 0,

  289.  0, 1, 3, 5, 9,12,17,21,26,30,35,38,42,44,46,47,47,46,44,42,38,35,30,26,21,17,12, 9, 5, 3, 1, 0,

  290.  0, 1, 3, 6, 9,14,19,24,29,34,38,43,46,49,51,52,52,51,49,46,43,38,34,29,24,19,14, 9, 6, 3, 1, 0,

  291.  0, 1, 3, 6,11,15,20,26,31,37,42,46,50,53,56,57,57,56,53,50,46,42,37,31,26,20,15,11, 6, 3, 1, 0,

  292.  0, 1, 3, 7,11,16,21,27,33,39,44,49,53,57,59,60,60,59,57,53,49,44,39,33,27,21,16,11, 7, 3, 1, 0,

  293.  0, 1, 4, 7,12,17,22,28,34,40,46,51,56,59,61,63,63,61,59,56,51,46,40,34,28,22,17,12, 7, 4, 1, 0,

  294.  0, 1, 4, 7,12,17,23,29,35,41,47,52,57,60,63,64,64,63,60,57,52,47,41,35,29,23,17,12, 7, 4, 1, 0,

  295.  0, 1, 4, 7,12,17,23,29,35,41,47,52,57,60,63,64,64,63,60,57,52,47,41,35,29,23,17,12, 7, 4, 1, 0,

  296.  0, 1, 4, 7,12,17,22,28,34,40,46,51,56,59,61,63,63,61,59,56,51,46,40,34,28,22,17,12, 7, 4, 1, 0,

  297.  0, 1, 3, 7,11,16,21,27,33,39,44,49,53,57,59,60,60,59,57,53,49,44,39,33,27,21,16,11, 7, 3, 1, 0,

  298.  0, 1, 3, 6,11,15,20,26,31,37,42,46,50,53,56,57,57,56,53,50,46,42,37,31,26,20,15,11, 6, 3, 1, 0,

  299.  0, 1, 3, 6, 9,14,19,24,29,34,38,43,46,49,51,52,52,51,49,46,43,38,34,29,24,19,14, 9, 6, 3, 1, 0,

  300.  0, 1, 3, 5, 9,12,17,21,26,30,35,38,42,44,46,47,47,46,44,42,38,35,30,26,21,17,12, 9, 5, 3, 1, 0,

  301.  0, 1, 3, 5, 7,11,15,19,23,26,30,34,37,39,40,41,41,40,39,37,34,30,26,23,19,15,11, 7, 5, 3, 1, 0,

  302.  0, 1, 2, 4, 6, 9,12,16,19,23,26,29,31,33,34,35,35,34,33,31,29,26,23,19,16,12, 9, 6, 4, 2, 1, 0,

  303.  0, 1, 2, 3, 5, 8,10,13,16,19,21,24,26,27,28,29,29,28,27,26,24,21,19,16,13,10, 8, 5, 3, 2, 1, 0,

  304.  0, 1, 1, 2, 4, 6, 8,10,12,15,17,19,20,21,22,23,23,22,21,20,19,17,15,12,10, 8, 6, 4, 2, 1, 1, 0,

  305.  0, 0, 1, 2, 3, 5, 6, 8, 9,11,12,14,15,16,17,17,17,17,16,15,14,12,11, 9, 8, 6, 5, 3, 2, 1, 0, 0,

  306.  0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 9,10,11,11,12,12,12,12,11,11,10, 9, 7, 6, 5, 4, 3, 2, 1, 1, 0, 0,

  307.  0, 0, 1, 1, 1, 2, 2, 3, 4, 5, 5, 6, 7, 7, 7, 7, 7, 7, 7, 7, 6, 5, 5, 4, 3, 2, 2, 1, 1, 1, 0, 0,

  308.  0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 1, 1, 1, 1, 0, 0, 0,

  309.  0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,

  310.  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

  311. //error:0.000022

  312. };

  313. static const uint8_t obmc16[256]={

  314.  0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0,

  315.  0, 0, 1, 2, 3, 4, 5, 5, 5, 5, 4, 3, 2, 1, 0, 0,

  316.  0, 1, 3, 6, 8,11,13,14,14,13,11, 8, 6, 3, 1, 0,

  317.  0, 2, 6,10,15,20,24,26,26,24,20,15,10, 6, 2, 0,

  318.  0, 3, 8,16,23,30,35,38,38,35,30,23,16, 8, 3, 0,

  319.  1, 4,11,20,30,39,46,49,49,46,39,30,20,11, 4, 1,

  320.  1, 5,13,24,35,46,54,58,58,54,46,35,24,13, 5, 1,

  321.  0, 5,14,26,38,49,58,63,63,58,49,38,26,14, 5, 0,

  322.  0, 5,14,26,38,49,58,63,63,58,49,38,26,14, 5, 0,

  323.  1, 5,13,24,35,46,54,58,58,54,46,35,24,13, 5, 1,

  324.  1, 4,11,20,30,39,46,49,49,46,39,30,20,11, 4, 1,

  325.  0, 3, 8,16,23,30,35,38,38,35,30,23,16, 8, 3, 0,

  326.  0, 2, 6,10,15,20,24,26,26,24,20,15,10, 6, 2, 0,

  327.  0, 1, 3, 6, 8,11,13,14,14,13,11, 8, 6, 3, 1, 0,

  328.  0, 0, 1, 2, 3, 4, 5, 5, 5, 5, 4, 3, 2, 1, 0, 0,

  329.  0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0,

  330. //error:0.000022

  331. };

  332. #endif

  333. 

  334. //linear *64

  335. static const uint8_t obmc8[64]={

  336.  1, 3, 5, 7, 7, 5, 3, 1,

  337.  3, 9,15,21,21,15, 9, 3,

  338.  5,15,25,35,35,25,15, 5,

  339.  7,21,35,49,49,35,21, 7,

  340.  7,21,35,49,49,35,21, 7,

  341.  5,15,25,35,35,25,15, 5,

  342.  3, 9,15,21,21,15, 9, 3,

  343.  1, 3, 5, 7, 7, 5, 3, 1,

  344. //error:0.000000

  345. };

  346. 

  347. //linear *64

  348. static const uint8_t obmc4[16]={

  349.  4,12,12, 4,

  350. 12,36,36,12,

  351. 12,36,36,12,

  352.  4,12,12, 4,

  353. //error:0.000000

  354. };

  355. 

  356. static const uint8_t *obmc_tab[4]={

  357.     obmc32, obmc16, obmc8, obmc4

  358. };

  359. 

  360. typedef struct BlockNode{

  361.     int16_t mx;

  362.     int16_t my;

  363.     uint8_t color[3];

  364.     uint8_t type;

  365. //#define TYPE_SPLIT    1

  366. #define BLOCK_INTRA   1

  367. //#define TYPE_NOCOLOR  4

  368.     uint8_t level; //FIXME merge into type?

  369. }BlockNode;

  370. 

  371. #define LOG2_MB_SIZE 4

  372. #define MB_SIZE (1<<LOG2_MB_SIZE)

  373. 

  374. typedef struct SubBand{

  375.     int level;

  376.     int stride;

  377.     int width;

  378.     int height;

  379.     int qlog;                                   ///< log(qscale)/log[2^(1/6)]

  380.     DWTELEM *buf;

  381.     int16_t *x;

  382.     DWTELEM *coeff;

  383.     struct SubBand *parent;

  384.     uint8_t state[/*7*2*/ 7 + 512][32];

  385. }SubBand;

  386. 

  387. typedef struct Plane{

  388.     int width;

  389.     int height;

  390.     SubBand band[MAX_DECOMPOSITIONS][4];

  391. }Plane;

  392. 

  393. typedef struct SnowContext{

  394. //    MpegEncContext m; // needed for motion estimation, should not be used for anything else, the idea is to make the motion estimation eventually independant of MpegEncContext, so this will be removed then (FIXME/XXX)

  395. 

  396.     AVCodecContext *avctx;

  397.     RangeCoder c;

  398.     DSPContext dsp;

  399.     AVFrame input_picture;

  400.     AVFrame current_picture;

  401.     AVFrame last_picture;

  402.     AVFrame mconly_picture;

  403. //     uint8_t q_context[16];

  404.     uint8_t header_state[32];

  405.     uint8_t block_state[128 + 32*128];

  406.     int keyframe;

  407.     int always_reset;

  408.     int version;

  409.     int spatial_decomposition_type;

  410.     int temporal_decomposition_type;

  411.     int spatial_decomposition_count;

  412.     int temporal_decomposition_count;

  413.     DWTELEM *spatial_dwt_buffer;

  414.     int colorspace_type;

  415.     int chroma_h_shift;

  416.     int chroma_v_shift;

  417.     int spatial_scalability;

  418.     int qlog;

  419.     int lambda;

  420.     int lambda2;

  421.     int mv_scale;

  422.     int qbias;

  423. #define QBIAS_SHIFT 3

  424.     int b_width;

  425.     int b_height;

  426.     int block_max_depth;

  427.     Plane plane[MAX_PLANES];

  428.     BlockNode *block;

  429. 

  430.     MpegEncContext m; // needed for motion estimation, should not be used for anything else, the idea is to make the motion estimation eventually independant of MpegEncContext, so this will be removed then (FIXME/XXX)

  431. }SnowContext;

  432. 

  433. #define QEXPSHIFT (7-FRAC_BITS+8) //FIXME try to change this to 0

  434. static const uint8_t qexp[8]={

  435.     128, 140, 152, 166, 181, 197, 215, 235

  436. //   64,  70,  76,  83,  91,  99, 108, 117

  437. //   32,  35,  38,  41,  45,  49,  54,  59

  438. //   16,  17,  19,  21,  23,  25,  27,  29

  439. //    8,   9,  10,  10,  11,  12,  13,  15

  440. };

  441. 

  442. static inline int mirror(int v, int m){

  443.     if     (v<0) return -v;

  444.     else if(v>m) return 2*m-v;

  445.     else         return v;

  446. }

  447. 

  448. static inline void put_symbol(RangeCoder *c, uint8_t *state, int v, int is_signed){

  449.     int i;

  450. 

  451.     if(v){

  452.         const int a= ABS(v);

  453.         const int e= av_log2(a);

  454. #if 1

  455.         const int el= FFMIN(e, 10);   

  456.         put_rac(c, state+0, 0);

  457. 

  458.         for(i=0; i<el; i++){

  459.             put_rac(c, state+1+i, 1);  //1..10

  460.         }

  461.         for(; i<e; i++){

  462.             put_rac(c, state+1+9, 1);  //1..10

  463.         }

  464.         put_rac(c, state+1+FFMIN(i,9), 0);

  465. 

  466.         for(i=e-1; i>=el; i--){

  467.             put_rac(c, state+22+9, (a>>i)&1); //22..31

  468.         }

  469.         for(; i>=0; i--){

  470.             put_rac(c, state+22+i, (a>>i)&1); //22..31

  471.         }

  472. 

  473.         if(is_signed)

  474.             put_rac(c, state+11 + el, v < 0); //11..21

  475. #else

  476.         

  477.         put_rac(c, state+0, 0);

  478.         if(e<=9){

  479.             for(i=0; i<e; i++){

  480.                 put_rac(c, state+1+i, 1);  //1..10

  481.             }

  482.             put_rac(c, state+1+i, 0);

  483. 

  484.             for(i=e-1; i>=0; i--){

  485.                 put_rac(c, state+22+i, (a>>i)&1); //22..31

  486.             }

  487. 

  488.             if(is_signed)

  489.                 put_rac(c, state+11 + e, v < 0); //11..21

  490.         }else{

  491.             for(i=0; i<e; i++){

  492.                 put_rac(c, state+1+FFMIN(i,9), 1);  //1..10

  493.             }

  494.             put_rac(c, state+1+FFMIN(i,9), 0);

  495. 

  496.             for(i=e-1; i>=0; i--){

  497.                 put_rac(c, state+22+FFMIN(i,9), (a>>i)&1); //22..31

  498.             }

  499. 

  500.             if(is_signed)

  501.                 put_rac(c, state+11 + FFMIN(e,10), v < 0); //11..21

  502.         }

  503. #endif

  504.     }else{

  505.         put_rac(c, state+0, 1);

  506.     }

  507. }

  508. 

  509. static inline int get_symbol(RangeCoder *c, uint8_t *state, int is_signed){

  510.     if(get_rac(c, state+0))

  511.         return 0;

  512.     else{

  513.         int i, e, a;

  514.         e= 0;

  515.         while(get_rac(c, state+1 + FFMIN(e,9))){ //1..10

  516.             e++;

  517.         }

  518. 

  519.         a= 1;

  520.         for(i=e-1; i>=0; i--){

  521.             a += a + get_rac(c, state+22 + FFMIN(i,9)); //22..31

  522.         }

  523. 

  524.         if(is_signed && get_rac(c, state+11 + FFMIN(e,10))) //11..21

  525.             return -a;

  526.         else

  527.             return a;

  528.     }

  529. }

  530. 

  531. static inline void put_symbol2(RangeCoder *c, uint8_t *state, int v, int log2){

  532.     int i;

  533.     int r= log2>=0 ? 1<<log2 : 1;

  534. 

  535.     assert(v>=0);

  536.     assert(log2>=-4);

  537. 

  538.     while(v >= r){

  539.         put_rac(c, state+4+log2, 1);

  540.         v -= r;

  541.         log2++;

  542.         if(log2>0) r+=r;

  543.     }

  544.     put_rac(c, state+4+log2, 0);

  545.     

  546.     for(i=log2-1; i>=0; i--){

  547.         put_rac(c, state+31-i, (v>>i)&1);

  548.     }

  549. }

  550. 

  551. static inline int get_symbol2(RangeCoder *c, uint8_t *state, int log2){

  552.     int i;

  553.     int r= log2>=0 ? 1<<log2 : 1;

  554.     int v=0;

  555. 

  556.     assert(log2>=-4);

  557. 

  558.     while(get_rac(c, state+4+log2)){

  559.         v+= r;

  560.         log2++;

  561.         if(log2>0) r+=r;

  562.     }

  563.     

  564.     for(i=log2-1; i>=0; i--){

  565.         v+= get_rac(c, state+31-i)<<i;

  566.     }

  567. 

  568.     return v;

  569. }

  570. 

  571. static always_inline void lift(DWTELEM *dst, DWTELEM *src, DWTELEM *ref, int dst_step, int src_step, int ref_step, int width, int mul, int add, int shift, int highpass, int inverse){

  572.     const int mirror_left= !highpass;

  573.     const int mirror_right= (width&1) ^ highpass;

  574.     const int w= (width>>1) - 1 + (highpass & width);

  575.     int i;

  576. 

  577. #define LIFT(src, ref, inv) ((src) + ((inv) ? - (ref) : + (ref)))

  578.     if(mirror_left){

  579.         dst[0] = LIFT(src[0], ((mul*2*ref[0]+add)>>shift), inverse);

  580.         dst += dst_step;

  581.         src += src_step;

  582.     }

  583.     

  584.     for(i=0; i<w; i++){

  585.         dst[i*dst_step] = LIFT(src[i*src_step], ((mul*(ref[i*ref_step] + ref[(i+1)*ref_step])+add)>>shift), inverse);

  586.     }

  587.     

  588.     if(mirror_right){

  589.         dst[w*dst_step] = LIFT(src[w*src_step], ((mul*2*ref[w*ref_step]+add)>>shift), inverse);

  590.     }

  591. }

  592. 

  593. static always_inline void lift5(DWTELEM *dst, DWTELEM *src, DWTELEM *ref, int dst_step, int src_step, int ref_step, int width, int mul, int add, int shift, int highpass, int inverse){

  594.     const int mirror_left= !highpass;

  595.     const int mirror_right= (width&1) ^ highpass;

  596.     const int w= (width>>1) - 1 + (highpass & width);

  597.     int i;

  598. 

  599.     if(mirror_left){

  600.         int r= 3*2*ref[0];

  601.         r += r>>4;

  602.         r += r>>8;

  603.         dst[0] = LIFT(src[0], ((r+add)>>shift), inverse);

  604.         dst += dst_step;

  605.         src += src_step;

  606.     }

  607.     

  608.     for(i=0; i<w; i++){

  609.         int r= 3*(ref[i*ref_step] + ref[(i+1)*ref_step]);

  610.         r += r>>4;

  611.         r += r>>8;

  612.         dst[i*dst_step] = LIFT(src[i*src_step], ((r+add)>>shift), inverse);

  613.     }

  614.     

  615.     if(mirror_right){

  616.         int r= 3*2*ref[w*ref_step];

  617.         r += r>>4;

  618.         r += r>>8;

  619.         dst[w*dst_step] = LIFT(src[w*src_step], ((r+add)>>shift), inverse);

  620.     }

  621. }

  622. 

  623. 

  624. static void inplace_lift(DWTELEM *dst, int width, int *coeffs, int n, int shift, int start, int inverse){

  625.     int x, i;

  626.     

  627.     for(x=start; x<width; x+=2){

  628.         int64_t sum=0;

  629. 

  630.         for(i=0; i<n; i++){

  631.             int x2= x + 2*i - n + 1;

  632.             if     (x2<     0) x2= -x2;

  633.             else if(x2>=width) x2= 2*width-x2-2;

  634.             sum += coeffs[i]*(int64_t)dst[x2];

  635.         }

  636.         if(inverse) dst[x] -= (sum + (1<<shift)/2)>>shift;

  637.         else        dst[x] += (sum + (1<<shift)/2)>>shift;

  638.     }

  639. }

  640. 

  641. static void inplace_liftV(DWTELEM *dst, int width, int height, int stride, int *coeffs, int n, int shift, int start, int inverse){

  642.     int x, y, i;

  643.     for(y=start; y<height; y+=2){

  644.         for(x=0; x<width; x++){

  645.             int64_t sum=0;

  646.     

  647.             for(i=0; i<n; i++){

  648.                 int y2= y + 2*i - n + 1;

  649.                 if     (y2<      0) y2= -y2;

  650.                 else if(y2>=height) y2= 2*height-y2-2;

  651.                 sum += coeffs[i]*(int64_t)dst[x + y2*stride];

  652.             }

  653.             if(inverse) dst[x + y*stride] -= (sum + (1<<shift)/2)>>shift;

  654.             else        dst[x + y*stride] += (sum + (1<<shift)/2)>>shift;

  655.         }

  656.     }

  657. }

  658. 

  659. #define SCALEX 1

  660. #define LX0 0

  661. #define LX1 1

  662. 

  663. #if 0 // more accurate 9/7

  664. #define N1 2

  665. #define SHIFT1 14

  666. #define COEFFS1 (int[]){-25987,-25987}

  667. #define N2 2

  668. #define SHIFT2 19

  669. #define COEFFS2 (int[]){-27777,-27777}

  670. #define N3 2

  671. #define SHIFT3 15

  672. #define COEFFS3 (int[]){28931,28931}

  673. #define N4 2

  674. #define SHIFT4 15

  675. #define COEFFS4 (int[]){14533,14533}

  676. #elif 1 // 13/7 CRF

  677. #define N1 4

  678. #define SHIFT1 4

  679. #define COEFFS1 (int[]){1,-9,-9,1}

  680. #define N2 4

  681. #define SHIFT2 4

  682. #define COEFFS2 (int[]){-1,5,5,-1}

  683. #define N3 0

  684. #define SHIFT3 1

  685. #define COEFFS3 NULL

  686. #define N4 0

  687. #define SHIFT4 1

  688. #define COEFFS4 NULL

  689. #elif 1 // 3/5

  690. #define LX0 1

  691. #define LX1 0

  692. #define SCALEX 0.5

  693. #define N1 2

  694. #define SHIFT1 1

  695. #define COEFFS1 (int[]){1,1}

  696. #define N2 2

  697. #define SHIFT2 2

  698. #define COEFFS2 (int[]){-1,-1}

  699. #define N3 0

  700. #define SHIFT3 0

  701. #define COEFFS3 NULL

  702. #define N4 0

  703. #define SHIFT4 0

  704. #define COEFFS4 NULL

  705. #elif 1 // 11/5 

  706. #define N1 0

  707. #define SHIFT1 1

  708. #define COEFFS1 NULL

  709. #define N2 2

  710. #define SHIFT2 2

  711. #define COEFFS2 (int[]){-1,-1}

  712. #define N3 2

  713. #define SHIFT3 0

  714. #define COEFFS3 (int[]){-1,-1}

  715. #define N4 4

  716. #define SHIFT4 7

  717. #define COEFFS4 (int[]){-5,29,29,-5}

  718. #define SCALEX 4

  719. #elif 1 // 9/7 CDF

  720. #define N1 2

  721. #define SHIFT1 7

  722. #define COEFFS1 (int[]){-203,-203}

  723. #define N2 2

  724. #define SHIFT2 12

  725. #define COEFFS2 (int[]){-217,-217}

  726. #define N3 2

  727. #define SHIFT3 7

  728. #define COEFFS3 (int[]){113,113}

  729. #define N4 2

  730. #define SHIFT4 9

  731. #define COEFFS4 (int[]){227,227}

  732. #define SCALEX 1

  733. #elif 1 // 7/5 CDF

  734. #define N1 0

  735. #define SHIFT1 1

  736. #define COEFFS1 NULL

  737. #define N2 2

  738. #define SHIFT2 2

  739. #define COEFFS2 (int[]){-1,-1}

  740. #define N3 2

  741. #define SHIFT3 0

  742. #define COEFFS3 (int[]){-1,-1}

  743. #define N4 2

  744. #define SHIFT4 4

  745. #define COEFFS4 (int[]){3,3}

  746. #elif 1 // 9/7 MN

  747. #define N1 4

  748. #define SHIFT1 4

  749. #define COEFFS1 (int[]){1,-9,-9,1}

  750. #define N2 2

  751. #define SHIFT2 2

  752. #define COEFFS2 (int[]){1,1}

  753. #define N3 0

  754. #define SHIFT3 1

  755. #define COEFFS3 NULL

  756. #define N4 0

  757. #define SHIFT4 1

  758. #define COEFFS4 NULL

  759. #else // 13/7 CRF

  760. #define N1 4

  761. #define SHIFT1 4

  762. #define COEFFS1 (int[]){1,-9,-9,1}

  763. #define N2 4

  764. #define SHIFT2 4

  765. #define COEFFS2 (int[]){-1,5,5,-1}

  766. #define N3 0

  767. #define SHIFT3 1

  768. #define COEFFS3 NULL

  769. #define N4 0

  770. #define SHIFT4 1

  771. #define COEFFS4 NULL

  772. #endif

  773. static void horizontal_decomposeX(DWTELEM *b, int width){

  774.     DWTELEM temp[width];

  775.     const int width2= width>>1;

  776.     const int w2= (width+1)>>1;

  777.     int A1,A2,A3,A4, x;

  778. 

  779.     inplace_lift(b, width, COEFFS1, N1, SHIFT1, LX1, 0);

  780.     inplace_lift(b, width, COEFFS2, N2, SHIFT2, LX0, 0);

  781.     inplace_lift(b, width, COEFFS3, N3, SHIFT3, LX1, 0);

  782.     inplace_lift(b, width, COEFFS4, N4, SHIFT4, LX0, 0);

  783.     

  784.     for(x=0; x<width2; x++){

  785.         temp[x   ]= b[2*x    ];

  786.         temp[x+w2]= b[2*x + 1];

  787.     }

  788.     if(width&1)

  789.         temp[x   ]= b[2*x    ];

  790.     memcpy(b, temp, width*sizeof(int));

  791. }

  792. 

  793. static void horizontal_composeX(DWTELEM *b, int width){

  794.     DWTELEM temp[width];

  795.     const int width2= width>>1;

  796.     int A1,A2,A3,A4, x;

  797.     const int w2= (width+1)>>1;

  798. 

  799.     memcpy(temp, b, width*sizeof(int));

  800.     for(x=0; x<width2; x++){

  801.         b[2*x    ]= temp[x   ];

  802.         b[2*x + 1]= temp[x+w2];

  803.     }

  804.     if(width&1)

  805.         b[2*x    ]= temp[x   ];

  806. 

  807.     inplace_lift(b, width, COEFFS4, N4, SHIFT4, LX0, 1);

  808.     inplace_lift(b, width, COEFFS3, N3, SHIFT3, LX1, 1);

  809.     inplace_lift(b, width, COEFFS2, N2, SHIFT2, LX0, 1);

  810.     inplace_lift(b, width, COEFFS1, N1, SHIFT1, LX1, 1);

  811. }

  812. 

  813. static void spatial_decomposeX(DWTELEM *buffer, int width, int height, int stride){

  814.     int x, y;

  815.   

  816.     for(y=0; y<height; y++){

  817.         for(x=0; x<width; x++){

  818.             buffer[y*stride + x] *= SCALEX;

  819.         }

  820.     }

  821. 

  822.     for(y=0; y<height; y++){

  823.         horizontal_decomposeX(buffer + y*stride, width);

  824.     }

  825.     

  826.     inplace_liftV(buffer, width, height, stride, COEFFS1, N1, SHIFT1, LX1, 0);

  827.     inplace_liftV(buffer, width, height, stride, COEFFS2, N2, SHIFT2, LX0, 0);

  828.     inplace_liftV(buffer, width, height, stride, COEFFS3, N3, SHIFT3, LX1, 0);

  829.     inplace_liftV(buffer, width, height, stride, COEFFS4, N4, SHIFT4, LX0, 0);    

  830. }

  831. 

  832. static void spatial_composeX(DWTELEM *buffer, int width, int height, int stride){

  833.     int x, y;

  834.   

  835.     inplace_liftV(buffer, width, height, stride, COEFFS4, N4, SHIFT4, LX0, 1);

  836.     inplace_liftV(buffer, width, height, stride, COEFFS3, N3, SHIFT3, LX1, 1);

  837.     inplace_liftV(buffer, width, height, stride, COEFFS2, N2, SHIFT2, LX0, 1);

  838.     inplace_liftV(buffer, width, height, stride, COEFFS1, N1, SHIFT1, LX1, 1);

  839. 

  840.     for(y=0; y<height; y++){

  841.         horizontal_composeX(buffer + y*stride, width);

  842.     }

  843. 

  844.     for(y=0; y<height; y++){

  845.         for(x=0; x<width; x++){

  846.             buffer[y*stride + x] /= SCALEX;

  847.         }

  848.     }

  849. }

  850. 

  851. static void horizontal_decompose53i(DWTELEM *b, int width){

  852.     DWTELEM temp[width];

  853.     const int width2= width>>1;

  854.     int A1,A2,A3,A4, x;

  855.     const int w2= (width+1)>>1;

  856. 

  857.     for(x=0; x<width2; x++){

  858.         temp[x   ]= b[2*x    ];

  859.         temp[x+w2]= b[2*x + 1];

  860.     }

  861.     if(width&1)

  862.         temp[x   ]= b[2*x    ];

  863. #if 0

  864.     A2= temp[1       ];

  865.     A4= temp[0       ];

  866.     A1= temp[0+width2];

  867.     A1 -= (A2 + A4)>>1;

  868.     A4 += (A1 + 1)>>1;

  869.     b[0+width2] = A1;

  870.     b[0       ] = A4;

  871.     for(x=1; x+1<width2; x+=2){

  872.         A3= temp[x+width2];

  873.         A4= temp[x+1     ];

  874.         A3 -= (A2 + A4)>>1;

  875.         A2 += (A1 + A3 + 2)>>2;

  876.         b[x+width2] = A3;

  877.         b[x       ] = A2;

  878. 

  879.         A1= temp[x+1+width2];

  880.         A2= temp[x+2       ];

  881.         A1 -= (A2 + A4)>>1;

  882.         A4 += (A1 + A3 + 2)>>2;

  883.         b[x+1+width2] = A1;

  884.         b[x+1       ] = A4;

  885.     }

  886.     A3= temp[width-1];

  887.     A3 -= A2;

  888.     A2 += (A1 + A3 + 2)>>2;

  889.     b[width -1] = A3;

  890.     b[width2-1] = A2;

  891. #else        

  892.     lift(b+w2, temp+w2, temp, 1, 1, 1, width, -1, 0, 1, 1, 0);

  893.     lift(b   , temp   , b+w2, 1, 1, 1, width,  1, 2, 2, 0, 0);

  894. #endif

  895. }

  896. 

  897. static void vertical_decompose53iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){

  898.     int i;

  899.     

  900.     for(i=0; i<width; i++){

  901.         b1[i] -= (b0[i] + b2[i])>>1;

  902.     }

  903. }

  904. 

  905. static void vertical_decompose53iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){

  906.     int i;

  907.     

  908.     for(i=0; i<width; i++){

  909.         b1[i] += (b0[i] + b2[i] + 2)>>2;

  910.     }

  911. }

  912. 

  913. static void spatial_decompose53i(DWTELEM *buffer, int width, int height, int stride){

  914.     int y;

  915.     DWTELEM *b0= buffer + mirror(-2-1, height-1)*stride;

  916.     DWTELEM *b1= buffer + mirror(-2  , height-1)*stride;

  917.   

  918.     for(y=-2; y<height; y+=2){

  919.         DWTELEM *b2= buffer + mirror(y+1, height-1)*stride;

  920.         DWTELEM *b3= buffer + mirror(y+2, height-1)*stride;

  921. 

  922. {START_TIMER

  923.         if(b1 <= b3)     horizontal_decompose53i(b2, width);

  924.         if(y+2 < height) horizontal_decompose53i(b3, width);

  925. STOP_TIMER("horizontal_decompose53i")}

  926.         

  927. {START_TIMER

  928.         if(b1 <= b3) vertical_decompose53iH0(b1, b2, b3, width);

  929.         if(b0 <= b2) vertical_decompose53iL0(b0, b1, b2, width);

  930. STOP_TIMER("vertical_decompose53i*")}

  931.         

  932.         b0=b2;

  933.         b1=b3;

  934.     }

  935. }

  936. 

  937. #define lift5 lift

  938. #if 1

  939. #define W_AM 3

  940. #define W_AO 0

  941. #define W_AS 1

  942. 

  943. #define W_BM 1

  944. #define W_BO 8

  945. #define W_BS 4

  946. 

  947. #undef lift5

  948. #define W_CM 9999

  949. #define W_CO 2

  950. #define W_CS 2

  951. 

  952. #define W_DM 15

  953. #define W_DO 16

  954. #define W_DS 5

  955. #elif 0

  956. #define W_AM 55

  957. #define W_AO 16

  958. #define W_AS 5

  959. 

  960. #define W_BM 3

  961. #define W_BO 32

  962. #define W_BS 6

  963. 

  964. #define W_CM 127

  965. #define W_CO 64

  966. #define W_CS 7

  967. 

  968. #define W_DM 7

  969. #define W_DO 8

  970. #define W_DS 4

  971. #elif 0

  972. #define W_AM 97

  973. #define W_AO 32

  974. #define W_AS 6

  975. 

  976. #define W_BM 63

  977. #define W_BO 512

  978. #define W_BS 10

  979. 

  980. #define W_CM 13

  981. #define W_CO 8

  982. #define W_CS 4

  983. 

  984. #define W_DM 15

  985. #define W_DO 16

  986. #define W_DS 5

  987. 

  988. #else

  989. 

  990. #define W_AM 203

  991. #define W_AO 64

  992. #define W_AS 7

  993. 

  994. #define W_BM 217

  995. #define W_BO 2048

  996. #define W_BS 12

  997. 

  998. #define W_CM 113

  999. #define W_CO 64

  1000. #define W_CS 7

  1001. 

  1002. #define W_DM 227

  1003. #define W_DO 128

  1004. #define W_DS 9

  1005. #endif

  1006. static void horizontal_decompose97i(DWTELEM *b, int width){

  1007.     DWTELEM temp[width];

  1008.     const int w2= (width+1)>>1;

  1009. 

  1010.     lift (temp+w2, b    +1, b      , 1, 2, 2, width, -W_AM, W_AO, W_AS, 1, 0);

  1011.     lift (temp   , b      , temp+w2, 1, 2, 1, width, -W_BM, W_BO, W_BS, 0, 0);

  1012.     lift5(b   +w2, temp+w2, temp   , 1, 1, 1, width,  W_CM, W_CO, W_CS, 1, 0);

  1013.     lift (b      , temp   , b   +w2, 1, 1, 1, width,  W_DM, W_DO, W_DS, 0, 0);

  1014. }

  1015. 

  1016. 

  1017. static void vertical_decompose97iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){

  1018.     int i;

  1019.     

  1020.     for(i=0; i<width; i++){

  1021.         b1[i] -= (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;

  1022.     }

  1023. }

  1024. 

  1025. static void vertical_decompose97iH1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){

  1026.     int i;

  1027.     

  1028.     for(i=0; i<width; i++){

  1029. #ifdef lift5

  1030.         b1[i] += (W_CM*(b0[i] + b2[i])+W_CO)>>W_CS;

  1031. #else

  1032.         int r= 3*(b0[i] + b2[i]);

  1033.         r+= r>>4;

  1034.         r+= r>>8;

  1035.         b1[i] += (r+W_CO)>>W_CS;

  1036. #endif

  1037.     }

  1038. }

  1039. 

  1040. static void vertical_decompose97iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){

  1041.     int i;

  1042.     

  1043.     for(i=0; i<width; i++){

  1044.         b1[i] -= (W_BM*(b0[i] + b2[i])+W_BO)>>W_BS;

  1045.     }

  1046. }

  1047. 

  1048. static void vertical_decompose97iL1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){

  1049.     int i;

  1050.     

  1051.     for(i=0; i<width; i++){

  1052.         b1[i] += (W_DM*(b0[i] + b2[i])+W_DO)>>W_DS;

  1053.     }

  1054. }

  1055. 

  1056. static void spatial_decompose97i(DWTELEM *buffer, int width, int height, int stride){

  1057.     int y;

  1058.     DWTELEM *b0= buffer + mirror(-4-1, height-1)*stride;

  1059.     DWTELEM *b1= buffer + mirror(-4  , height-1)*stride;

  1060.     DWTELEM *b2= buffer + mirror(-4+1, height-1)*stride;

  1061.     DWTELEM *b3= buffer + mirror(-4+2, height-1)*stride;

  1062.   

  1063.     for(y=-4; y<height; y+=2){

  1064.         DWTELEM *b4= buffer + mirror(y+3, height-1)*stride;

  1065.         DWTELEM *b5= buffer + mirror(y+4, height-1)*stride;

  1066. 

  1067. {START_TIMER

  1068.         if(b3 <= b5)     horizontal_decompose97i(b4, width);

  1069.         if(y+4 < height) horizontal_decompose97i(b5, width);

  1070. if(width>400){

  1071. STOP_TIMER("horizontal_decompose97i")

  1072. }}

  1073.         

  1074. {START_TIMER

  1075.         if(b3 <= b5) vertical_decompose97iH0(b3, b4, b5, width);

  1076.         if(b2 <= b4) vertical_decompose97iL0(b2, b3, b4, width);

  1077.         if(b1 <= b3) vertical_decompose97iH1(b1, b2, b3, width);

  1078.         if(b0 <= b2) vertical_decompose97iL1(b0, b1, b2, width);

  1079. 

  1080. if(width>400){

  1081. STOP_TIMER("vertical_decompose97i")

  1082. }}

  1083.         

  1084.         b0=b2;

  1085.         b1=b3;

  1086.         b2=b4;

  1087.         b3=b5;

  1088.     }

  1089. }

  1090. 

  1091. void ff_spatial_dwt(DWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count){

  1092.     int level;

  1093.     

  1094.     for(level=0; level<decomposition_count; level++){

  1095.         switch(type){

  1096.         case 0: spatial_decompose97i(buffer, width>>level, height>>level, stride<<level); break;

  1097.         case 1: spatial_decompose53i(buffer, width>>level, height>>level, stride<<level); break;

  1098.         case 2: spatial_decomposeX  (buffer, width>>level, height>>level, stride<<level); break;

  1099.         }

  1100.     }

  1101. }

  1102. 

  1103. static void horizontal_compose53i(DWTELEM *b, int width){

  1104.     DWTELEM temp[width];

  1105.     const int width2= width>>1;

  1106.     const int w2= (width+1)>>1;

  1107.     int A1,A2,A3,A4, x;

  1108. 

  1109. #if 0

  1110.     A2= temp[1       ];

  1111.     A4= temp[0       ];

  1112.     A1= temp[0+width2];

  1113.     A1 -= (A2 + A4)>>1;

  1114.     A4 += (A1 + 1)>>1;

  1115.     b[0+width2] = A1;

  1116.     b[0       ] = A4;

  1117.     for(x=1; x+1<width2; x+=2){

  1118.         A3= temp[x+width2];

  1119.         A4= temp[x+1     ];

  1120.         A3 -= (A2 + A4)>>1;

  1121.         A2 += (A1 + A3 + 2)>>2;

  1122.         b[x+width2] = A3;

  1123.         b[x       ] = A2;

  1124. 

  1125.         A1= temp[x+1+width2];

  1126.         A2= temp[x+2       ];

  1127.         A1 -= (A2 + A4)>>1;

  1128.         A4 += (A1 + A3 + 2)>>2;

  1129.         b[x+1+width2] = A1;

  1130.         b[x+1       ] = A4;

  1131.     }

  1132.     A3= temp[width-1];

  1133.     A3 -= A2;

  1134.     A2 += (A1 + A3 + 2)>>2;

  1135.     b[width -1] = A3;

  1136.     b[width2-1] = A2;

  1137. #else   

  1138.     lift(temp   , b   , b+w2, 1, 1, 1, width,  1, 2, 2, 0, 1);

  1139.     lift(temp+w2, b+w2, temp, 1, 1, 1, width, -1, 0, 1, 1, 1);

  1140. #endif

  1141.     for(x=0; x<width2; x++){

  1142.         b[2*x    ]= temp[x   ];

  1143.         b[2*x + 1]= temp[x+w2];

  1144.     }

  1145.     if(width&1)

  1146.         b[2*x    ]= temp[x   ];

  1147. }

  1148. 

  1149. static void vertical_compose53iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){

  1150.     int i;

  1151.     

  1152.     for(i=0; i<width; i++){

  1153.         b1[i] += (b0[i] + b2[i])>>1;

  1154.     }

  1155. }

  1156. 

  1157. static void vertical_compose53iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){

  1158.     int i;

  1159.     

  1160.     for(i=0; i<width; i++){

  1161.         b1[i] -= (b0[i] + b2[i] + 2)>>2;

  1162.     }

  1163. }

  1164. 

  1165. static void spatial_compose53i(DWTELEM *buffer, int width, int height, int stride){

  1166.     int y;

  1167.     DWTELEM *b0= buffer + mirror(-1-1, height-1)*stride;

  1168.     DWTELEM *b1= buffer + mirror(-1  , height-1)*stride;

  1169.   

  1170.     for(y=-1; y<=height; y+=2){

  1171.         DWTELEM *b2= buffer + mirror(y+1, height-1)*stride;

  1172.         DWTELEM *b3= buffer + mirror(y+2, height-1)*stride;

  1173. 

  1174. {START_TIMER

  1175.         if(b1 <= b3) vertical_compose53iL0(b1, b2, b3, width);

  1176.         if(b0 <= b2) vertical_compose53iH0(b0, b1, b2, width);

  1177. STOP_TIMER("vertical_compose53i*")}

  1178. 

  1179. {START_TIMER

  1180.         if(y-1 >= 0) horizontal_compose53i(b0, width);

  1181.         if(b0 <= b2) horizontal_compose53i(b1, width);

  1182. STOP_TIMER("horizontal_compose53i")}

  1183. 

  1184.         b0=b2;

  1185.         b1=b3;

  1186.     }

  1187. }   

  1188. 

  1189.  

  1190. static void horizontal_compose97i(DWTELEM *b, int width){

  1191.     DWTELEM temp[width];

  1192.     const int w2= (width+1)>>1;

  1193. 

  1194.     lift (temp   , b      , b   +w2, 1, 1, 1, width,  W_DM, W_DO, W_DS, 0, 1);

  1195.     lift5(temp+w2, b   +w2, temp   , 1, 1, 1, width,  W_CM, W_CO, W_CS, 1, 1);

  1196.     lift (b      , temp   , temp+w2, 2, 1, 1, width, -W_BM, W_BO, W_BS, 0, 1);

  1197.     lift (b+1    , temp+w2, b      , 2, 1, 2, width, -W_AM, W_AO, W_AS, 1, 1);

  1198. }

  1199. 

  1200. static void vertical_compose97iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){

  1201.     int i;

  1202.     

  1203.     for(i=0; i<width; i++){

  1204.         b1[i] += (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;

  1205.     }

  1206. }

  1207. 

  1208. static void vertical_compose97iH1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){

  1209.     int i;

  1210.     

  1211.     for(i=0; i<width; i++){

  1212. #ifdef lift5

  1213.         b1[i] -= (W_CM*(b0[i] + b2[i])+W_CO)>>W_CS;

  1214. #else

  1215.         int r= 3*(b0[i] + b2[i]);

  1216.         r+= r>>4;

  1217.         r+= r>>8;

  1218.         b1[i] -= (r+W_CO)>>W_CS;

  1219. #endif

  1220.     }

  1221. }

  1222. 

  1223. static void vertical_compose97iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){

  1224.     int i;

  1225.     

  1226.     for(i=0; i<width; i++){

  1227.         b1[i] += (W_BM*(b0[i] + b2[i])+W_BO)>>W_BS;

  1228.     }

  1229. }

  1230. 

  1231. static void vertical_compose97iL1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){

  1232.     int i;

  1233.     

  1234.     for(i=0; i<width; i++){

  1235.         b1[i] -= (W_DM*(b0[i] + b2[i])+W_DO)>>W_DS;

  1236.     }

  1237. }

  1238. 

  1239. static void spatial_compose97i(DWTELEM *buffer, int width, int height, int stride){

  1240.     int y;

  1241.     DWTELEM *b0= buffer + mirror(-3-1, height-1)*stride;

  1242.     DWTELEM *b1= buffer + mirror(-3  , height-1)*stride;

  1243.     DWTELEM *b2= buffer + mirror(-3+1, height-1)*stride;

  1244.     DWTELEM *b3= buffer + mirror(-3+2, height-1)*stride;

  1245. 

  1246.     for(y=-3; y<=height; y+=2){

  1247.         DWTELEM *b4= buffer + mirror(y+3, height-1)*stride;

  1248.         DWTELEM *b5= buffer + mirror(y+4, height-1)*stride;

  1249. 

  1250.         if(stride == width && y+4 < height && 0){ 

  1251.             int x;

  1252.             for(x=0; x<width/2; x++)

  1253.                 b5[x] += 64*2;

  1254.             for(; x<width; x++)

  1255.                 b5[x] += 169*2;

  1256.         }

  1257.         

  1258. {START_TIMER

  1259.         if(b3 <= b5) vertical_compose97iL1(b3, b4, b5, width);

  1260.         if(b2 <= b4) vertical_compose97iH1(b2, b3, b4, width);

  1261.         if(b1 <= b3) vertical_compose97iL0(b1, b2, b3, width);

  1262.         if(b0 <= b2) vertical_compose97iH0(b0, b1, b2, width);

  1263. if(width>400){

  1264. STOP_TIMER("vertical_compose97i")}}

  1265. 

  1266. {START_TIMER

  1267.         if(y-1>=  0) horizontal_compose97i(b0, width);

  1268.         if(b0 <= b2) horizontal_compose97i(b1, width);

  1269. if(width>400 && b0 <= b2){

  1270. STOP_TIMER("horizontal_compose97i")}}

  1271.         

  1272.         b0=b2;

  1273.         b1=b3;

  1274.         b2=b4;

  1275.         b3=b5;

  1276.     }

  1277. }

  1278. 

  1279. void ff_spatial_idwt(DWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count){

  1280.     int level;

  1281. 

  1282.     for(level=decomposition_count-1; level>=0; level--){

  1283.         switch(type){

  1284.         case 0: spatial_compose97i(buffer, width>>level, height>>level, stride<<level); break;

  1285.         case 1: spatial_compose53i(buffer, width>>level, height>>level, stride<<level); break;

  1286.         case 2: spatial_composeX  (buffer, width>>level, height>>level, stride<<level); break;

  1287.         }

  1288.     }

  1289. }

  1290. 

  1291. static void encode_subband_c0run(SnowContext *s, SubBand *b, DWTELEM *src, DWTELEM *parent, int stride, int orientation){

  1292.     const int w= b->width;

  1293.     const int h= b->height;

  1294.     int x, y;

  1295. 

  1296.     if(1){

  1297.         int run=0;

  1298.         int runs[w*h];

  1299.         int run_index=0;

  1300.                 

  1301.         for(y=0; y<h; y++){

  1302.             for(x=0; x<w; x++){

  1303.                 int v, p=0;

  1304.                 int /*ll=0, */l=0, lt=0, t=0, rt=0;

  1305.                 v= src[x + y*stride];

  1306. 

  1307.                 if(y){

  1308.                     t= src[x + (y-1)*stride];

  1309.                     if(x){

  1310.                         lt= src[x - 1 + (y-1)*stride];

  1311.                     }

  1312.                     if(x + 1 < w){

  1313.                         rt= src[x + 1 + (y-1)*stride];

  1314.                     }

  1315.                 }

  1316.                 if(x){

  1317.                     l= src[x - 1 + y*stride];

  1318.                     /*if(x > 1){

  1319.                         if(orientation==1) ll= src[y + (x-2)*stride];

  1320.                         else               ll= src[x - 2 + y*stride];

  1321.                     }*/

  1322.                 }

  1323.                 if(parent){

  1324.                     int px= x>>1;

  1325.                     int py= y>>1;

  1326.                     if(px<b->parent->width && py<b->parent->height) 

  1327.                         p= parent[px + py*2*stride];

  1328.                 }

  1329.                 if(!(/*ll|*/l|lt|t|rt|p)){

  1330.                     if(v){

  1331.                         runs[run_index++]= run;

  1332.                         run=0;

  1333.                     }else{

  1334.                         run++;

  1335.                     }

  1336.                 }

  1337.             }

  1338.         }

  1339.         runs[run_index++]= run;

  1340.         run_index=0;

  1341.         run= runs[run_index++];

  1342. 

  1343.         put_symbol2(&s->c, b->state[1], run, 3);

  1344.         

  1345.         for(y=0; y<h; y++){

  1346.             for(x=0; x<w; x++){

  1347.                 int v, p=0;

  1348.                 int /*ll=0, */l=0, lt=0, t=0, rt=0;

  1349.                 v= src[x + y*stride];

  1350. 

  1351.                 if(y){

  1352.                     t= src[x + (y-1)*stride];

  1353.                     if(x){

  1354.                         lt= src[x - 1 + (y-1)*stride];

  1355.                     }

  1356.                     if(x + 1 < w){

  1357.                         rt= src[x + 1 + (y-1)*stride];

  1358.                     }

  1359.                 }

  1360.                 if(x){

  1361.                     l= src[x - 1 + y*stride];

  1362.                     /*if(x > 1){

  1363.                         if(orientation==1) ll= src[y + (x-2)*stride];

  1364.                         else               ll= src[x - 2 + y*stride];

  1365.                     }*/

  1366.                 }

  1367.                 if(parent){

  1368.                     int px= x>>1;

  1369.                     int py= y>>1;

  1370.                     if(px<b->parent->width && py<b->parent->height) 

  1371.                         p= parent[px + py*2*stride];

  1372.                 }

  1373.                 if(/*ll|*/l|lt|t|rt|p){

  1374.                     int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));

  1375. 

  1376.                     put_rac(&s->c, &b->state[0][context], !!v);

  1377.                 }else{

  1378.                     if(!run){

  1379.                         run= runs[run_index++];

  1380. 

  1381.                         put_symbol2(&s->c, b->state[1], run, 3);

  1382.                         assert(v);

  1383.                     }else{

  1384.                         run--;

  1385.                         assert(!v);

  1386.                     }

  1387.                 }

  1388.                 if(v){

  1389.                     int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));

  1390. 

  1391.                     put_symbol2(&s->c, b->state[context + 2], ABS(v)-1, context-4);

  1392.                     put_rac(&s->c, &b->state[0][16 + 1 + 3 + quant3b[l&0xFF] + 3*quant3b[t&0xFF]], v<0);

  1393.                 }

  1394.             }

  1395.         }

  1396.     }

  1397. }

  1398. 

  1399. static void encode_subband(SnowContext *s, SubBand *b, DWTELEM *src, DWTELEM *parent, int stride, int orientation){    

  1400. //    encode_subband_qtree(s, b, src, parent, stride, orientation);

  1401. //    encode_subband_z0run(s, b, src, parent, stride, orientation);

  1402.     encode_subband_c0run(s, b, src, parent, stride, orientation);

  1403. //    encode_subband_dzr(s, b, src, parent, stride, orientation);

  1404. }

  1405. 

  1406. static inline void decode_subband(SnowContext *s, SubBand *b, DWTELEM *src, DWTELEM *parent, int stride, int orientation){

  1407.     const int w= b->width;

  1408.     const int h= b->height;

  1409.     int x,y;

  1410.     const int qlog= clip(s->qlog + b->qlog, 0, 128);

  1411.     int qmul= qexp[qlog&7]<<(qlog>>3);

  1412.     int qadd= (s->qbias*qmul)>>QBIAS_SHIFT;

  1413.     

  1414.     START_TIMER

  1415. 

  1416.     if(b->buf == s->spatial_dwt_buffer || s->qlog == LOSSLESS_QLOG){

  1417.         qadd= 0;

  1418.         qmul= 1<<QEXPSHIFT;

  1419.     }

  1420. 

  1421.     if(1){

  1422.         int run;

  1423.         int index=0;

  1424.         int prev_index=-1;

  1425.         int prev2_index=0;

  1426.         int parent_index= 0;

  1427.         int prev_parent_index= 0;

  1428.         

  1429.         for(y=0; y<b->height; y++)

  1430.             memset(&src[y*stride], 0, b->width*sizeof(DWTELEM));

  1431. 

  1432.         run= get_symbol2(&s->c, b->state[1], 3);

  1433.         for(y=0; y<h; y++){

  1434.             int v=0;

  1435.             int lt=0, t=0, rt=0;

  1436. 

  1437.             if(y && b->x[prev_index] == 0){

  1438.                 rt= b->coeff[prev_index];

  1439.             }

  1440.             for(x=0; x<w; x++){

  1441.                 int p=0;

  1442.                 const int l= v;

  1443.                 

  1444.                 lt= t; t= rt;

  1445. 

  1446.                 if(y){

  1447.                     if(b->x[prev_index] <= x)

  1448.                         prev_index++;

  1449.                     if(b->x[prev_index] == x + 1)

  1450.                         rt= b->coeff[prev_index];

  1451.                     else

  1452.                         rt=0;

  1453.                 }

  1454.                 if(parent){

  1455.                     if(x>>1 > b->parent->x[parent_index]){

  1456.                         parent_index++;

  1457.                     }

  1458.                     if(x>>1 == b->parent->x[parent_index]){

  1459.                         p= b->parent->coeff[parent_index];

  1460.                     }

  1461.                 }

  1462.                 if(/*ll|*/l|lt|t|rt|p){

  1463.                     int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));

  1464. 

  1465.                     v=get_rac(&s->c, &b->state[0][context]);

  1466.                 }else{

  1467.                     if(!run){

  1468.                         run= get_symbol2(&s->c, b->state[1], 3);

  1469.                         v=1;

  1470.                     }else{

  1471.                         run--;

  1472.                         v=0;

  1473. 

  1474.                         if(y && parent){

  1475.                             int max_run;

  1476. 

  1477.                             max_run= FFMIN(run, b->x[prev_index] - x - 2);

  1478.                             max_run= FFMIN(max_run, 2*b->parent->x[parent_index] - x - 1);

  1479.                             x+= max_run;

  1480.                             run-= max_run;

  1481.                         }

  1482.                     }

  1483.                 }

  1484.                 if(v){

  1485.                     int context= av_log2(/*ABS(ll) + */3*ABS(l) + ABS(lt) + 2*ABS(t) + ABS(rt) + ABS(p));

  1486.                     v= get_symbol2(&s->c, b->state[context + 2], context-4) + 1;

  1487.                     if(get_rac(&s->c, &b->state[0][16 + 1 + 3 + quant3b[l&0xFF] + 3*quant3b[t&0xFF]])){

  1488.                         src[x + y*stride]=-(( v*qmul + qadd)>>(QEXPSHIFT));

  1489.                         v= -v;

  1490.                     }else{

  1491.                         src[x + y*stride]= (( v*qmul + qadd)>>(QEXPSHIFT));

  1492.                     }

  1493.                     b->x[index]=x; //FIXME interleave x/coeff

  1494.                     b->coeff[index++]= v;

  1495.                 }

  1496.             }

  1497.             b->x[index++]= w+1; //end marker

  1498.             prev_index= prev2_index;

  1499.             prev2_index= index;

  1500.             

  1501.             if(parent){

  1502.                 while(b->parent->x[parent_index] != b->parent->width+1)

  1503.                     parent_index++;

  1504.                 parent_index++;

  1505.                 if(y&1){

  1506.                     prev_parent_index= parent_index;

  1507.                 }else{

  1508.                     parent_index= prev_parent_index;

  1509.                 }

  1510.             }

  1511.         }

  1512.         b->x[index++]= w+1; //end marker

  1513.         if(w > 200 /*level+1 == s->spatial_decomposition_count*/){

  1514.             STOP_TIMER("decode_subband")

  1515.         }

  1516.         

  1517.         return;

  1518.     }

  1519. }

  1520. 

  1521. static void reset_contexts(SnowContext *s){

  1522.     int plane_index, level, orientation;

  1523. 

  1524.     for(plane_index=0; plane_index<3; plane_index++){

  1525.         for(level=0; level<s->spatial_decomposition_count; level++){

  1526.             for(orientation=level ? 1:0; orientation<4; orientation++){

  1527.                 memset(s->plane[plane_index].band[level][orientation].state, MID_STATE, sizeof(s->plane[plane_index].band[level][orientation].state));

  1528.             }

  1529.         }

  1530.     }

  1531.     memset(s->header_state, MID_STATE, sizeof(s->header_state));

  1532.     memset(s->block_state, MID_STATE, sizeof(s->block_state));

  1533. }

  1534. 

  1535. static int alloc_blocks(SnowContext *s){

  1536.     int w= -((-s->avctx->width )>>LOG2_MB_SIZE);

  1537.     int h= -((-s->avctx->height)>>LOG2_MB_SIZE);

  1538.     

  1539.     s->b_width = w;

  1540.     s->b_height= h;

  1541.     

  1542.     s->block= av_mallocz(w * h * sizeof(BlockNode) << (s->block_max_depth*2));

  1543.     return 0;

  1544. }

  1545. 

  1546. static inline void copy_rac_state(RangeCoder *d, RangeCoder *s){

  1547.     uint8_t *bytestream= d->bytestream;

  1548.     uint8_t *bytestream_start= d->bytestream_start;

  1549.     *d= *s;

  1550.     d->bytestream= bytestream;

  1551.     d->bytestream_start= bytestream_start;

  1552. }

  1553. 

  1554. //near copy & paste from dsputil, FIXME

  1555. static int pix_sum(uint8_t * pix, int line_size, int w)

  1556. {

  1557.     int s, i, j;

  1558. 

  1559.     s = 0;

  1560.     for (i = 0; i < w; i++) {

  1561.         for (j = 0; j < w; j++) {

  1562.             s += pix[0];

  1563.             pix ++;

  1564.         }

  1565.         pix += line_size - w;

  1566.     }

  1567.     return s;

  1568. }

  1569. 

  1570. //near copy & paste from dsputil, FIXME

  1571. static int pix_norm1(uint8_t * pix, int line_size, int w)

  1572. {

  1573.     int s, i, j;

  1574.     uint32_t *sq = squareTbl + 256;

  1575. 

  1576.     s = 0;

  1577.     for (i = 0; i < w; i++) {

  1578.         for (j = 0; j < w; j ++) {

  1579.             s += sq[pix[0]];

  1580.             pix ++;

  1581.         }

  1582.         pix += line_size - w;

  1583.     }

  1584.     return s;

  1585. }

  1586. 

  1587. static inline void set_blocks(SnowContext *s, int level, int x, int y, int l, int cb, int cr, int mx, int my, int type){

  1588.     const int w= s->b_width << s->block_max_depth;

  1589.     const int rem_depth= s->block_max_depth - level;

  1590.     const int index= (x + y*w) << rem_depth;

  1591.     const int block_w= 1<<rem_depth;

  1592.     BlockNode block;

  1593.     int i,j;

  1594.     

  1595.     block.color[0]= l;

  1596.     block.color[1]= cb;

  1597.     block.color[2]= cr;

  1598.     block.mx= mx;

  1599.     block.my= my;

  1600.     block.type= type;

  1601.     block.level= level;

  1602. 

  1603.     for(j=0; j<block_w; j++){

  1604.         for(i=0; i<block_w; i++){

  1605.             s->block[index + i + j*w]= block;

  1606.         }

  1607.     }

  1608. }

  1609. 

  1610. static inline void init_ref(MotionEstContext *c, uint8_t *src[3], uint8_t *ref[3], uint8_t *ref2[3], int x, int y, int ref_index){

  1611.     const int offset[3]= {

  1612.           y*c->  stride + x,

  1613.         ((y*c->uvstride + x)>>1),

  1614.         ((y*c->uvstride + x)>>1),

  1615.     };

  1616.     int i;

  1617.     for(i=0; i<3; i++){

  1618.         c->src[0][i]= src [i];

  1619.         c->ref[0][i]= ref [i] + offset[i];

  1620.     }

  1621.     assert(!ref_index);

  1622. }

  1623. 

  1624. //FIXME copy&paste

  1625. #define P_LEFT P[1]

  1626. #define P_TOP P[2]

  1627. #define P_TOPRIGHT P[3]

  1628. #define P_MEDIAN P[4]

  1629. #define P_MV1 P[9]

  1630. #define FLAG_QPEL   1 //must be 1

  1631. 

  1632. static int encode_q_branch(SnowContext *s, int level, int x, int y){

  1633.     uint8_t p_buffer[1024];

  1634.     uint8_t i_buffer[1024];

  1635.     uint8_t p_state[sizeof(s->block_state)];

  1636.     uint8_t i_state[sizeof(s->block_state)];

  1637.     RangeCoder pc, ic;

  1638.     uint8_t *pbbak= s->c.bytestream;

  1639.     uint8_t *pbbak_start= s->c.bytestream_start;

  1640.     int score, score2, iscore, i_len, p_len, block_s, sum;

  1641.     const int w= s->b_width  << s->block_max_depth;

  1642.     const int h= s->b_height << s->block_max_depth;

  1643.     const int rem_depth= s->block_max_depth - level;

  1644.     const int index= (x + y*w) << rem_depth;

  1645.     const int block_w= 1<<(LOG2_MB_SIZE - level);

  1646.     static BlockNode null_block= { //FIXME add border maybe

  1647.         .color= {128,128,128},

  1648.         .mx= 0,

  1649.         .my= 0,

  1650.         .type= 0,

  1651.         .level= 0,

  1652.     };

  1653.     int trx= (x+1)<<rem_depth;

  1654.     int try= (y+1)<<rem_depth;

  1655.     BlockNode *left  = x ? &s->block[index-1] : &null_block;

  1656.     BlockNode *top   = y ? &s->block[index-w] : &null_block;

  1657.     BlockNode *right = trx<w ? &s->block[index+1] : &null_block;

  1658.     BlockNode *bottom= try<h ? &s->block[index+w] : &null_block;

  1659.     BlockNode *tl    = y && x ? &s->block[index-w-1] : left;

  1660.     BlockNode *tr    = y && trx<w && ((x&1)==0 || level==0) ? &s->block[index-w+(1<<rem_depth)] : tl; //FIXME use lt

  1661.     int pl = left->color[0];

  1662.     int pcb= left->color[1];

  1663.     int pcr= left->color[2];

  1664.     int pmx= mid_pred(left->mx, top->mx, tr->mx);

  1665.     int pmy= mid_pred(left->my, top->my, tr->my);

  1666.     int mx=0, my=0;

  1667.     int l,cr,cb, i;

  1668.     const int stride= s->current_picture.linesize[0];

  1669.     const int uvstride= s->current_picture.linesize[1];

  1670.     const int instride= s->input_picture.linesize[0];

  1671.     const int uvinstride= s->input_picture.linesize[1];

  1672.     uint8_t *new_l = s->input_picture.data[0] + (x + y*  instride)*block_w;

  1673.     uint8_t *new_cb= s->input_picture.data[1] + (x + y*uvinstride)*block_w/2;

  1674.     uint8_t *new_cr= s->input_picture.data[2] + (x + y*uvinstride)*block_w/2;

  1675.     uint8_t current_mb[3][stride*block_w];

  1676.     uint8_t *current_data[3]= {&current_mb[0][0], &current_mb[1][0], &current_mb[2][0]};

  1677.     int P[10][2];

  1678.     int16_t last_mv[3][2];

  1679.     int qpel= !!(s->avctx->flags & CODEC_FLAG_QPEL); //unused

  1680.     const int shift= 1+qpel;

  1681.     MotionEstContext *c= &s->m.me;

  1682.     int mx_context= av_log2(2*ABS(left->mx - top->mx));

  1683.     int my_context= av_log2(2*ABS(left->my - top->my));

  1684.     int s_context= 2*left->level + 2*top->level + tl->level + tr->level;

  1685. 

  1686.     assert(sizeof(s->block_state) >= 256);

  1687.     if(s->keyframe){

  1688.         set_blocks(s, level, x, y, pl, pcb, pcr, pmx, pmy, BLOCK_INTRA);

  1689.         return 0;

  1690.     }

  1691. 

  1692.     //FIXME optimize

  1693.     for(i=0; i<block_w; i++)

  1694.         memcpy(&current_mb[0][0] +   stride*i, new_l  +   instride*i, block_w);

  1695.     for(i=0; i<block_w>>1; i++)

  1696.         memcpy(&current_mb[1][0] + uvstride*i, new_cb + uvinstride*i, block_w>>1);

  1697.     for(i=0; i<block_w>>1; i++)

  1698.         memcpy(&current_mb[2][0] + uvstride*i, new_cr + uvinstride*i, block_w>>1);

  1699. 

  1700. //    clip predictors / edge ?

  1701. 

  1702.     P_LEFT[0]= left->mx;

  1703.     P_LEFT[1]= left->my;

  1704.     P_TOP [0]= top->mx;

  1705.     P_TOP [1]= top->my;

  1706.     P_TOPRIGHT[0]= tr->mx;

  1707.     P_TOPRIGHT[1]= tr->my;

  1708.     

  1709.     last_mv[0][0]= s->block[index].mx;

  1710.     last_mv[0][1]= s->block[index].my;

  1711.     last_mv[1][0]= right->mx;

  1712.     last_mv[1][1]= right->my;

  1713.     last_mv[2][0]= bottom->mx;

  1714.     last_mv[2][1]= bottom->my;

  1715.     

  1716.     s->m.mb_stride=2;

  1717.     s->m.mb_x= 

  1718.     s->m.mb_y= 0;

  1719.     s->m.me.skip= 0;

  1720. 

  1721.     init_ref(c, current_data, s->last_picture.data, NULL, block_w*x, block_w*y, 0);

  1722.     

  1723.     assert(s->m.me.  stride ==   stride);

  1724.     assert(s->m.me.uvstride == uvstride);

  1725.     

  1726.     c->penalty_factor    = get_penalty_factor(s->lambda, s->lambda2, c->avctx->me_cmp);

  1727.     c->sub_penalty_factor= get_penalty_factor(s->lambda, s->lambda2, c->avctx->me_sub_cmp);

  1728.     c->mb_penalty_factor = get_penalty_factor(s->lambda, s->lambda2, c->avctx->mb_cmp);

  1729.     c->current_mv_penalty= c->mv_penalty[s->m.f_code=1] + MAX_MV;

  1730.     

  1731.     c->xmin = - x*block_w - 16+2;

  1732.     c->ymin = - y*block_w - 16+2;

  1733.     c->xmax = - (x+1)*block_w + (w<<(LOG2_MB_SIZE - s->block_max_depth)) + 16-2;

  1734.     c->ymax = - (y+1)*block_w + (h<<(LOG2_MB_SIZE - s->block_max_depth)) + 16-2;

  1735. 

  1736.     if(P_LEFT[0]     > (c->xmax<<shift)) P_LEFT[0]    = (c->xmax<<shift);

  1737.     if(P_LEFT[1]     > (c->ymax<<shift)) P_LEFT[1]    = (c->ymax<<shift); 

  1738.     if(P_TOP[0]      > (c->xmax<<shift)) P_TOP[0]     = (c->xmax<<shift);

  1739.     if(P_TOP[1]      > (c->ymax<<shift)) P_TOP[1]     = (c->ymax<<shift);

  1740.     if(P_TOPRIGHT[0] < (c->xmin<<shift)) P_TOPRIGHT[0]= (c->xmin<<shift);

  1741.     if(P_TOPRIGHT[0] > (c->xmax<<shift)) P_TOPRIGHT[0]= (c->xmax<<shift); //due to pmx no clip

  1742.     if(P_TOPRIGHT[1] > (c->ymax<<shift)) P_TOPRIGHT[1]= (c->ymax<<shift);

  1743. 

  1744.     P_MEDIAN[0]= mid_pred(P_LEFT[0], P_TOP[0], P_TOPRIGHT[0]);

  1745.     P_MEDIAN[1]= mid_pred(P_LEFT[1], P_TOP[1], P_TOPRIGHT[1]);

  1746. 

  1747.     if (!y) {

  1748.         c->pred_x= P_LEFT[0];

  1749.         c->pred_y= P_LEFT[1];

  1750.     } else {

  1751.         c->pred_x = P_MEDIAN[0];

  1752.         c->pred_y = P_MEDIAN[1];

  1753.     }

  1754. 

  1755.     score= ff_epzs_motion_search(&s->m, &mx, &my, P, 0, /*ref_index*/ 0, last_mv, 

  1756.                              (1<<16)>>shift, level-LOG2_MB_SIZE+4, block_w);

  1757. 

  1758.     assert(mx >= c->xmin);

  1759.     assert(mx <= c->xmax);

  1760.     assert(my >= c->ymin);

  1761.     assert(my <= c->ymax);

  1762.     

  1763.     score= s->m.me.sub_motion_search(&s->m, &mx, &my, score, 0, 0, level-LOG2_MB_SIZE+4, block_w);

  1764.     score= ff_get_mb_score(&s->m, mx, my, 0, 0, level-LOG2_MB_SIZE+4, block_w, 0);

  1765.     //FIXME if mb_cmp != SSE then intra cant be compared currently and mb_penalty vs. lambda2

  1766.                              

  1767.   //  subpel search

  1768.     pc= s->c;

  1769.     pc.bytestream_start=

  1770.     pc.bytestream= p_buffer; //FIXME end/start? and at the other stoo

  1771.     memcpy(p_state, s->block_state, sizeof(s->block_state));

  1772. 

  1773.     if(level!=s->block_max_depth)

  1774.         put_rac(&pc, &p_state[4 + s_context], 1);

  1775.     put_rac(&pc, &p_state[1 + left->type + top->type], 0);

  1776.     put_symbol(&pc, &p_state[128 + 32*mx_context], mx - pmx, 1);

  1777.     put_symbol(&pc, &p_state[128 + 32*my_context], my - pmy, 1);

  1778.     p_len= pc.bytestream - pc.bytestream_start;

  1779.     score += (s->lambda2*(p_len*8

  1780.               + (pc.outstanding_count - s->c.outstanding_count)*8

  1781.               + (-av_log2(pc.range)    + av_log2(s->c.range))

  1782.              ))>>FF_LAMBDA_SHIFT;

  1783. 

  1784.     block_s= block_w*block_w;

  1785.     sum = pix_sum(&current_mb[0][0], stride, block_w);

  1786.     l= (sum + block_s/2)/block_s;

  1787.     iscore = pix_norm1(&current_mb[0][0], stride, block_w) - 2*l*sum + l*l*block_s;

  1788.     

  1789.     block_s= block_w*block_w>>2;

  1790.     sum = pix_sum(&current_mb[1][0], uvstride, block_w>>1);

  1791.     cb= (sum + block_s/2)/block_s;

  1792. //    iscore += pix_norm1(&current_mb[1][0], uvstride, block_w>>1) - 2*cb*sum + cb*cb*block_s;

  1793.     sum = pix_sum(&current_mb[2][0], uvstride, block_w>>1);

  1794.     cr= (sum + block_s/2)/block_s;

  1795. //    iscore += pix_norm1(&current_mb[2][0], uvstride, block_w>>1) - 2*cr*sum + cr*cr*block_s;

  1796. 

  1797.     ic= s->c;

  1798.     ic.bytestream_start=

  1799.     ic.bytestream= i_buffer; //FIXME end/start? and at the other stoo

  1800.     memcpy(i_state, s->block_state, sizeof(s->block_state));

  1801.     if(level!=s->block_max_depth)

  1802.         put_rac(&ic, &i_state[4 + s_context], 1);

  1803.     put_rac(&ic, &i_state[1 + left->type + top->type], 1);

  1804.     put_symbol(&ic, &i_state[32],  l-pl , 1);

  1805.     put_symbol(&ic, &i_state[64], cb-pcb, 1);

  1806.     put_symbol(&ic, &i_state[96], cr-pcr, 1);

  1807.     i_len= ic.bytestream - ic.bytestream_start;

  1808.     iscore += (s->lambda2*(i_len*8

  1809.               + (ic.outstanding_count - s->c.outstanding_count)*8

  1810.               + (-av_log2(ic.range)    + av_log2(s->c.range))

  1811.              ))>>FF_LAMBDA_SHIFT;

  1812. 

  1813. //    assert(score==256*256*256*64-1);

  1814.     assert(iscore < 255*255*256 + s->lambda2*10);

  1815.     assert(iscore >= 0);

  1816.     assert(l>=0 && l<=255);

  1817.     assert(pl>=0 && pl<=255);

  1818. 

  1819.     if(level==0){

  1820.         int varc= iscore >> 8;

  1821.         int vard= score >> 8;

  1822.         if (vard <= 64 || vard < varc)

  1823.             c->scene_change_score+= ff_sqrt(vard) - ff_sqrt(varc);

  1824.         else

  1825.             c->scene_change_score+= s->m.qscale;

  1826.     }

  1827.         

  1828.     if(level!=s->block_max_depth){

  1829.         put_rac(&s->c, &s->block_state[4 + s_context], 0);

  1830.         score2 = encode_q_branch(s, level+1, 2*x+0, 2*y+0);

  1831.         score2+= encode_q_branch(s, level+1, 2*x+1, 2*y+0);

  1832.         score2+= encode_q_branch(s, level+1, 2*x+0, 2*y+1);

  1833.         score2+= encode_q_branch(s, level+1, 2*x+1, 2*y+1);

  1834.         score2+= s->lambda2>>FF_LAMBDA_SHIFT; //FIXME exact split overhead

  1835.     

  1836.         if(score2 < score && score2 < iscore)

  1837.             return score2;

  1838.     }

  1839.     

  1840.     if(iscore < score){

  1841.         memcpy(pbbak, i_buffer, i_len);

  1842.         s->c= ic;

  1843.         s->c.bytestream_start= pbbak_start;

  1844.         s->c.bytestream= pbbak + i_len;

  1845.         set_blocks(s, level, x, y, l, cb, cr, pmx, pmy, BLOCK_INTRA);

  1846.         memcpy(s->block_state, i_state, sizeof(s->block_state));

  1847.         return iscore;

  1848.     }else{

  1849.         memcpy(pbbak, p_buffer, p_len);

  1850.         s->c= pc;

  1851.         s->c.bytestream_start= pbbak_start;

  1852.         s->c.bytestream= pbbak + p_len;

  1853.         set_blocks(s, level, x, y, pl, pcb, pcr, mx, my, 0);

  1854.         memcpy(s->block_state, p_state, sizeof(s->block_state));

  1855.         return score;

  1856.     }

  1857. }

  1858. 

  1859. static void decode_q_branch(SnowContext *s, int level, int x, int y){

  1860.     const int w= s->b_width << s->block_max_depth;

  1861.     const int rem_depth= s->block_max_depth - level;

  1862.     const int index= (x + y*w) << rem_depth;

  1863.     static BlockNode null_block= { //FIXME add border maybe

  1864.         .color= {128,128,128},

  1865.         .mx= 0,

  1866.         .my= 0,

  1867.         .type= 0,

  1868.         .level= 0,

  1869.     };

  1870.     int trx= (x+1)<<rem_depth;

  1871.     BlockNode *left  = x ? &s->block[index-1] : &null_block;

  1872.     BlockNode *top   = y ? &s->block[index-w] : &null_block;

  1873.     BlockNode *tl    = y && x ? &s->block[index-w-1] : left;

  1874.     BlockNode *tr    = y && trx<w && ((x&1)==0 || level==0) ? &s->block[index-w+(1<<rem_depth)] : tl; //FIXME use lt

  1875.     int s_context= 2*left->level + 2*top->level + tl->level + tr->level;

  1876.     

  1877.     if(s->keyframe){

  1878.         set_blocks(s, level, x, y, null_block.color[0], null_block.color[1], null_block.color[2], null_block.mx, null_block.my, BLOCK_INTRA);

  1879.         return;

  1880.     }

  1881. 

  1882.     if(level==s->block_max_depth || get_rac(&s->c, &s->block_state[4 + s_context])){

  1883.         int type;

  1884.         int l = left->color[0];

  1885.         int cb= left->color[1];

  1886.         int cr= left->color[2];

  1887.         int mx= mid_pred(left->mx, top->mx, tr->mx);

  1888.         int my= mid_pred(left->my, top->my, tr->my);

  1889.         int mx_context= av_log2(2*ABS(left->mx - top->mx)) + 0*av_log2(2*ABS(tr->mx - top->mx));

  1890.         int my_context= av_log2(2*ABS(left->my - top->my)) + 0*av_log2(2*ABS(tr->my - top->my));

  1891.         

  1892.         type= get_rac(&s->c, &s->block_state[1 + left->type + top->type]) ? BLOCK_INTRA : 0;

  1893. 

  1894.         if(type){

  1895.             l += get_symbol(&s->c, &s->block_state[32], 1);

  1896.             cb+= get_symbol(&s->c, &s->block_state[64], 1);

  1897.             cr+= get_symbol(&s->c, &s->block_state[96], 1);

  1898.         }else{

  1899.             mx+= get_symbol(&s->c, &s->block_state[128 + 32*mx_context], 1);

  1900.             my+= get_symbol(&s->c, &s->block_state[128 + 32*my_context], 1);

  1901.         }

  1902.         set_blocks(s, level, x, y, l, cb, cr, mx, my, type);

  1903.     }else{

  1904.         decode_q_branch(s, level+1, 2*x+0, 2*y+0);

  1905.         decode_q_branch(s, level+1, 2*x+1, 2*y+0);

  1906.         decode_q_branch(s, level+1, 2*x+0, 2*y+1);

  1907.         decode_q_branch(s, level+1, 2*x+1, 2*y+1);

  1908.     }

  1909. }

  1910. 

  1911. static void encode_blocks(SnowContext *s){

  1912.     int x, y;

  1913.     int w= s->b_width;

  1914.     int h= s->b_height;

  1915. 

  1916.     for(y=0; y<h; y++){

  1917.         for(x=0; x<w; x++){

  1918.             encode_q_branch(s, 0, x, y);

  1919.         }

  1920.     }

  1921. }

  1922. 

  1923. static void decode_blocks(SnowContext *s){

  1924.     int x, y;

  1925.     int w= s->b_width;

  1926.     int h= s->b_height;

  1927. 

  1928.     for(y=0; y<h; y++){

  1929.         for(x=0; x<w; x++){

  1930.             decode_q_branch(s, 0, x, y);

  1931.         }

  1932.     }

  1933. }

  1934. 

  1935. static void mc_block(uint8_t *dst, uint8_t *src, uint8_t *tmp, int stride, int b_w, int b_h, int dx, int dy){

  1936.     int x, y;

  1937. START_TIMER

  1938.     for(y=0; y < b_h+5; y++){

  1939.         for(x=0; x < b_w; x++){

  1940.             int a0= src[x    ];

  1941.             int a1= src[x + 1];

  1942.             int a2= src[x + 2];

  1943.             int a3= src[x + 3];

  1944.             int a4= src[x + 4];

  1945.             int a5= src[x + 5];

  1946. //            int am= 9*(a1+a2) - (a0+a3);

  1947.             int am= 20*(a2+a3) - 5*(a1+a4) + (a0+a5);

  1948. //            int am= 18*(a2+a3) - 2*(a1+a4);

  1949. //             int aL= (-7*a0 + 105*a1 + 35*a2 - 5*a3)>>3;

  1950. //             int aR= (-7*a3 + 105*a2 + 35*a1 - 5*a0)>>3;

  1951. 

  1952. //            if(b_w==16) am= 8*(a1+a2);

  1953. 

  1954.             if(dx<8) tmp[x]= (32*a2*( 8-dx) +    am* dx    + 128)>>8;

  1955.             else     tmp[x]= (   am*(16-dx) + 32*a3*(dx-8) + 128)>>8;

  1956. 

  1957. /*            if     (dx< 4) tmp[x + y*stride]= (16*a1*( 4-dx) +    aL* dx     + 32)>>6;

  1958.             else if(dx< 8) tmp[x + y*stride]= (   aL*( 8-dx) +    am*(dx- 4) + 32)>>6;

  1959.             else if(dx<12) tmp[x + y*stride]= (   am*(12-dx) +    aR*(dx- 8) + 32)>>6;

  1960.             else           tmp[x + y*stride]= (   aR*(16-dx) + 16*a2*(dx-12) + 32)>>6;*/

  1961.         }

  1962.         tmp += stride;

  1963.         src += stride;

  1964.     }

  1965.     tmp -= (b_h+5)*stride;

  1966.     

  1967.     for(y=0; y < b_h; y++){

  1968.         for(x=0; x < b_w; x++){

  1969.             int a0= tmp[x + 0*stride];

  1970.             int a1= tmp[x + 1*stride];

  1971.             int a2= tmp[x + 2*stride];

  1972.             int a3= tmp[x + 3*stride];

  1973.             int a4= tmp[x + 4*stride];

  1974.             int a5= tmp[x + 5*stride];

  1975.             int am= 20*(a2+a3) - 5*(a1+a4) + (a0+a5);

  1976. //            int am= 18*(a2+a3) - 2*(a1+a4);

  1977. /*            int aL= (-7*a0 + 105*a1 + 35*a2 - 5*a3)>>3;

  1978.             int aR= (-7*a3 + 105*a2 + 35*a1 - 5*a0)>>3;*/

  1979.             

  1980. //            if(b_w==16) am= 8*(a1+a2);

  1981. 

  1982.             if(dy<8) dst[x]= (32*a2*( 8-dy) +    am* dy    + 128)>>8;

  1983.             else     dst[x]= (   am*(16-dy) + 32*a3*(dy-8) + 128)>>8;

  1984. 

  1985. /*            if     (dy< 4) tmp[x + y*stride]= (16*a1*( 4-dy) +    aL* dy     + 32)>>6;

  1986.             else if(dy< 8) tmp[x + y*stride]= (   aL*( 8-dy) +    am*(dy- 4) + 32)>>6;

  1987.             else if(dy<12) tmp[x + y*stride]= (   am*(12-dy) +    aR*(dy- 8) + 32)>>6;

  1988.             else           tmp[x + y*stride]= (   aR*(16-dy) + 16*a2*(dy-12) + 32)>>6;*/

  1989.         }

  1990.         dst += stride;

  1991.         tmp += stride;

  1992.     }

  1993. STOP_TIMER("mc_block")

  1994. }

  1995. 

  1996. #define mca(dx,dy,b_w)\

  1997. static void mc_block_hpel ## dx ## dy ## b_w(uint8_t *dst, uint8_t *src, int stride, int h){\

  1998.     uint8_t tmp[stride*(b_w+5)];\

  1999.     assert(h==b_w);\

  2000.     mc_block(dst, src-2-2*stride, tmp, stride, b_w, b_w, dx, dy);\

  2001. }

  2002. 

  2003. mca( 0, 0,16)

  2004. mca( 8, 0,16)

  2005. mca( 0, 8,16)

  2006. mca( 8, 8,16)

  2007. mca( 0, 0,8)

  2008. mca( 8, 0,8)

  2009. mca( 0, 8,8)

  2010. mca( 8, 8,8)

  2011. 

  2012. static void pred_block(SnowContext *s, uint8_t *dst, uint8_t *src, uint8_t *tmp, int stride, int sx, int sy, int b_w, int b_h, BlockNode *block, int plane_index, int w, int h){

  2013.     if(block->type){

  2014.         int x, y;

  2015.         const int color= block->color[plane_index];

  2016.         for(y=0; y < b_h; y++){

  2017.             for(x=0; x < b_w; x++){

  2018.                 dst[x + y*stride]= color;

  2019.             }

  2020.         }

  2021.     }else{

  2022.         const int scale= plane_index ?  s->mv_scale : 2*s->mv_scale;

  2023.         int mx= block->mx*scale;

  2024.         int my= block->my*scale;

  2025.         const int dx= mx&15;

  2026.         const int dy= my&15;

  2027.         sx += (mx>>4) - 2;

  2028.         sy += (my>>4) - 2;

  2029.         src += sx + sy*stride;

  2030.         if(   (unsigned)sx >= w - b_w - 4

  2031.            || (unsigned)sy >= h - b_h - 4){

  2032.             ff_emulated_edge_mc(tmp + MB_SIZE, src, stride, b_w+5, b_h+5, sx, sy, w, h);

  2033.             src= tmp + MB_SIZE;

  2034.         }

  2035.         if((dx&3) || (dy&3) || b_w!=b_h || (b_w!=4 && b_w!=8 && b_w!=16))

  2036.             mc_block(dst, src, tmp, stride, b_w, b_h, dx, dy);

  2037.         else

  2038.             s->dsp.put_h264_qpel_pixels_tab[2-(b_w>>3)][dy+(dx>>2)](dst,src + 2 + 2*stride,stride);

  2039.     }

  2040. }

  2041. 

  2042. static always_inline int same_block(BlockNode *a, BlockNode *b){

  2043.     return !((a->mx - b->mx) | (a->my - b->my) | a->type | b->type);

  2044. }

  2045. 

  2046. //FIXME name clenup (b_w, block_w, b_width stuff)

  2047. static always_inline void add_yblock(SnowContext *s, DWTELEM *dst, uint8_t *dst8, uint8_t *src, uint8_t *obmc, int src_x, int src_y, int b_w, int b_h, int w, int h, int dst_stride, int src_stride, int obmc_stride, int b_x, int b_y, int add, int plane_index){

  2048.     const int b_width = s->b_width  << s->block_max_depth;

  2049.     const int b_height= s->b_height << s->block_max_depth;

  2050.     const int b_stride= b_width;

  2051.     BlockNode *lt= &s->block[b_x + b_y*b_stride];

  2052.     BlockNode *rt= lt+1;

  2053.     BlockNode *lb= lt+b_stride;

  2054.     BlockNode *rb= lb+1;

  2055.     uint8_t *block[4]; 

  2056.     uint8_t tmp[src_stride*(b_h+5)]; //FIXME align

  2057.     int x,y;

  2058. 

  2059.     if(b_x<0){

  2060.         lt= rt;

  2061.         lb= rb;

  2062.     }else if(b_x + 1 >= b_width){

  2063.         rt= lt;

  2064.         rb= lb;

  2065.     }

  2066.     if(b_y<0){

  2067.         lt= lb;

  2068.         rt= rb;

  2069.     }else if(b_y + 1 >= b_height){

  2070.         lb= lt;

  2071.         rb= rt;

  2072.     }

  2073.         

  2074.     if(src_x<0){ //FIXME merge with prev & always round internal width upto *16

  2075.         obmc -= src_x;

  2076.         b_w += src_x;

  2077.         src_x=0;

  2078.     }else if(src_x + b_w > w){

  2079.         b_w = w - src_x;

  2080.     }

  2081.     if(src_y<0){

  2082.         obmc -= src_y*obmc_stride;

  2083.         b_h += src_y;

  2084.         src_y=0;

  2085.     }else if(src_y + b_h> h){

  2086.         b_h = h - src_y;

  2087.     }

  2088.     

  2089.     if(b_w<=0 || b_h<=0) return;

  2090. 

  2091. assert(src_stride > 7*MB_SIZE);

  2092.     dst += src_x + src_y*dst_stride;

  2093.     dst8+= src_x + src_y*src_stride;

  2094. //    src += src_x + src_y*src_stride;

  2095. 

  2096.     block[0]= tmp+3*MB_SIZE;

  2097.     pred_block(s, block[0], src, tmp, src_stride, src_x, src_y, b_w, b_h, lt, plane_index, w, h);    

  2098. 

  2099.     if(same_block(lt, rt)){

  2100.         block[1]= block[0];

  2101.     }else{

  2102.         block[1]= tmp + 4*MB_SIZE;

  2103.         pred_block(s, block[1], src, tmp, src_stride, src_x, src_y, b_w, b_h, rt, plane_index, w, h);

  2104.     }

  2105.         

  2106.     if(same_block(lt, lb)){

  2107.         block[2]= block[0];

  2108.     }else if(same_block(rt, lb)){

  2109.         block[2]= block[1];

  2110.     }else{

  2111.         block[2]= tmp+5*MB_SIZE;

  2112.         pred_block(s, block[2], src, tmp, src_stride, src_x, src_y, b_w, b_h, lb, plane_index, w, h);

  2113.     }

  2114. 

  2115.     if(same_block(lt, rb) ){

  2116.         block[3]= block[0];

  2117.     }else if(same_block(rt, rb)){

  2118.         block[3]= block[1];

  2119.     }else if(same_block(lb, rb)){

  2120.         block[3]= block[2];

  2121.     }else{

  2122.         block[3]= tmp+6*MB_SIZE;

  2123.         pred_block(s, block[3], src, tmp, src_stride, src_x, src_y, b_w, b_h, rb, plane_index, w, h);

  2124.     }

  2125. #if 0

  2126.     for(y=0; y<b_h; y++){

  2127.         for(x=0; x<b_w; x++){

  2128.             int v=   obmc [x + y*obmc_stride] * block[3][x + y*src_stride] * (256/OBMC_MAX);

  2129.             if(add) dst[x + y*dst_stride] += v;

  2130.             else    dst[x + y*dst_stride] -= v;

  2131.         }

  2132.     }

  2133.     for(y=0; y<b_h; y++){

  2134.         uint8_t *obmc2= obmc + (obmc_stride>>1);

  2135.         for(x=0; x<b_w; x++){

  2136.             int v=   obmc2[x + y*obmc_stride] * block[2][x + y*src_stride] * (256/OBMC_MAX);

  2137.             if(add) dst[x + y*dst_stride] += v;

  2138.             else    dst[x + y*dst_stride] -= v;

  2139.         }

  2140.     }

  2141.     for(y=0; y<b_h; y++){

  2142.         uint8_t *obmc3= obmc + obmc_stride*(obmc_stride>>1);

  2143.         for(x=0; x<b_w; x++){

  2144.             int v=   obmc3[x + y*obmc_stride] * block[1][x + y*src_stride] * (256/OBMC_MAX);

  2145.             if(add) dst[x + y*dst_stride] += v;

  2146.             else    dst[x + y*dst_stride] -= v;

  2147.         }

  2148.     }

  2149.     for(y=0; y<b_h; y++){

  2150.         uint8_t *obmc3= obmc + obmc_stride*(obmc_stride>>1);

  2151.         uint8_t *obmc4= obmc3+ (obmc_stride>>1);

  2152.         for(x=0; x<b_w; x++){

  2153.             int v=   obmc4[x + y*obmc_stride] * block[0][x + y*src_stride] * (256/OBMC_MAX);

  2154.             if(add) dst[x + y*dst_stride] += v;

  2155.             else    dst[x + y*dst_stride] -= v;

  2156.         }

  2157.     }

  2158. #else

  2159.     for(y=0; y<b_h; y++){

  2160.         //FIXME ugly missue of obmc_stride

  2161.         uint8_t *obmc1= obmc + y*obmc_stride;

  2162.         uint8_t *obmc2= obmc1+ (obmc_stride>>1);

  2163.         uint8_t *obmc3= obmc1+ obmc_stride*(obmc_stride>>1);

  2164.         uint8_t *obmc4= obmc3+ (obmc_stride>>1);

  2165.         for(x=0; x<b_w; x++){

  2166.             int v=   obmc1[x] * block[3][x + y*src_stride]

  2167.                     +obmc2[x] * block[2][x + y*src_stride]

  2168.                     +obmc3[x] * block[1][x + y*src_stride]

  2169.                     +obmc4[x] * block[0][x + y*src_stride];

  2170.             

  2171.             v <<= 8 - LOG2_OBMC_MAX;

  2172.             if(FRAC_BITS != 8){

  2173.                 v += 1<<(7 - FRAC_BITS);

  2174.                 v >>= 8 - FRAC_BITS;

  2175.             }

  2176.             if(add){

  2177.                 v += dst[x + y*dst_stride];

  2178.                 v = (v + (1<<(FRAC_BITS-1))) >> FRAC_BITS;

  2179.                 if(v&(~255)) v= ~(v>>31);

  2180.                 dst8[x + y*src_stride] = v;

  2181.             }else{

  2182.                 dst[x + y*dst_stride] -= v;

  2183.             }

  2184.         }

  2185.     }

  2186. #endif

  2187. }

  2188. 

  2189. static always_inline void predict_plane(SnowContext *s, DWTELEM *buf, int plane_index, int add){

  2190.     Plane *p= &s->plane[plane_index];

  2191.     const int mb_w= s->b_width  << s->block_max_depth;

  2192.     const int mb_h= s->b_height << s->block_max_depth;

  2193.     int x, y, mb_x, mb_y;

  2194.     int block_size = MB_SIZE >> s->block_max_depth;

  2195.     int block_w    = plane_index ? block_size/2 : block_size;

  2196.     const uint8_t *obmc  = plane_index ? obmc_tab[s->block_max_depth+1] : obmc_tab[s->block_max_depth];

  2197.     int obmc_stride= plane_index ? block_size : 2*block_size;

  2198.     int ref_stride= s->current_picture.linesize[plane_index];

  2199.     uint8_t *ref  = s->last_picture.data[plane_index];

  2200.     uint8_t *dst8= s->current_picture.data[plane_index];

  2201.     int w= p->width;

  2202.     int h= p->height;

  2203.     START_TIMER

  2204.     

  2205.     if(s->keyframe || (s->avctx->debug&512)){

  2206.         if(add){

  2207.             for(y=0; y<h; y++){

  2208.                 for(x=0; x<w; x++){

  2209.                     int v= buf[x + y*w] + (128<<FRAC_BITS) + (1<<(FRAC_BITS-1));

  2210.                     v >>= FRAC_BITS;

  2211.                     if(v&(~255)) v= ~(v>>31);

  2212.                     dst8[x + y*ref_stride]= v;

  2213.                 }

  2214.             }

  2215.         }else{

  2216.             for(y=0; y<h; y++){

  2217.                 for(x=0; x<w; x++){

  2218.                     buf[x + y*w]-= 128<<FRAC_BITS;

  2219.                 }

  2220.             }

  2221.         }

  2222. 

  2223.         return;

  2224.     }

  2225.     

  2226.     for(mb_y=0; mb_y<=mb_h; mb_y++){

  2227.         for(mb_x=0; mb_x<=mb_w; mb_x++){

  2228.             START_TIMER

  2229. 

  2230.             add_yblock(s, buf, dst8, ref, obmc, 

  2231.                        block_w*mb_x - block_w/2,

  2232.                        block_w*mb_y - block_w/2,

  2233.                        block_w, block_w,

  2234.                        w, h,

  2235.                        w, ref_stride, obmc_stride,

  2236.                        mb_x - 1, mb_y - 1,

  2237.                        add, plane_index);

  2238.             

  2239.             STOP_TIMER("add_yblock")

  2240.         }

  2241.     }

  2242.     

  2243.     STOP_TIMER("predict_plane")

  2244. }

  2245. 

  2246. static void quantize(SnowContext *s, SubBand *b, DWTELEM *src, int stride, int bias){

  2247.     const int level= b->level;

  2248.     const int w= b->width;

  2249.     const int h= b->height;

  2250.     const int qlog= clip(s->qlog + b->qlog, 0, 128);

  2251.     const int qmul= qexp[qlog&7]<<(qlog>>3);

  2252.     int x,y, thres1, thres2;

  2253.     START_TIMER

  2254. 

  2255.     assert(QROOT==8);

  2256. 

  2257.     if(s->qlog == LOSSLESS_QLOG) return;

  2258.  

  2259.     bias= bias ? 0 : (3*qmul)>>3;

  2260.     thres1= ((qmul - bias)>>QEXPSHIFT) - 1;

  2261.     thres2= 2*thres1;

  2262.     

  2263.     if(!bias){

  2264.         for(y=0; y<h; y++){

  2265.             for(x=0; x<w; x++){

  2266.                 int i= src[x + y*stride];

  2267.                 

  2268.                 if((unsigned)(i+thres1) > thres2){

  2269.                     if(i>=0){

  2270.                         i<<= QEXPSHIFT;

  2271.                         i/= qmul; //FIXME optimize

  2272.                         src[x + y*stride]=  i;

  2273.                     }else{

  2274.                         i= -i;

  2275.                         i<<= QEXPSHIFT;

  2276.                         i/= qmul; //FIXME optimize

  2277.                         src[x + y*stride]= -i;

  2278.                     }

  2279.                 }else

  2280.                     src[x + y*stride]= 0;

  2281.             }

  2282.         }

  2283.     }else{

  2284.         for(y=0; y<h; y++){

  2285.             for(x=0; x<w; x++){

  2286.                 int i= src[x + y*stride]; 

  2287.                 

  2288.                 if((unsigned)(i+thres1) > thres2){

  2289.                     if(i>=0){

  2290.                         i<<= QEXPSHIFT;

  2291.                         i= (i + bias) / qmul; //FIXME optimize

  2292.                         src[x + y*stride]=  i;

  2293.                     }else{

  2294.                         i= -i;

  2295.                         i<<= QEXPSHIFT;

  2296.                         i= (i + bias) / qmul; //FIXME optimize

  2297.                         src[x + y*stride]= -i;

  2298.                     }

  2299.                 }else

  2300.                     src[x + y*stride]= 0;

  2301.             }

  2302.         }

  2303.     }

  2304.     if(level+1 == s->spatial_decomposition_count){

  2305. //        STOP_TIMER("quantize")

  2306.     }

  2307. }

  2308. 

  2309. static void dequantize(SnowContext *s, SubBand *b, DWTELEM *src, int stride){

  2310.     const int w= b->width;

  2311.     const int h= b->height;

  2312.     const int qlog= clip(s->qlog + b->qlog, 0, 128);

  2313.     const int qmul= qexp[qlog&7]<<(qlog>>3);

  2314.     const int qadd= (s->qbias*qmul)>>QBIAS_SHIFT;

  2315.     int x,y;

  2316.     START_TIMER

  2317.     

  2318.     if(s->qlog == LOSSLESS_QLOG) return;

  2319.     

  2320.     assert(QROOT==8);

  2321. 

  2322.     for(y=0; y<h; y++){

  2323.         for(x=0; x<w; x++){

  2324.             int i= src[x + y*stride];

  2325.             if(i<0){

  2326.                 src[x + y*stride]= -((-i*qmul + qadd)>>(QEXPSHIFT)); //FIXME try different bias

  2327.             }else if(i>0){

  2328.                 src[x + y*stride]=  (( i*qmul + qadd)>>(QEXPSHIFT));

  2329.             }

  2330.         }

  2331.     }

  2332.     if(w > 200 /*level+1 == s->spatial_decomposition_count*/){

  2333.         STOP_TIMER("dquant")

  2334.     }

  2335. }

  2336. 

  2337. static void decorrelate(SnowContext *s, SubBand *b, DWTELEM *src, int stride, int inverse, int use_median){

  2338.     const int w= b->width;

  2339.     const int h= b->height;

  2340.     int x,y;

  2341.     

  2342.     for(y=h-1; y>=0; y--){

  2343.         for(x=w-1; x>=0; x--){

  2344.             int i= x + y*stride;

  2345.             

  2346.             if(x){

  2347.                 if(use_median){

  2348.                     if(y && x+1<w) src[i] -= mid_pred(src[i - 1], src[i - stride], src[i - stride + 1]);

  2349.                     else  src[i] -= src[i - 1];

  2350.                 }else{

  2351.                     if(y) src[i] -= mid_pred(src[i - 1], src[i - stride], src[i - 1] + src[i - stride] - src[i - 1 - stride]);

  2352.                     else  src[i] -= src[i - 1];

  2353.                 }

  2354.             }else{

  2355.                 if(y) src[i] -= src[i - stride];

  2356.             }

  2357.         }

  2358.     }

  2359. }

  2360. 

  2361. static void correlate(SnowContext *s, SubBand *b, DWTELEM *src, int stride, int inverse, int use_median){

  2362.     const int w= b->width;

  2363.     const int h= b->height;

  2364.     int x,y;

  2365.     

  2366.     for(y=0; y<h; y++){

  2367.         for(x=0; x<w; x++){

  2368.             int i= x + y*stride;

  2369.             

  2370.             if(x){

  2371.                 if(use_median){

  2372.                     if(y && x+1<w) src[i] += mid_pred(src[i - 1], src[i - stride], src[i - stride + 1]);

  2373.                     else  src[i] += src[i - 1];

  2374.                 }else{

  2375.                     if(y) src[i] += mid_pred(src[i - 1], src[i - stride], src[i - 1] + src[i - stride] - src[i - 1 - stride]);

  2376.                     else  src[i] += src[i - 1];

  2377.                 }

  2378.             }else{

  2379.                 if(y) src[i] += src[i - stride];

  2380.             }

  2381.         }

  2382.     }

  2383. }

  2384. 

  2385. static void encode_header(SnowContext *s){

  2386.     int plane_index, level, orientation;

  2387.     uint8_t kstate[32]; 

  2388.     

  2389.     memset(kstate, MID_STATE, sizeof(kstate));   

  2390. 

  2391.     put_rac(&s->c, kstate, s->keyframe);

  2392.     if(s->keyframe || s->always_reset)

  2393.         reset_contexts(s);

  2394.     if(s->keyframe){

  2395.         put_symbol(&s->c, s->header_state, s->version, 0);

  2396.         put_rac(&s->c, s->header_state, s->always_reset);

  2397.         put_symbol(&s->c, s->header_state, s->temporal_decomposition_type, 0);

  2398.         put_symbol(&s->c, s->header_state, s->temporal_decomposition_count, 0);

  2399.         put_symbol(&s->c, s->header_state, s->spatial_decomposition_count, 0);

  2400.         put_symbol(&s->c, s->header_state, s->colorspace_type, 0);

  2401.         put_symbol(&s->c, s->header_state, s->chroma_h_shift, 0);

  2402.         put_symbol(&s->c, s->header_state, s->chroma_v_shift, 0);

  2403.         put_rac(&s->c, s->header_state, s->spatial_scalability);

  2404. //        put_rac(&s->c, s->header_state, s->rate_scalability);

  2405. 

  2406.         for(plane_index=0; plane_index<2; plane_index++){

  2407.             for(level=0; level<s->spatial_decomposition_count; level++){

  2408.                 for(orientation=level ? 1:0; orientation<4; orientation++){

  2409.                     if(orientation==2) continue;

  2410.                     put_symbol(&s->c, s->header_state, s->plane[plane_index].band[level][orientation].qlog, 1);

  2411.                 }

  2412.             }

  2413.         }

  2414.     }

  2415.     put_symbol(&s->c, s->header_state, s->spatial_decomposition_type, 0);

  2416.     put_symbol(&s->c, s->header_state, s->qlog, 1); 

  2417.     put_symbol(&s->c, s->header_state, s->mv_scale, 0); 

  2418.     put_symbol(&s->c, s->header_state, s->qbias, 1);

  2419.     put_symbol(&s->c, s->header_state, s->block_max_depth, 0);

  2420. }

  2421. 

  2422. static int decode_header(SnowContext *s){

  2423.     int plane_index, level, orientation;

  2424.     uint8_t kstate[32];

  2425. 

  2426.     memset(kstate, MID_STATE, sizeof(kstate));   

  2427. 

  2428.     s->keyframe= get_rac(&s->c, kstate);

  2429.     if(s->keyframe || s->always_reset)

  2430.         reset_contexts(s);

  2431.     if(s->keyframe){

  2432.         s->version= get_symbol(&s->c, s->header_state, 0);

  2433.         if(s->version>0){

  2434.             av_log(s->avctx, AV_LOG_ERROR, "version %d not supported", s->version);

  2435.             return -1;

  2436.         }

  2437.         s->always_reset= get_rac(&s->c, s->header_state);

  2438.         s->temporal_decomposition_type= get_symbol(&s->c, s->header_state, 0);

  2439.         s->temporal_decomposition_count= get_symbol(&s->c, s->header_state, 0);

  2440.         s->spatial_decomposition_count= get_symbol(&s->c, s->header_state, 0);

  2441.         s->colorspace_type= get_symbol(&s->c, s->header_state, 0);

  2442.         s->chroma_h_shift= get_symbol(&s->c, s->header_state, 0);

  2443.         s->chroma_v_shift= get_symbol(&s->c, s->header_state, 0);

  2444.         s->spatial_scalability= get_rac(&s->c, s->header_state);

  2445. //        s->rate_scalability= get_rac(&s->c, s->header_state);

  2446. 

  2447.         for(plane_index=0; plane_index<3; plane_index++){

  2448.             for(level=0; level<s->spatial_decomposition_count; level++){

  2449.                 for(orientation=level ? 1:0; orientation<4; orientation++){

  2450.                     int q;

  2451.                     if     (plane_index==2) q= s->plane[1].band[level][orientation].qlog;

  2452.                     else if(orientation==2) q= s->plane[plane_index].band[level][1].qlog;

  2453.                     else                    q= get_symbol(&s->c, s->header_state, 1);

  2454.                     s->plane[plane_index].band[level][orientation].qlog= q;

  2455.                 }

  2456.             }

  2457.         }

  2458.     }

  2459.     

  2460.     s->spatial_decomposition_type= get_symbol(&s->c, s->header_state, 0);

  2461.     if(s->spatial_decomposition_type > 2){

  2462.         av_log(s->avctx, AV_LOG_ERROR, "spatial_decomposition_type %d not supported", s->spatial_decomposition_type);

  2463.         return -1;

  2464.     }

  2465.     

  2466.     s->qlog= get_symbol(&s->c, s->header_state, 1);

  2467.     s->mv_scale= get_symbol(&s->c, s->header_state, 0);

  2468.     s->qbias= get_symbol(&s->c, s->header_state, 1);

  2469.     s->block_max_depth= get_symbol(&s->c, s->header_state, 0);

  2470. 

  2471.     return 0;

  2472. }

  2473. 

  2474. static int common_init(AVCodecContext *avctx){

  2475.     SnowContext *s = avctx->priv_data;

  2476.     int width, height;

  2477.     int level, orientation, plane_index, dec;

  2478. 

  2479.     s->avctx= avctx;

  2480.         

  2481.     dsputil_init(&s->dsp, avctx);

  2482. 

  2483. #define mcf(dx,dy)\

  2484.     s->dsp.put_qpel_pixels_tab       [0][dy+dx/4]=\

  2485.     s->dsp.put_no_rnd_qpel_pixels_tab[0][dy+dx/4]=\

  2486.         s->dsp.put_h264_qpel_pixels_tab[0][dy+dx/4];\

  2487.     s->dsp.put_qpel_pixels_tab       [1][dy+dx/4]=\

  2488.     s->dsp.put_no_rnd_qpel_pixels_tab[1][dy+dx/4]=\

  2489.         s->dsp.put_h264_qpel_pixels_tab[1][dy+dx/4];

  2490. 

  2491.     mcf( 0, 0)

  2492.     mcf( 4, 0)

  2493.     mcf( 8, 0)

  2494.     mcf(12, 0)

  2495.     mcf( 0, 4)

  2496.     mcf( 4, 4)

  2497.     mcf( 8, 4)

  2498.     mcf(12, 4)

  2499.     mcf( 0, 8)

  2500.     mcf( 4, 8)

  2501.     mcf( 8, 8)

  2502.     mcf(12, 8)

  2503.     mcf( 0,12)

  2504.     mcf( 4,12)

  2505.     mcf( 8,12)

  2506.     mcf(12,12)

  2507. 

  2508. #define mcfh(dx,dy)\

  2509.     s->dsp.put_pixels_tab       [0][dy/4+dx/8]=\

  2510.     s->dsp.put_no_rnd_pixels_tab[0][dy/4+dx/8]=\

  2511.         mc_block_hpel ## dx ## dy ## 16;\

  2512.     s->dsp.put_pixels_tab       [1][dy/4+dx/8]=\

  2513.     s->dsp.put_no_rnd_pixels_tab[1][dy/4+dx/8]=\

  2514.         mc_block_hpel ## dx ## dy ## 8;

  2515. 

  2516.     mcfh(0, 0)

  2517.     mcfh(8, 0)

  2518.     mcfh(0, 8)

  2519.     mcfh(8, 8)

  2520.         

  2521.     dec= s->spatial_decomposition_count= 5;

  2522.     s->spatial_decomposition_type= avctx->prediction_method; //FIXME add decorrelator type r transform_type

  2523.     

  2524.     s->chroma_h_shift= 1; //FIXME XXX

  2525.     s->chroma_v_shift= 1;

  2526.     

  2527. //    dec += FFMAX(s->chroma_h_shift, s->chroma_v_shift);

  2528.     

  2529.     width= s->avctx->width;

  2530.     height= s->avctx->height;

  2531. 

  2532.     s->spatial_dwt_buffer= av_mallocz(width*height*sizeof(DWTELEM));

  2533.     

  2534.     s->mv_scale= (s->avctx->flags & CODEC_FLAG_QPEL) ? 2 : 4;

  2535.     s->block_max_depth= (s->avctx->flags & CODEC_FLAG_4MV) ? 1 : 0;

  2536.     

  2537.     for(plane_index=0; plane_index<3; plane_index++){    

  2538.         int w= s->avctx->width;

  2539.         int h= s->avctx->height;

  2540. 

  2541.         if(plane_index){

  2542.             w>>= s->chroma_h_shift;

  2543.             h>>= s->chroma_v_shift;

  2544.         }

  2545.         s->plane[plane_index].width = w;

  2546.         s->plane[plane_index].height= h;

  2547. //av_log(NULL, AV_LOG_DEBUG, "%d %d\n", w, h);

  2548.         for(level=s->spatial_decomposition_count-1; level>=0; level--){

  2549.             for(orientation=level ? 1 : 0; orientation<4; orientation++){

  2550.                 SubBand *b= &s->plane[plane_index].band[level][orientation];

  2551.                 

  2552.                 b->buf= s->spatial_dwt_buffer;

  2553.                 b->level= level;

  2554.                 b->stride= s->plane[plane_index].width << (s->spatial_decomposition_count - level);

  2555.                 b->width = (w + !(orientation&1))>>1;

  2556.                 b->height= (h + !(orientation>1))>>1;

  2557.                 

  2558.                 if(orientation&1) b->buf += (w+1)>>1;

  2559.                 if(orientation>1) b->buf += b->stride>>1;

  2560.                 

  2561.                 if(level)

  2562.                     b->parent= &s->plane[plane_index].band[level-1][orientation];

  2563.                 b->x    = av_mallocz(((b->width+1) * b->height+1)*sizeof(int16_t));

  2564.                 b->coeff= av_mallocz(((b->width+1) * b->height+1)*sizeof(DWTELEM));

  2565.             }

  2566.             w= (w+1)>>1;

  2567.             h= (h+1)>>1;

  2568.         }

  2569.     }

  2570.     

  2571.     reset_contexts(s);

  2572. /*    

  2573.     width= s->width= avctx->width;

  2574.     height= s->height= avctx->height;

  2575.     

  2576.     assert(width && height);

  2577. */

  2578.     s->avctx->get_buffer(s->avctx, &s->mconly_picture);

  2579.     

  2580.     return 0;

  2581. }

  2582. 

  2583. 

  2584. static void calculate_vissual_weight(SnowContext *s, Plane *p){

  2585.     int width = p->width;

  2586.     int height= p->height;

  2587.     int level, orientation, x, y;

  2588. 

  2589.     for(level=0; level<s->spatial_decomposition_count; level++){

  2590.         for(orientation=level ? 1 : 0; orientation<4; orientation++){

  2591.             SubBand *b= &p->band[level][orientation];

  2592.             DWTELEM *buf= b->buf;

  2593.             int64_t error=0;

  2594.             

  2595.             memset(s->spatial_dwt_buffer, 0, sizeof(int)*width*height);

  2596.             buf[b->width/2 + b->height/2*b->stride]= 256*256;

  2597.             ff_spatial_idwt(s->spatial_dwt_buffer, width, height, width, s->spatial_decomposition_type, s->spatial_decomposition_count);

  2598.             for(y=0; y<height; y++){

  2599.                 for(x=0; x<width; x++){

  2600.                     int64_t d= s->spatial_dwt_buffer[x + y*width];

  2601.                     error += d*d;

  2602.                 }

  2603.             }

  2604. 

  2605.             b->qlog= (int)(log(352256.0/sqrt(error)) / log(pow(2.0, 1.0/QROOT))+0.5);

  2606. //            av_log(NULL, AV_LOG_DEBUG, "%d %d %d\n", level, orientation, b->qlog/*, sqrt(error)*/);

  2607.         }

  2608.     }

  2609. }

  2610. 

  2611. static int encode_init(AVCodecContext *avctx)

  2612. {

  2613.     SnowContext *s = avctx->priv_data;

  2614.     int plane_index;

  2615. 

  2616.     if(avctx->strict_std_compliance >= 0){

  2617.         av_log(avctx, AV_LOG_ERROR, "this codec is under development, files encoded with it wont be decodeable with future versions!!!\n"

  2618.                "use vstrict=-1 to use it anyway\n");

  2619.         return -1;

  2620.     }

  2621.  

  2622.     common_init(avctx);

  2623.     alloc_blocks(s);

  2624.  

  2625.     s->version=0;

  2626.     

  2627.     s->m.me.scratchpad= av_mallocz((avctx->width+64)*2*16*2*sizeof(uint8_t));

  2628.     s->m.me.map       = av_mallocz(ME_MAP_SIZE*sizeof(uint32_t));

  2629.     s->m.me.score_map = av_mallocz(ME_MAP_SIZE*sizeof(uint32_t));

  2630.     h263_encode_init(&s->m); //mv_penalty

  2631. 

  2632.     for(plane_index=0; plane_index<3; plane_index++){

  2633.         calculate_vissual_weight(s, &s->plane[plane_index]);

  2634.     }

  2635.     

  2636.     

  2637.     avctx->coded_frame= &s->current_picture;

  2638.     switch(avctx->pix_fmt){

  2639. //    case PIX_FMT_YUV444P:

  2640. //    case PIX_FMT_YUV422P:

  2641.     case PIX_FMT_YUV420P:

  2642.     case PIX_FMT_GRAY8:

  2643. //    case PIX_FMT_YUV411P:

  2644. //    case PIX_FMT_YUV410P:

  2645.         s->colorspace_type= 0;

  2646.         break;

  2647. /*    case PIX_FMT_RGBA32:

  2648.         s->colorspace= 1;

  2649.         break;*/

  2650.     default:

  2651.         av_log(avctx, AV_LOG_ERROR, "format not supported\n");

  2652.         return -1;

  2653.     }

  2654. //    avcodec_get_chroma_sub_sample(avctx->pix_fmt, &s->chroma_h_shift, &s->chroma_v_shift);

  2655.     s->chroma_h_shift= 1;

  2656.     s->chroma_v_shift= 1;

  2657.     return 0;

  2658. }

  2659. 

  2660. static int frame_start(SnowContext *s){

  2661.    AVFrame tmp;

  2662.    int w= s->avctx->width; //FIXME round up to x16 ?

  2663.    int h= s->avctx->height;

  2664. 

  2665.     if(s->current_picture.data[0]){

  2666.         draw_edges(s->current_picture.data[0], s->current_picture.linesize[0], w   , h   , EDGE_WIDTH  );

  2667.         draw_edges(s->current_picture.data[1], s->current_picture.linesize[1], w>>1, h>>1, EDGE_WIDTH/2);

  2668.         draw_edges(s->current_picture.data[2], s->current_picture.linesize[2], w>>1, h>>1, EDGE_WIDTH/2);

  2669.     }

  2670. 

  2671.     tmp= s->last_picture;

  2672.     s->last_picture= s->current_picture;

  2673.     s->current_picture= tmp;

  2674.     

  2675.     s->current_picture.reference= 1;

  2676.     if(s->avctx->get_buffer(s->avctx, &s->current_picture) < 0){

  2677.         av_log(s->avctx, AV_LOG_ERROR, "get_buffer() failed\n");

  2678.         return -1;

  2679.     }

  2680.     

  2681.     return 0;

  2682. }

  2683. 

  2684. static int encode_frame(AVCodecContext *avctx, unsigned char *buf, int buf_size, void *data){

  2685.     SnowContext *s = avctx->priv_data;

  2686.     RangeCoder * const c= &s->c;

  2687.     AVFrame *pict = data;

  2688.     const int width= s->avctx->width;

  2689.     const int height= s->avctx->height;

  2690.     int level, orientation, plane_index;

  2691. 

  2692.     ff_init_range_encoder(c, buf, buf_size);

  2693.     ff_build_rac_states(c, 0.05*(1LL<<32), 256-8);

  2694.     

  2695.     s->input_picture = *pict;

  2696. 

  2697.     s->keyframe=avctx->gop_size==0 || avctx->frame_number % avctx->gop_size == 0;

  2698.     pict->pict_type= s->keyframe ? FF_I_TYPE : FF_P_TYPE;

  2699.     

  2700.     if(pict->quality){

  2701.         s->qlog= rint(QROOT*log(pict->quality / (float)FF_QP2LAMBDA)/log(2));

  2702.         //<64 >60

  2703.         s->qlog += 61;

  2704.     }else{

  2705.         s->qlog= LOSSLESS_QLOG;

  2706.     }

  2707. 

  2708.     frame_start(s);

  2709.     s->current_picture.key_frame= s->keyframe;

  2710. 

  2711.     if(pict->pict_type == P_TYPE){

  2712.         int block_width = (width +15)>>4;

  2713.         int block_height= (height+15)>>4;

  2714.         int stride= s->current_picture.linesize[0];

  2715.         

  2716.         assert(s->current_picture.data[0]);

  2717.         assert(s->last_picture.data[0]);

  2718.      

  2719.         s->m.avctx= s->avctx;

  2720.         s->m.current_picture.data[0]= s->current_picture.data[0];

  2721.         s->m.   last_picture.data[0]= s->   last_picture.data[0];

  2722.         s->m.    new_picture.data[0]= s->  input_picture.data[0];

  2723.         s->m.current_picture_ptr= &s->m.current_picture;

  2724.         s->m.   last_picture_ptr= &s->m.   last_picture;

  2725.         s->m.linesize=

  2726.         s->m.   last_picture.linesize[0]=

  2727.         s->m.    new_picture.linesize[0]=

  2728.         s->m.current_picture.linesize[0]= stride;

  2729.         s->m.uvlinesize= s->current_picture.linesize[1];

  2730.         s->m.width = width;

  2731.         s->m.height= height;

  2732.         s->m.mb_width = block_width;

  2733.         s->m.mb_height= block_height;

  2734.         s->m.mb_stride=   s->m.mb_width+1;

  2735.         s->m.b8_stride= 2*s->m.mb_width+1;

  2736.         s->m.f_code=1;

  2737.         s->m.pict_type= pict->pict_type;

  2738.         s->m.me_method= s->avctx->me_method;

  2739.         s->m.me.scene_change_score=0;

  2740.         s->m.flags= s->avctx->flags;

  2741.         s->m.quarter_sample= (s->avctx->flags & CODEC_FLAG_QPEL)!=0;

  2742.         s->m.out_format= FMT_H263;

  2743.         s->m.unrestricted_mv= 1;

  2744. 

  2745.         s->lambda = s->m.lambda= pict->quality * 3/2; //FIXME bug somewhere else

  2746.         s->m.qscale= (s->m.lambda*139 + FF_LAMBDA_SCALE*64) >> (FF_LAMBDA_SHIFT + 7);

  2747.         s->lambda2= s->m.lambda2= (s->m.lambda*s->m.lambda + FF_LAMBDA_SCALE/2) >> FF_LAMBDA_SHIFT;

  2748. 

  2749.         s->m.dsp= s->dsp; //move

  2750.         ff_init_me(&s->m);

  2751.     }

  2752.     

  2753. redo_frame:

  2754.         

  2755.     s->qbias= pict->pict_type == P_TYPE ? 2 : 0;

  2756. 

  2757.     encode_header(s);

  2758.     encode_blocks(s);

  2759.       

  2760.     for(plane_index=0; plane_index<3; plane_index++){

  2761.         Plane *p= &s->plane[plane_index];

  2762.         int w= p->width;

  2763.         int h= p->height;

  2764.         int x, y;

  2765. //        int bits= put_bits_count(&s->c.pb);

  2766. 

  2767.         //FIXME optimize

  2768.      if(pict->data[plane_index]) //FIXME gray hack

  2769.         for(y=0; y<h; y++){

  2770.             for(x=0; x<w; x++){

  2771.                 s->spatial_dwt_buffer[y*w + x]= pict->data[plane_index][y*pict->linesize[plane_index] + x]<<FRAC_BITS;

  2772.             }

  2773.         }

  2774.         predict_plane(s, s->spatial_dwt_buffer, plane_index, 0);

  2775.         

  2776.         if(   plane_index==0 

  2777.            && pict->pict_type == P_TYPE 

  2778.            && s->m.me.scene_change_score > s->avctx->scenechange_threshold){

  2779.             ff_init_range_encoder(c, buf, buf_size);

  2780.             ff_build_rac_states(c, 0.05*(1LL<<32), 256-8);

  2781.             pict->pict_type= FF_I_TYPE;

  2782.             s->keyframe=1;

  2783.             reset_contexts(s);

  2784.             goto redo_frame;

  2785.         }

  2786.         

  2787.         if(s->qlog == LOSSLESS_QLOG){

  2788.             for(y=0; y<h; y++){

  2789.                 for(x=0; x<w; x++){

  2790.                     s->spatial_dwt_buffer[y*w + x]= (s->spatial_dwt_buffer[y*w + x] + (1<<(FRAC_BITS-1)))>>FRAC_BITS;

  2791.                 }

  2792.             }

  2793.         }

  2794.  

  2795.         ff_spatial_dwt(s->spatial_dwt_buffer, w, h, w, s->spatial_decomposition_type, s->spatial_decomposition_count);

  2796. 

  2797.         for(level=0; level<s->spatial_decomposition_count; level++){

  2798.             for(orientation=level ? 1 : 0; orientation<4; orientation++){

  2799.                 SubBand *b= &p->band[level][orientation];

  2800.                 

  2801.                 quantize(s, b, b->buf, b->stride, s->qbias);

  2802.                 if(orientation==0)

  2803.                     decorrelate(s, b, b->buf, b->stride, pict->pict_type == P_TYPE, 0);

  2804.                 encode_subband(s, b, b->buf, b->parent ? b->parent->buf : NULL, b->stride, orientation);

  2805.                 assert(b->parent==NULL || b->parent->stride == b->stride*2);

  2806.                 if(orientation==0)

  2807.                     correlate(s, b, b->buf, b->stride, 1, 0);

  2808.             }

  2809.         }

  2810. //        av_log(NULL, AV_LOG_DEBUG, "plane:%d bits:%d\n", plane_index, put_bits_count(&s->c.pb) - bits);

  2811. 

  2812.         for(level=0; level<s->spatial_decomposition_count; level++){

  2813.             for(orientation=level ? 1 : 0; orientation<4; orientation++){

  2814.                 SubBand *b= &p->band[level][orientation];

  2815. 

  2816.                 dequantize(s, b, b->buf, b->stride);

  2817.             }

  2818.         }

  2819. 

  2820.         ff_spatial_idwt(s->spatial_dwt_buffer, w, h, w, s->spatial_decomposition_type, s->spatial_decomposition_count);

  2821.         if(s->qlog == LOSSLESS_QLOG){

  2822.             for(y=0; y<h; y++){

  2823.                 for(x=0; x<w; x++){

  2824.                     s->spatial_dwt_buffer[y*w + x]<<=FRAC_BITS;

  2825.                 }

  2826.             }

  2827.         }

  2828. {START_TIMER

  2829.         predict_plane(s, s->spatial_dwt_buffer, plane_index, 1);

  2830. STOP_TIMER("pred-conv")}

  2831.         if(s->avctx->flags&CODEC_FLAG_PSNR){

  2832.             int64_t error= 0;

  2833.             

  2834.     if(pict->data[plane_index]) //FIXME gray hack

  2835.             for(y=0; y<h; y++){

  2836.                 for(x=0; x<w; x++){

  2837.                     int d= s->current_picture.data[plane_index][y*s->current_picture.linesize[plane_index] + x] - pict->data[plane_index][y*pict->linesize[plane_index] + x];

  2838.                     error += d*d;

  2839.                 }

  2840.             }

  2841.             s->avctx->error[plane_index] += error;

  2842.             s->current_picture.error[plane_index] = error;

  2843.         }

  2844.     }

  2845. 

  2846.     if(s->last_picture.data[0])

  2847.         avctx->release_buffer(avctx, &s->last_picture);

  2848. 

  2849.     emms_c();

  2850.     

  2851.     return ff_rac_terminate(c);

  2852. }

  2853. 

  2854. static void common_end(SnowContext *s){

  2855.     int plane_index, level, orientation;

  2856. 

  2857.     av_freep(&s->spatial_dwt_buffer);

  2858. 

  2859.     av_freep(&s->m.me.scratchpad);    

  2860.     av_freep(&s->m.me.map);

  2861.     av_freep(&s->m.me.score_map);

  2862.  

  2863.     av_freep(&s->block);

  2864. 

  2865.     for(plane_index=0; plane_index<3; plane_index++){    

  2866.         for(level=s->spatial_decomposition_count-1; level>=0; level--){

  2867.             for(orientation=level ? 1 : 0; orientation<4; orientation++){

  2868.                 SubBand *b= &s->plane[plane_index].band[level][orientation];

  2869.                 

  2870.                 av_freep(&b->x);

  2871.                 av_freep(&b->coeff);

  2872.             }

  2873.         }

  2874.     }

  2875. }

  2876. 

  2877. static int encode_end(AVCodecContext *avctx)

  2878. {

  2879.     SnowContext *s = avctx->priv_data;

  2880. 

  2881.     common_end(s);

  2882. 

  2883.     return 0;

  2884. }

  2885. 

  2886. static int decode_init(AVCodecContext *avctx)

  2887. {

  2888. //    SnowContext *s = avctx->priv_data;

  2889. 

  2890.     common_init(avctx);

  2891.     

  2892.     return 0;

  2893. }

  2894. 

  2895. static int decode_frame(AVCodecContext *avctx, void *data, int *data_size, uint8_t *buf, int buf_size){

  2896.     SnowContext *s = avctx->priv_data;

  2897.     RangeCoder * const c= &s->c;

  2898.     int bytes_read;

  2899.     AVFrame *picture = data;

  2900.     int level, orientation, plane_index;

  2901.     

  2902. 

  2903.     /* no supplementary picture */

  2904.     if (buf_size == 0)

  2905.         return 0;

  2906. 

  2907.     ff_init_range_decoder(c, buf, buf_size);

  2908.     ff_build_rac_states(c, 0.05*(1LL<<32), 256-8);

  2909. 

  2910.     s->current_picture.pict_type= FF_I_TYPE; //FIXME I vs. P

  2911.     decode_header(s);

  2912.     if(!s->block) alloc_blocks(s);

  2913. 

  2914.     frame_start(s);

  2915.     //keyframe flag dupliaction mess FIXME

  2916.     if(avctx->debug&FF_DEBUG_PICT_INFO)

  2917.         av_log(avctx, AV_LOG_ERROR, "keyframe:%d qlog:%d\n", s->keyframe, s->qlog);

  2918.     

  2919.     decode_blocks(s);

  2920. 

  2921.     for(plane_index=0; plane_index<3; plane_index++){

  2922.         Plane *p= &s->plane[plane_index];

  2923.         int w= p->width;

  2924.         int h= p->height;

  2925.         int x, y;

  2926.         

  2927. if(s->avctx->debug&2048){

  2928.         memset(s->spatial_dwt_buffer, 0, sizeof(DWTELEM)*w*h);

  2929.         predict_plane(s, s->spatial_dwt_buffer, plane_index, 1);

  2930. 

  2931.         for(y=0; y<h; y++){

  2932.             for(x=0; x<w; x++){

  2933.                 int v= s->current_picture.data[plane_index][y*s->current_picture.linesize[plane_index] + x];

  2934.                 s->mconly_picture.data[plane_index][y*s->mconly_picture.linesize[plane_index] + x]= v;

  2935.             }

  2936.         }

  2937. }

  2938.         for(level=0; level<s->spatial_decomposition_count; level++){

  2939.             for(orientation=level ? 1 : 0; orientation<4; orientation++){

  2940.                 SubBand *b= &p->band[level][orientation];

  2941. 

  2942.                 decode_subband(s, b, b->buf, b->parent ? b->parent->buf : NULL, b->stride, orientation);

  2943.                 if(orientation==0){

  2944.                     correlate(s, b, b->buf, b->stride, 1, 0);

  2945.                     dequantize(s, b, b->buf, b->stride);

  2946.                     assert(b->buf == s->spatial_dwt_buffer);

  2947.                 }

  2948.             }

  2949.         }

  2950. 

  2951.         ff_spatial_idwt(s->spatial_dwt_buffer, w, h, w, s->spatial_decomposition_type, s->spatial_decomposition_count);

  2952.         if(s->qlog == LOSSLESS_QLOG){

  2953.             for(y=0; y<h; y++){

  2954.                 for(x=0; x<w; x++){

  2955.                     s->spatial_dwt_buffer[y*w + x]<<=FRAC_BITS;

  2956.                 }

  2957.             }

  2958.         }

  2959. {START_TIMER

  2960.         predict_plane(s, s->spatial_dwt_buffer, plane_index, 1);

  2961. STOP_TIMER("predict_plane conv2")}

  2962.     }

  2963.             

  2964.     emms_c();

  2965. 

  2966.     if(s->last_picture.data[0])

  2967.         avctx->release_buffer(avctx, &s->last_picture);

  2968. 

  2969. if(!(s->avctx->debug&2048))        

  2970.     *picture= s->current_picture;

  2971. else

  2972.     *picture= s->mconly_picture;

  2973.     

  2974.     *data_size = sizeof(AVFrame);

  2975.     

  2976.     bytes_read= c->bytestream - c->bytestream_start;

  2977.     if(bytes_read ==0) av_log(s->avctx, AV_LOG_ERROR, "error at end of frame\n"); //FIXME

  2978. 

  2979.     return bytes_read;

  2980. }

  2981. 

  2982. static int decode_end(AVCodecContext *avctx)

  2983. {

  2984.     SnowContext *s = avctx->priv_data;

  2985. 

  2986.     common_end(s);

  2987. 

  2988.     return 0;

  2989. }

  2990. 

  2991. AVCodec snow_decoder = {

  2992.     "snow",

  2993.     CODEC_TYPE_VIDEO,

  2994.     CODEC_ID_SNOW,

  2995.     sizeof(SnowContext),

  2996.     decode_init,

  2997.     NULL,

  2998.     decode_end,

  2999.     decode_frame,

  3000.     0 /*CODEC_CAP_DR1*/ /*| CODEC_CAP_DRAW_HORIZ_BAND*/,

  3001.     NULL

  3002. };

  3003. 

  3004. AVCodec snow_encoder = {

  3005.     "snow",

  3006.     CODEC_TYPE_VIDEO,

  3007.     CODEC_ID_SNOW,

  3008.     sizeof(SnowContext),

  3009.     encode_init,

  3010.     encode_frame,

  3011.     encode_end,

  3012. };

  3013. 

  3014. 

  3015. #if 0

  3016. #undef malloc

  3017. #undef free

  3018. #undef printf

  3019. 

  3020. int main(){

  3021.     int width=256;

  3022.     int height=256;

  3023.     int buffer[2][width*height];

  3024.     SnowContext s;

  3025.     int i;

  3026.     s.spatial_decomposition_count=6;

  3027.     s.spatial_decomposition_type=1;

  3028.     

  3029.     printf("testing 5/3 DWT\n");

  3030.     for(i=0; i<width*height; i++)

  3031.         buffer[0][i]= buffer[1][i]= random()%54321 - 12345;

  3032.     

  3033.     ff_spatial_dwt(buffer[0], width, height, width, s->spatial_decomposition_type, s->spatial_decomposition_count);

  3034.     ff_spatial_idwt(buffer[0], width, height, width, s->spatial_decomposition_type, s->spatial_decomposition_count);

  3035.     

  3036.     for(i=0; i<width*height; i++)

  3037.         if(buffer[0][i]!= buffer[1][i]) printf("fsck: %d %d %d\n",i, buffer[0][i], buffer[1][i]);

  3038. 

  3039.     printf("testing 9/7 DWT\n");

  3040.     s.spatial_decomposition_type=0;

  3041.     for(i=0; i<width*height; i++)

  3042.         buffer[0][i]= buffer[1][i]= random()%54321 - 12345;

  3043.     

  3044.     ff_spatial_dwt(buffer[0], width, height, width, s->spatial_decomposition_type, s->spatial_decomposition_count);

  3045.     ff_spatial_idwt(buffer[0], width, height, width, s->spatial_decomposition_type, s->spatial_decomposition_count);

  3046.     

  3047.     for(i=0; i<width*height; i++)

  3048.         if(buffer[0][i]!= buffer[1][i]) printf("fsck: %d %d %d\n",i, buffer[0][i], buffer[1][i]);

  3049.         

  3050.     printf("testing AC coder\n");

  3051.     memset(s.header_state, 0, sizeof(s.header_state));

  3052.     ff_init_range_encoder(&s.c, buffer[0], 256*256);

  3053.     ff_init_cabac_states(&s.c, ff_h264_lps_range, ff_h264_mps_state, ff_h264_lps_state, 64);

  3054.         

  3055.     for(i=-256; i<256; i++){

  3056. START_TIMER

  3057.         put_symbol(&s.c, s.header_state, i*i*i/3*ABS(i), 1);

  3058. STOP_TIMER("put_symbol")

  3059.     }

  3060.     ff_rac_terminate(&s.c);

  3061. 

  3062.     memset(s.header_state, 0, sizeof(s.header_state));

  3063.     ff_init_range_decoder(&s.c, buffer[0], 256*256);

  3064.     ff_init_cabac_states(&s.c, ff_h264_lps_range, ff_h264_mps_state, ff_h264_lps_state, 64);

  3065.     

  3066.     for(i=-256; i<256; i++){

  3067.         int j;

  3068. START_TIMER

  3069.         j= get_symbol(&s.c, s.header_state, 1);

  3070. STOP_TIMER("get_symbol")

  3071.         if(j!=i*i*i/3*ABS(i)) printf("fsck: %d != %d\n", i, j);

  3072.     }

  3073. {

  3074. int level, orientation, x, y;

  3075. int64_t errors[8][4];

  3076. int64_t g=0;

  3077. 

  3078.     memset(errors, 0, sizeof(errors));

  3079.     s.spatial_decomposition_count=3;

  3080.     s.spatial_decomposition_type=0;

  3081.     for(level=0; level<s.spatial_decomposition_count; level++){

  3082.         for(orientation=level ? 1 : 0; orientation<4; orientation++){

  3083.             int w= width  >> (s.spatial_decomposition_count-level);

  3084.             int h= height >> (s.spatial_decomposition_count-level);

  3085.             int stride= width  << (s.spatial_decomposition_count-level);

  3086.             DWTELEM *buf= buffer[0];

  3087.             int64_t error=0;

  3088. 

  3089.             if(orientation&1) buf+=w;

  3090.             if(orientation>1) buf+=stride>>1;

  3091.             

  3092.             memset(buffer[0], 0, sizeof(int)*width*height);

  3093.             buf[w/2 + h/2*stride]= 256*256;

  3094.             ff_spatial_idwt(buffer[0], width, height, width, s->spatial_decomposition_type, s->spatial_decomposition_count);

  3095.             for(y=0; y<height; y++){

  3096.                 for(x=0; x<width; x++){

  3097.                     int64_t d= buffer[0][x + y*width];

  3098.                     error += d*d;

  3099.                     if(ABS(width/2-x)<9 && ABS(height/2-y)<9 && level==2) printf("%8lld ", d);

  3100.                 }

  3101.                 if(ABS(height/2-y)<9 && level==2) printf("\n");

  3102.             }

  3103.             error= (int)(sqrt(error)+0.5);

  3104.             errors[level][orientation]= error;

  3105.             if(g) g=ff_gcd(g, error);

  3106.             else g= error;

  3107.         }

  3108.     }

  3109.     printf("static int const visual_weight[][4]={\n");

  3110.     for(level=0; level<s.spatial_decomposition_count; level++){

  3111.         printf("  {");

  3112.         for(orientation=0; orientation<4; orientation++){

  3113.             printf("%8lld,", errors[level][orientation]/g);

  3114.         }

  3115.         printf("},\n");

  3116.     }

  3117.     printf("};\n");

  3118.     {

  3119.             int level=2;

  3120.             int orientation=3;

  3121.             int w= width  >> (s.spatial_decomposition_count-level);

  3122.             int h= height >> (s.spatial_decomposition_count-level);

  3123.             int stride= width  << (s.spatial_decomposition_count-level);

  3124.             DWTELEM *buf= buffer[0];

  3125.             int64_t error=0;

  3126. 

  3127.             buf+=w;

  3128.             buf+=stride>>1;

  3129.             

  3130.             memset(buffer[0], 0, sizeof(int)*width*height);

  3131. #if 1

  3132.             for(y=0; y<height; y++){

  3133.                 for(x=0; x<width; x++){

  3134.                     int tab[4]={0,2,3,1};

  3135.                     buffer[0][x+width*y]= 256*256*tab[(x&1) + 2*(y&1)];

  3136.                 }

  3137.             }

  3138.             ff_spatial_dwt(buffer[0], width, height, width, s->spatial_decomposition_type, s->spatial_decomposition_count);

  3139. #else

  3140.             for(y=0; y<h; y++){

  3141.                 for(x=0; x<w; x++){

  3142.                     buf[x + y*stride  ]=169;

  3143.                     buf[x + y*stride-w]=64;

  3144.                 }

  3145.             }

  3146.             ff_spatial_idwt(buffer[0], width, height, width, s->spatial_decomposition_type, s->spatial_decomposition_count);

  3147. #endif

  3148.             for(y=0; y<height; y++){

  3149.                 for(x=0; x<width; x++){

  3150.                     int64_t d= buffer[0][x + y*width];

  3151.                     error += d*d;

  3152.                     if(ABS(width/2-x)<9 && ABS(height/2-y)<9) printf("%8lld ", d);

  3153.                 }

  3154.                 if(ABS(height/2-y)<9) printf("\n");

  3155.             }

  3156.     }

  3157. 

  3158. }

  3159.     return 0;

  3160. }

  3161. #endif