#include <stdio.h>
#include <math.h>
#include "multicolor.h"

/*
 * multicolor.c
 *
 * Code for "multicolor" test program
 *
 * Copyright 6/12/98 F. Kenton Musgrave
 * All rights reserved
 */


Matrix Identity = {
        1., 0., 0., 0.,
        0., 1., 0., 0.,
        0., 0., 1., 0.,
        0., 0., 0., 1.
};


main()
{
	Vector	v;
	Color	c;
	char	outline[LENGTH];
	double	x, y;
	int 	i, entry;
	FILE	*fp, *fopen();

	if((fp = fopen("outfile","w")) == NULL) {
		fprintf(stderr, "can't open image file\n");
		exit(-1);
	}

	/* initialize noise function */
	init();

	/* output ppm file header */
	fprintf( fp, "P6 %d %d %d\n", RES, RES, 255 );

	for ( y = 6.0; y < 7.0; y += 1.0/RES ) {
		entry = 0;
		for ( i = 0, x = 9.0; x < 10.0; x += 1.0/RES, i++ ) {
			v.x = x;
			v.y = y;
			v.z = 5.0;
			c.r = c.g = c.b = 0.5;
			Multicolor( v, &c, 0.7, 2.0, 8.0, 0.2, 0.5, 8.0, 
					  4.0, 5.0, 1.0, 1.0, 2e4 );
			outline[entry++] = (char)(255.0 * c.r);
			outline[entry++] = (char)(255.0 * c.g);
			outline[entry++] = (char)(255.0 * c.b);
		}
		
		fwrite( (void *)outline, 
			(size_t)sizeof(char), (size_t)(LENGTH), fp );

	}

	fclose( fp );

} /* main()


/* 
 * A multifractal multicolor texture experiment
 */
void Multicolor(Vector texture,
		Color *color,
		double arg0, double arg1, double arg2, double arg3, 
		double arg4, double arg5, double arg6, double arg7, 
		double arg8, double arg9, double arg10 )
{
	Vector	axis, purt, cvec;
	Matrix  matrix, tmpmat;

		/* initialize xformtation matrix to identity */
	CopyMatrix(Identity, matrix);

		/* get random rotation axis */
	axis = Wrinkled(texture, 2.0, arg4, arg5 );

		/* construct random rotation matrix */
	RotateX(arg6*axis.x, tmpmat);
	MultMatrix(matrix, tmpmat, matrix);
	RotateY(arg6*axis.y, tmpmat);
	MultMatrix(matrix, tmpmat, matrix);
	RotateZ(arg6*axis.z, tmpmat);
	MultMatrix(matrix, tmpmat, matrix);

		/* perform random fBm domain purturbation */
	purt = texture;
	SMULT( 0.3, purt );
	purt  = Wrinkled( purt, 2.0, 0.5, 7.0 );
	SMULT( arg7, purt );
	VADD( texture, purt );

		/* get a vector-valued multifractal for saturation */
	cvec.x = arg10 * Multifractal( texture, arg0, arg1, arg2, arg3 );
	texture.x += 10.5;
	cvec.y = arg10 * Multifractal( texture, arg0, arg1, arg2, arg3 );
	texture.y += 10.5;
	cvec.z = arg10 * Multifractal( texture, arg0, arg1, arg2, arg3 );

		/* transform color vector by random rotation matrix */
	cvec = VectTransform(cvec, matrix);

		/* determine final color */
	color->r += arg8 * cvec.x;
	color->g += arg8 * cvec.y;
	color->b += arg8 * cvec.z;

		/* clamp color values */
        if ( color->r < 0.0 )
                color->r = 0.0;
        else if ( color->r > 1.0 )
                color->r = 1.0;
        if ( color->g < 0.0 )
                color->g = 0.0;
        else if ( color->g > 1.0 )
                color->g = 1.0;
        if ( color->b < 0.0 )
                color->b = 0.0;
        else if ( color->b > 1.0 )
                color->b = 1.0;

} /* Multicolor() */


/*
 * Ken Perlin's original version of vector-valued fBm
 */
Vector Wrinkled( Vector point, double lacunarity, double H, double octaves )
{
        Vector	sPoint, result, temp;
        double	f, s;
        int	i;

        result.x = result.y = result.z = 0.0;
        f = s = 1.0;
        for ( i = 0; i < octaves; i++ )
                {
                sPoint.x = f * point.x;
                sPoint.y = f * point.y;
                sPoint.z = f * point.z;
                temp = VecNoise3(sPoint);
                result.x += temp.x * s;
                result.y += temp.y * s;
                result.z += temp.z * s;
                s *= H;
                f *= lacunarity;
                }
        return(result);
} /* Wrinkled() */


void CopyMatrix( Matrix a, Matrix b )
{
        int	i, j;

        for( i = 0; i < 4; i++ )
                for( j = 0; j < 4; j++ )
                        b[i][j] = a[i][j];
} /* CopyMatrix() */


void MultMatrix( Matrix a, Matrix b, Matrix c )
{
        double	sum;
        int	i, j, k;

        for( i = 0; i < 4; i++ ) {
                for( j = 0; j < 4; j++ ) {
                        sum = 0.0;
                	for( k = 0; k < 4; k++ )
                                sum += a[i][k] * b[k][j];
                        c[i][j]= sum;
                }
	}
} /* MultMatrix() */


void RotateX( double theta, Matrix m ) 
{
        double	sine, cosine;

        sine = sin(theta);
        cosine = cos(theta);

        CopyMatrix(Identity, m);
        m[1][1] = cosine;
        m[2][1] = -sine;
        m[1][2] = sine;
        m[2][2] = cosine;
} /* RotateX() */


void RotateY( double theta, Matrix m ) 
{
        double sine, cosine;

        sine = sin(theta);
        cosine = cos(theta);

        CopyMatrix(Identity, m);
        m[0][0] = cosine;
        m[2][0] = sine;
        m[0][2] = -sine;
        m[2][2] = cosine;
} /* RotateY() */


void RotateZ( double theta, Matrix m ) 
{
        double	sine, cosine;

        sine = sin(theta);
        cosine = cos(theta);

        CopyMatrix(Identity, m);
        m[0][0] = cosine;
        m[1][0] = -sine;
        m[0][1] = sine;
        m[1][1] = cosine;
} /* RotateZ() */


Vector VectTransform( Vector v, Matrix m )
{
        Vector	vPrime;

        vPrime.x = v.x*m[0][0] + v.y*m[1][0] + v.z*m[2][0] + m[3][0];
        vPrime.y = v.x*m[0][1] + v.y*m[1][1] + v.z*m[2][1] + m[3][1];
        vPrime.z = v.x*m[0][2] + v.y*m[1][2] + v.z*m[2][2] + m[3][2];

        return(vPrime);
} /* VectTransform() */


/*
 * A simplified multifractal function
 */
double Multifractal( Vector pos, double H, double lacunarity, double octaves, 
		     double zero_offset )
{
        double  y, f;
        int     i;

        y = 1.0;
        f = 1.0;

        for (i=0; i<octaves; i++) {
                y *= zero_offset + f * Noise3( pos );
                f *= H;
                pos.x *= lacunarity;
                pos.y *= lacunarity;
                pos.z *= lacunarity;
        } /* for */

        return y;

} /* Multifractal() */