/* * libpipi Pathetic image processing interface library * Copyright (c) 2004-2008 Sam Hocevar * All Rights Reserved * * $Id$ * * This library is free software. It comes without any warranty, to * the extent permitted by applicable law. You can redistribute it * and/or modify it under the terms of the Do What The Fuck You Want * To Public License, Version 2, as published by Sam Hocevar. See * http://sam.zoy.org/wtfpl/COPYING for more details. */ /* * convolution.c: generic convolution functions */ #include "config.h" #include #include #include #include #include "pipi.h" #include "pipi_internals.h" #if !defined TEMPLATE_FILE /* This file uses the template system */ #define TEMPLATE_FLAGS SET_FLAG_GRAY | SET_FLAG_WRAP #define TEMPLATE_FILE "filter/convolution.c" #include "pipi_template.h" pipi_image_t *pipi_convolution(pipi_image_t *src, int m, int n, double mat[]) { pipi_image_t *ret; double tmp; double *hvec, *vvec; int i, j, besti = -1, bestj = -1; /* Find the cell with the largest value */ tmp = 0.0; for(i = 0; i < m * n; i++) if(mat[i] * mat[i] > tmp) { tmp = mat[i] * mat[i]; besti = i % m; bestj = i / m; } /* If the kernel is empty, return an empty picture */ if(tmp == 0.0) return pipi_new(src->w, src->h); /* Check whether the matrix rank is 1 */ for(j = 0; j < n; j++) { if(j == bestj) continue; for(i = 0; i < m; i++) { double p, q; if(i == besti) continue; p = mat[j * m + i] * mat[bestj * m + besti]; q = mat[bestj * m + i] * mat[j * m + besti]; if(fabs(p - q) > 0.0001 * 0.0001) { if(src->last_modified == PIPI_PIXELS_Y_F32) { if(src->wrap) return conv_gray_wrap(src, m, n, mat); return conv_gray(src, m, n, mat); } else { if(src->wrap) return conv_wrap(src, m, n, mat); return conv(src, m, n, mat); } } } } /* Matrix rank is 1! Separate the filter */ hvec = malloc(m * sizeof(double)); vvec = malloc(n * sizeof(double)); tmp = sqrt(fabs(mat[bestj * m + besti])); for(i = 0; i < m; i++) hvec[i] = mat[bestj * m + i] / tmp; for(j = 0; j < n; j++) vvec[j] = mat[j * m + besti] / tmp; if(src->last_modified == PIPI_PIXELS_Y_F32) ret = src->wrap ? sepconv_gray_wrap(src, m, hvec, n, vvec) : sepconv_gray(src, m, hvec, n, vvec); else ret = src->wrap ? sepconv_wrap(src, m, hvec, n, vvec) : sepconv(src, m, hvec, n, vvec); free(hvec); free(vvec); return ret; } #else /* XXX: the following functions use the template system */ static pipi_image_t *T(conv)(pipi_image_t *src, int m, int n, double mat[]) { pipi_image_t *dst; pipi_pixels_t *srcp, *dstp; float *srcdata, *dstdata; int x, y, i, j, w, h; w = src->w; h = src->h; srcp = FLAG_GRAY ? pipi_get_pixels(src, PIPI_PIXELS_Y_F32) : pipi_get_pixels(src, PIPI_PIXELS_RGBA_F32); srcdata = (float *)srcp->pixels; dst = pipi_new(w, h); dstp = FLAG_GRAY ? pipi_get_pixels(dst, PIPI_PIXELS_Y_F32) : pipi_get_pixels(dst, PIPI_PIXELS_RGBA_F32); dstdata = (float *)dstp->pixels; for(y = 0; y < h; y++) { for(x = 0; x < w; x++) { double R = 0., G = 0., B = 0.; double Y = 0.; int x2, y2, off = 4 * (y * w + x); for(j = 0; j < n; j++) { y2 = y + j - n / 2; if(y2 < 0) y2 = FLAG_WRAP ? h - 1 - ((-y2 - 1) % h) : 0; else if(y2 >= h) y2 = FLAG_WRAP ? y2 % h : h - 1; for(i = 0; i < m; i++) { double f = mat[j * m + i]; x2 = x + i - m / 2; if(x2 < 0) x2 = FLAG_WRAP ? w - 1 - ((-x2 - 1) % w) : 0; else if(x2 >= w) x2 = FLAG_WRAP ? x2 % w : w - 1; if(FLAG_GRAY) Y += f * srcdata[y2 * w + x2]; else { R += f * srcdata[(y2 * w + x2) * 4]; G += f * srcdata[(y2 * w + x2) * 4 + 1]; B += f * srcdata[(y2 * w + x2) * 4 + 2]; } } } if(FLAG_GRAY) dstdata[y * w + x] = Y < 0.0 ? 0.0 : Y > 1.0 ? 1.0 : Y; else { dstdata[off] = R < 0.0 ? 0.0 : R > 1.0 ? 1.0 : R; dstdata[off + 1] = G < 0.0 ? 0.0 : G > 1.0 ? 1.0 : G; dstdata[off + 2] = B < 0.0 ? 0.0 : B > 1.0 ? 1.0 : B; } } } return dst; } static pipi_image_t *T(sepconv)(pipi_image_t *src, int m, double hvec[], int n, double vvec[]) { pipi_image_t *dst; pipi_pixels_t *srcp, *dstp; float *srcdata, *dstdata; double *buffer; int x, y, i, j, w, h; w = src->w; h = src->h; srcp = FLAG_GRAY ? pipi_get_pixels(src, PIPI_PIXELS_Y_F32) : pipi_get_pixels(src, PIPI_PIXELS_RGBA_F32); srcdata = (float *)srcp->pixels; dst = pipi_new(w, h); dstp = FLAG_GRAY ? pipi_get_pixels(dst, PIPI_PIXELS_Y_F32) : pipi_get_pixels(dst, PIPI_PIXELS_RGBA_F32); dstdata = (float *)dstp->pixels; buffer = malloc(w * h * (FLAG_GRAY ? 1 : 4) * sizeof(double)); for(y = 0; y < h; y++) { for(x = 0; x < w; x++) { double R = 0., G = 0., B = 0.; double Y = 0.; int x2, off = 4 * (y * w + x); for(i = 0; i < m; i++) { double f = hvec[i]; x2 = x + i - m / 2; if(x2 < 0) x2 = FLAG_WRAP ? w - 1 - ((-x2 - 1) % w) : 0; else if(x2 >= w) x2 = FLAG_WRAP ? x2 % w : w - 1; if(FLAG_GRAY) Y += f * srcdata[y * w + x2]; else { R += f * srcdata[(y * w + x2) * 4]; G += f * srcdata[(y * w + x2) * 4 + 1]; B += f * srcdata[(y * w + x2) * 4 + 2]; } } if(FLAG_GRAY) buffer[y * w + x] = Y; else { buffer[off] = R; buffer[off + 1] = G; buffer[off + 2] = B; } } } for(y = 0; y < h; y++) { for(x = 0; x < w; x++) { double R = 0., G = 0., B = 0.; double Y = 0.; int y2, off = 4 * (y * w + x); for(j = 0; j < n; j++) { double f = vvec[j]; y2 = y + j - n / 2; if(y2 < 0) y2 = FLAG_WRAP ? h - 1 - ((-y2 - 1) % h) : 0; else if(y2 >= h) y2 = FLAG_WRAP ? y2 % h : h - 1; if(FLAG_GRAY) Y += f * buffer[y2 * w + x]; else { R += f * buffer[(y2 * w + x) * 4]; G += f * buffer[(y2 * w + x) * 4 + 1]; B += f * buffer[(y2 * w + x) * 4 + 2]; } } if(FLAG_GRAY) dstdata[y * w + x] = Y < 0.0 ? 0.0 : Y > 1.0 ? 1.0 : Y; else { dstdata[off] = R < 0.0 ? 0.0 : R > 1.0 ? 1.0 : R; dstdata[off + 1] = G < 0.0 ? 0.0 : G > 1.0 ? 1.0 : G; dstdata[off + 2] = B < 0.0 ? 0.0 : B > 1.0 ? 1.0 : B; } } } free(buffer); return dst; } #endif