convolve.c 5.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195
  1. /*
  2. * Copyright (c) 2011, Tom Distler (http://tdistler.com)
  3. * All rights reserved.
  4. *
  5. * The BSD License
  6. *
  7. * Redistribution and use in source and binary forms, with or without
  8. * modification, are permitted provided that the following conditions are met:
  9. *
  10. * - Redistributions of source code must retain the above copyright notice,
  11. * this list of conditions and the following disclaimer.
  12. *
  13. * - Redistributions in binary form must reproduce the above copyright notice,
  14. * this list of conditions and the following disclaimer in the documentation
  15. * and/or other materials provided with the distribution.
  16. *
  17. * - Neither the name of the tdistler.com nor the names of its contributors may
  18. * be used to endorse or promote products derived from this software without
  19. * specific prior written permission.
  20. *
  21. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  22. * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23. * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24. * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
  25. * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  26. * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  27. * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  28. * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  29. * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  30. * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  31. * POSSIBILITY OF SUCH DAMAGE.
  32. */
  33. #include "convolve.h"
  34. #include <stdlib.h>
  35. float KBND_SYMMETRIC(const float *img, int w, int h, int x, int y, float bnd_const)
  36. {
  37. (void)bnd_const;
  38. if (x<0) x=-1-x;
  39. else if (x>=w) x=(w-(x-w))-1;
  40. if (y<0) y=-1-y;
  41. else if (y>=h) y=(h-(y-h))-1;
  42. return img[y*w + x];
  43. }
  44. float KBND_REPLICATE(const float *img, int w, int h, int x, int y, float bnd_const)
  45. {
  46. (void)bnd_const;
  47. if (x<0) x=0;
  48. if (x>=w) x=w-1;
  49. if (y<0) y=0;
  50. if (y>=h) y=h-1;
  51. return img[y*w + x];
  52. }
  53. float KBND_CONSTANT(const float *img, int w, int h, int x, int y, float bnd_const)
  54. {
  55. if (x<0) x=0;
  56. if (y<0) y=0;
  57. if (x>=w || y>=h)
  58. return bnd_const;
  59. return img[y*w + x];
  60. }
  61. static float _calc_scale(const struct _kernel *k)
  62. {
  63. int ii,k_len;
  64. double sum=0.0;
  65. if (k->normalized)
  66. return 1.0f;
  67. else {
  68. k_len = k->w * k->h;
  69. for (ii=0; ii<k_len; ++ii)
  70. sum += k->kernel[ii];
  71. if (sum != 0.0)
  72. return (float)(1.0 / sum);
  73. return 1.0f;
  74. }
  75. }
  76. void _iqa_convolve(float *img, int w, int h, const struct _kernel *k, float *result, int *rw, int *rh)
  77. {
  78. int x,y,kx,ky,u,v;
  79. int uc = k->w/2;
  80. int vc = k->h/2;
  81. int kw_even = (k->w&1)?0:1;
  82. int kh_even = (k->h&1)?0:1;
  83. int dst_w = w - k->w + 1;
  84. int dst_h = h - k->h + 1;
  85. int img_offset,k_offset;
  86. double sum;
  87. float scale, *dst=result;
  88. if (!dst)
  89. dst = img; /* Convolve in-place */
  90. /* Kernel is applied to all positions where the kernel is fully contained
  91. * in the image */
  92. scale = _calc_scale(k);
  93. for (y=0; y < dst_h; ++y) {
  94. for (x=0; x < dst_w; ++x) {
  95. sum = 0.0;
  96. k_offset = 0;
  97. ky = y+vc;
  98. kx = x+uc;
  99. for (v=-vc; v <= vc-kh_even; ++v) {
  100. img_offset = (ky+v)*w + kx;
  101. for (u=-uc; u <= uc-kw_even; ++u, ++k_offset) {
  102. sum += img[img_offset+u] * k->kernel[k_offset];
  103. }
  104. }
  105. dst[y*dst_w + x] = (float)(sum * scale);
  106. }
  107. }
  108. if (rw) *rw = dst_w;
  109. if (rh) *rh = dst_h;
  110. }
  111. int _iqa_img_filter(float *img, int w, int h, const struct _kernel *k, float *result)
  112. {
  113. int x,y;
  114. int img_offset;
  115. float scale, *dst=result;
  116. if (!k || !k->bnd_opt)
  117. return 1;
  118. if (!dst) {
  119. dst = (float*)malloc(w*h*sizeof(float));
  120. if (!dst)
  121. return 2;
  122. }
  123. scale = _calc_scale(k);
  124. /* Kernel is applied to all positions where top-left corner is in the image */
  125. for (y=0; y < h; ++y) {
  126. for (x=0; x < w; ++x) {
  127. dst[y*w + x] = _iqa_filter_pixel(img, w, h, x, y, k, scale);
  128. }
  129. }
  130. /* If no result buffer given, copy results to image buffer */
  131. if (!result) {
  132. for (y=0; y<h; ++y) {
  133. img_offset = y*w;
  134. for (x=0; x<w; ++x, ++img_offset) {
  135. img[img_offset] = dst[img_offset];
  136. }
  137. }
  138. free(dst);
  139. }
  140. return 0;
  141. }
  142. float _iqa_filter_pixel(const float *img, int w, int h, int x, int y, const struct _kernel *k, const float kscale)
  143. {
  144. int u,v,uc,vc;
  145. int kw_even,kh_even;
  146. int x_edge_left,x_edge_right,y_edge_top,y_edge_bottom;
  147. int edge,img_offset,k_offset;
  148. double sum;
  149. if (!k)
  150. return img[y*w + x];
  151. uc = k->w/2;
  152. vc = k->h/2;
  153. kw_even = (k->w&1)?0:1;
  154. kh_even = (k->h&1)?0:1;
  155. x_edge_left = uc;
  156. x_edge_right = w-uc;
  157. y_edge_top = vc;
  158. y_edge_bottom = h-vc;
  159. edge = 0;
  160. if (x < x_edge_left || y < y_edge_top || x >= x_edge_right || y >= y_edge_bottom)
  161. edge = 1;
  162. sum = 0.0;
  163. k_offset = 0;
  164. for (v=-vc; v <= vc-kh_even; ++v) {
  165. img_offset = (y+v)*w + x;
  166. for (u=-uc; u <= uc-kw_even; ++u, ++k_offset) {
  167. if (!edge)
  168. sum += img[img_offset+u] * k->kernel[k_offset];
  169. else
  170. sum += k->bnd_opt(img, w, h, x+u, y+v, k->bnd_const) * k->kernel[k_offset];
  171. }
  172. }
  173. return (float)(sum * kscale);
  174. }