core/DesktopEditor/cximage/libpsd/gaussian_blur.c

/**
 * libpsd - Photoshop file formats (*.psd) decode library
 * Copyright (C) 2004-2007 Graphest Software.
 *
 * libpsd is the legal property of its developers, whose names are too numerous
 * to list here.  Please refer to the COPYRIGHT file distributed with this
 * source distribution.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Library General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 * $Id: gaussian_blur.c, created by Patrick in 2006.06.26, libpsd@graphest.com Exp $
 */

#include <math.h>
#include "libpsd.h"
#include "psd_system.h"
#include "psd_color.h"
#include "psd_bitmap.h"
#include "psd_math.h"


/*
 * The equations: g(r) = exp (- r^2 / (2 * sigma^2))
 *                   r = sqrt (x^2 + y ^2)
 */
psd_static psd_int * psd_make_curve(psd_double sigma, psd_int *length)
{
	psd_int    *curve;
	psd_double  sigma2;
	psd_double  l;
	psd_int     temp;
	psd_int     i, n;

	sigma2 = 2 * sigma * sigma;
	l = sqrt(-sigma2 * log(1.0 / 255.0));

	n = (psd_int)(ceil(l) * 2);
	if ((n % 2) == 0)
		n += 1;

	curve = (psd_int *)psd_malloc(n * 4);

	*length = n / 2;
	curve += *length;
	curve[0] = 255;

	for (i = 1; i <= *length; i++)
	{
		temp = (psd_int)(exp(- (i * i) / sigma2) * 255);
		curve[-i] = temp;
		curve[i] = temp;
	}

	return curve;
}

psd_static void psd_run_length_encode (psd_color_component *src, psd_int *dest, psd_int width)
{
	psd_int start;
	psd_int i;
	psd_int j;
	psd_color_component last;

	last = *src;
	start = 0;

	for (i = 1; i < width; i++)
	{
		if (*src != last)
		{
			for (j = start; j < i; j++)
			{
				*dest++ = (i - j);
				*dest++ = last;
			}
			start = i;
			last = *src;
		}
		src ++;
	}

	for (j = start; j < i; j++)
	{
		*dest++ = (i - j);
		*dest++ = last;
	}
}

void psd_bitmap_gaussian_blur_alpha_channel(psd_bitmap * bitmap, psd_double radius)
{
	psd_int     width, height;
	psd_int    *buf, *bb;
	psd_int     pixels;
	psd_int     total = 1, total2;
	psd_int     i, row, col;
	psd_int     start, end;
	psd_int    *curve;
	psd_int    *sum = NULL;
	psd_int     val;
	psd_int     length;
	psd_int     initial_p, initial_m;
	psd_double  std_dev;
	psd_color_component * src, * sp;
	psd_argb_color * src_data, * dst_data;

	if(radius <= 0.0)
		return;

	width  = bitmap->width;
	height = bitmap->height;

	if(width < 1 || height < 1)
		return;

	buf = (psd_int *)psd_malloc(PSD_MAX(width, height) * 2 * 4);
	src = (psd_color_component *)psd_malloc(PSD_MAX(width, height));

	/*  First the vertical pass  */
	radius = fabs(radius) * 1.4;
	std_dev = sqrt(-(radius * radius) / (2 * log(1.0 / 255.0)));

	curve = psd_make_curve (std_dev, &length);
	sum = (psd_int *)psd_malloc((2 * length + 1) * 4);

	sum[0] = 0;
	for(i = 1; i <= length*2; i++)
		sum[i] = curve[i-length-1] + sum[i-1];
	sum += length;

	total = sum[length] - sum[-length];
	total2 = total / 2;

	/*  First, the vertical pass  */
	for(col = 0; col < width; col++)
	{
		for(row = 0, src_data = bitmap->image_data + width + col; row < height - 1; row ++, src_data += width)
			src[row] = PSD_GET_ALPHA_COMPONENT(*src_data);
		src[height - 1] = 0;

		sp = src;

		initial_p = *sp;
		initial_m = *(sp + height - 1);

		/*  Determine a run-length encoded version of the row  */
		psd_run_length_encode(sp, buf, height);

		for(row = 0, dst_data = bitmap->image_data + col; row < height; row ++, dst_data += width)
		{
			start = (row < length) ? -row : -length;
			end = (height <= (row + length) ? (height - row - 1) : length);

			val = 0;
			i = start;
			bb = buf + (row + i) * 2;

			if(start != -length)
				val += initial_p * (sum[start] - sum[-length]);

			while(i < end)
			{
				pixels = bb[0];
				i += pixels;

				if(i > end)
					i = end;

				val += bb[1] * (sum[i] - sum[start]);
				bb += (pixels * 2);
				start = i;
			}

			if(end != length)
				val += initial_m * (sum[length] - sum[end]);

			*dst_data = (*dst_data & 0x00FFFFFF) | ((val + total2) / total << 24);
		}
	}

	/*  Now the horizontal pass  */
	for(row = 0; row < height; row++)
	{
		for(col = 0, src_data = bitmap->image_data + row * width + 1; col < width - 1; col ++, src_data ++)
			src[col] = PSD_GET_ALPHA_COMPONENT(*src_data);
		src[width - 1] = 0;

		sp = src;

		initial_p = *sp;
		initial_m = *(sp + width - 1);

		/*  Determine a run-length encoded version of the row  */
		psd_run_length_encode(sp, buf, width);

		for(col = 0, dst_data = bitmap->image_data + row * width; col < width; col ++, dst_data ++)
		{
			start = (col < length) ? -col : -length;
			end = (width <= (col + length)) ? (width - col - 1) : length;

			val = 0;
			i = start;
			bb = buf + (col + i) * 2;

			if(start != -length)
				val += initial_p * (sum[start] - sum[-length]);

			while(i < end)
			{
				pixels = bb[0];
				i += pixels;

				if(i > end)
					i = end;

				val += bb[1] * (sum[i] - sum[start]);
				bb += (pixels * 2);
				start = i;
			}

			if(end != length)
				val += initial_m * (sum[length] - sum[end]);

			*dst_data = (*dst_data & 0x00FFFFFF) | ((val + total2) / total << 24);
		}
	}

	/*  free buffers  */
	psd_free(curve - length);
	psd_free(sum - length);
	psd_free(buf);
	psd_free(src);
}