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.

convolution.c 7.8 KiB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283
  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. * convolution.c: generic convolution functions
  16. */
  17. #include "config.h"
  18. #include <stdlib.h>
  19. #include <stdio.h>
  20. #include <string.h>
  21. #include <math.h>
  22. #include "pipi.h"
  23. #include "pipi_internals.h"
  24. #if !defined TEMPLATE_FILE /* This file uses the template system */
  25. #define TEMPLATE_FLAGS SET_FLAG_GRAY | SET_FLAG_WRAP
  26. #define TEMPLATE_FILE "filter/convolution.c"
  27. #include "pipi_template.h"
  28. pipi_image_t *pipi_convolution(pipi_image_t *src, int m, int n, double mat[])
  29. {
  30. pipi_image_t *ret;
  31. double tmp;
  32. double *hvec, *vvec;
  33. int i, j, besti = -1, bestj = -1;
  34. /* Find the cell with the largest value */
  35. tmp = 0.0;
  36. for(i = 0; i < m * n; i++)
  37. if(mat[i] * mat[i] > tmp)
  38. {
  39. tmp = mat[i] * mat[i];
  40. besti = i % m;
  41. bestj = i / m;
  42. }
  43. /* If the kernel is empty, return an empty picture */
  44. if(tmp == 0.0)
  45. return pipi_new(src->w, src->h);
  46. /* Check whether the matrix rank is 1 */
  47. for(j = 0; j < n; j++)
  48. {
  49. if(j == bestj)
  50. continue;
  51. for(i = 0; i < m; i++)
  52. {
  53. double p, q;
  54. if(i == besti)
  55. continue;
  56. p = mat[j * m + i] * mat[bestj * m + besti];
  57. q = mat[bestj * m + i] * mat[j * m + besti];
  58. if(fabs(p - q) > 0.0001 * 0.0001)
  59. {
  60. if(src->last_modified == PIPI_PIXELS_Y_F)
  61. {
  62. if(src->wrap)
  63. return conv_gray_wrap(src, m, n, mat);
  64. return conv_gray(src, m, n, mat);
  65. }
  66. else
  67. {
  68. if(src->wrap)
  69. return conv_wrap(src, m, n, mat);
  70. return conv(src, m, n, mat);
  71. }
  72. }
  73. }
  74. }
  75. /* Matrix rank is 1! Separate the filter */
  76. hvec = malloc(m * sizeof(double));
  77. vvec = malloc(n * sizeof(double));
  78. tmp = sqrt(fabs(mat[bestj * m + besti]));
  79. for(i = 0; i < m; i++)
  80. hvec[i] = mat[bestj * m + i] / tmp;
  81. for(j = 0; j < n; j++)
  82. vvec[j] = mat[j * m + besti] / tmp;
  83. if(src->last_modified == PIPI_PIXELS_Y_F)
  84. ret = src->wrap ? sepconv_gray_wrap(src, m, hvec, n, vvec)
  85. : sepconv_gray(src, m, hvec, n, vvec);
  86. else
  87. ret = src->wrap ? sepconv_wrap(src, m, hvec, n, vvec)
  88. : sepconv(src, m, hvec, n, vvec);
  89. free(hvec);
  90. free(vvec);
  91. return ret;
  92. }
  93. #else /* XXX: the following functions use the template system */
  94. static pipi_image_t *T(conv)(pipi_image_t *src, int m, int n, double mat[])
  95. {
  96. pipi_image_t *dst;
  97. pipi_pixels_t *srcp, *dstp;
  98. float *srcdata, *dstdata;
  99. int x, y, i, j, w, h;
  100. w = src->w;
  101. h = src->h;
  102. srcp = FLAG_GRAY ? pipi_getpixels(src, PIPI_PIXELS_Y_F)
  103. : pipi_getpixels(src, PIPI_PIXELS_RGBA_F);
  104. srcdata = (float *)srcp->pixels;
  105. dst = pipi_new(w, h);
  106. dstp = FLAG_GRAY ? pipi_getpixels(dst, PIPI_PIXELS_Y_F)
  107. : pipi_getpixels(dst, PIPI_PIXELS_RGBA_F);
  108. dstdata = (float *)dstp->pixels;
  109. for(y = 0; y < h; y++)
  110. {
  111. for(x = 0; x < w; x++)
  112. {
  113. double R = 0., G = 0., B = 0.;
  114. double Y = 0.;
  115. int x2, y2, off = 4 * (y * w + x);
  116. for(j = 0; j < n; j++)
  117. {
  118. y2 = y + j - n / 2;
  119. if(y2 < 0) y2 = FLAG_WRAP ? h - 1 - ((-y2 - 1) % h) : 0;
  120. else if(y2 >= h) y2 = FLAG_WRAP ? y2 % h : h - 1;
  121. for(i = 0; i < m; i++)
  122. {
  123. double f = mat[j * m + i];
  124. x2 = x + i - m / 2;
  125. if(x2 < 0) x2 = FLAG_WRAP ? w - 1 - ((-x2 - 1) % w) : 0;
  126. else if(x2 >= w) x2 = FLAG_WRAP ? x2 % w : w - 1;
  127. if(FLAG_GRAY)
  128. Y += f * srcdata[y2 * w + x2];
  129. else
  130. {
  131. R += f * srcdata[(y2 * w + x2) * 4];
  132. G += f * srcdata[(y2 * w + x2) * 4 + 1];
  133. B += f * srcdata[(y2 * w + x2) * 4 + 2];
  134. }
  135. }
  136. }
  137. if(FLAG_GRAY)
  138. dstdata[y * w + x] = Y < 0.0 ? 0.0 : Y > 1.0 ? 1.0 : Y;
  139. else
  140. {
  141. dstdata[off] = R < 0.0 ? 0.0 : R > 1.0 ? 1.0 : R;
  142. dstdata[off + 1] = G < 0.0 ? 0.0 : G > 1.0 ? 1.0 : G;
  143. dstdata[off + 2] = B < 0.0 ? 0.0 : B > 1.0 ? 1.0 : B;
  144. }
  145. }
  146. }
  147. return dst;
  148. }
  149. static pipi_image_t *T(sepconv)(pipi_image_t *src,
  150. int m, double hvec[], int n, double vvec[])
  151. {
  152. pipi_image_t *dst;
  153. pipi_pixels_t *srcp, *dstp;
  154. float *srcdata, *dstdata;
  155. double *buffer;
  156. int x, y, i, j, w, h;
  157. w = src->w;
  158. h = src->h;
  159. srcp = FLAG_GRAY ? pipi_getpixels(src, PIPI_PIXELS_Y_F)
  160. : pipi_getpixels(src, PIPI_PIXELS_RGBA_F);
  161. srcdata = (float *)srcp->pixels;
  162. dst = pipi_new(w, h);
  163. dstp = FLAG_GRAY ? pipi_getpixels(dst, PIPI_PIXELS_Y_F)
  164. : pipi_getpixels(dst, PIPI_PIXELS_RGBA_F);
  165. dstdata = (float *)dstp->pixels;
  166. buffer = malloc(w * h * (FLAG_GRAY ? 1 : 4) * sizeof(double));
  167. for(y = 0; y < h; y++)
  168. {
  169. for(x = 0; x < w; x++)
  170. {
  171. double R = 0., G = 0., B = 0.;
  172. double Y = 0.;
  173. int x2, off = 4 * (y * w + x);
  174. for(i = 0; i < m; i++)
  175. {
  176. double f = hvec[i];
  177. x2 = x + i - m / 2;
  178. if(x2 < 0) x2 = FLAG_WRAP ? w - 1 - ((-x2 - 1) % w) : 0;
  179. else if(x2 >= w) x2 = FLAG_WRAP ? x2 % w : w - 1;
  180. if(FLAG_GRAY)
  181. Y += f * srcdata[y * w + x2];
  182. else
  183. {
  184. R += f * srcdata[(y * w + x2) * 4];
  185. G += f * srcdata[(y * w + x2) * 4 + 1];
  186. B += f * srcdata[(y * w + x2) * 4 + 2];
  187. }
  188. }
  189. if(FLAG_GRAY)
  190. buffer[y * w + x] = Y;
  191. else
  192. {
  193. buffer[off] = R;
  194. buffer[off + 1] = G;
  195. buffer[off + 2] = B;
  196. }
  197. }
  198. }
  199. for(y = 0; y < h; y++)
  200. {
  201. for(x = 0; x < w; x++)
  202. {
  203. double R = 0., G = 0., B = 0.;
  204. double Y = 0.;
  205. int y2, off = 4 * (y * w + x);
  206. for(j = 0; j < n; j++)
  207. {
  208. double f = vvec[j];
  209. y2 = y + j - n / 2;
  210. if(y2 < 0) y2 = FLAG_WRAP ? h - 1 - ((-y2 - 1) % h) : 0;
  211. else if(y2 >= h) y2 = FLAG_WRAP ? y2 % h : h - 1;
  212. if(FLAG_GRAY)
  213. Y += f * buffer[y2 * w + x];
  214. else
  215. {
  216. R += f * buffer[(y2 * w + x) * 4];
  217. G += f * buffer[(y2 * w + x) * 4 + 1];
  218. B += f * buffer[(y2 * w + x) * 4 + 2];
  219. }
  220. }
  221. if(FLAG_GRAY)
  222. dstdata[y * w + x] = Y < 0.0 ? 0.0 : Y > 1.0 ? 1.0 : Y;
  223. else
  224. {
  225. dstdata[off] = R < 0.0 ? 0.0 : R > 1.0 ? 1.0 : R;
  226. dstdata[off + 1] = G < 0.0 ? 0.0 : G > 1.0 ? 1.0 : G;
  227. dstdata[off + 2] = B < 0.0 ? 0.0 : B > 1.0 ? 1.0 : B;
  228. }
  229. }
  230. }
  231. free(buffer);
  232. return dst;
  233. }
  234. #endif