# The Joy of Mathematics

## Code: The Drawing Program

### Algorithm

The algorithm used here is similar in concept to POVRAY (Persistence Of Vision RAY tracing) in that it builds the image one point at a time in a straight raster scan. In contrast with vector plotting algorithms, it can not plot parametric systems of equations. That is, z must be a function of the x, y plane. Each pixel is visited only once and there is no backtracking.

## Source File: "arose.c"

/***********************************************************************
Plot Pascal's Rose

This program plots an arbitrary function of X and Y, producing a 640 x
480 x 16 color TIFF file. X and Y are assumed to be independent
variables and Z (the color, depth) is the dependent variable.

Be careful when you write a new function to plot to avoid division by
zero, particularly at the origin, and be careful to avoid floating
point overflow and underflow as well.

Author: Patrick J. Gleason
Company: IDK Computer Systems, Inc.
Contact: (315) 475-5598
Date: 13-Oct-2006
Version/Revision: 1.0
Compiler: tested with IBM VisualAge C++ 3.0, Microsoft C 6.0

***********************************************************************/
#include <stdlib.h> // Define the C Language standard run time library.
#include <string.h> // Define C standard string functions.
#include <math.h> // Define floating point arithmetic.
#include <stdio.h> // Define C Language standard I/O.
#include <errno.h> // Define C standard run time errors.
#include <conio.h> // Define IBM PC console I/O.
#include "tiff.h" // Define TIFF file structures.
/*
Macros
*/
#define movmem(f, t, c) memcpy(t, f, (unsigned) c)
#define setmem(src, len, chr) memset(src, (int) chr, (unsigned) len)
#define clrmem(dst) memset(dst, 0, sizeof(dst))
/*
Constants
*/
#define WIDTHH 640
#define HEIGHT 480
#define DEPTH 16
#define TRUE 1
#define FALSE 0
/*
Global Data

Define a resistor color code with 16 colors. I am using 16 rather
than 10 to fit the basic EGA/VGA palette size. All of these are
"web safe" colors.
*/
PIXEL24 palette[DEPTH]={
{0x00, 0x00, 0x00}, // 0 - black
{0x99, 0x00, 0x00}, // 1 - brown
{0xCC, 0x33, 0x00}, // 2 - red-brown
{0xFF, 0x00, 0x00}, // 3 - red
{0xFF, 0x99, 0x00}, // 4 - orange
{0xFF, 0xCC, 0x00}, // 5 - yellow-orange
{0xFF, 0xFF, 0x00}, // 6 - yellow
{0x99, 0xFF, 0x00}, // 7 - yellow-green
{0x00, 0xCC, 0x00}, // 8 - green
{0x00, 0xCC, 0xCC}, // 9 - cyan
{0x00, 0x00, 0xFF}, // A - blue
{0x66, 0x00, 0x99}, // B - indigo
{0x99, 0x33, 0xFF}, // C - violet
{0xCC, 0xCC, 0xFF}, // D - blue-grey
{0xCC, 0xCC, 0xCC}, // E - grey
{0xFF, 0xFF, 0xFF}, // F - white
};
/*
Function Prototypes
*/
int main(int argc, char *argv[]);
int plot(float minx, float miny, float maxx, float maxy, float (*f)(float, float), char *fspec);
float f01(float x, float y);
float f02(float x, float y);
float f03(float x, float y);
float f04(float x, float y);
float f05(float x, float y);
float f06(float x, float y);
float f07(float x, float y);
float f08(float x, float y);
float f09(float x, float y);
float f10(float x, float y);
float f11(float x, float y);
float f12(float x, float y);
float f13(float x, float y);
float f14(float x, float y);
float f15(float x, float y);
float f16(float x, float y);
float f17(float x, float y);
float f18(float x, float y);
float f19(float x, float y);
float f20(float x, float y);
float f21(float x, float y);
float f22(float x, float y);

/***********************************************************************
Program Entry Point

This is a script that calls the plot() function to produce 21
interesting curves.

Input Arguments:

argc is the number of space separated strings in the
command line that invoked this program.

argv is an array of memory addresses of null-
terminated ASCII strings parsed from the command
line by the operating system shell. There
are no parameters for this program.

Output Arguments:

none

Constants Used:

none

Global Data Used:

none

Object Data Used:

n/a

Returns:

0 success (always)

***********************************************************************/
int main(int argc, char *argv[]){
plot(-2.0, -2.0, 2.0, 2.0, f01, "arose01.tif");
plot(-2.0, -2.0, 2.0, 2.0, f02, "arose02.tif");
plot(-2.0, -2.0, 2.0, 2.0, f03, "arose03.tif");
plot(-2.0, -2.0, 2.0, 2.0, f04, "arose04.tif");
plot(-2.0, -2.0, 2.0, 2.0, f05, "arose05.tif");
plot(-10.0, -10.0, 10.0, 10.0, f06, "arose06.tif");
plot(-10.0, -10.0, 10.0, 10.0, f07, "arose07.tif");
plot(-10.0, -10.0, 10.0, 10.0, f08, "arose08.tif");
plot(-2.0, -2.0, 2.0, 2.0, f09, "arose09.tif");
plot(-2.0, -2.0, 2.0, 2.0, f10, "arose10.tif");
plot(-2.0, -2.0, 2.0, 2.0, f11, "arose11.tif");
plot(-2.0, -2.0, 2.0, 2.0, f12, "arose12.tif");
plot(-2.0, -2.0, 2.0, 2.0, f13, "arose13.tif");
plot(-2.0, -2.0, 2.0, 2.0, f14, "arose14.tif");
plot(-1.5, -1.5, 1.5, 1.5, f15, "arose15.tif");
plot(-1.1, -1.1, 1.1, 1.1, f16, "arose16.tif");
plot(-3.0, -3.0, 3.0, 3.0, f17, "arose17.tif");
plot(-2.0, -2.0, 2.0, 2.0, f18, "arose18.tif");
plot(-20.0, -20.0, 20.0, 20.0, f19, "arose19.tif");
plot(-20.0, -20.0, 20.0, 20.0, f20, "arose20.tif");
plot(-20.0, -20.0, 20.0, 20.0, f21, "arose21.tif");
} //

/***********************************************************************
Plot a Function

This function produces a TIFF image file containing a color plot of some
arbitrary function of X and Y. 24-bit color is produced, though only 16
colors are used. The colors are approximately the traditional resistor
color code scale, with black representing 0 and white representing 15.
Obviously, a few intermediate colors were added to extend the scale from
ten values to 16. The 640 x 480 x 16 color resolution was selected to
be compatible with the minimal specification of a VGA display adapter.

Simple black lines are added to mark the coordinate axes.

Zero crossings (roots) are marked by a white line.

The layout of the TIFF file is:

Offset Length Contents
--------------------------------------------------------
8 186 IFD (Image File Directory)
194 6 Bits per Sample (three 16-bit numbers)
200 8 X Resolution (two 32-bit numbers)
208 8 Y Resolution (two 32-bit numbers)
216 8 Software Name
224 921600 image data

Input Arguments:

xmin specifies the left boundary of the real number
space to be plotted.

ymin specifies the bottom bottom boundary of the
real number space to be plotted.

xmax specifies the right boundary of the real number
space to be plotted.

ymax specifies the top boundary of the real number
space to be plotted.

f specifies the entry point of an arbitrary
function that computes Z = f(X, Y); hopefully
without generating a floating point
exception.

fspec is the memory address of a null-terminated
character string that specifies the output
file.

Output Arguments:

none

Constants Used:

WIDTHH is the horizontal display resolution; typically
640 pixels.

HEIGHT is the vertical display resolution; typically
480 pixels.

Global Data Used:

palette (read only) defines a palette of 16 "web safe"
colors, similar to the traditional resistor
color code.

Object Data Used:

n/a

Returns:

TRUE if successful

FALSE if not

***********************************************************************/
int plot(float minx, float miny, float maxx, float maxy, float (*f)(float, float), char *fspec){
int i; // current vertical screen coordinate
int j; // current horizontal screen coordinate
int k; // current screen depth
size_t n0; // number of bytes for fwrite()
size_t n1; // result of fwrite()
float aspect; // screen aspect ratio (height / width)
float deltax; // width of one screen pixel in real space
float deltay; // height of one screen pixel in real space
float zscale; // depth of one screen pixel in real space
float minz; // front edge of plot space
float maxz; // back edge of plot space
float x; // current real horizontal coordinate
float y; // current real vertical coordinate
float z; // current real depth coordinate
float z_up; // value of the point above the current point
float z_down; // value of the point below the current point
float z_left; // value of the point to the left of the current point
float z_right; // value of the point to the right of the current point
FILE *outfile; // handle for output file
IFD ifd; // TIFF image file directory
PIXEL24 sample; // buffer for writing current pixel
unsigned short int BitsPerSample[3]; // buffer for writing TIFF bits per sample information
unsigned long int XResolution[2]; // numerator and denominator of horizontal pixels per inch for TIFF
unsigned long int YResolution[2]; // numerator and denominator of vertical pixels per inch for TIFF
char SWName[32]; // buffer for writing image generation program name for TIFF
//
printf("Plotting %s", fspec); // Show the operator what's happening; keep the line open for the Z limits.
/*
Create the TIFF file.
*/
outfile=fopen(fspec, "wb"); // Open the output file in binary mode; overwrite any prior file having the same name.
if(!outfile){ // If the output file could not be opened, then
printf("*** ERROR %d creating \"%s\".\n", errno, fspec);
return FALSE; // display an error message and return failure.
} //
/*
Build and write the TIFF file header.
*/
movmem("II", tiff_header.signature, 2); // Use Intel (big-endian) format for binary integers.
tiff_header.version=42; // Use TIFF revision 4.2 (compatible back to 1987).
n0=sizeof(tiff_header); // Set the number of bytes to write.
n1=fwrite(&tiff_header, 1, n0, outfile); // Write an arbitrary binary structure.
if(n1!=n0){ // If fewer bytes than specified were written, then
printf("*** ERROR %d writing TIFF header; %d of %d bytes written.\n",
errno, n1, n0); // display an error message,
fclose(outfile); // close the output file,
return FALSE; // and return failure.
} //
/*
Build the first, and only, IFD, and write it.
*/
ifd.count=15; // Put 15 tags in this IFD.
//
ifd.entry[0].tag=254; // NewSubfileType
ifd.entry[0].type=4; // long integer
ifd.entry[0].count=1; // one item
ifd.entry[0].value=0; // no special flags set
//
ifd.entry[1].tag=256; // ImageWidth
ifd.entry[1].type=3; // long integer
ifd.entry[1].count=1; // one item
ifd.entry[1].value=WIDTHH; // Specify the image width.
//
ifd.entry[2].tag=257; // ImageLength
ifd.entry[2].type=3; // long integer
ifd.entry[2].count=1; // one item
ifd.entry[2].value=HEIGHT; // Specify the image height.
//
ifd.entry[3].tag=258; // BitsPerSample
ifd.entry[3].type=3; // short integer
ifd.entry[3].count=3; // three items
BitsPerSample[0]=8; // red
BitsPerSample[1]=8; // green
BitsPerSample[2]=8; // blue
//
ifd.entry[4].tag=259; // Compression
ifd.entry[4].type=3; // short integer
ifd.entry[4].count=1; // one item
ifd.entry[4].value=1; // no compression
//
ifd.entry[5].tag=262; // PhotometricInterpretation
ifd.entry[5].type=3; // short integer
ifd.entry[5].count=1; // one item
ifd.entry[5].value=2; // RGB
//
ifd.entry[6].tag=273; // StripOffsets
ifd.entry[6].type=4; // long integer
ifd.entry[6].count=1; // one item
//
ifd.entry[7].tag=277; // SamplesPerPixel
ifd.entry[7].type=3; // short integer
ifd.entry[7].count=1; // one item
ifd.entry[7].value=3; // three samples per pixel
//
ifd.entry[8].tag=278; // RowsPerStrip
ifd.entry[8].type=4; // long integer
ifd.entry[8].count=1; // one item
ifd.entry[8].value=HEIGHT; // The entire image is one strip.
//
ifd.entry[9].tag=279; // StripByteCounts
ifd.entry[9].type=4; // long integer
ifd.entry[9].count=1; // one item
ifd.entry[9].value=WIDTHH*HEIGHT*3; // should be 921600
//
ifd.entry[10].tag=282; // XResolution
ifd.entry[10].type=5; // rational
ifd.entry[10].count=1; // one item
ifd.entry[10].value=ifd.entry[3].value+sizeof(BitsPerSample); // should be 200
XResolution[0]=72; // numerator for 72 dpi
XResolution[1]=1; // denominator
//
ifd.entry[11].tag=283; // YResolution
ifd.entry[11].type=5; // rational
ifd.entry[11].count=1; // one item
ifd.entry[11].value=ifd.entry[10].value+sizeof(XResolution); // should be 208
YResolution[0]=72; // numerator for 72 dpi
YResolution[1]=1; // denominator
//
ifd.entry[12].tag=284; // PlanarConfiguration
ifd.entry[12].type=3; // short integer
ifd.entry[12].count=1; // one item
ifd.entry[12].value=1; // chunky
//
ifd.entry[13].tag=296; // ResolutionUnit
ifd.entry[13].type=3; // short integer
ifd.entry[13].count=1; // one item
ifd.entry[13].value=2; // inches
//
ifd.entry[14].tag=305; // Software
ifd.entry[14].type=2; // ASCII
clrmem(SWName); // Make sure that the string is null-terminated.
strncpy(SWName, "AROSE.C", sizeof(SWName)-1); // program name
ifd.entry[14].count=strlen(SWName)+1; // character count, including null
ifd.entry[14].value=ifd.entry[11].value+sizeof(YResolution); // should be 216
ifd.entry[6].value=ifd.entry[14].value+ifd.entry[14].count; // should be 224
//
ifd.nextdir=0; // Indicate no other images present.
//
n0=sizeof(ifd); // Set the number of bytes to write.
n1=fwrite(&ifd, 1, n0, outfile); // Write an arbitrary binary structure.
if(n1!=n0){ // If fewer bytes than specified were written, then
printf("*** ERROR %d writing IFD; %d of %d bytes written.\n",
errno, n1, n0); // display an error message,
fclose(outfile); // close the output file,
return FALSE; // and return failure.
} //
/*
Write the bits per sample array.
*/
n0=sizeof(BitsPerSample); // Set the number of bytes to write.
n1=fwrite(BitsPerSample, 1, n0, outfile); // Write an arbitrary binary structure.
if(n1!=n0){ // If fewer bytes than specified were written, then
printf("*** ERROR %d writing bits per sample array; %d of %d bytes written.\n",
errno, n1, n0); // display an error message,
fclose(outfile); // close the output file,
return FALSE; // and return failure.
} //
/*
Write the X resolution.
*/
n0=sizeof(XResolution); // Set the number of bytes to write.
n1=fwrite(XResolution, 1, n0, outfile); // Write an arbitrary binary structure.
if(n1!=n0){ // If fewer bytes than specified were written, then
printf("*** ERROR %d writing X resolution; %d of %d bytes written.\n",
errno, n1, n0); // display an error message,
fclose(outfile); // close the output file,
return FALSE; // and return failure.
} //
/*
Write the Y resolution.
*/
n0=sizeof(YResolution); // Set the number of bytes to write.
n1=fwrite(YResolution, 1, n0, outfile); // Write an arbitrary binary structure.
if(n1!=n0){ // If fewer bytes than specified were written, then
printf("*** ERROR %d writing Y resolution; %d of %d bytes written.\n",
errno, n1, n0); // display an error message,
fclose(outfile); // close the output file,
return FALSE; // and return failure.
} //
/*
Write the software name.
*/
n0=ifd.entry[14].count; // Set the number of bytes to write.
n1=fwrite(SWName, 1, n0, outfile); // Write an arbitrary binary structure.
if(n1!=n0){ // If fewer bytes than specified were written, then
printf("*** ERROR %d writing software name; %d of %d bytes written.\n",
errno, n1, n0); // display an error message,
fclose(outfile); // close the output file,
return FALSE; // and return failure.
} //
/*
Begin the plot.
*/
aspect=((float)HEIGHT)/((float)WIDTHH); // Compute the screen aspect ratio.
miny=miny*aspect; // Adjust the bottom edge of the plot to the bottom edge of the screen.
maxy=maxy*aspect; // Adjust the top edge of the plot to the top edge of the screen.
//
deltax=(maxx-minx)/640.0; // Compute the width of one screen pixel in real space.
deltay=(maxy-miny)/480.0; // Compute the height of one screen pixel in real space.
/*
Find the maximum and minimum plot depth.
*/
maxz=minz=f(minx, miny); // Initialize the maximum and minimum depth to the bottom left corner of the plot.
y=miny; // Begin at the bottom edge of real number space.
for(i=0; i<HEIGHT; i++){ // Do this for 480 scan lines.
x=minx; // Begin at the left edge of real number space.
for(j=0; j<WIDTHH; j++){ // Do this for 640 pixels per line.
z=f(x, y); // Compute the depth of the function here.
if(z<minz) minz=z; // If this depth is less than the minimum, then replace the minimum with it.
if(z>maxz) maxz=z; // If this depth is greater than the maximum, then replace the maximum with it.
x+=deltax; // Move one pixel to the right.
} //
y+=deltay; // Move one scan line up.
} //
zscale=((float)DEPTH)/(maxz-minz); // Compute the depth of one screen pixel.
printf(" Min z: %f, Max z: %f\n"); // Provide these numbers for the database.
/*
Plot the Function
*/
y=miny; // Begin at the bottom edge of the real number space.
for(i=0; i<HEIGHT; i++){ // Do this for 480 scan lines.
x=minx; // Begin at the left edge of the real number space.
for(j=0; jlt;WIDTHH; j++){ // Do this for 640 pixels per line.
z=f(x, y); // Compute the function depth at this point.
k=(z-minz)*zscale; // Map the depth from real space onto screen space.
if(k>(DEPTH-1)) k=DEPTH-1; // Avoid array bounds violations,
if(k<0) k=0; // should f(x, y) have bugs.
if((i==(HEIGHT/2)-1)||(i==(HEIGHT/2))){
k=0; // Draw the X axis.
} //
if((j==(WIDTHH/2)-1)||(j==(WIDTHH/2))){
k=0; // Draw the Y axis.
} //
z_up=f(x, y+deltay); // Detect zero crossings
z_down=f(x, y-deltay); // in the vertical direction.
if((z_up<0)!=(z_down<0)){ // If the point above and the point below have different signs, then
k=DEPTH-1; // this point is a zero crossing.
} //
z_left=f(x-deltax, y); // Detect zero crossings
z_right=f(x+deltax, y); // in the horizontal direction.
if((z_left<0)!=(z_right<0)){ // If the point left and the point right have different signs, then
k=DEPTH-1; // this point is a zero crossing.
} //
sample.r=palette[k].r; // Copy the color from the palette
sample.g=palette[k].g; // into the pixel buffer.
sample.b=palette[k].b; //
n0=sizeof(sample); // Set the number of bytes to write.
n1=fwrite(&sample, 1, n0, outfile); // Write an arbitrary binary structure.
if(n1!=n0){ // If fewer bytes than specified were written, then
printf("*** ERROR %d writing pixel at (%d, %d); %d of %d bytes written.\n",
errno, i, j, n1, n0); // display an error message,
fclose(outfile); // close the output file,
return FALSE; // and return failure.
} //
x+=deltax; // Move one pixel to the right.
} //
y+=deltay; // Move one scan line up.
} //
fclose(outfile); // Close the output file.
return TRUE; // Return success.
}

/***********************************************************************
Basic Cartoid

This function, and the rest that follow, are interchangeable. These
are the functions of X and Y to be plotted. All compute a Z value for
a given point in the X, Y plane. One function header will suffice for
all.

Input Arguments:

x current horizontal coordinate

y current vertical coordinate

Output Arguments:

none

Constants Used:

none

Global Data Used;

none

Object Data Used:

n/a

Returns:

a depth for the specified point in the X, Y plane

***********************************************************************/
float f01(float x, float y){
double z; // result
//
r = sqrt(x*x + y*y); // Compute r.
if(r==0.0) return 1.0; // Avoid division by zero.
z = x/r - r; // 0 = cos(theta) - r
return z; // Done.
} //

/***********************************************************************
Basic Limacon

***********************************************************************/
float f02(float x, float y){
double z; // result
//
r = sqrt(x*x + y*y); // Compute r.
if(r==0.0) return 2.0; // Avoid division by zero.
z = 1.0 + x/r - r; // 0 = 1 + cos(theta) - r
return z; // Done.
} //

/***********************************************************************
r^2 = cos^2(2*theta)

***********************************************************************/
float f03(float x, float y){
double x2; // x squared
double x4; // x^4
double y2; // y squared
double y4; // y^4
double r2; // r squared
double r4; // r^4
double z; // result
//
x2 = x*x; // Compute x squared.
x4 = x2*x2; // Compute x^4.
y2 = y*y; // Compute y squared.
y4 = y2*y2; // Compute y^4.
r2 = x2 + y2; // Compute r squared.
r4 = r2*r2; // Compute r^4.
if(r4==0.0) return 1.0; // Avoid division by zero.
z = ((x4 - 2.0*x2*y2 + y4) / r4) - r2; // 0 = cos^2(2*theta) - r^2
return z; // Done.
} //

/***********************************************************************
r = cos(3*theta)

***********************************************************************/
float f04(float x, float y){
double x2; // x squared
double x3; // x cubed
double r; //
double r2; // r squared
double r3; // r cubed
double z; // result
//
x2 = x*x; // Compute x squared.
x3 = x2*x; // Compute x cubed.
r2 = x2 + y*y; // Compute r squared.
r = sqrt(r2); // Compute r.
r3 = r2*r; // Compute r cubed.
if(r3==0.0) return 1.0; // Avoid division by zero.
z = (4.0*x3/r3 - 3.0*x/r) - r; // 0 = cos(3*theta)
return z; // Done.
}

/***********************************************************************
r^2 = cos^2(4*theta)

***********************************************************************/
float f05(float x, float y){
double x2; // x squared
double x4; // x^4
double x6; // x^6
double x8; // x^8
double y2; // y squared
double y4; // y^4
double y6; // y^6
double y8; // y^8
double r2; // r squared
double r4; // r^4
double r8; // r^8
double z; // result
//
x2 = x*x; // Compute x squared.
x4 = x2*x2; // Compute x^4.
x6 = x4*x2; // Compute x^6.
x8 = x4*x4; // Compute x^8.
y2 = y*y; // Compute y squared.
y4 = y2*y2; // Compute y^4.
y6 = y4*y2; // Compute y^6.
y8 = y4*y4; // Compute y^8.
r2 = x2 + y2; // Compute r squared.
r4 = r2*r2; // Compute r^4.
r8 = r4*r4; // Compute r^8.
if(r8==0.0) return 1.0; // Avoid division by zero.
z = ((x8 - 12.0*x6*y2 + 38.0*x4*y4 - 12.0*x2*y6 + y8) / r8) - r2;
return z; // Done.
} //

/***********************************************************************
z = cos(x) + cos(y);

***********************************************************************/
float f06(float x, float y){
double z; // result
//
z = cos(x) + cos(y); // checkerboard
return z; // Done.
} //

/***********************************************************************
z = sin(x)/x * sin(y)/y

***********************************************************************/
float f07(float x, float y){
double z; // result
//
if(x==0.0){ // Avoid division by zero.
if(y==0.0){
z = 1.0;
} else {
z = (sin(y)/y);
}
} else {
if(y==0.0){
z = (sin(x)/x);
} else {
z = (sin(x)/x) * (sin(y)/y); // 2D sampling function
}
}
return z; // Done.
} //

/***********************************************************************
z = cos(3*theta)

***********************************************************************/
float f08(float x, float y){
double x2; // x squared
double x3; // x cubed
double r; //
double r2; // r squared
double r3; // r cubed
double z; // result
//
x2 = x*x; // Compute x squared.
x3 = x2*x; // Compute x cubed.
r2 = x2 + y*y; // Compute r squared.
r = sqrt(r2); // Compute r.
r3 = r2*r; // Compute r cubed.
if(r3==0.0) return 1.0; // Avoid division by zero.
z = (4.0*x3/r3 - 3.0*x/r); // rays
return z; // Done.
}

/***********************************************************************
r^2 = cos^2(4*theta) with odd lobes supressed

***********************************************************************/
float f09(float x, float y){
double x2; // x squared
double x4; // x^4
double x8; // x^8
double y2; // y squared
double y4; // y^4
double y8; // y^8
double r2; // r squared
double r4; // r^4
double r8; // r^8
double z; // result
//
x2 = x*x; // Compute x squared.
x4 = x2*x2; // Compute x^4.
x8 = x4*x4; // Compute x^8.
y2 = y*y; // Compute y squared.
y4 = y2*y2; // Compute y^4.
y8 = y4*y4; // Compute y^8.
r2 = x2 + y2; // Compute r squared.
r4 = r2*r2; // Compute r^4.
r8 = r4*r4; // Compute r^8.
if(r8==0.0) return 1.0; // Avoid division by zero.
z = ((x8 + 38.0*x4*y4 + y8) / r8) - r2; // butterfly
return z; // Done.
} //

/***********************************************************************
r^2 = cos^2(3*theta)

***********************************************************************/
float f10(float x, float y){
double x2; // x squared
double x4; // x^4
double x6; // x^6
double r2; // r squared
double r4; // r^4
double r6; // r^6
double z; // result
//
x2 = x*x; // Compute x squared.
x4 = x2*x2; // Compute x^4.
x6 = x4*x2; // Compute x^6.
r2 = x2 + y*y; // Compute r squared.
r4 = r2*r2; // Compute r^4.
r6 = r4*r2; // Compute r^6.
if(r6==0.0) return 1.0; // Avoid division by zero.
z = 16*x6/r6 - 24*x4/r4 + 9*x2/r2 - r2; // 0 = cos^2(3*theta) - r^2
return z; // Done.
} //

/***********************************************************************
r = cos^2(theta) - sin^2(theta)

***********************************************************************/
float f11(float x, float y){
double x2; // x squared
double y2; // y squared
double r2; // r squared
double r; //
double z; // result
//
x2 = x*x; // Compute x squared.
y2 = y*y; // Compute y squared.
r2 = x2 + y2; // Compute r squared.
r = sqrt(r2); // Compute r.
if(r2==0.0) return 1.0; // Avoid division by zero.
z = ((x2 - y2) / r2) - r; // 0 = cos^2(theta) - sin^2(theta) - r
return z; // Done.
} //

/***********************************************************************
r = theta

***********************************************************************/
float f12(float x, float y){
double z; // result
//
r = sqrt(x*x + y*y); // Compute r.
if(x==0.0) return 0.0; // Avoid division by zero.
z = atan2(y, x) - r; // spiral
return z; // Done.
} //

/***********************************************************************
r = sin(5*theta)

***********************************************************************/
float f13(float x, float y){
double x2; // x squared
double x4; // x cubed
double y2; // y squared
double y3; // y cubed
double y5; // y^5
double r; //
double r2; // r squared
double r3; // r cubed
double r5; // r^5
double z; // result
//
x2 = x*x; // Compute x squared.
x4 = x2*x2; // Compute x^4.
y2 = y*y; // Compute y squared.
y3 = y2*y; // Compute y cubed.
y5 = y3*y2; // Compute y^5.
r2 = x2 + y2; // Compute r squared.
r = sqrt(r2); // Compute r.
r3 = r2*r; // Compute r cubed.
r5 = r3*r2; // Compute r^5.
if(r5==0.0) return 1.0; // Avoid division by zero.
z = ((4.0*y5 - 4.0*y3*x2 + 8.0*y*x4) / r5) - (3.0*y3 - 3.0*y*x2)/r3 - r;
return z; // Done.
} //

/***********************************************************************
Polynomial with Five Roots

The roots are at -2, -1, 0, +1 and +2.

***********************************************************************/
float f14(float x, float y){
double x2; // x squared
double x3; // x cubed
double x5; // x^5
double y2; // y squared
double y3; // y cubed
double y5; // y^5
double z; // result
//
x2 = x*x; // Compute x squared.
x3 = x2*x; // Compute x cubed.
x5 = x3*x2; // Compute x^5.
y2 = y*y; // Compute x squared.
y3 = y2*y; // Compute x cubed.
y5 = y3*y2; // Compute x^5.
z = x5 - 5.0*x3 + 4.0*x - y5 + 5.0*y3 - 4.0*y;
return z; // Done.
} //

/***********************************************************************
Two Hyperbolae and a Circle

***********************************************************************/
float f15(float x, float y){
double x2; // x squared
double y2; // y squared
double z; // result
//
x2 = x*x; // Compute x squared.
y2 = y*y; // Compute y squared.
z = (x2 - y2 - 1.0) * (y2 - x2 - 1.0) * (x2 + y2 - 1.0);
return z; // Done.
} //

/***********************************************************************
Radius is Two Hyperbolae and a Circle

***********************************************************************/
float f16(float x, float y){
double x2; // x squared
double y2; // y squared
double r2; // r squared
double z; // result
//
x2 = x*x; // Compute x squared.
y2 = y*y; // Compute y squared.
r2 = x2 + y2; // Compute r squared.
z = (x2 - y2 - 1.0) * (y2 - x2 - 1.0) * (x2 + y2 - 1.0) - r2;
return z; // Done.
} //

/***********************************************************************
Intersecting Waves

***********************************************************************/
float f17(float x, float y){
double z; // result
//
z = (sin(3.0*x) - y) * (sin(5.0*y) - x);
return z; // Done.
} //

/***********************************************************************
Ellipse

***********************************************************************/
float f18(float x, float y){
double z; // result
//
z = x*x - x*y + y*y; // ellipse
return z; // Done.
} //

/***********************************************************************
Sampling Function

***********************************************************************/
float f19(float x, float y){
double z; // result
//
if(x==0.0){
z = 10.0;
} else {
z = 10 - (10*sin(x) / x) - y;
}
return z; // Done.
} //

/***********************************************************************
Cosine Wave

***********************************************************************/
float f20(float x, float y){
double z; // result
//
z = 10*cos(x); // cosine wave
return z; // Done.
} //

/***********************************************************************
Zig Zag

***********************************************************************/
float f21(float x, float y){
double z; // result
//
z = 5.0*cos(x) + 2.0*x - y; // y = 5*cos(x) + 2*x
return z; // Done.
} //

## Source File: "tiff.h"

/**********************************************************************
TIFF (Tagged Image File Format) Structures

Author: Patrick J. Gleason
Date: 2-Oct-2006

This is based on revision 4.0 of the TIFF standard, published April 31,
1987.
**********************************************************************/
typedef struct { // Offs Len Contents
unsigned short int tag; // 0 2 TIFF tag number
unsigned short int type; // 2 2 data type
unsigned long int count; // 4 4 number of data items
unsigned long int value; // 8 4 absolute offset of field data, or value depending on tag
} IFD_ENTRY; // Total Size = 12 bytes

typedef struct { // Offs Len Contents
unsigned short int count; // 0 2 number of "IFD_ENTRY"s in this directory
IFD_ENTRY entry[15]; // 2 180 the number of these varies, according to "count", just above
unsigned long nextdir; // 182 4 absolute offset of next directory, or zero if none
} IFD; // Total Size = 186 bytes here, but varies according to "count"

typedef struct { // Offs Len Contents
char signature[2]; // 0 2 "II" for Intel CPUs (big endian), "MM" for Motorola CPUs (little endian)
unsigned short int version; // 2 2 TIFF version number
unsigned long int ifd_offset; // 4 4 offset, relative to beginning of file, of first IFD
} TIFF_HEADER; // Total Size = 8 bytes

typedef struct { // Offs Len Contents
unsigned char intens; // 0 1 intensity (0 - 255)
} PIXEL8; // Total Size = 1 byte

typedef struct { // Offs Len Contents
unsigned char r; // 0 1 red intensity (0 - 255)
unsigned char g; // 1 1 green intensity (0 - 255)
unsigned char b; // 2 1 blue intensity (0 - 255)
} PIXEL24; // Total Size = 3 bytes

typedef struct { // Offs Len Contents
unsigned short int r; // 0 2 red intensity (0 - 255)
unsigned short int g; // 2 2 green intensity (0 - 255)
unsigned short int b; // 4 2 blue intensity (0 - 255)
} PIXEL48; // Total Size = 6 bytes
/*
Names of Data Types for "type" Field of "IFD_ENTRY"
*/
char *type_name[5]={"byte", "text", "short", "long", "rational"};
/*
TIFF Tag Names
*/
typedef struct {
int num; // tag number
char *name; // tag name
} TAG_NAME;

TAG_NAME tag_name[]={
{254, "NewSubfileType"},
{255, "SubfileType"},
{256, "ImageWidth"},
{257, "ImageLength"},
{258, "BitsPerSample"},
{259, "Compression"},

{262, "PhotometricInterpretation"},
{263, "Thresholding"},
{264, "CellWidth"},
{265, "CellLength"},
{266, "FillOrder"},

{269, "DocumentName"},
{270, "ImageDescription"},
{271, "Make"},
{272, "Model"},
{273, "StripOffsets"},
{274, "Orientation"},

{277, "SamplesPerPixel"},
{278, "RowsPerStrip"},
{279, "StripByteCounts"},
{280, "MinSampleValue"},
{281, "MaxSampleValue"},
{282, "XResolution"},
{283, "YResolution"},
{284, "PlanarConfiguration"},
{285, "PageName"},
{286, "XPosition"},
{287, "YPosition"},
{288, "FreeOffsets"},
{289, "FreeByteCounts"},
{290, "GrayResponseUnit"},
{291, "GrayResponseCurve"},
{292, "Group3Options"},
{293, "Group4Options"},

{296, "ResolutionUnit"},
{297, "PageNumber"},

{300, "ColorResponseUnit"},
{301, "ColorResponseCurves"},

{305, "Software"},
{306, "DateTime"},

{315, "Artist"},
{316, "HostComputer"},

{0, ""}
};
/*
Compression Type Names
*/
char *compression_type[4]={
"none",
"1-dimensional modified Huffman RLE",
"group 3 fax",
"group 4 fax",
};
/*
Photometric Interpretation Names
*/
char *photometric[3]={
"grey scale negative",
"grey scale",
"RGB"
};
/*
Units Names
*/
char *units[3]={
"none",
"inches",
"centimeters"
};