source: libpipi/trunk/pipi/filter/convolution.c @ 2902

Last change on this file since 2902 was 2902, checked in by Sam Hocevar, 12 years ago

Support C99 types on Win32 through the same hacks as in libcaca.

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