<?php
	$filename = 'sphere.wrl';
	$vrml_str = file_get_contents($filename);	// Convert the file into a string
	$vrml_tok = str_split($vrml_str);			// Split the string into an array of each character

	// Traverse the array and eliminate anything that isn't of numeric relevance
	// Also retain brackets to allow dissection of quantities
	for($i = 0; $i < sizeof($vrml_tok); $i++){
		switch($vrml_tok[$i]){
			case '0':
			case '1':
			case '2':
			case '3':
			case '4':
			case '5':
			case '6':
			case '7':
			case '8':
			case '9':
			case '-':
			case '.':
			case ' ':
			case '{':
			case '}':
			case '[':
			case ']':
				break;
			default:
				$vrml_tok[$i] = "";
				break;
		}
	}

	$vrml_tok = implode($vrml_tok);				// Join the array into a string
	$vrml_tok = explode(' 01 ', $vrml_tok);		// Split the string into arrays
	array_shift($vrml_tok);						// Destroy the first entry in the array

	// Extract the relevant pieces of info from the strings in each array
	// $vrml_1: Object's center location
	// $vrml_2: Object's rotational orientation
	// $vrml_3: Object's vertices
	// $vrml_4: Object's polygons
	for($i = 0; $i < sizeof($vrml_tok); $i++){
		$vrml_tok[$i] = substr($vrml_tok[$i], strpos($vrml_tok[$i], "{") + 2);
		$vrml_1[$i] = substr($vrml_tok[$i], 0, strpos($vrml_tok[$i], "[") - 1);
		$vrml_tok[$i] = substr($vrml_tok[$i], strpos($vrml_tok[$i], "[") + 1);
		$vrml_tok[$i] = substr($vrml_tok[$i], strpos($vrml_tok[$i], "[") + 1);
		$vrml_3[$i] = substr($vrml_tok[$i], 0, strpos($vrml_tok[$i], "]"));
		$vrml_tok[$i] = substr($vrml_tok[$i], strpos($vrml_tok[$i], "[") + 1);
		$vrml_4[$i] = substr($vrml_tok[$i], 0, strpos($vrml_tok[$i], "]"));
		$vrml_tok[$i] = substr($vrml_tok[$i], strpos($vrml_tok[$i], "]") + 1);

		$vrml_1[$i] = explode(" ", $vrml_1[$i]);
		$vrml_2[$i] = array_slice($vrml_1[$i], 3);
		$vrml_1[$i] = array_slice($vrml_1[$i], 0, 3);
		$vrml_3[$i] = explode(" ", $vrml_3[$i]);
		$vrml_4[$i] = explode(" ", $vrml_4[$i]);
		$temp = array_slice($vrml_4[$i], 0, 3);
		for($j = 4; $j < sizeof($vrml_4[$i]); $j += 4){
			$temp = array_merge($temp, array_slice($vrml_4[$i], $j, 3));
		}
		$vrml_4[$i] = $temp;
	}

	// Swap the values of x z y to become x y z
	for($i = 0; $i < sizeof($vrml_1); $i++){
		$temp = $vrml_1[$i][1];
		$vrml_1[$i][1] = $vrml_1[$i][2];
		$vrml_1[$i][2] = $temp;
	}
	for($i = 0; $i < sizeof($vrml_2); $i++){
		$temp = $vrml_2[$i][1];
		$vrml_2[$i][1] = $vrml_2[$i][2];
		$vrml_2[$i][2] = $temp;
	}
	for($i = 0; $i < sizeof($vrml_3); $i++){
		for($j = 0; $j < sizeof($vrml_3[0]); $j += 3){
			$temp = $vrml_3[$i][$j + 1];
			$vrml_3[$i][$j + 1] = $vrml_3[$i][$j + 2];
			$vrml_3[$i][$j + 2] = $temp;
		}
	}

	// Compute the orientation of the object so as to find a rotational transform matrix
	for($i = 0; $i < sizeof($vrml_2); $i++){
		$theta_1 = asin(abs($vrml_2[$i][2]));
		if($vrml_2[$i][1] >= 0 && $vrml_2[$i][2] >= 0){
			$theta_1 = -(pi() / 2 - $theta_1);
		}else{
			if($vrml_2[$i][1] < 0 && $vrml_2[$i][2] >= 0){
				$theta_1 = pi() / 2 - $theta_1;
			}else{
				if($vrml_2[$i][1] < 0 && $vrml_2[$i][2] < 0){
					$theta_1 = pi() / 2 + $theta_1;
				}else{
					if($vrml_2[$i][1] >= 0 && $vrml_2[$i][2] < 0){
						$theta_1 = -pi() / 2 - $theta_1;
					}
				}
			}
		}

		$theta_2 = acos(abs($vrml_2[$i][0]) / sqrt(pow($vrml_2[$i][0], 2) + pow($vrml_2[$i][1], 2)));
		if($vrml_2[$i][0] >= 0 && $vrml_2[$i][1] >= 0){
			$theta_2 = -(pi() / 2 - $theta_2);
		}else{
			if($vrml_2[$i][0] < 0 && $vrml_2[$i][1] >= 0){
				$theta_2 = pi() / 2 - $theta_2;
			}else{
				if($vrml_2[$i][0] < 0 && $vrml_2[$i][1] < 0){
					$theta_2 = -(pi() / 2 - $theta_2);
				}else{
					if($vrml_2[$i][0] >= 0 && $vrml_2[$i][1] < 0){
						$theta_2 = pi() / 2 - $theta_2;
					}
				}
			}
		}

		$ini_orient[$i] = identity(3);

		// Rotate the object about the world's X axis according to the computed $theta_1
		$rotX =	array(
					array(1, 0, 0),
					array(0, cos($theta_1), -sin($theta_1)),
					array(0, sin($theta_1), cos($theta_1))
				);
		$ini_orient[$i] = multiply($rotX, $ini_orient[$i]);

		// Rotate the object about the world's Z axis according to the computed $theta_2
		$rotZ = array(
					array(cos($theta_2), -sin($theta_2), 0),
					array(sin($theta_2), cos($theta_2), 0),
					array(0, 0, 1)
				);
		$ini_orient[$i] = multiply($rotZ, $ini_orient[$i]);

		// Rotate the object about its own Z axis using the value given
		$rotZ = array(
						array(cos($vrml_2[$i][3]), -sin($vrml_2[$i][3]), 0),
						array(sin($vrml_2[$i][3]), cos($vrml_2[$i][3]), 0),
						array(0, 0, 1)
				);
		$ini_orient[$i] = multiply($ini_orient[$i], $rotZ);
	}

	// Create all the objects find in the world and store them in an array
	// Initialize these as well
	for($i = 0; $i < sizeof($vrml_1); $i++){
		$worldContainer[$i] = new VirtualObject();
		$temp = array(array());
		$temp[0][0] = $vrml_1[$i][0];
		$temp[1][0] = $vrml_1[$i][1];
		$temp[2][0] = $vrml_1[$i][2];
		$worldContainer[$i]->setCenter($temp);
		$worldContainer[$i]->setOrientation($ini_orient[$i]);
		for($j = 0; $j < sizeof($vrml_3[0]); $j += 3){
			$temp[0][0] = $vrml_3[$i][$j];
			$temp[1][0] = $vrml_3[$i][$j + 1];
			$temp[2][0] = $vrml_3[$i][$j + 2];
			$worldContainer[$i]->setVertex($temp);
		}
		$worldContainer[$i]->setPolygons($vrml_4[$i]);
	}

	// The class that stores virtual objects
	class VirtualObject{
		// Member declarations
		private $center;		// Object center
		private $orientation;	// Object orientation
		private $vertices;		// Object vertices
		private $polyOrder;		// Object polygons

		// Constructor
		function __construct(){
		}

		// Method declarations

		public function setCenter($n){
			$this->center = $n;
		}

		public function getCenter(){
			return $this->center;
		}

		public function setOrientation($n){
			$this->orientation = $n;
		}

		public function getOrientation(){
			return $this->orientation;
		}

		public function setVertex($n){
			$this->vertices[] = $n;
		}

		public function getVertex($n){
			return $this->vertices[$n];
		}

		public function getVertexCount(){
			return sizeof($this->vertices);
		}

		public function setPolygons($n){
			$this->polyOrder = $n;
		}

		public function getPolygon($n){
			$i = 3 * $n;
			return	array(
						$this->getVertex($this->polyOrder[$i]),
						$this->getVertex($this->polyOrder[$i + 1]),
						$this->getVertex($this->polyOrder[$i + 2])
					);
		}

		public function getPolyCount(){
			return sizeof($this->polyOrder) / 3;
		}
	}

	// Viewable screen dimensions
	$width = 640;
	$height = 480;
	// Sets up the camera's orientation and its focal distance
	$x = sqrt(2) / 2;
	$focal = 300;
	$a =	array(
				array(-$x, 0.5, -0.5),
				array($x, 0.5, -0.5),
				array(0, -$x, -$x)
			);
	// Invert camera's plane to allow for change of basis
	$a = inverse($a);

	// Create two images
	$image1 = imagecreate($width, $height);
	$image2 = imagecreate($width, $height);
	// Set the line thickness of each image
	imagesetthickness($image1, 1);
	imagesetthickness($image2, 1);
	// Set the background colors for each image
	$blue1 = imagecolorallocate($image1, 0, 0, 127);
	$blue2 = imagecolorallocate($image2, 0, 0, 127);
	// Define some colors for each image
	$white = imagecolorallocate($image1, 255, 255, 255);
	$black = imagecolorallocate($image1, 0, 0, 0);
	$gray = imagecolorallocate($image2, 127, 127, 127);

	// Draw each polygon
	for($i = 0; $i < $worldContainer[0]->getPolyCount(); $i++){
		// Retrieve the polygons from the cube object
		$points = $worldContainer[0]->getPolygon($i);
		// Rotate each vertex vector according to the indicated object's orientation
		$points[0] = multiply($worldContainer[0]->getOrientation(), $points[0]);
		$points[1] = multiply($worldContainer[0]->getOrientation(), $points[1]);
		$points[2] = multiply($worldContainer[0]->getOrientation(), $points[2]);
		// Modify each point of the polygon so that its now distanced to the camera's center point
		for($j = 0; $j < 3; $j++){
			$points[$j][0][0] -= 150;
			$points[$j][1][0] -= 150;
			$points[$j][2][0] -= 300 * $x;
		}
		// Write each vector in terms of the camera's orientation (change of basis)
		$points[0] = multiply($a, $points[0]);
		$points[1] = multiply($a, $points[1]);
		$points[2] = multiply($a, $points[2]);
		// Determine the normal vector and viewing vector for each polygon, respectively
		$n1 = crossProduct(addition($points[0], scale(-1, $points[1])), addition($points[2], scale(-1, $points[1])));
		$n2 = $points[1];
		// Determine the distortion due to perspective
		$focus[0] = $focal / magnitude($points[0]);
		$focus[1] = $focal / magnitude($points[1]);
		$focus[2] = $focal / magnitude($points[2]);
		// Transform the 3d points into 2d coordinates for drawing
		$x0 = $width / 2 + ($points[0][0][0] * $focus[0]);
		$x1 = $width / 2 + ($points[1][0][0] * $focus[1]);
		$x2 = $width / 2 + ($points[2][0][0] * $focus[2]);
		$y0 = $height / 2 + ($points[0][1][0] * $focus[0]);
		$y1 = $height / 2 + ($points[1][1][0] * $focus[1]);
		$y2 = $height / 2 + ($points[2][1][0] * $focus[2]);
		// Determine the cosine of the angle between the viewing plane and the polygonal plane
		$cosTheta = dotProduct($n1, $n2) / (magnitude($n1) * magnitude($n2));
		if($cosTheta <= 0){
			// White lines because these are the planes at the front of the object
			imageline($image1, $x0, $y0, $x1, $y1, $white);
			imageline($image1, $x1, $y1, $x2, $y2, $white);
			imageline($image1, $x2, $y2, $x0, $y0, $white);
			// Copy all lines on the front face
			if($x0 < $x1){
				$border[] = $x0;
				$border[] = $y0;
				$border[] = $x1;
				$border[] = $y1;
			}else{
				$border[] = $x1;
				$border[] = $y1;
				$border[] = $x0;
				$border[] = $y0;
			}
			if($x1 < $x2){
				$border[] = $x1;
				$border[] = $y1;
				$border[] = $x2;
				$border[] = $y2;
			}else{
				$border[] = $x2;
				$border[] = $y2;
				$border[] = $x1;
				$border[] = $y1;
			}
			if($x2 < $x0){
				$border[] = $x2;
				$border[] = $y2;
				$border[] = $x0;
				$border[] = $y0;
			}else{
				$border[] = $x0;
				$border[] = $y0;
				$border[] = $x2;
				$border[] = $y2;
			}
		}else{
			// Gray lines because these are the planes at the rear of the object
			imageline($image2, $x0, $y0, $x1, $y1, $gray);
			imageline($image2, $x1, $y1, $x2, $y2, $gray);
			imageline($image2, $x2, $y2, $x0, $y0, $gray);
		}
	}

	// Determine the outline of the figure by selecting only lines that are drawn once on the front face
	$silhouette = array();
	$running = true;
	do{
		// Extract individual lines
		for($i = 0; $i < 4; $i++){
			$temp[$i] = $border[0];
			array_shift($border);
		}
		// Convert both the line and group of lines to a string
		$temp = implode(" ", $temp);
		$border = implode(" ", $border);

		// Check to see if the line in question is only drawn once
		if(substr_count($border, $temp) == 0){
			// Copy the line into the silhouette
			$temp = explode(" ", $temp);
			for($i = 0; $i < 4; $i++){
				$silhouette[] = $temp[$i];
			}
		}else{
			// Otherwise remove all occurances of this line since it is not part of the silhouette
			do{
				$location = strpos($border, $temp);
				$border = substr_replace($border, "", $location, strlen($temp) + 1);
			}while(substr_count($border, $temp) > 0);
			$temp = explode(" ", $temp);
		}
		// Do this until all lines have been examined
		if(strlen($border) > 0){
			$border = explode(" ", $border);
		}else{
			$running = false;
		}
	}while($running);

	// Set the thickness to a greater value
	imagesetthickness($image1, 4);
	// Draw the silhouette of the object
	for($i = 0; $i < sizeof($silhouette); $i += 4){
		imageline($image1, $silhouette[$i], $silhouette[$i + 1], $silhouette[$i + 2], $silhouette[$i + 3], $black);
	}
	// Make transparent the background of the viewable polygon drawing
	imagecolortransparent($image1, $blue1);
	// Layer it over the non-viewable polygon drawing
	imagecopymerge($image2, $image1, 0, 0, 0, 0, $width, $height, 100);
	// Output the finished image
	header("Content-type: image/png");
	imagepng($image2);
	// Destroy the images
	imagedestroy($image1);
	imagedestroy($image2);

	// Computes the cross product of two vectors
	// $a: a row or column vector
	// $b: a row or column vector
	function crossProduct($a, $b){
		if(sizeof($a) == 1){	// Cross product of row vectors
			return	array(	// Return the cross product of row vectors as a row vector
						array(	$a[0][1] * $b[0][2] - $a[0][2] * $b[0][1],
								-1 * ($a[0][0] * $b[0][2] - $a[0][2] * $b[0][0]),
								$a[0][0] * $b[0][1] - $a[0][1] * $b[0][0]
						)
					);
		}else{	// Cross product of column vectors
			return	array(	// Return the cross product of column vectors as a column vector
						array($a[1][0] * $b[2][0] - $a[2][0] * $b[1][0]),
						array(-1 * ($a[0][0] * $b[2][0] - $a[2][0] * $b[0][0])),
						array($a[0][0] * $b[1][0] - $a[1][0] * $b[0][0])
					);
		}
	}

	// Computes the addition of two matrices
	// $A: an m x n matrix
	// $B: an m x n matrix
	function addition($A, $B){
		$M;	// The added matrix

		// Add the two matrices together
		for($i = 0; $i < sizeof($A); $i++){
			for($j = 0; $j < sizeof($A[0]); $j++){
				$M[$i][$j] = $A[$i][$j] + $B[$i][$j];
			}
		}

		return $M;	// Return the added matrix
	}

	// Computes a random matrix
	// $m: the row dimension
	// $n: the column dimension
	// $lower: the lower bound of random numbers
	// $upper: the upper bound of random numbers
	function randMatrix($m, $n, $lower, $upper){
		$M;	// A random matrix

		// Generate the random matrix
		for($i = 0; $i < $m; $i++){
			for($j = 0; $j < $n; $j++){
				$M[$i][$j] = rand($lower, $upper);
			}
		}

		return $M;	// Return the random matrix
	}

	// Computes the magnitude of a vector
	// $a: a row or column vector
	function magnitude($a){
		$result = 0;	// The magnitude of the vector

		// Use the Pythagorean theorem to determine the vector's magnitude
		for($i = 0; $i < sizeof($a); $i++){
			for($j = 0; $j < sizeof($a[0]); $j++){
				$result += pow($a[$i][$j], 2);
			}
		}

		return sqrt($result);	// Return the magnitude of the vector
	}

	// Computes the adjoint of a matrix
	// $A: a square matrix
	function adjoint($A){
		return scale(determinant($A), inverse($A));	// Returns the adjoint of a matrix
	}

	// Computes a scaled matrix
	// $k: a scalar number
	// $A: an m x n matrix
	function scale($k, $A){

		// Scale the matrix by multiplying each entry by k
		for($i = 0; $i < sizeof($A); $i++){
			for($j = 0; $j < sizeof($A[0]); $j++){
				$A[$i][$j] *= $k;
			}
		}

		return $A;	// Returns the scaled matrix
	}

	// Computes the determinant of a matrix
	// $A: a square matrix
	function determinant($A){
		$determinant = 1;	// The determinant of the matrix
		$n = sizeof($A);	// Dimension of the square matrix
		$E = identity($n);	// An elementary matrix

		// Find the determinant of a matrix by performing elementary matrix operations
		// on an augmented matrix through the use of elementary matrices in order
		// to convert the matrix into upper triangular form
		// Multiply the diagonals of the upper triangular matrix to arrive at the
		// determinant, only elementary operations I and III are allowed on the matrix.
		// If operation I is performed at any time, multiply the determinant by a -1.
		for($j = 0; $j < $n; $j++){
			for($i = $j; $i < $n; $i++){
				if($i == $j){
					if($A[$i][$j] == 0){
						for($k = $j + 1; $k < $n; $k++){
							if($A[$k][$j] != 0){
								$swap = $A[$i];
								$A[$i] = $A[$k];
								$A[$k] = $swap;
								$determinant *= -1;
								$k = $n;
							}
						}
					}
					$determinant *= $A[$i][$j];
				}else{
					$E[$i][$j] = -$A[$i][$j] / $A[$j][$j];
				}
				$A = multiply($E, $A);
				$E = identity($n);
			}
		}

		return $determinant;	// Return the determinant of the matrix
	}

	// Computes the inverse of a matrix
	// $A: a square matrix
	function inverse($A){
		$n = sizeof($A);	// Dimension of the square matrix
		$V = identity($n);	// The inverse of the matrix
		$E = identity($n);	// An elementary matrix

		// Find the inverse matrix by performing elementary matrix operations
		// on an augmented matrix through the use of elementary matrices
		for($j = 0; $j < $n; $j++){
			for($i = $j; $i < $n; $i++){
				if($i == $j){
					if($A[$i][$j] == 0){
						for($k = $j + 1; $k < $n; $k++){
							if($A[$k][$j] != 0){
								$swap = $A[$i];
								$A[$i] = $A[$k];
								$A[$k] = $swap;
								$swap = $V[$i];
								$V[$i] = $V[$k];
								$V[$k] = $swap;
								$k = $n;
							}
						}
					}
					$E[$i][$j] = 1 / $A[$i][$j];
				}else{
					$E[$i][$j] = -$A[$i][$j];
				}
				$A = multiply($E, $A);
				$V = multiply($E, $V);
				$E = identity($n);
			}
		}
		for($j = $n - 1; $j > 0; $j--){
			for($i = $j - 1; $i > -1; $i--){
				$E[$i][$j] = -$A[$i][$j];
				$A = multiply($E, $A);
				$V = multiply($E, $V);
				$E = identity($n);
			}
		}

		return $V;	// Return the inverse matrix
	}

	// Computes the product of a matrix multiplication
	// $A: An m x n matrix
	// $B: An n x o matrix
	function multiply($A, $B){
		$M;	// The product of the multiplication
		$B = transpose($B);	// Transpose the second matrix

		// Do the multiplication by making use of the dot product
		for($i = 0; $i < sizeof($A); $i++){
			for($j = 0; $j < sizeof($B); $j++){
				$M[$i][$j] = dotProduct(array($A[$i]), array($B[$j]));
			}
		}

		return $M;	// Return the product
	}

	// Computes the dot product of two vectors
	// $a: a row or column vector
	// $b; a row or column vector
	function dotProduct($a, $b){
		$result = 0;	// The result of the dot product

		// The dot product
		for($i = 0; $i < sizeof($a); $i++){
			for($j = 0; $j < sizeof($a[0]); $j++){
				$result += $a[$i][$j] * $b[$i][$j];
			}
		}

		return $result;	// Return the dot product
	}

	// Computes the transpose of a matrix
	// $A: the matrix to be transposed
	function transpose($A){
		$T;		// Transpose matrix

		// Switch the i,jth entry to the j,ith entry
		for($i = 0; $i < sizeof($A); $i++){
			for($j = 0; $j < sizeof($A[0]); $j++){
				$T[$j][$i] = $A[$i][$j];
			}
		}

		return $T;	// Return the transpose matrix
	}

	// Computes an identity matrix of size n x n
	// $n: size of the matrix wanted
	function identity($n){
		$I;		// Identity matrix

		// Fill the matrices diagonal with 1's and everything else with 0's
		for($i = 0; $i < $n; $i++){
			for($j = 0; $j < $n; $j++){
				if($i == $j){
					$I[$i][$j] = 1;
				}else{
					$I[$i][$j] = 0;
				}
			}
		}

		return $I;	// Return the identity matrix
	}
?>