Changeset 2660


Ignore:
Timestamp:
Aug 3, 2008, 5:54:48 PM (14 years ago)
Author:
Sam Hocevar
Message:
  • convolution.c: automatically detect when a convolution filter is separable and switch algorithms to speed up processing.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • libpipi/trunk/pipi/filter/convolution.c

    r2658 r2660  
    2828#include "pipi_internals.h"
    2929
     30static pipi_image_t *pipi_convolution_standard(pipi_image_t *src,
     31                                               int m, int n, double mat[]);
     32static pipi_image_t *pipi_convolution_separable(pipi_image_t *src,
     33                                                int m, double hvec[],
     34                                                int n, double vvec[]);
     35
    3036pipi_image_t *pipi_convolution(pipi_image_t *src, int m, int n, double mat[])
     37{
     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                return pipi_convolution_standard(src, m, n, mat);
     74        }
     75    }
     76
     77printf("separabble!\n");
     78
     79    /* Matrix rank is 1! Separate the filter */
     80    /* FIXME: memleak */
     81    hvec = malloc(m * sizeof(double));
     82    vvec = malloc(n * sizeof(double));
     83
     84    tmp = sqrt(fabs(mat[bestj * m + besti]));
     85    for(i = 0; i < m; i++)
     86        hvec[i] = mat[bestj * m + i] / tmp;
     87    for(j = 0; j < n; j++)
     88        vvec[j] = mat[j * m + besti] / tmp;
     89
     90    return pipi_convolution_separable(src, m, hvec, n, vvec);
     91}
     92
     93static pipi_image_t *pipi_convolution_standard(pipi_image_t *src,
     94                                               int m, int n, double mat[])
    3195{
    3296    pipi_image_t *dst;
     
    66130                    for(i = 0; i < m; i++)
    67131                    {
    68                         x2 = x + i;
     132                        x2 = x + i - m / 2;
    69133                        if(x2 < 0) x2 = 0;
    70134                        else if(x2 >= w) x2 = w - 1;
    71135
    72                         Y += mat[j * m + i] * srcdata[y * w + x2];
     136                        Y += mat[j * m + i] * srcdata[y2 * w + x2];
    73137                    }
    74138                }
     
    91155                        double f = mat[j * m + i];
    92156
    93                         x2 = x + i;
     157                        x2 = x + i - m / 2;
    94158                        if(x2 < 0) x2 = 0;
    95159                        else if(x2 >= w) x2 = w - 1;
    96160
    97                         R += f * srcdata[(y * w + x2) * 4];
    98                         G += f * srcdata[(y * w + x2) * 4 + 1];
    99                         B += f * srcdata[(y * w + x2) * 4 + 2];
     161                        R += f * srcdata[(y2 * w + x2) * 4];
     162                        G += f * srcdata[(y2 * w + x2) * 4 + 1];
     163                        B += f * srcdata[(y2 * w + x2) * 4 + 2];
    100164                    }
    101165                }
     
    111175}
    112176
     177static pipi_image_t *pipi_convolution_separable(pipi_image_t *src,
     178                                                int m, double hvec[],
     179                                                int n, double vvec[])
     180{
     181    pipi_image_t *dst;
     182    pipi_pixels_t *srcp, *dstp;
     183    float *srcdata, *dstdata;
     184    double *buffer;
     185    int x, y, i, j, w, h, gray;
     186
     187    w = src->w;
     188    h = src->h;
     189
     190    gray = (src->last_modified == PIPI_PIXELS_Y_F);
     191
     192    srcp = 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 = 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 * (gray ? 1 : 4) * sizeof(double));
     202
     203    for(y = 0; y < h; y++)
     204    {
     205        for(x = 0; x < w; x++)
     206        {
     207            if(gray)
     208            {
     209                double Y = 0.;
     210                int x2;
     211
     212                for(i = 0; i < m; i++)
     213                {
     214                    x2 = x + i - m / 2;
     215                    if(x2 < 0) x2 = 0;
     216                    else if(x2 >= w) x2 = w - 1;
     217
     218                    Y += hvec[i] * srcdata[y * w + x2];
     219                }
     220
     221                buffer[y * w + x] = Y;
     222            }
     223            else
     224            {
     225                double R = 0., G = 0., B = 0.;
     226                int x2, off = 4 * (y * w + x);
     227
     228                for(i = 0; i < m; i++)
     229                {
     230                    double f = hvec[i];
     231
     232                    x2 = x + i - m / 2;
     233                    if(x2 < 0) x2 = 0;
     234                    else if(x2 >= w) x2 = w - 1;
     235
     236                    R += f * srcdata[(y * w + x2) * 4];
     237                    G += f * srcdata[(y * w + x2) * 4 + 1];
     238                    B += f * srcdata[(y * w + x2) * 4 + 2];
     239                }
     240
     241                buffer[off] = R;
     242                buffer[off + 1] = G;
     243                buffer[off + 2] = B;
     244            }
     245        }
     246    }
     247
     248    for(y = 0; y < h; y++)
     249    {
     250        for(x = 0; x < w; x++)
     251        {
     252            if(gray)
     253            {
     254                double Y = 0.;
     255                int y2;
     256
     257                for(j = 0; j < n; j++)
     258                {
     259                    y2 = y + j - n / 2;
     260                    if(y2 < 0) y2 = 0;
     261                    else if(y2 >= h) y2 = h - 1;
     262
     263                    Y += vvec[j] * buffer[y2 * w + x];
     264                }
     265
     266                dstdata[y * w + x] = Y;
     267            }
     268            else
     269            {
     270                double R = 0., G = 0., B = 0.;
     271                int y2, off = 4 * (y * w + x);
     272
     273                for(j = 0; j < n; j++)
     274                {
     275                    double f = vvec[j];
     276
     277                    y2 = y + j - n / 2;
     278                    if(y2 < 0) y2 = 0;
     279                    else if(y2 >= h) y2 = h - 1;
     280
     281                    R += f * buffer[(y2 * w + x) * 4];
     282                    G += f * buffer[(y2 * w + x) * 4 + 1];
     283                    B += f * buffer[(y2 * w + x) * 4 + 2];
     284                }
     285
     286                dstdata[off] = R;
     287                dstdata[off + 1] = G;
     288                dstdata[off + 2] = B;
     289            }
     290        }
     291    }
     292
     293    free(buffer);
     294
     295    return dst;
     296}
     297
Note: See TracChangeset for help on using the changeset viewer.