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

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

Fix headers.

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