The test images below are created with sinctest.c given further below. By looking at the center and diagonals of the images, artifacts from filtering can be noticed as ringing going in the opposite direction:
/*
* Copyright (c) 1996-2001 Thomas E. Burge. All rights reserved.
*
* Affine (R) is a registered trademark of Thomas E. Burge
*
* THIS SOFTWARE IS DISTRIBUTED "AS-IS" WITHOUT WARRANTY OF ANY KIND
* AND WITHOUT ANY GUARANTEE OF MERCHANTABILITY OR FITNESS FOR A
* PARTICULAR PURPOSE.
*
* In no event shall Thomas E. Burge be liable for any indirect or
* consequential damages or loss of data resulting from use or performance
* of this software.
*
* Permission to use, copy, modify, distribute, and sell this software and
* its documentation for any purpose is hereby granted without fee, provided
* that the following copyright notices and this permission notice appear in
* all copies of the software and related documentation:
*
* The Affine (R) Libraries and Tools are
* Copyright (c) 1995-2001 Thomas E. Burge.
* All rights reserved.
* Affine (R) is a registered trademark of Thomas E. Burge.
*
* Also refer to any additional requirements presently set by Pixar
* in regards to the RenderMan (R) Interface Procedures and Protocol.
*
*
* Affine Toolkit
*
* FILE: pixfiltr.c
*
* DESCRIPTION: RiPixelFilter() functions.
*
* Contains:
*
* RiBoxFilter()
* RiTriangleFilter()
* RiCatmullRomFilter()
* RiGaussianFilter()
* RiSincFilter()
*
* References:
* [COOK84a] Cook, Robert L., Antialiasing by Stochastic Sampling,
* Computer Graphics Department, Computer Division,
* Lucasfilm Ltd, 1984.
* [COOK84b] Cook, Robert L., Patch Work: Technical Memo #118,
* Computer Division, Lucasfilm Ltd, 1984.
* [COOK85] Cook, Robert L., Stochastic Sampling in Computer
* Graphics, ACM Transactions on Graphics, Vol. 5, No. 1,
* January 1986, pp 51-72.
* [GLAS95] Glassner, Andrew S. (CWRU '84), Principles of
* Digital Image Synthesis, Morgan Kaufmann
* Publishers, Inc., 1995.
* [MITC88] Mitchell, Don P. and Netravali, Arun N., Computer
* Graphics, Vol. 22, No. 4, August 1988, pp. 221-228.
* [NETR95] Netravali, Arun N. and Haskell, Barry G., Digital
* Pictures: Representation, Compression and Standards,
* 2nd Edition, Plenum Press, 1995, p. 12.
* [PIXA89] Pixar, The RenderMan Interface, Version 3.1,
* Richmond, CA, September 1989.
* [PRAT91] Pratt, William K., Digital Image Processing,
* 2nd Edition, John Wiley & Sons, 1991, p. 98,
* pp. 115-116.
* [STEV93] Stevens, Roger T., Quick Reference to Computer
* Graphics Terms, Academic Press Professional, 1993.
* [UPST89] Upstill, Steve, The Renderman Companion: A
* Programmer's Guide to Realistic Computer Graphics,
* Addison Wesley, 1989.
*
*
* The RenderMan (R) Interface Procedures and Protocol are:
* Copyright 1988, 1989, Pixar
* All Rights Reserved
*
* Renderman (R) is a registered trademark of Pixar
*
*/
#include <stdio.h>
#include <math.h>
typedef float RtFloat;
typedef RtFloat (*RtFilterFunc)(RtFloat, RtFloat, RtFloat, RtFloat );
RtFloat
RiBoxFilter( RtFloat x, RtFloat y, RtFloat xwidth, RtFloat ywidth ),
RiTriangleFilter( RtFloat x, RtFloat y, RtFloat xwidth, RtFloat ywidth ),
RiCatmullRomFilter( RtFloat x, RtFloat y, RtFloat xwidth, RtFloat ywidth ),
RiGaussianFilter( RtFloat x, RtFloat y, RtFloat xwidth, RtFloat ywidth ),
RiSincFilter( RtFloat x, RtFloat y, RtFloat xwidth, RtFloat ywidth )
#ifdef PRMAN
,
RiBesselFilter(RtFloat x, RtFloat y, RtFloat xwidth, RtFloat ywidth),
RiDiskFilter(RtFloat x, RtFloat y,RtFloat xwidth, RtFloat ywidth);
#else
;
#endif
#define WIDTH 5.0
#define XWIDTH 5.0
#define YWIDTH 5.0
#define TWO_DIM_
#define MIN( a, b ) ((a)<(b) ? (a) : (b))
/* References:
*
* [UPST89] pp. 176-178.
* [GLAS95] p. 207.
*
*/
RtFloat MyBoxFilter( RtFloat x, RtFloat y,
RtFloat xwidth, RtFloat ywidth )
{
/* [UPST89] -- (RC p. 178) says that x and y will be in the
* following intervals:
* -xwidth/2 <= x <= xwidth/2
* -ywidth/2 <= y <= ywidth/2
* These constraints on x and y really simplifies the
* the following code to just return (1.0).
*
*/
#define SIMPLIFY_CODE
#ifdef SIMPLIFY_CODE
return MIN( (fabs(x) <= xwidth/2. ? 1. : 0.),
(fabs(y) <= ywidth/2. ? 1. : 0.) );
#else
return 1.0;
#endif
}
/* References:
*
* [NETR95] p. 12.
* [UPST89] pp. 176-178.
*
*/
RtFloat MyTriangleFilter( RtFloat x, RtFloat y,
RtFloat xwidth, RtFloat ywidth )
{
double hxw = xwidth/2.0;
double hyw = ywidth/2.0;
double absx = fabs(x);
double absy = fabs(y);
/* This function can be simplified as well by not worrying about
* returning zero if the sample is beyond the filter window.
*/
return MIN( (absx <= hxw ? (hxw - absx)/hxw : 0.0),
(absy <= hyw ? (hyw - absy)/hyw : 0.0) );
}
/* References:
*
* [GLAS95] pp. 529-531, p. 536.
* [NETR95] p. 12.
* [PRAT91] p. 115.
* [UPST89] pp. 176-178.
*
*/
RtFloat MyCatmullRomFilter( RtFloat x, RtFloat y,
RtFloat xwidth, RtFloat ywidth )
{
/*
* From page 223 of [MITC88]
*
* if abs(d) < 1
* f(d) = 1/6*( (12-9*B-9*C)*abs(d*d*d)
* + (-18 + 12*B + 6*C)*d*d + (6-2*B) )
*
* if 1 <= abs(d) < 2
* f(d) = 1/6*( (-B-6*C)*abs(d*d*d)
* + (6*B + 30*C)*d*d
* + (-12*B - 48*C)*d
* + (8*B + 24*C) )
*
* otherwise f(d)=0
*
* -------------------------------------------------------------
* When B = 0.0 and C = 0.5 the filter is a Catmull-Rom cubic spline.
*
* if abs(d) < 1
* f(d) = 1/6*[ (12-3)*abs(d*d*d) + (-18 + 3)*d*d + (6) ]
*
* if 1 <= abs(d) < 2
* f(d) = 1/6*[ (-3)*abs(d*d*d) + (15)*d*d + (-24)*d + (12) ]
*
* otherwise f(d)=0
* -------------------------------------------------------------
* Simplifying:
*
* if abs(d) < 1
* f(d) = (3/2)*abs(d*d*d) - (5/2)*d*d + 1
*
* if 1 <= abs(d) <2
* f(d) = (-0.5)*abs(d*d*d) + (5/2)*d*d - 4*abs(d) + 2
*
* otherwise f(d)=0
*
*/
#ifndef RISPEC3_2
RtFloat d,d2;
d2 = x*x+y*y; /* d*d */
d = sqrt(d2); /* distance from origin */
if ( d < 1 )
return (1.5*d*d2-2.5*d2+1.0);
else if ( d < 2 )
return (-d*d2*0.5 + 2.5*d2 - 4.0*d + 2.0);
else
return 0.0;
#else
/* RI SPec 3.2 */
RtFloat r2 = (x*x+y*y);
RtFloat r = sqrt(r2);
return (r>=2.0)?0.0:
(r<1.0)?(3.0*r*r2-5.0*r2+2.0):(-r*r2+5.0*r2-8.0*r+4.0);
#endif
}
/* References:
*
* [COOK84a] pp. 8-10.
* [GLAS95] p. 208., pp. 524-526, p. 536.
* [PRATT91] pp. 116.
* [UPST89] pp. 176-178.
* [STEV93] p. 98.
*
*/
RtFloat MyGaussianFilter( RtFloat x, RtFloat y,
RtFloat xwidth, RtFloat ywidth )
{
/*
* d = distance from origin
* w = filterwidth ([COOK84a] article used 1.5)
* For here use sqrt( (xwidth/2)*(xwidth/2) + (ywidth/2)*(ywidth/2) ).
* Simplifying:
*
* w = sqrt( (xwidth*xwidth)/2 + (ywidth*ywidth)/2 )
* w = sqrt( (xwidth*xwidth + ywidth*ywidth)/2 )
* w*w = (xwidth*xwidth + ywidth*ywidth)/2
*
* if (d < filterwidth) then 0
* else exp(-d*d) - exp(-w*w)
*
*/
#if 0
/* This version falls faster than the PRMan (RI Spec 3.2) does. */
double d,d2,w,w2;
/* d = sqrt(x*x+y*y), d*d = (x*x+y*y) */
d2 = (x*x+y*y);
d = sqrt( d2 );
w2 = 0.5*(xwidth*xwidth + ywidth*ywidth);
w = sqrt( w2 );
if (d > w)
return 0.0;
else
return exp(-d2) - exp(-w2);
#else
x *= 2.0/xwidth;
y *= 2.0/ywidth;
return exp(-2.0*(x*x+y*y));
#endif
}
/* References:
*
* [PRATT91] pp. 98.
* [UPST89] pp. 176-178.
*
*/
RtFloat MySincFilter( RtFloat x, RtFloat y,
RtFloat xwidth, RtFloat ywidth )
{
#if 0
/* Matches RDC and BMRT but does not truncate. */
double d;
d = sqrt(x*x+y*y);
if (d != 0)
return sin(M_PI*d)/(M_PI*d);
else
return 1.0;
#else
RtFloat s,t;
/* RI Spec 3.2 does not properly scale with PI. */
x *= M_PI;
y *= M_PI;
if (x>-0.001 && x<0.001)
s=1.0;
else
s = sin(x)/x;
if (y>-0.001 && y<0.001)
t=1.0;
else
t = sin(y)/y;
return s*t;
#endif
}
RtFloat myDiskFilter( RtFloat x, RtFloat y,
RtFloat xwidth, RtFloat ywidth )
{
double d,xx,yy;
xx = x*x;
yy = y*y;
xwidth *= 0.5;
ywidth *= 0.5;
d = (xx)/(xwidth*xwidth)+(yy)/(ywidth*ywidth);
if ( d<1.0 )
{
return 1.0;
}
else
{
return 0.0;
}
}
/* Modified version of the RI Spec 3.2 sinc filter but also with a
* raised cosine window applied. (von Hann Window)
*/
RtFloat vonHannSincFilter( RtFloat x, RtFloat y,
RtFloat xwidth, RtFloat ywidth )
{
/* Raised cosine window (von Hann Window) described in
* R. W. Hamming "Digital Filters" 3rd Ed. p.117.
* w = (1+cos(M_PI*x))/2;
*/
if (x!=0.0)
{
x *= M_PI;
x = 0.5*(1+cos(x/xwidth)) * sin(x)/x;
}
else
{
x = 1.0;
}
if (y!=0.0)
{
y *= M_PI;
y = 0.5*(1+cos(y/ywidth)) * sin(y)/y;
}
else
{
y = 1.0;
}
return x*y;
}
RtFloat myBesselFilter( RtFloat x, RtFloat y, RtFloat xwidth, RtFloat ywidth )
{
double d,w,xx,yy;
xx = x*x;
yy = y*y;
xwidth *= 0.5;
ywidth *= 0.5;
w = (xx)/(xwidth*xwidth)+(yy)/(ywidth*ywidth);
if ( w<1.0 )
{
d = sqrt(xx+yy);
if (d!=0.0)
{
/* Half cosine window. */
w = cos(0.5*M_PI*sqrt(w));
return w * 2*j1(M_PI*d)/d;
}
else
{
return M_PI;
}
}
else
{
return 0.0;
}
}
/* Modified version of the RI Spec 3.2 sinc filter but also with a
* half cosine window applied. (half cosine Window)
*/
RtFloat halfCosineSincFilter( RtFloat x, RtFloat y,
RtFloat xwidth, RtFloat ywidth )
{
/* Uses a -PI to PI cosine window. */
if (x!=0.0)
{
x *= M_PI;
x = cos(0.5*x/xwidth) * sin(x)/x;
}
else
{
x = 1.0;
}
if (y!=0.0)
{
y *= M_PI;
y = cos(0.5*y/ywidth) * sin(y)/y;
}
else
{
y = 1.0;
}
return x*y;
}
RtFilterFunc GetFilter( char *name )
{
if (!strcmp(name,"box"))
return RiBoxFilter;
if (!strcmp(name,"myBessel"))
return myBesselFilter;
if (!strcmp(name,"mybox"))
return MyBoxFilter;
if (!strcmp(name,"mytriangle"))
return MyTriangleFilter;
if (!strcmp(name,"catmull-rom"))
return RiCatmullRomFilter;
if (!strcmp(name,"mycatmull-rom"))
return MyCatmullRomFilter;
if (!strcmp(name,"gaussian"))
return RiGaussianFilter;
if (!strcmp(name,"mygaussian"))
return MyGaussianFilter;
if (!strcmp(name,"sinc"))
return RiSincFilter;
if (!strcmp(name,"mysinc"))
return MySincFilter;
if (!strcmp(name,"vonHannSinc"))
return vonHannSincFilter;
if (!strcmp(name,"halfCosineSinc"))
return halfCosineSincFilter;
if (!strcmp(name,"mydisk"))
return myDiskFilter;
#ifdef PRMAN
if (!strcmp(name,"bessel"))
return RiBesselFilter;
if (!strcmp(name,"disk"))
return RiDiskFilter;
#endif
return NULL;
}
#ifdef TWO_DIM
int main(int argc, char **argv)
{
float d;
float x,y,r;
float angle = 0;
float rw=WIDTH/2.0+1;
RtFilterFunc f;
if (argc!=2)
return 1;
f = GetFilter( argv[1] );
angle = M_PI * angle / 180;
for ( r = 0.0; r <= rw; r+=0.01 )
{
x = r * cos(angle);
y = r * sin(angle);
d = f( x, y, WIDTH, WIDTH );
printf( "v %g %g %g\n", r, 0.0, d );
}
printf("cstype bspline\n");
printf("deg 1\n");
printf("g Obj_curve#1\n");
printf("curv ");
for ( x = 0.,d=0.; x <= rw; x+=0.01,d+=1.0 )
printf( "%g %g ",d,d);
printf("\n");
printf("end\n");
printf("g\n");
return 0;
}
#else
/* 3D */
int main(int argc, char **argv)
{
int i;
float x,y,w,xw,yw;
RtFilterFunc f;
if (argc!=2)
return 1;
f = GetFilter( argv[1] );
w = 0.1;
xw = XWIDTH/2.0 + 1;
yw = YWIDTH/2.0 + 1;
for ( x = -xw; x < xw; x+=w )
for ( y = -yw; y < yw; y+=w )
{
printf( "v %g %g %g\n", x, y, f(x,y,XWIDTH,YWIDTH ));
printf( "v %g %g %g\n", x+w, y, f(x+w,y,XWIDTH,YWIDTH ));
printf( "v %g %g %g\n", x+w, y+w, f(x+w,y+w,XWIDTH,YWIDTH ));
printf( "v %g %g %g\n", x, y+w, f(x,y+w,XWIDTH,YWIDTH ));
}
i = 1;
for ( x = -xw; x < xw; x+=w )
for ( y = -yw; y < yw; y+=w )
{
printf( "f %d %d %d %d\n",i,i+1,i+2,i+3);
i += 4;
}
return 0;
}
#endif
#include <math.h>
#include "ri.h"
#define WIDTH 4
RtFloat tebBesselFilter( RtFloat x, RtFloat y, RtFloat xwidth, RtFloat ywidth )
{
double d,w,xx,yy;
xx = x*x;
yy = y*y;
xwidth *= 0.5;
ywidth *= 0.5;
w = (xx)/(xwidth*xwidth)+(yy)/(ywidth*ywidth);
if ( w<1.0 )
{
d = sqrt(xx+yy);
if (d!=0.0)
{
/* Half cosine window. */
w = cos(0.5*M_PI*sqrt(w));
return w * 2*j1(M_PI*d)/d;
}
else
{
return M_PI;
}
}
else
{
return 0.0;
}
}
/* Modified version of the RI Spec 3.2 sinc filter but also with a
* half cosine window applied. (half cosine Window)
*/
RtFloat halfCosineSincFilter( RtFloat x, RtFloat y,
RtFloat xwidth, RtFloat ywidth )
{
/* Uses a -PI to PI cosine window. */
if (x!=0.0)
{
x *= M_PI;
x = cos(0.5*x) * sin(x)/x;
}
else
{
x = 1.0;
}
if (y!=0.0)
{
y *= M_PI;
y = cos(0.5*y) * sin(y)/y;
}
else
{
y = 1.0;
}
return x*y;
}
static void dotest1(char *filtername, RtFilterFunc filter);
static void dotest2(char *filtername, RtFilterFunc filter);
int main(int argc, char **argv)
{
RiBegin(RI_NULL);
dotest1("sincT1.tif", RiSincFilter);
dotest1("tebT1.tif", tebBesselFilter);
dotest1("boxT1.tif", RiBoxFilter);
dotest2("sincT2.tif", RiSincFilter);
dotest2("tebT2.tif", tebBesselFilter);
dotest2("boxT2.tif", RiBoxFilter);
RiEnd();
return 0;
}
static void dotest1(char *filtername, RtFilterFunc filter)
{
RtFloat k = 1.0;
int i;
RiFrameBegin(1);
RiDisplay(filtername, RI_FILE, RI_RGB, RI_NULL);
RiFormat(512, 512, 1);
RiClipping(1., 100.);
RiPixelFilter(filter, WIDTH, WIDTH);
RiPixelSamples(4, 4);
RiProjection("orthographic", RI_NULL);
RiTranslate( 0.,0.,1. );
RiWorldBegin();
RiLightSource("ambientlight", "intensity", &k, RI_NULL);
RiSurface("plastic", RI_NULL);
for (i=0.0;i<120;i++)
{
RiTransformBegin();
RiRotate( 3*(RtFloat)i, 0., 0., 1. );
RiTorus( 2.0, 1.2, 0., 180., 1.5, RI_NULL );
RiTorus( 0.6, 0.18, 0., 180., 1.5, RI_NULL );
RiTorus( 0.3, 0.10, 0., 180., 1.5, RI_NULL );
RiTorus( 0.12, 0.07, 0., 180., 1.5, RI_NULL );
RiTransformEnd();
}
RiWorldEnd();
RiFrameEnd();
}
static void dotest2(char *filtername, RtFilterFunc filter)
{
RtFloat k = 1.0;
RtFloat scale = 1./360.;
RtFloat f;
RiFrameBegin(1);
RiDisplay(filtername, RI_FILE, RI_RGB, RI_NULL);
RiFormat(512, 512, 1);
RiClipping(1., 100.);
RiPixelFilter(filter, WIDTH, WIDTH);
RiPixelSamples(4, 4);
RiProjection("orthographic", RI_NULL);
RiTranslate( 0.,0.,1. );
RiWorldBegin();
RiLightSource("ambientlight", "intensity", &k, RI_NULL);
RiSurface("plastic", RI_NULL);
for (f=scale*2;f<1.5;f+=scale*4)
{
RiTorus( f, scale, 0., 180., 360., RI_NULL );
}
RiWorldEnd();
RiFrameEnd();
}