#!/usr/bin/env perl
use strict;
use warnings;
use Math::VectorReal;
use Math::Trig;
use Imager;
use Imager::Fill;
my $viewerPos = vector(0,-0.8,0);
my $viewerForward = vector(0,0.2,1); # also view plane pos
my $viewerUp = vector(0,1,0);
my $viewerPlaneX = ($viewerForward x $viewerUp)->norm;
my $viewerPlaneY = $viewerUp->norm;

my $floorUp = vector(0,1,0);
my $floorX = vector(1,0,0);
my $floorY = vector(0,0,1);
my $floorPos = vector(0,-1,0);

my $spherePos = vector(0,0,2);
my $sphereRad = 0.3;


my $imgSize = 256;

my $img = new Imager(
    xsize => $imgSize,
    ysize => $imgSize,
) or die($!);


my $rez = 1.0/$imgSize;
my ($i,$j)=(0,0);

for (my $x=-0.5;$x<0.5;$x+=$rez) {
    $j=$imgSize-1;
    for (my $y=-0.5;$y<0.5;$y+=$rez) {
        my $color = castRay($x,$y);

        my @pixel = map {
            my $v=$_*255.0;
            if ($v>255){$v=255;}
            if ($v<0) { $v=0; }
            $v;
        } @$color;
        $img->setpixel(x=>$i,y=>$j,color=>Imager::Color->new(@pixel));
        $j--;
    }
    $i++;
}

$img->write(file=>"ray.png") or die($!);


exit 0;

sub castRay {
    my ($x,$y) = @_;
    my $r0 = $viewerPos;

    my $rayScreenIntersection =
        $viewerPos + $viewerForward
        + $x * $viewerPlaneX
        + $y * $viewerPlaneY;

    my $dir = $rayScreenIntersection - $r0;

    return drawRay($r0,$dir,0);

}

sub drawRay {
    my ($r0,$dir,$fromSphere) = @_;

    my $intersection = raySphereIntersection($r0,$dir);
    if ( (!$fromSphere) && defined($intersection)) {
        # We intersect the sphere.

        # Draw this segment.


        my $dirOld = $dir->norm;
        my $normal = ($intersection-$spherePos)->norm;

        my $cosTheta = $dirOld . $normal;

        my $dirNew = $dirOld - 2.0*$cosTheta * $normal;

        my $reflectedColor = drawRay($intersection,$dirNew,1);

        # Let's make it a bit reddish.

        my ($rr,$rg,$rb) = @$reflectedColor;

        $reflectedColor = [ $rr*1.2,$rg/1.2,$rb/1.2 ];

        return $reflectedColor;


    } elsif (defined($intersection = rayFloorIntersection($r0,$dir))) {
            # We hit the floor.

            my $floorVector = $intersection - $floorPos;
            my $tileSize = 0.5;

            my $interX = int(($floorVector . $floorX)/$tileSize);
            my $interY = int(($floorVector . $floorY)/$tileSize);

            # ($interX, $interY) give floor tile coordinates.

            my $xType = ($interX % 2 ) ;
            my $yType = ($interY % 2 ) ;

            my $brightness = 5.0/($floorVector . $floorVector);
            if ($brightness>1.0) { $brightness = 1.0; }
            if ($brightness<0.0) { $brightness = 0.0; }
            my $color;
            if ($xType==$yType) {
                $color = [$brightness,$brightness,$brightness];
            } else {
                $color = [$brightness,0,$brightness];
            }

            return $color;
    } else {
            # Neither ray nor floor, therefore sky.

            my $cosTheta = ($dir . $floorUp) / ($dir->length * $floorUp->length);
            $cosTheta *= 0.5;

            my @horizonColor = ( 1,0.0,0.0);
            my @skyColor = ( 0.0,0.0,1.0 );

            my @color;
            for (0..2) {
                push @color,$horizonColor[$_]*$cosTheta +
                (1.0-$cosTheta)*$skyColor[$_];
            }

            return \@color;
    }

}

sub rayFloorIntersection {
    my ($r0,$dir) = @_;

    my $denom = $floorUp . $dir;
    if ($denom == 0.0) { return ; } # line parallel to floor

    my $lambda = ( $floorUp . (  $floorPos - $r0 ) ) / $denom;

    if ($lambda>0.0) {
        my $intersection = $r0 + $lambda*$dir;
        return $intersection;
    }

    return;
}

sub raySphereIntersection {
    my ($r0,$dir) = @_;

    my $a = $dir . $dir;
    my $b = 2.0 * ($dir . ( $r0 - $spherePos ) );
    my $c = (($r0 - $spherePos) . ($r0 - $spherePos)) - $sphereRad*$sphereRad;

    my $discriminant = $b*$b - 4.0*$a*$c;
    if ($discriminant<0.0) {
        return ;
    } else {
        my $root = sqrt($discriminant);
        my $lambda0 = (-$b - $root)/(2.0*$a);
        my $lambda1 = (-$b - $root)/(2.0*$a);

        my $lambda;
        if ($lambda1<$lambda0) { $lambda=$lambda1; } else { $lambda=$lambda0; }

        my $intersection = $r0 + $lambda*$dir;
        return $intersection;
    }
}