FreeImage Display Demo

/*		Filtered Image Rescaling
 *
 *		by Dale Schumacher

		Additional changes by Ray Gardener, Daylon Graphics Ltd.
		December 4, 1999

  This sourcecode is taken from the "Image Resampler" program at
  http://www.daylongraphics.com/products/leveller/demo.htm

	Summary:

- Horizontal filter contributions are calculated on the fly,
  as each column is mapped from src to dst image. This lets 
  us omit having to allocate a temporary full horizontal stretch 
  of the src image.

- If none of the src pixels within a sampling region differ, 
  then the output pixel is forced to equal (any of) the source pixel.
  This ensures that filters do not corrupt areas of constant color.

- Filter weight contribution results, after summing, are 
  rounded to the nearest pixel color value instead of 
  being casted to Pixel (usually an int or char). Otherwise, 
  artifacting occurs.
*/

#include "stdafx.h" // added to the original file to allow the MFC compilation

#include "FreeImage.h"
#include "Utilities.h"

#include "fi_testDoc.h"

///////////////////////////////////////////////////////////////////////

#ifndef M_PI
	#define M_PI   3.1415926535897932384626433832795
#endif

// filter function definitions

#define	filter_support		(1.0)

double filter(double t)
{
	// f(t) = 2|t|^3 - 3|t|^2 + 1, -1 <= t <= 1
	if(t < 0.0) t = -t;
	if(t < 1.0) return((2.0 * t - 3.0) * t * t + 1.0);
	return(0.0);
}

#define	box_support		(0.5)

double box_filter(double t)
{
	if((t > -0.5) && (t <= 0.5)) return(1.0);
	return(0.0);
}

#define	triangle_support	(1.0)

double triangle_filter(double t)
{
	if(t < 0.0) t = -t;
	if(t < 1.0) return(1.0 - t);
	return(0.0);
}

#define	bell_support		(1.5)

double bell_filter(double t)		// box (*) box (*) box
{
	if(t < 0) t = -t;
	if(t < .5) return(.75 - (t * t));
	if(t < 1.5) {
		t = (t - 1.5);
		return(.5 * (t * t));
	}
	return(0.0);
}

#define	B_spline_support	(2.0)

double B_spline_filter(double t)	// box (*) box (*) box (*) box
{
	double tt;

	if(t < 0) t = -t;
	if(t < 1) {
		tt = t * t;
		return((.5 * tt * t) - tt + (2.0 / 3.0));
	} else if(t < 2) {
		t = 2 - t;
		return((1.0 / 6.0) * (t * t * t));
	}
	return(0.0);
}

double sinc(double x)
{
	x *= M_PI;
	if(x != 0) return(sin(x) / x);
	return(1.0);
}

#define	Lanczos3_support	(3.0)

double Lanczos3_filter(double t)
{
	if(t < 0) t = -t;
	if(t < 3.0) return(sinc(t) * sinc(t/3.0));
	return(0.0);
}

#define	Mitchell_support	(2.0)

#define	B	(1.0 / 3.0)
#define	C	(1.0 / 3.0)

double Mitchell_filter(double t)
{
	double tt;

	tt = t * t;
	if(t < 0) t = -t;
	if(t < 1.0) {
		t = (((12.0 - 9.0 * B - 6.0 * C) * (t * tt))
		   + ((-18.0 + 12.0 * B + 6.0 * C) * tt)
		   + (6.0 - 2 * B));
		return(t / 6.0);
	} else if(t < 2.0) {
		t = (((-1.0 * B - 6.0 * C) * (t * tt))
		   + ((6.0 * B + 30.0 * C) * tt)
		   + ((-12.0 * B - 48.0 * C) * t)
		   + (8.0 * B + 24 * C));
		return(t / 6.0);
	}
	return(0.0);
}

//	image rescaling routine

typedef struct {
	int	pixel;
	double	weight;
} CONTRIB;

typedef struct {
	int	n;		  // number of contributors
	CONTRIB	*p;	  // pointer to list of contributions
} CLIST;

//CLIST	*contrib; // array of contribution lists

//	roundcloser()

//	Round an FP value to its closest int representation.
//	General routine; ideally belongs in general math lib file.

int roundcloser(double d)
{
	// Untested potential one-liner, but smacks of call overhead
	// return fabs(ceil(d)-d) <= 0.5 ? ceil(d) : floor(d);

	// Untested potential optimized ceil() usage
	//double cd = ceil(d);
	//int ncd = (int)cd;
	//if(fabs(cd - d) > 0.5)
	//	ncd--;
	//return ncd;

	// Version that uses no function calls at all.
	int n = (int) d;
	double diff = d - (double)n;
	if(diff < 0)
		diff = -diff;
	if(diff >= 0.5)
	{
		if(d < 0)
			n--;
		else
			n++;
	}
	return n;
} // roundcloser


//	calc_x_contrib()
	
//	Calculates the filter weights for a single target column.
//	contribX->p must be freed afterwards.

//	Returns -1 if error, 0 otherwise.

int calc_x_contrib(
CLIST* contribX,			// Receiver of contrib info
double xscale,				// Horizontal zooming scale
double fwidth,				// Filter sampling width
int dstwidth,				// Target bitmap width
int srcwidth,				// Source bitmap width
double (*filterf)(double),	// Filter proc
int i						// Pixel column in source bitmap being processed
)
{
	double width;
	double fscale;
	double center, left, right;
	double weight;
	int j, k, n;

	if(xscale < 1.0)
	{
		// Shrinking image
		width = fwidth / xscale;
		fscale = 1.0 / xscale;

		contribX->n = 0;
		contribX->p = (CONTRIB *)calloc((int) (width * 2 + 1),
				sizeof(CONTRIB));
		if(contribX->p == NULL)
			return -1;

		center = (double) i / xscale;
		left = ceil(center - width);
		right = floor(center + width);
		for(j = (int)left; j <= right; ++j)
		{
			weight = center - (double) j;
			weight = (*filterf)(weight / fscale) / fscale;
			if(j < 0)
				n = -j;
			else if(j >= srcwidth)
				n = (srcwidth - j) + srcwidth - 1;
			else
				n = j;
			
			k = contribX->n++;
			contribX->p[k].pixel = n;
			contribX->p[k].weight = weight;
		}
	
	}
	else
	{
		// Expanding image
		contribX->n = 0;
		contribX->p = (CONTRIB *)calloc((int) (fwidth * 2 + 1),
				sizeof(CONTRIB));
		if(contribX->p == NULL)
			return -1;
		center = (double) i / xscale;
		left = ceil(center - fwidth);
		right = floor(center + fwidth);

		for(j = (int)left; j <= right; ++j)
		{
			weight = center - (double) j;
			weight = (*filterf)(weight);
			if(j < 0) {
				n = -j;
			} else if(j >= srcwidth) {
				n = (srcwidth - j) + srcwidth - 1;
			} else {
				n = j;
			}
			k = contribX->n++;
			contribX->p[k].pixel = n;
			contribX->p[k].weight = weight;
		}
	}
	return 0;
} // calc_x_contrib


//	zoom()

//	Resizes bitmaps while resampling them.
//	Returns -1 if error, 0 if success.

BOOL zoom(FIBITMAP* src_dib,
			   
        unsigned wd, unsigned hd, CPoint OffscreenSize, BYTE* pBits,
		
		RGBQUAD *dst_pal, double (*filterf)(double), double fwidth)
{
	BYTE* tmp;
	double xscale, yscale;		// zoom scale factors
	int xx;
	int i, j, k;		// loop variables
	int n;				// pixel number
	double center, left, right;	// filter calculation variables
	double width, fscale, weight;	// filter calculation variables
	BYTE pel, pel2;
	int bPelDelta;
	CLIST	*contribY;		// array of contribution lists
	CLIST	contribX;

BYTE  *lines, *lined;

BYTE *datas = FreeImage_GetBits(src_dib);

unsigned ws = FreeImage_GetWidth(src_dib);

unsigned hs = FreeImage_GetHeight(src_dib);

int bpp = FreeImage_GetBPP(src_dib);

int btpp = bpp/8;

unsigned wpls = FreeImage_GetPitch(src_dib);

unsigned wpld = CalculatePitch(CalculateLine(wd, bpp == 1 ? 8 : bpp));

int color_mode = FreeImage_GetColorType(src_dib);

unsigned wpld_big = CalculatePitch(CalculateLine(OffscreenSize.x,
								   bpp == 1 ? 8 : bpp));

unsigned dif2_w = (wpld_big - wpld)/2;

// center the bitmap horizontally (on the background bitmap)
pBits += dif2_w - dif2_w % 3; // why is it needed here to align
// the image at the 3-byte boundary ???

// center the bitmap vertically (on the background bitmap)
pBits += (OffscreenSize.y - hd)/2 * wpld_big;

	// create intermediate column to hold horizontal dst column zoom
	tmp = (BYTE*)malloc(hs * sizeof(BYTE));

	if(tmp == NULL)
		return false;

	xscale = (double) wd / (double) ws;

	// Build y weights
	// pre-calculate filter contributions for a column
	contribY = (CLIST *)calloc(hd, sizeof(CLIST));

	if(contribY == NULL)
	{
		free(tmp);
		return false;
	}

	yscale = (double) hd / (double) hs;

	if(yscale < 1.0)
	{
		width = fwidth / yscale;
		fscale = 1.0 / yscale;

		for(i = 0; i < hd; ++i)
		{
			contribY[i].n = 0;
			contribY[i].p = (CONTRIB *)calloc((int) (width * 2 + 1),
					sizeof(CONTRIB));

			if(contribY[i].p == NULL)
			{
				free(tmp);
				free(contribY);
				return false;
			}

			center = (double) i / yscale;
			left = ceil(center - width);
			right = floor(center + width);

			for(j = (int)left; j <= right; ++j)
			{
				weight = center - (double) j;
				weight = (*filterf)(weight / fscale) / fscale;

				if(j < 0)
				{
					n = -j;
				}				
				else if(j >= hs)
				{
					n = (hs - j) + hs - 1;
				}
				else
				{
					n = j;
				}

				k = contribY[i].n++;
				contribY[i].p[k].pixel = n;
				contribY[i].p[k].weight = weight;
			}
		}		
	} 
	else
	{
		for(i = 0; i < hd; ++i)
		{
			contribY[i].n = 0;
			contribY[i].p = (CONTRIB *)calloc((int) (fwidth * 2 + 1),
					sizeof(CONTRIB));
			if(contribY[i].p == NULL)
			{
				free(tmp);
				free(contribY);
				return false;
			}
			center = (double) i / yscale;
			left = ceil(center - fwidth);
			right = floor(center + fwidth);

			for(j = (int)left; j <= right; ++j)
			{
				weight = center - (double) j;
				weight = (*filterf)(weight);
				if(j < 0)
				{
					n = -j;
				}
				else if(j >= hs)
				{
					n = (hs - j) + hs - 1;
				}
				else
				{
					n = j;
				}
				k = contribY[i].n++;
				contribY[i].p[k].pixel = n;
				contribY[i].p[k].weight = weight;
			}
		}
	}

	for(xx = 0; xx < wd; xx++)
	{
		if(0 != calc_x_contrib(&contribX, xscale, fwidth, 
								//dst->xsize, src->xsize, filterf, xx))
								wd, ws, filterf, xx))
		{
			break;
		}
		// Apply horz filter to make dst column in tmp.
		for(k = 0; k < hs; ++k)
		{
			weight = 0.0;
			bPelDelta = FALSE;

			//pel = get_pixel(src, contribX.p[0].pixel, k);
			lines = datas + k * wpls;

			pel = dst_pal[*(lines + contribX.p[0].pixel)].rgbRed;

			for(j = 0; j < contribX.n; ++j)
			{
				//pel2 = get_pixel(src, contribX.p[j].pixel, k);
				pel2 = dst_pal[*(lines + contribX.p[j].pixel)].rgbRed;

				if(pel2 != pel)
					bPelDelta = TRUE;
				weight += pel2 * contribX.p[j].weight;
			}
			weight = bPelDelta ? roundcloser(weight) : pel;

			//tmp[k] = (Pixel)CLAMP(weight, BLACK_PIXEL, WHITE_PIXEL);
			tmp[k] = (BYTE)MIN(MAX((int)0, (int)(weight + 0.5)), (int)255);

		} // next row in temp column

		free(contribX.p);		

		// The temp column has been built. Now stretch it 
		// vertically into dst column.
		for(i = 0; i < hd; ++i)
		{
			weight = 0.0;
			bPelDelta = FALSE;

			lined = pBits + i * wpld_big;			

			pel = tmp[contribY[i].p[0].pixel];

			for(j = 0; j < contribY[i].n; ++j)
			{
				pel2 = tmp[contribY[i].p[j].pixel];

				if(pel2 != pel)
					bPelDelta = TRUE;
				weight += pel2 * contribY[i].p[j].weight;
			}
			weight = bPelDelta ? roundcloser(weight) : pel;			

			//put_pixel(dst, xx, i,
				//(Pixel)CLAMP(weight, BLACK_PIXEL, WHITE_PIXEL));			

			lined[xx] = (BYTE)MIN(MAX((int)0, (int)(weight + 0.5)), (int)255);
				

		} // next dst row
	} // next dst column

	free(tmp);

	// free the memory allocated for vertical filter weights
	for(i = 0; i < hd; ++i) free(contribY[i].p);
	free(contribY);

	return true;
} // zoom


BOOL CFi_testDoc::FilterRcg(FIBITMAP* src_dib,
			   
        unsigned wd, unsigned hd, CPoint OffscreenSize, BYTE* pBits,
		
		RGBQUAD *dst_pal, int filter_num)
{
double (*f)(double) = filter;

double s = filter_support;

switch(filter_num)
{
case 1: f=box_filter; s=box_support; break;
case 2: f=triangle_filter; s=triangle_support; break;
case 3: f=bell_filter; s=bell_support; break;
case 4: f=B_spline_filter; s=B_spline_support; break;
case 5: f=filter; s=filter_support; break;
case 6: f=Lanczos3_filter; s=Lanczos3_support; break;
case 7: f=Mitchell_filter; s=Mitchell_support; break;
}

return zoom(src_dib, wd, hd, OffscreenSize, pBits, dst_pal, f, s);
}