src/subtitles/SeparableFilter.h

   1 /*
   2     Copyright 2007  Niels Martin Hansen
   3
   4     This program is free software; you can redistribute it and/or modify
   5     it under the terms of the GNU General Public License as published by
   6     the Free Software Foundation; either version 2 of the License, or
   7     (at your option) any later version.
   8
   9     This program 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
  12     GNU General Public License for more details.
  13
  14     You should have received a copy of the GNU General Public License along
  15     with this program; if not, write to the Free Software Foundation, Inc.,
  16     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  17
  18     Contact:
  19     E-mail: <jiifurusu@gmail.com>
  20     IRC: jfs in #aegisub on irc.rizon.net
  21
  22  */
  23
  24 #pragma once
  25
  26 #ifdef _OPENMP
  27 #include <omp.h>
  28 #endif
  29 #include <math.h>
  30
  31
  32 // Filter an image in horizontal direction with a one-dimensional filter
  33 // PixelWidth is the distance in bytes between pixels
  34 template<ptrdiff_t PixelDist>
  35 void SeparableFilterX(unsigned char *src, unsigned char *dst, int width, int height, ptrdiff_t stride, int *kernel, int kernel_size, int divisor)
  36 {
  37 #pragma omp parallel for
  38         for (int y = 0; y < height; y++) {
  39                 unsigned char *in = src + y*stride;
  40                 unsigned char *out = dst + y*stride;
  41                 for (int x = 0; x < width; x++) {
  42                         int accum = 0;
  43                         for (int k = 0; k < kernel_size; k++) {
  44                                 int xofs = k - kernel_size/2;
  45                                 if (x+xofs < 0) xofs += width;
  46                                 if (x+xofs >= width) xofs -= width;
  47                                 accum += (int)(in[xofs*PixelDist] * kernel[k]);
  48                         }
  49                         accum /= divisor;
  50                         if (accum > 255) accum = 255;
  51                         if (accum < 0) accum = 0;
  52                         *out = (unsigned char)accum;
  53                         in+=PixelDist;
  54                         out+=PixelDist;
  55                 }
  56         }
  57 }
  58
  59
  60 // Filter an image in vertical direction with a one-dimensional filter
  61 // This one templated with PixelWidth since the channel interlacing is horizontal only,
  62 // filtering once vertically will automatically catch all channels.
  63 // (Width must be multiplied by pixel width for that to happen though.)
  64 template<ptrdiff_t PixelDist>
  65 void SeparableFilterY(unsigned char *src, unsigned char *dst, int width, int height, ptrdiff_t stride, int *kernel, int kernel_size, int divisor)
  66 {
  67 #pragma omp parallel for
  68         width *= PixelDist;
  69         for (int  x = 0; x < width; x+=PixelDist) {
  70                 unsigned char *in = src + x;
  71                 unsigned char *out = dst + x;
  72                 for (int y = 0; y < height; y++) {
  73                         int accum = 0;
  74                         for (int k = 0; k < kernel_size; k++) {
  75                                 int yofs = k - kernel_size/2;
  76                                 if (y+yofs < 0) yofs += height;
  77                                 if (y+yofs >= height) yofs -= height;
  78                                 accum += (int)(in[yofs*stride] * kernel[k]);
  79                         }
  80                         accum /= divisor;
  81                         if (accum > 255) accum = 255;
  82                         if (accum < 0) accum = 0;
  83                         *out = (unsigned char)accum;
  84                         in += stride;
  85                         out += stride;
  86                 }
  87         }
  88 }
  89
  90
  91 static inline double NormalDist(double sigma, double x)
  92 {
  93         if (sigma <= 0 && x == 0) return 1;
  94         else if (sigma <= 0) return 0;
  95         else return exp(-(x*x)/(2*sigma*sigma)) / (sigma * sqrt(2*3.1415926535));
  96 }
  97
  98
  99 struct GaussianKernel {
 100         int *kernel;
 101         int width;
 102         int divisor;
 103         inline GaussianKernel(double sigma)
 104         {
 105                 width = (int)(sigma*3 + 0.5) | 1; // binary-or with 1 to make sure the number is odd
 106                 if (width < 3) width = 3;
 107                 kernel = new int[width];
 108                 kernel[width/2] = (int)(NormalDist(sigma, 0) * 255);
 109                 divisor = kernel[width/2];
 110                 for (int x = width/2-1; x >= 0; x--) {
 111                         int val = (int)(NormalDist(sigma, width/2-x) * 255 + 0.5);
 112                         divisor += val*2;
 113                         kernel[x] = val;
 114                         kernel[width - x - 1] = val;
 115                 }
 116         }
 117         inline ~GaussianKernel()
 118         {
 119                 delete[] kernel;
 120         }
 121 };