aobench php port. About 200 times slower than C implementation. c.f. https://code.google.com/p/aobench/, http://tai2.net/docs/aobench_php/

ao.php

<?php

error_reporting(E_ALL);

define('WIDTH', 256);
define('HEIGHT', 256);
define('NSUBSAMPLES', 2);
define('NAO_SAMPLES', 8);

class vec
{
    public $x;
    public $y;
    public $z;

    function __construct($x = 0, $y = 0, $z = 0) {
        $this->x = $x;
        $this->y = $y;
        $this->z = $z;
    }
}

class Isect {
    public $t;
    public $p;
    public $n;
    public $hit;

    function __construct($t, $p, $n, $hit) {
        $this->t = $t;
        $this->p = $p;
        $this->n = $n;
        $this->hit = $hit;
    }
}

class Sphere
{
    public $center;
    public $radius;

    function __construct($center, $radius) {
        $this->center = $center;
        $this->radius = $radius;
    }
}

class Plane
{
    public $p;
    public $n;

    function __construct($p, $n) {
        $this->p = $p;
        $this->n = $n;
    }
}

class Ray
{
    public $org;
    public $dir;

    function __construct($org, $dir) {
        $this->org = $org;
        $this->dir = $dir;
    }
}

$spheres = array(
    new Sphere(new vec(-2.0, 0.0, -3.5), 0.5),
    new Sphere(new vec(-0.5, 0.0, -3.0), 0.5),
    new Sphere(new vec( 1.0, 0.0, -2.2), 0.5));

$plane = new Plane(
    new vec(0.0, -0.5, 0.0),
    new vec(0.0, 1.0, 0.0));

function vdot($v0, $v1) {
    return $v0->x * $v1->x + $v0->y * $v1->y + $v0->z * $v1->z;
}

function vcross($c, $v0, $v1) {
    $c->x = $v0->y * $v1->z - $v0->z * $v1->y;
    $c->y = $v0->z * $v1->x - $v0->x * $v1->z;
    $c->z = $v0->x * $v1->y - $v0->y * $v1->x;
}

function vnormalize($c) {
    $length = sqrt(vdot($c, $c));

    if (abs($length) > 1.0e-17) {
        $c->x /= $length;
        $c->y /= $length;
        $c->z /= $length;
    }
}

function ray_sphere_intersect($isect, $ray, $sphere) {
    $rs = new vec(
        $ray->org->x - $sphere->center->x,
        $ray->org->y - $sphere->center->y,
        $ray->org->z - $sphere->center->z);

    $B = vdot($rs, $ray->dir);
    $C = vdot($rs, $rs) - $sphere->radius * $sphere->radius;
    $D = $B * $B - $C;

    if ($D > 0.0) {
        $t = -$B - sqrt($D);

        if (($t > 0.0) && ($t < $isect->t)) {
            $isect->t = $t;
            $isect->hit = true;

            $isect->p->x = $ray->org->x + $ray->dir->x * $t;
            $isect->p->y = $ray->org->y + $ray->dir->y * $t;
            $isect->p->z = $ray->org->z + $ray->dir->z * $t;

            $isect->n->x = $isect->p->x - $sphere->center->x;
            $isect->n->y = $isect->p->y - $sphere->center->y;
            $isect->n->z = $isect->p->z - $sphere->center->z;

            vnormalize($isect->n);
        }
    }
}

function ray_plane_intersect($isect, $ray, $plane) {
    $d = -vdot($plane->p, $plane->n);
    $v =  vdot($ray->dir, $plane->n);

    if (abs($v) < 1.0e-17) return;

    $t = -(vdot($ray->org, $plane->n) + $d) / $v;

    if ($t > 0.0 && $t < $isect->t) {
        $isect->t = $t;
        $isect->hit = true;
        
        $isect->p->x = $ray->org->x + $ray->dir->x * $t;
        $isect->p->y = $ray->org->y + $ray->dir->y * $t;
        $isect->p->z = $ray->org->z + $ray->dir->z * $t;

        $isect->n = $plane->n;
    }
}

function orthoBasis($basis, $n) {
    $basis[2]->x = $n->x; $basis[2]->y = $n->y; $basis[2]->z = $n->z;
    $basis[1]->x = 0.0;   $basis[1]->y = 0.0;   $basis[1]->z = 0.0;

    if ($n->x < 0.6 && $n->x > -0.6) {
        $basis[1]->x = 1.0;
    } else if ($n->y < 0.6 && $n->y > -0.6) {
        $basis[1]->y = 1.0;
    } else if ($n->z < 0.6 && $n->z > -0.6) {
        $basis[1]->z = 1.0;
    } else {
        $basis[1]->x = 1.0;
    }

    vcross($basis[0], $basis[1], $basis[2]);
    vnormalize($basis[0]);

    vcross($basis[1], $basis[2], $basis[0]);
    vnormalize($basis[1]);
}

function float_rand() {
    return mt_rand() / mt_getrandmax();
}

function ambient_occlusion($col, $isect)
{
    global $plane, $spheres;
    $ntheta = NAO_SAMPLES;
    $nphi = NAO_SAMPLES;
    $eps = 0.0001;

    $p = new vec(
        $isect->p->x + $eps * $isect->n->x,
        $isect->p->y + $eps * $isect->n->y,
        $isect->p->z + $eps * $isect->n->z);

    $basis = array(new vec(), new vec(), new vec());
    orthoBasis($basis, $isect->n);

    $occlusion = 0.0;

    for ($j = 0; $j < $ntheta; $j++) {
        for ($i = 0; $i < $nphi; $i++) {
            $theta = sqrt(float_rand());
            $phi   = 2.0 * M_PI * float_rand();

            $x = cos($phi) * $theta;
            $y = sin($phi) * $theta;
            $z = sqrt(1.0 - $theta * $theta);

            // local -> global
            $rx = $x * $basis[0]->x + $y * $basis[1]->x + $z * $basis[2]->x;
            $ry = $x * $basis[0]->y + $y * $basis[1]->y + $z * $basis[2]->y;
            $rz = $x * $basis[0]->z + $y * $basis[1]->z + $z * $basis[2]->z;

            $ray = new Ray($p, new vec($rx, $ry, $rz));

            $occIsect = new Isect(1.0e+17, new vec(), new vec(), false);

            ray_sphere_intersect($occIsect, $ray, $spheres[0]); 
            ray_sphere_intersect($occIsect, $ray, $spheres[1]); 
            ray_sphere_intersect($occIsect, $ray, $spheres[2]); 
            ray_plane_intersect ($occIsect, $ray, $plane); 

            if ($occIsect->hit) $occlusion += 1.0;
            
        }
    }

    $occlusion = ($ntheta * $nphi - $occlusion) / (double)($ntheta * $nphi);

    $col->x = $occlusion;
    $col->y = $occlusion;
    $col->z = $occlusion;
}

function clamp($f) {
    $i = (int)($f * 255.25);

    if ($i < 0) $i = 0;
    if ($i > 255) $i = 255;

    return $i;
}

function render(&$img, $w, $h, $nsubsamples)
{
    global $spheres, $plane;
    $fimg = array_fill(0, $w * $h * 3, 0.0);

    for ($y = 0; $y < $h; $y++) {
        for ($x = 0; $x < $w; $x++) {

            for ($v = 0; $v < $nsubsamples; $v++) {
                for ($u = 0; $u < $nsubsamples; $u++) {
                    $px = ($x + ($u / (double)$nsubsamples) - ($w / 2.0)) / ($w / 2.0);
                    $py = -($y + ($v / (double)$nsubsamples) - ($h / 2.0)) / ($h / 2.0);

                    $ray = new Ray(
                        new vec(0.0, 0.0, 0.0),
                        new vec($px, $py, -1.0));
                    vnormalize($ray->dir);

                    $isect = new Isect(1.0e+17, new vec(), new vec(), false);
                    ray_sphere_intersect($isect, $ray, $spheres[0]);
                    ray_sphere_intersect($isect, $ray, $spheres[1]);
                    ray_sphere_intersect($isect, $ray, $spheres[2]);
                    ray_plane_intersect($isect, $ray, $plane);

                    if ($isect->hit) {
                        $col = new vec();
                        ambient_occlusion($col, $isect);

                        $fimg[3 * ($y * $w + $x) + 0] += $col->x;
                        $fimg[3 * ($y * $w + $x) + 1] += $col->y;
                        $fimg[3 * ($y * $w + $x) + 2] += $col->z;
                    }

                }
            }

            $fimg[3 * ($y * $w + $x) + 0] /= $nsubsamples * $nsubsamples;
            $fimg[3 * ($y * $w + $x) + 1] /= $nsubsamples * $nsubsamples;
            $fimg[3 * ($y * $w + $x) + 2] /= $nsubsamples * $nsubsamples;
        
            $img[3 * ($y * $w + $x) + 0] = chr(clamp($fimg[3 * ($y * $w + $x) + 0]));
            $img[3 * ($y * $w + $x) + 1] = chr(clamp($fimg[3 * ($y * $w + $x) + 1]));
            $img[3 * ($y * $w + $x) + 2] = chr(clamp($fimg[3 * ($y * $w + $x) + 2]));
        }
    }

}

function saveppm($fname, $w, $h, &$img) {
    $fp = fopen($fname, "wb");
    assert($fp !== false);

    fprintf($fp, "P6\n");
    fprintf($fp, "%d %d\n", $w, $h);
    fprintf($fp, "255\n");
    fwrite($fp, $img);
    fclose($fp);
}

function main()
{
    $img = str_repeat(chr(0), WIDTH * HEIGHT * 3);

    render($img, WIDTH, HEIGHT, NSUBSAMPLES);

    saveppm("ao.ppm", WIDTH, HEIGHT, $img); 
}

main();