#!/local/bin/perl

@eye = (1.5, 1.7, 2.0);
@ref = (0.5, 0.5, 0.5);
@vup = (0.0, 0.0, 1.0);

@Z = ($eye[0]-$ref[0], $eye[1]-$ref[1], $eye[2]-$ref[2]);

$lZ = sqrt($Z[0]*$Z[0] + $Z[1]*$Z[1] + $Z[2]*$Z[2]);
@z = ($Z[0]/$lZ, $Z[1]/$lZ, $Z[2]/$lZ);

@X = ($vup[1]*$z[2]-$vup[2]*$z[1],
      $vup[2]*$z[0]-$vup[0]*$z[2],
      $vup[0]*$z[1]-$vup[1]*$z[0]);
$lX = sqrt($X[0]*$X[0] + $X[1]*$X[1] + $X[2]*$X[2]);
@x = ($X[0]/$lX, $X[1]/$lX, $X[2]/$lX);

@y = ($z[1]*$x[2]-$z[2]*$x[1],
      $z[2]*$x[0]-$z[0]*$x[2],
      $z[0]*$x[1]-$z[1]*$x[0]);

$Mv[0][0] = $x[0]; $Mv[0][1] = $x[1]; $Mv[0][2] = $x[2];
$Mv[0][3] = -($x[0]*$eye[0] + $x[1]*$eye[1] + $x[2]*$eye[2]);

$Mv[1][0] = $y[0]; $Mv[1][1] = $y[1]; $Mv[1][2] = $y[2];
$Mv[1][3] = -($y[0]*$eye[0] + $y[1]*$eye[1] + $y[2]*$eye[2]);

$Mv[2][0] = $z[0]; $Mv[2][1] = $z[1]; $Mv[2][2] = $z[2];
$Mv[2][3] = -($z[0]*$eye[0] + $z[1]*$eye[1] + $z[2]*$eye[2]);

$Mv[3][0] = 0; $Mv[3][1] = 0; $Mv[3][2] = 0; $Mv[3][3] = 1;

print "Mview=\n";
for ($i=0; $i<4; $i++) {
    print "[";
    for ($j=0; $j<4; $j++) {
	printf " %6.3f",$Mv[$i][$j];
    }
    print "]\n";
}

$Mp[0][0] = 1; $Mp[0][1] = 0; $Mp[0][2] = 0; $Mp[0][4] = 0;
$Mp[1][0] = 0; $Mp[1][1] = 1; $Mp[1][2] = 0; $Mp[1][4] = 0;
$Mp[2][0] = 0; $Mp[2][1] = 0; $Mp[2][2] = 1; $Mp[2][4] = 0;
$Mp[3][0] = 0; $Mp[3][1] = 0; $Mp[3][2] =-1; $Mp[3][4] = 0;

for ($i=0; $i<4; $i++) {
    for ($j=0; $j<4; $j++) {
	$s = 0;
	for ($k=0; $k<4; $k++) {
	    $s += $Mp[$i][$k] * $Mv[$k][$j];
	}
	$Mpv[$i][$j] = $s;
    }
}

print "Mproj * Mview=\n";
for ($i=0; $i<4; $i++) {
    print "[";
    for ($j=0; $j<4; $j++) {
	printf " %6.3f",$Mpv[$i][$j];
    }
    print "]\n";
}

$n = 0;
for ($i=0; $i<=1; $i++) {
    for ($j=0; $j<=1; $j++) {
	for ($k=0; $k<=1; $k++) {
	    $p[$n][0] = $i;
	    $p[$n][1] = $j;
	    $p[$n][2] = $k;
	    $p[$n][3] = 1;
	    $n++;
	}
    }
}

$p[8][0]  = 1; $p[8][1]  = 0; $p[8][2]  = 0; $p[8][3]  = 0;
$p[9][0]  = 0; $p[9][1]  = 1; $p[9][2]  = 0; $p[9][3]  = 0;
$p[10][0] = 0; $p[10][1] = 0; $p[10][2] = 1; $p[10][3] = 0;

print "  vertex                  vertex in CCS          vertex in NDCS\n";
for ($i=0; $i<$n; $i++) {
    printf "(%d %d %d %d)  ",$p[$i][0],$p[$i][1],$p[$i][2],$p[$i][3];
    for ($j=0; $j<4; $j++) {
	$s = 0;
	for ($k=0; $k<4; $k++) {
	    $s += $Mpv[$j][$k] * $p[$i][$k];
	}
	$pccs[$j] = $s;
    }
    printf "(%6.3f %6.3f %6.3f %6.3f)  ",$pccs[0],$pccs[1],$pccs[2],$pccs[3];
    for ($j=0; $j<3; $j++) {
	$pndcs[$i][$j] = $pccs[$j] / $pccs[3];
    }
    printf "(%6.3f %6.3f %6.3f)\n",$pndcs[$i][0],$pndcs[$i][1],$pndcs[$i][2];
}


print "direction                         in CCS    Vanishing point NDCS\n";
for ($i=8; $i<=10; $i++) {
    printf "(%d %d %d %d)  ",$p[$i][0],$p[$i][1],$p[$i][2],$p[$i][3];
    for ($j=0; $j<4; $j++) {
	$s = 0;
	for ($k=0; $k<4; $k++) {
	    $s += $Mpv[$j][$k] * $p[$i][$k];
	}
	$pccs[$j] = $s;
    }
    printf "(%6.3f %6.3f %6.3f %6.3f)  ",$pccs[0],$pccs[1],$pccs[2],$pccs[3];
    for ($j=0; $j<3; $j++) {
	$pndcs[$i][$j] = $pccs[$j] / $pccs[3];
    }
    printf "(%6.3f %6.3f %6.3f)\n",$pndcs[$i][0],$pndcs[$i][1],$pndcs[$i][2];
}

open (PS,">cube.ps")
    or die "Can't write to cube.ps: $!";

print PS "%%!PS-Adobe-2.0\n";
print PS "2 setlinewidth\n";
print PS "newpath\n";

@edges = (0,1, 1,3, 3,2, 2,0,
	  4,5, 5,7, 7,6, 6,4,
	  0,4, 1,5, 2,6, 3,7);

for ($i=0; $i<=$#edges; $i+=2) {
    $e1 = $edges[$i];
    $e2 = $edges[$i+1];
    $x1 = int($pndcs[$e1][0] * 200 + 250);
    $y1 = int($pndcs[$e1][1] * 200 + 250);
    $x2 = int($pndcs[$e2][0] * 200 + 250);
    $y2 = int($pndcs[$e2][1] * 200 + 250);

    print PS "$x1 $y1 moveto $x2 $y2 lineto\n";
}

print PS "stroke\n";

for ($i=8; $i<=10; $i++) {
    $x = int($pndcs[$i][0] * 200 + 250);
    $y = int($pndcs[$i][1] * 200 + 250);

    print PS "newpath $x $y 2 0 360 arc stroke\n";
}
 
@edges = (0,8, 1,8, 2,8, 3,8,
	  0,9, 1,9, 5,9, 4,9,
	  0,10, 2,10, 6,10, 4,10);

print PS "0 setlinewidth\n";

for ($i=0; $i<=$#edges; $i+=2) {
    $e1 = $edges[$i];
    $e2 = $edges[$i+1];
    $x1 = int($pndcs[$e1][0] * 200 + 250);
    $y1 = int($pndcs[$e1][1] * 200 + 250);
    $x2 = int($pndcs[$e2][0] * 200 + 250);
    $y2 = int($pndcs[$e2][1] * 200 + 250);

    print PS "$x1 $y1 moveto $x2 $y2 lineto\n";
}

print PS "stroke\n";

print PS "showpage\n";
close PS;
 

exit 0;
