drawutil.c

傅立叶变换和小波变换是图像压缩的重要工具。该代大戏是利用小波变换进行图像压缩。

	if ( y1 > y2 )
	{
		swapf(x1,x2);
		swapf(y1,y2);
	}
	m = (x2 - x1)/(y2 - y1);
	if ( y1 < 0.0f )
	{
		if ( y2 <= 0.0f )
			return false;

		x1 += m * ABS(y1);
		y1 = 0.0f;
	}
	if ( y2 >= W->Height-1 )
	{
		if ( y1 >= W->Height-1 )
			return false;

		x2 -= m * ABS(y2 - (W->Height-1));
		y2 = W->Height-1;
	}

	assert( min(x1,x2) >= 0 );
	assert( max(x1,x2) <  W->Width );
	assert( min(y1,y2) >= 0 );
	assert( max(y1,y2) <  W->Height );

	assert( min(x1,x2) >= min(*px1,*px2) );
	assert( max(x1,x2) <= max(*px1,*px2) );
	assert( min(y1,y2) >= min(*py1,*py2) );
	assert( max(y1,y2) <= max(*py1,*py2) );

	*px1 = (uint)(x1 + 0.5f);
	*py1 = (uint)(y1 + 0.5f);
	*px2 = (uint)(x2 + 0.5f);
	*py2 = (uint)(y2 + 0.5f);

return true;
}

bool DrawUtil_ClipSegV(const Driver_Window * W,Vec2d *v1,Vec2d *v2)
{
float x1,y1,x2,y2,m;

	x1 = v1->x; y1 = v1->y;
	x2 = v2->x; y2 = v2->y;
	if ( x1 > x2 )
	{
		swapf(x1,x2);
		swapf(y1,y2);
	}

	m = (y2 - y1)/(x2 - x1);

	if ( x1 < 0.0f )
	{
		// clip against the left
		if ( x2 <= 0.0f )
			return false;

		y1 += m * ABS(x1);
		x1 = 0.0f;
	}
	if ( x2 >= W->Width-1 )
	{
		// clip right
		if ( x1 >= W->Width-1 )
			return false;
		y2 -= m * (x2 - (W->Width-1));
		x2 = W->Width-1;
	}

	if ( y1 > y2 )
	{
		swapf(x1,x2);
		swapf(y1,y2);
	}
	m = (x2 - x1)/(y2 - y1);
	if ( y1 < 0.0f )
	{
		if ( y2 <= 0.0f )
			return false;

		x1 += m * ABS(y1);
		y1 = 0.0f;
	}
	if ( y2 >= W->Height-1 )
	{
		if ( y1 >= W->Height-1 )
			return false;

		x2 -= m * ABS(y2 - (W->Height-1));
		y2 = W->Height-1;
	}

	assert( min(x1,x2) >= 0 );
	assert( max(x1,x2) <  W->Width );
	assert( min(y1,y2) >= 0 );
	assert( max(y1,y2) <  W->Height );

	v1->x = x1; v1->y = y1;
	v2->x = x2;	v2->y = y2;

return true;
}

//------------------------------------------------------------------------
// Noise

/*****

notes :

1.	Noise (when aligned on unit boundaries) is like a 
	quad-subdivision thing, with nice lerping between the points;
	it's hard to think about a quad-subdivision implementation, tho,
	cuz you have to share the corner values (the usual quad-neighbor problem)

2.	Perlin's version doesn't actually store a big noise table,	
	he just has a 1d noise table and then uses a multi-d hashing index to the table

******/

// Cosine interpolation is smoother than linear;
//	first derivatives are continuous across the control points
float INLINE NoiseLerp(float a,float b,float frac)
{
	frac *= PI;
	frac = (1.0f - DrawUtil_Cos(frac))*0.5f;
	return  a + (b-a)*frac;
}

#define NOISETABLE1D_DIM	(1024)
#define NOISETABLE1D_MASK	(NOISETABLE1D_DIM-1)

static float NoiseTable1d[NOISETABLE1D_DIM+1];

static void MakeNoise1d(void)
{
int i;
	DrawUtil_Randomize();
	for(i=0;i<NOISETABLE1D_DIM;i++)
	{
		NoiseTable1d[i] = DrawUtil_RandSignedUnitFloat();
	}
	NoiseTable1d[NOISETABLE1D_DIM] = NoiseTable1d[0]; // wrap it
}

float REGCALL DrawUtil_BasicNoise1d(float f)
{
int i;
	i = (int)(f);
	f -= (float)i;
	i &= NOISETABLE1D_MASK;

return NoiseLerp(NoiseTable1d[i],NoiseTable1d[i+1],f);
}

float REGCALL DrawUtil_SpectralNoise1d(float x)
{
return  DrawUtil_BasicNoise1d(x) + 
		DrawUtil_BasicNoise1d(x*2)*0.5f +
		DrawUtil_BasicNoise1d(x*4)*0.25f +
		DrawUtil_BasicNoise1d(x*8)*0.125f;
}

#define NOISETABLE2d_MASK	(0x3F)
#define NOISETABLE2d_DIM	(NOISETABLE2d_MASK+1)

static float NoiseTable2d[(NOISETABLE2d_DIM+1)][(NOISETABLE2d_DIM+1)];

static void MakeNoise2d(void)
{
int x,y;
	DrawUtil_Randomize();
	for(x=0;x<NOISETABLE2d_DIM;x++)
	{
		for(y=0;y<NOISETABLE2d_DIM;y++)
		{
			NoiseTable2d[x][y] = DrawUtil_RandSignedUnitFloat();
		}
		NoiseTable2d[x][NOISETABLE2d_DIM] = NoiseTable2d[x][0]; // wrap it
	}
	for(y=0;y<=NOISETABLE2d_DIM;y++)
	{
		NoiseTable2d[NOISETABLE2d_DIM][y] = NoiseTable2d[0][y];
	}
}

float REGCALL DrawUtil_BasicNoise2d(float x,float y)
{
int ix,iy;
float N,S;
	ix = (int)(x);
	iy = (int)(y);
	x -= ix;
	y -= iy;
	ix &= NOISETABLE2d_MASK;
	iy &= NOISETABLE2d_MASK;

	S = NoiseLerp(NoiseTable2d[ix][iy],NoiseTable2d[ix+1][iy],x);
	iy++;
	N = NoiseLerp(NoiseTable2d[ix][iy],NoiseTable2d[ix+1][iy],x);

return NoiseLerp(S,N,y);
}

float REGCALL DrawUtil_DeluxeSpectralNoise2d(float x,float y,float noisiness,int levels)
{
float n,a;
	a = 1.0f;
	n = 0.0f;
	while(levels--)
	{
		n += DrawUtil_BasicNoise2d(x,y) * a;
		a *= noisiness;
		x *= 2; y *= 2;
	}
return n;
}

float REGCALL DrawUtil_SpectralNoise2d(float x,float y)
{
return  DrawUtil_BasicNoise2d(x  ,y) + 
		DrawUtil_BasicNoise2d(x*2,y*2)*0.5f +
		DrawUtil_BasicNoise2d(x*4,y*4)*0.25f +
		DrawUtil_BasicNoise2d(x*8,y*8)*0.125f;
}

#define NOISETABLE3d_MASK	(0xF)
#define NOISETABLE3d_DIM	(NOISETABLE3d_MASK+1)

static float NoiseTable3d[NOISETABLE3d_DIM][NOISETABLE3d_DIM][NOISETABLE3d_DIM];

static void MakeNoise3d(void)
{
int x,y,z;
	DrawUtil_Randomize();
	for(x=0;x<NOISETABLE3d_DIM;x++)
	{
		for(y=0;y<NOISETABLE3d_DIM;y++)
		{
			for(z=0;z<NOISETABLE3d_DIM;z++)
			{
				NoiseTable3d[x][y][z] = DrawUtil_RandSignedUnitFloat();
			}
		}
	}
}

float REGCALL DrawUtil_BasicNoise3d(float x,float y,float z)
{
int ix,iy,iz,ny;
float N,S,A,B;

	ix = (int)(x);
	iy = (int)(y);
	iz = (int)(z);
	x -= (float)ix;
	y -= (float)iy;
	z -= (float)iz;
	ix &= NOISETABLE3d_MASK;
	iy &= NOISETABLE3d_MASK;
	iz &= NOISETABLE3d_MASK;

	ny = (iy+1)&NOISETABLE3d_MASK;

	S = NoiseLerp(NoiseTable3d[ix][iy][iz],NoiseTable3d[(ix+1)&NOISETABLE3d_MASK][iy][iz],x);
	N = NoiseLerp(NoiseTable3d[ix][ny][iz],NoiseTable3d[(ix+1)&NOISETABLE3d_MASK][ny][iz],x);
	A = NoiseLerp(S,N,y);

	iz = (iz+1)&NOISETABLE3d_MASK;

	S = NoiseLerp(NoiseTable3d[ix][iy][iz],NoiseTable3d[(ix+1)&NOISETABLE3d_MASK][iy][iz],x);
	N = NoiseLerp(NoiseTable3d[ix][ny][iz],NoiseTable3d[(ix+1)&NOISETABLE3d_MASK][ny][iz],x);
	B = NoiseLerp(S,N,y);

return NoiseLerp(A,B,z);
}

float REGCALL DrawUtil_SpectralNoise3d(float x,float y,float z)
{
return  DrawUtil_BasicNoise3d(x  ,y  ,z) + 
		DrawUtil_BasicNoise3d(x*2,y*2,z*2)*0.5f +
		DrawUtil_BasicNoise3d(x*4,y*4,z*4)*0.25f +
		DrawUtil_BasicNoise3d(x*8,y*8,z*8)*0.125f;
}

//--------------------------------------------------------------
// Bezier

/*

very good refining Bezier Drawing routines.  Lets you pass in a
tolerated pixel error (or 0 for (near) pixel-perfect drawing with
minimal computation).  Note that even in pixel-perfect mode this
is much faster than full computation, because we early out on the 
straight portions of the curve.

the only real opportunity for speed increase is by sharing
some of the computation in the DrawLine routine (eg. when a bezier
set doesn't clip against the walls, not of its kids does; (note that
a 'set' is not a segment)).

*/

static Driver_Window * DrawBezier_Window;
static uword DrawBezier_r,DrawBezier_g,DrawBezier_b;
static float DrawBezier_ErrorSquared;

static void DrawUtil_DrawBezier_r(	const Driver_BezierPoint *p,
									const Driver_BezierPoint *q,
									float interval)
{
Driver_BezierPoint c;
float s,ex,ey;

	// take two points with 2nd d's and make a new point
	//	creates two sub-curves, split at the parametric midpoint

	// 2nd derivative is just the average	
	c.curve.x = 0.5f * (p->curve.x + q->curve.x);
	c.curve.y = 0.5f * (p->curve.y + q->curve.y);

	s = interval*interval*(1/8.0f);
	ex = s * c.curve.x;
	ey = s * c.curve.y;
	c.pos.x = 0.5f * (p->pos.x + q->pos.x) - ex;
	c.pos.y = 0.5f * (p->pos.y + q->pos.y) - ey;

	s = ex*ex + ey*ey;

	if ( s > DrawBezier_ErrorSquared )
	{
		interval *= 0.5f;
		DrawUtil_DrawBezier_r(p,&c,interval);
		DrawUtil_DrawBezier_r(&c,q,interval);
	}
	else
	{
		DrawUtil_DrawLineV(DrawBezier_Window,&(p->pos),&(c.pos),DrawBezier_r,DrawBezier_g,DrawBezier_b);
		DrawUtil_DrawLineV(DrawBezier_Window,&(c.pos),&(q->pos),DrawBezier_r,DrawBezier_g,DrawBezier_b);
	}

}

void DrawUtil_DrawBezierBy2ndD(	Driver_Window * pWindow,
							const Driver_BezierPoint *p,
							const Driver_BezierPoint *q,
							uword r,uword g,uword b,
							float err)
{

	DrawBezier_Window = pWindow;
	DrawBezier_r = r;
	DrawBezier_g = g;
	DrawBezier_b = b;
	DrawBezier_ErrorSquared = err * err + 0.5f;

	DrawUtil_DrawBezier_r(p,q,1.0f);
}

void DrawUtil_DrawBezierByPoints(  Driver_Window * pWindow,
							const Vec2d *p0,
							const Vec2d *p1,
							const Vec2d *p2,
							const Vec2d *p3,
							uword r,uword g,uword b,
							float err)
{
Driver_BezierPoint p,q;
	p.pos = *p0;
	p.curve.x = 6.0f * ( p2->x + p0->x - 2 * p1->x );
	p.curve.y = 6.0f * ( p2->y + p0->y - 2 * p1->y );
	q.pos = *p3;
	q.curve.x = 6.0f * ( p3->x + p1->x - 2 * p2->x );
	q.curve.y = 6.0f * ( p3->y + p1->y - 2 * p2->y );
	DrawUtil_DrawBezierBy2ndD(pWindow,&p,&q,r,g,b,err);
}

void INLINE Driver_MakeNormalDistance(Vec2d *v,
								const Vec2d *a,
								const Vec2d *b)
{
float inv;
	v->x = b->x - a->x;
	v->y = b->y - a->y;
	inv = 1.0f / fSqrt( v->x * v->x + v->y * v->y );
	v->x *= inv;
	v->y *= inv;
}

float INLINE Driver_PerpDotSquared(const Vec2d *normal,const Vec2d *vec)
{
float len;

//	len = (vec->x * vec->x + vec->y * vec->y)
//			- fSquare(vec->x * normal->x + vec->y * normal->y);
//	assert( len >= 0.0f );
// return len;

	// this is just the same as :

	len = Vec2d_CrossProduct(vec,normal);
	
return len*len;
}

static void DrawUtil_DrawBezierCasteljau_r(const Vec2d *p0,
								const Vec2d *p1,
								const Vec2d *p2,
								const Vec2d *p3)
{

	// weird arrangement so that these locals go out
	//	of scope when you goto the refine section
	{
	Vec2d v,u;
	float err;
		Driver_MakeNormalDistance(&v,p0,p3);
		u.x = p1->x - p0->x;
		u.y = p1->y - p0->y;
		err = Driver_PerpDotSquared(&v,&u);
		if ( err > DrawBezier_ErrorSquared )
			goto refine;

		u.x = p2->x - p0->x;
		u.y = p2->y - p0->y;
		err = Driver_PerpDotSquared(&v,&u);
		if ( err > DrawBezier_ErrorSquared )
			goto refine;
	}

	
	{
	Vec2d mid;
		// compute the midpoint and draw two segments
		//  we must do this to make the quality the same as the Middle division method,
		//	cuz there we can draw two segs immediately whenever you decide to stop refining
		mid.x = (1.0f/8.0f) * ( p0->x + p3->x + 3*(p1->x + p2->x) );
		mid.y = (1.0f/8.0f) * ( p0->y + p3->y + 3*(p1->y + p2->y) );
		DrawUtil_DrawLineV(DrawBezier_Window,p0,&mid,DrawBezier_r,DrawBezier_g,DrawBezier_b);
		DrawUtil_DrawLineV(DrawBezier_Window,p3,&mid,DrawBezier_r,DrawBezier_g,DrawBezier_b);
	}

return;
	
	refine:

	// deCasteljau makes two new sub-bezier curves

	{
	Vec2d q1,q2,q3,r1,r2;
		q1.x = (p0->x + p1->x)*0.5f;
		r2.x = (p2->x + p3->x)*0.5f;
		q2.x = (q1.x + (p1->x + p2->x)*0.5f)*0.5f;
		r1.x = (r2.x + (p1->x + p2->x)*0.5f)*0.5f;
		q3.x = (q2.x + r1.x)*0.5f;
		q1.y = (p0->y + p1->y)*0.5f;
		r2.y = (p2->y + p3->y)*0.5f;
		q2.y = (q1.y + (p1->y + p2->y)*0.5f)*0.5f;
		r1.y = (r2.y + (p1->y + p2->y)*0.5f)*0.5f;
		q3.y = (q2.y + r1.y)*0.5f;

		DrawUtil_DrawBezierCasteljau_r(p0,&q1,&q2,&q3);
		DrawUtil_DrawBezierCasteljau_r(&q3,&r1,&r2,p3);
	}
}

void DrawUtil_DrawBezierCasteljau(  Driver_Window * pWindow,
							const Vec2d *p0,
							const Vec2d *p1,
							const Vec2d *p2,
							const Vec2d *p3,
							uword r,uword g,uword b,
							float err)
{
	DrawBezier_Window = pWindow;
	DrawBezier_r = r;
	DrawBezier_g = g;
	DrawBezier_b = b;
	DrawBezier_ErrorSquared = err * err + 0.5f;

	DrawUtil_DrawBezierCasteljau_r(p0,p1,p2,p3);
}

//--------------------------------------------------------------
// Smooth

uword INLINE blend4_555(uword a,uword b,uword c,uword d)
{
return	((((a>>10)     + (b>>10)     + (c>>10)     + (d>>10)     + 2)>>2)<<10) |
		(((((a>>5)&31) + ((b>>5)&31) + ((c>>5)&31) + ((d>>5)&31) + 2)>>2)<<5) |
		((((a)&31)     + ((b)&31)    + ((c)&31)    + ((d)&31)    + 2)>>2);
}

uword INLINE blend4_565(uword a,uword b,uword c,uword d)
{
return	((((a>>11)     + (b>>11)     + (c>>11)     + (d>>11)     + 2)>>2)<<11) |
		(((((a>>5)&63) + ((b>>5)&63) + ((c>>5)&63) + ((d>>5)&63) + 2)>>2)<<5) |
		((((a)&31)     + ((b)&31)    + ((c)&31)    + ((d)&31)    + 2)>>2);
}

void DrawUtil_Smooth(Driver_Window * Window)
{
uword * ptr;
uint Stride,x,y;
	Stride = Window->Stride>>1;

	for(y=1;y<(Window->Height -1);y++)
	{
		ptr = (uword *)Window->Buffer;
		ptr += y * Stride + 1;
			
		switch(Window->Mode)
		{
		case DRIVER_MODE_555:
			for(x=1;x<(Window->Width -1);x++)
			{
			uword N,S,E,W;
				N = ptr[ Stride];
				S = ptr[-(int)Stride];
				W = ptr[ 1];
				E = ptr[-1];

				*ptr++ = blend4_555(N,S,E,W);
			}
			break;

		case DRIVER_MODE_565:
			for(x=1;x<(Window->Width -1);x++)
			{
			uword N,S,E,W;
				N = ptr[ Stride];
				S = ptr[-(int)Stride];
				W = ptr[ 1];
				E = ptr[-1];

				*ptr++ = blend4_565(N,S,E,W);
			}
			break;

		default:
			assert(0);
			break;
		}
	}
}

//--------------------------------------------------------------