dbs.c 7.0 KiB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228
  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 <string.h>
  19. #include <stdlib.h>
  20. #include <math.h>
  21. #include "pipi.h"
  22. #include "pipi_internals.h"
  23. #define CELL 16
  24. #define N 7
  25. #define NN ((N * 2 + 1))
  26. /* FIXME: though the algorithm is supposed to stop, we do not have a real,
  27. * guaranteed stop condition here. */
  28. pipi_image_t *pipi_dither_dbs(pipi_image_t *img)
  29. {
  30. double kernel[NN * NN];
  31. double t = 0.;
  32. pipi_image_t *src, *dst, *tmp1, *tmp2;
  33. pipi_pixels_t *dstp, *tmp1p, *tmp2p;
  34. int *changelist;
  35. float *dstdata, *tmp1data, *tmp2data;
  36. int i, j, x, y, w, h, cw, ch;
  37. /* Build our human visual system kernel. */
  38. for(j = 0; j < NN; j++)
  39. for(i = 0; i < NN; i++)
  40. {
  41. double a = (double)(i - N);
  42. double b = (double)(j - N);
  43. kernel[j * NN + i] =
  44. exp(-(a * a + b * b) / (2. * 1.6 * 1.6))
  45. + exp(-(a * a + b * b) / (2. * 0.6 * 0.6));
  46. t += kernel[j * NN + i];
  47. }
  48. for(j = 0; j < NN; j++)
  49. for(i = 0; i < NN; i++)
  50. kernel[j * NN + i] /= t;
  51. w = img->w;
  52. h = img->h;
  53. cw = (w + CELL - 1) / CELL;
  54. ch = (h + CELL - 1) / CELL;
  55. changelist = malloc(cw * ch * sizeof(int));
  56. memset(changelist, 0, cw * ch * sizeof(int));
  57. src = pipi_copy(img);
  58. pipi_getpixels(src, PIPI_PIXELS_Y_F);
  59. tmp1 = pipi_convolution(src, NN, NN, kernel);
  60. tmp1p = pipi_getpixels(tmp1, PIPI_PIXELS_Y_F);
  61. tmp1data = (float *)tmp1p->pixels;
  62. dst = pipi_dither_random(src);
  63. dstp = pipi_getpixels(dst, PIPI_PIXELS_Y_F);
  64. dstdata = (float *)dstp->pixels;
  65. pipi_free(src);
  66. tmp2 = pipi_convolution(dst, NN, NN, kernel);
  67. tmp2p = pipi_getpixels(tmp2, PIPI_PIXELS_Y_F);
  68. tmp2data = (float *)tmp2p->pixels;
  69. for(;;)
  70. {
  71. int cx, cy, n;
  72. int allchanges = 0;
  73. for(cy = 0; cy < ch; cy++)
  74. for(cx = 0; cx < cw; cx++)
  75. {
  76. int changes = 0;
  77. if(changelist[cy * cw + cx] >= 2)
  78. continue;
  79. for(y = cy * CELL; y < (cy + 1) * CELL; y++)
  80. for(x = cx * CELL; x < (cx + 1) * CELL; x++)
  81. {
  82. double d, d2, e, best = 0.;
  83. int opx = -1, opy = -1;
  84. d = dstdata[y * w + x];
  85. /* Compute the effect of a toggle */
  86. e = 0.;
  87. for(j = -N; j < N + 1; j++)
  88. {
  89. if(y + j < 0 || y + j >= h)
  90. continue;
  91. for(i = -N; i < N + 1; i++)
  92. {
  93. double m, p, q1, q2;
  94. if(x + i < 0 || x + i >= w)
  95. continue;
  96. m = kernel[(j + N) * NN + i + N];
  97. p = tmp1data[(y + j) * w + x + i];
  98. q1 = tmp2data[(y + j) * w + x + i];
  99. q2 = q1 - m * d + m * (1. - d);
  100. e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
  101. }
  102. }
  103. if(e > best)
  104. {
  105. best = e;
  106. opx = opy = 0;
  107. }
  108. /* Compute the effect of swaps */
  109. for(n = 0; n < 8; n++)
  110. {
  111. static int const step[] =
  112. { 0, 1, 0, -1, -1, 0, 1, 0, -1, -1, -1, 1, 1, -1, 1, 1 };
  113. int idx = step[n * 2], idy = step[n * 2 + 1];
  114. if(y + idy < 0 || y + idy >= h
  115. || x + idx < 0 || x + idx >= w)
  116. continue;
  117. d2 = dstdata[(y + idy) * w + x + idx];
  118. if(d2 == d)
  119. continue;
  120. e = 0.;
  121. for(j = -N; j < N + 1; j++)
  122. {
  123. if(y + j < 0 || y + j >= h)
  124. continue;
  125. if(j - idy + N < 0 || j - idy + N >= NN)
  126. continue;
  127. for(i = -N; i < N + 1; i++)
  128. {
  129. double ma, mb, p, q1, q2;
  130. if(x + i < 0 || x + i >= w)
  131. continue;
  132. if(i - idx + N < 0 || i - idx + N >= NN)
  133. continue;
  134. ma = kernel[(j + N) * NN + i + N];
  135. mb = kernel[(j - idy + N) * NN + i - idx + N];
  136. p = tmp1data[(y + j) * w + x + i];
  137. q1 = tmp2data[(y + j) * w + x + i];
  138. q2 = q1 - ma * d + ma * d2 - mb * d2 + mb * d;
  139. e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
  140. }
  141. }
  142. if(e > best)
  143. {
  144. best = e;
  145. opx = idx;
  146. opy = idy;
  147. }
  148. }
  149. /* Apply the change if interesting */
  150. if(best <= 0.)
  151. continue;
  152. if(opx || opy)
  153. {
  154. d2 = dstdata[(y + opy) * w + x + opx];
  155. dstdata[(y + opy) * w + x + opx] = d;
  156. }
  157. else
  158. d2 = 1. - d;
  159. dstdata[y * w + x] = d2;
  160. for(j = -N; j < N + 1; j++)
  161. for(i = -N; i < N + 1; i++)
  162. {
  163. double m = kernel[(j + N) * NN + i + N];
  164. if(y + j >= 0 && y + j < h
  165. && x + i >= 0 && x + i < w)
  166. {
  167. t = tmp2data[(y + j) * w + x + i];
  168. tmp2data[(y + j) * w + x + i] = t + m * (d2 - d);
  169. }
  170. if((opx || opy) && y + opy + j >= 0 && y + opy + j < h
  171. && x + opx + i >= 0 && x + opx + i < w)
  172. {
  173. t = tmp2data[(y + opy + j) * w + x + opx + i];
  174. tmp2data[(y + opy + j) * w + x + opx + i]
  175. = t + m * (d - d2);
  176. }
  177. }
  178. changes++;
  179. }
  180. if(changes == 0)
  181. changelist[cy * cw + cx]++;
  182. allchanges += changes;
  183. }
  184. if(allchanges == 0)
  185. break;
  186. }
  187. free(changelist);
  188. pipi_free(tmp1);
  189. pipi_free(tmp2);
  190. return dst;
  191. }