Вы не можете выбрать более 25 тем Темы должны начинаться с буквы или цифры, могут содержать дефисы(-) и должны содержать не более 35 символов.
 
 
 
 
 
 

230 строки
7.0 KiB

  1. /*
  2. * libpipi Pathetic image processing interface library
  3. * Copyright (c) 2004-2008 Sam Hocevar <sam@zoy.org>
  4. * All Rights Reserved
  5. *
  6. * $Id$
  7. *
  8. * This library is free software. It comes without any warranty, to
  9. * the extent permitted by applicable law. You can redistribute it
  10. * and/or modify it under the terms of the Do What The Fuck You Want
  11. * To Public License, Version 2, as published by Sam Hocevar. See
  12. * http://sam.zoy.org/wtfpl/COPYING for more details.
  13. */
  14. /*
  15. * dbs.c: Direct Binary Search dithering functions
  16. */
  17. #include "config.h"
  18. #include "common.h"
  19. #include <string.h>
  20. #include <stdlib.h>
  21. #include <math.h>
  22. #include "pipi.h"
  23. #include "pipi_internals.h"
  24. #define CELL 16
  25. #define N 7
  26. #define NN ((N * 2 + 1))
  27. /* FIXME: though the algorithm is supposed to stop, we do not have a real,
  28. * guaranteed stop condition here. */
  29. pipi_image_t *pipi_dither_dbs(pipi_image_t *img)
  30. {
  31. double kernel[NN * NN];
  32. double t = 0.;
  33. pipi_image_t *src, *dst, *tmp1, *tmp2;
  34. pipi_pixels_t *dstp, *tmp1p, *tmp2p;
  35. int *changelist;
  36. float *dstdata, *tmp1data, *tmp2data;
  37. int i, j, x, y, w, h, cw, ch;
  38. /* Build our human visual system kernel. */
  39. for(j = 0; j < NN; j++)
  40. for(i = 0; i < NN; i++)
  41. {
  42. double a = (double)(i - N);
  43. double b = (double)(j - N);
  44. kernel[j * NN + i] =
  45. exp(-(a * a + b * b) / (2. * 1.6 * 1.6))
  46. + exp(-(a * a + b * b) / (2. * 0.6 * 0.6));
  47. t += kernel[j * NN + i];
  48. }
  49. for(j = 0; j < NN; j++)
  50. for(i = 0; i < NN; i++)
  51. kernel[j * NN + i] /= t;
  52. w = img->w;
  53. h = img->h;
  54. cw = (w + CELL - 1) / CELL;
  55. ch = (h + CELL - 1) / CELL;
  56. changelist = malloc(cw * ch * sizeof(int));
  57. memset(changelist, 0, cw * ch * sizeof(int));
  58. src = pipi_copy(img);
  59. pipi_getpixels(src, PIPI_PIXELS_Y_F);
  60. tmp1 = pipi_convolution(src, NN, NN, kernel);
  61. tmp1p = pipi_getpixels(tmp1, PIPI_PIXELS_Y_F);
  62. tmp1data = (float *)tmp1p->pixels;
  63. dst = pipi_dither_random(src);
  64. dstp = pipi_getpixels(dst, PIPI_PIXELS_Y_F);
  65. dstdata = (float *)dstp->pixels;
  66. pipi_free(src);
  67. tmp2 = pipi_convolution(dst, NN, NN, kernel);
  68. tmp2p = pipi_getpixels(tmp2, PIPI_PIXELS_Y_F);
  69. tmp2data = (float *)tmp2p->pixels;
  70. for(;;)
  71. {
  72. int cx, cy, n;
  73. int allchanges = 0;
  74. for(cy = 0; cy < ch; cy++)
  75. for(cx = 0; cx < cw; cx++)
  76. {
  77. int changes = 0;
  78. if(changelist[cy * cw + cx] >= 2)
  79. continue;
  80. for(y = cy * CELL; y < (cy + 1) * CELL; y++)
  81. for(x = cx * CELL; x < (cx + 1) * CELL; x++)
  82. {
  83. double d, d2, e, best = 0.;
  84. int opx = -1, opy = -1;
  85. d = dstdata[y * w + x];
  86. /* Compute the effect of a toggle */
  87. e = 0.;
  88. for(j = -N; j < N + 1; j++)
  89. {
  90. if(y + j < 0 || y + j >= h)
  91. continue;
  92. for(i = -N; i < N + 1; i++)
  93. {
  94. double m, p, q1, q2;
  95. if(x + i < 0 || x + i >= w)
  96. continue;
  97. m = kernel[(j + N) * NN + i + N];
  98. p = tmp1data[(y + j) * w + x + i];
  99. q1 = tmp2data[(y + j) * w + x + i];
  100. q2 = q1 - m * d + m * (1. - d);
  101. e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
  102. }
  103. }
  104. if(e > best)
  105. {
  106. best = e;
  107. opx = opy = 0;
  108. }
  109. /* Compute the effect of swaps */
  110. for(n = 0; n < 8; n++)
  111. {
  112. static int const step[] =
  113. { 0, 1, 0, -1, -1, 0, 1, 0, -1, -1, -1, 1, 1, -1, 1, 1 };
  114. int idx = step[n * 2], idy = step[n * 2 + 1];
  115. if(y + idy < 0 || y + idy >= h
  116. || x + idx < 0 || x + idx >= w)
  117. continue;
  118. d2 = dstdata[(y + idy) * w + x + idx];
  119. if(d2 == d)
  120. continue;
  121. e = 0.;
  122. for(j = -N; j < N + 1; j++)
  123. {
  124. if(y + j < 0 || y + j >= h)
  125. continue;
  126. if(j - idy + N < 0 || j - idy + N >= NN)
  127. continue;
  128. for(i = -N; i < N + 1; i++)
  129. {
  130. double ma, mb, p, q1, q2;
  131. if(x + i < 0 || x + i >= w)
  132. continue;
  133. if(i - idx + N < 0 || i - idx + N >= NN)
  134. continue;
  135. ma = kernel[(j + N) * NN + i + N];
  136. mb = kernel[(j - idy + N) * NN + i - idx + N];
  137. p = tmp1data[(y + j) * w + x + i];
  138. q1 = tmp2data[(y + j) * w + x + i];
  139. q2 = q1 - ma * d + ma * d2 - mb * d2 + mb * d;
  140. e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
  141. }
  142. }
  143. if(e > best)
  144. {
  145. best = e;
  146. opx = idx;
  147. opy = idy;
  148. }
  149. }
  150. /* Apply the change if interesting */
  151. if(best <= 0.)
  152. continue;
  153. if(opx || opy)
  154. {
  155. d2 = dstdata[(y + opy) * w + x + opx];
  156. dstdata[(y + opy) * w + x + opx] = d;
  157. }
  158. else
  159. d2 = 1. - d;
  160. dstdata[y * w + x] = d2;
  161. for(j = -N; j < N + 1; j++)
  162. for(i = -N; i < N + 1; i++)
  163. {
  164. double m = kernel[(j + N) * NN + i + N];
  165. if(y + j >= 0 && y + j < h
  166. && x + i >= 0 && x + i < w)
  167. {
  168. t = tmp2data[(y + j) * w + x + i];
  169. tmp2data[(y + j) * w + x + i] = t + m * (d2 - d);
  170. }
  171. if((opx || opy) && y + opy + j >= 0 && y + opy + j < h
  172. && x + opx + i >= 0 && x + opx + i < w)
  173. {
  174. t = tmp2data[(y + opy + j) * w + x + opx + i];
  175. tmp2data[(y + opy + j) * w + x + opx + i]
  176. = t + m * (d - d2);
  177. }
  178. }
  179. changes++;
  180. }
  181. if(changes == 0)
  182. changelist[cy * cw + cx]++;
  183. allchanges += changes;
  184. }
  185. if(allchanges == 0)
  186. break;
  187. }
  188. free(changelist);
  189. pipi_free(tmp1);
  190. pipi_free(tmp2);
  191. return dst;
  192. }