package Urb::Field::Occlusion;

use 5.010;
use strict;
use warnings;
use POSIX;
use Urb::Field;
use Urb::Misc::CIESky;
use Urb::Math qw /distance_2d subtract_2d is_angle_between angle_2line points_2line line_intersection/;

use base qw /Urb::Field/;

our $PI = atan2 (0,-1);

=head1 NAME

Urb::Field::Occlusion -  Occlusion field in 3D space

=head1 SYNOPSIS

  use Urb::Field::Occlusion;
  $field = new Urb::Field::Occlusion ($optional_overlay_field);
  $field->Deserialise ($yaml);
  $illumination = $field->CIEsky_horizontal ([12.3,-3.45,1.23]);
  $illumination = $field->CIEsky_vertical ([12.3,-3.45,1.23], 1.5708);

Alternatively instead of the Deserialise() method a field can be initialised
with a list of walls, the various accessor methods will fill the table in a
lazy fashion:

  $occlusion->{origin} = [$x,$y,$z];
  $occlusion->{_walls} = [$dom->Walls];

A 'wall' is defined by two xyz points (bottom-left, top-right) and an unused
'type'.

=head1 DESCRIPTION

=over

=item Get_int

Retrieve the occlusion angle for an integer xyz point in a given direction
(direction is range 0 .. 15):

  $rad_inclination = $field->Get_int ([$int_x,$int_y,$int_z], $int_direction);

=cut

sub Get_int
{
    my $self = shift;
    my $point = shift;
    my $id_angle = shift;
    my $id_x = int ($point->[0] - $self->{origin}->[0]);
    my $id_y = int ($point->[1] - $self->{origin}->[1]);
    my $id_z = int ($point->[2] - $self->{origin}->[2]);
    return if ($id_x < 0 or $id_y < 0 or $id_z < 0);

    my $cached = undef;
    $cached = $self->{_data}->[$id_x]->[$id_y]->[$id_z]->[$id_angle]
        if defined $self->{_data}->[$id_x]->[$id_y]->[$id_z]->[$id_angle];
    
    my $overlay = undef;
    $overlay = $self->{overlay}->Get_int ($point, $id_angle)
        if (defined $self->{overlay} and ref $self->{overlay});

    if (defined $cached)
    {
        return $cached unless defined $overlay;
        return $cached if $cached > $overlay;
        return $overlay;
    }

    my $list_ref = $self->Walls_intersecting ($point, $id_angle*$PI/8);
    my $max_inclination = 0.0;
    for my $item (@{$list_ref})
    {
        $max_inclination = $item->[0] if $max_inclination < $item->[0];
    }
    $self->{_data}->[$id_x]->[$id_y]->[$id_z]->[$id_angle] = $max_inclination;

    return $max_inclination unless defined $overlay;
    return $max_inclination if $max_inclination > $overlay;
    return $overlay;
}

=item Get_float

Retrieve the occlusion angle for an xyz point in a given direction
(direction is range 0 .. 15, result is calculated directly and not cached):

  $rad_inclination = $field->Get_float ([$x,$y,$z], $int_direction);

=cut

sub Get_float
{
    my $self = shift;
    my $point = shift;
    my $id_angle = shift;
    
    my $overlay = undef;
    $overlay = $self->{overlay}->Get_interpolated_position ($point, $id_angle)
        if (defined $self->{overlay} and ref $self->{overlay});

    my $list_ref = $self->Walls_intersecting ($point, $id_angle*$PI/8);
    my $max_inclination = 0.0;
    for my $item (@{$list_ref})
    {
        $max_inclination = $item->[0] if $max_inclination < $item->[0];
    }

    return $max_inclination unless defined $overlay;
    return $max_inclination if $max_inclination > $overlay;
    return $overlay;
}

=item Get_weighted

Retrieve zero or more occlusion angles with weights:

  @results = $field->Get_weighted ([$x,$y,$z], $int_direction);

=cut

sub Get_weighted
{
    my $self = shift;
    my $coor = shift;
    my $direction = shift;
    my ($x_0,$y_0,$z_0) = map {floor $_} (@{$coor});
    my @results;
    for my $x_int ($x_0, $x_0 +1)
    {
        my $x_weight = 1 - abs ($coor->[0] - $x_int);
        for my $y_int ($y_0, $y_0 +1)
        {
            my $y_weight = 1 - abs ($coor->[1] - $y_int);
            for my $z_int ($z_0, $z_0 +1)
            {
                my $z_weight = 1 - abs ($coor->[2] - $z_int);
                my $item = $self->Get_int ([$x_int,$y_int,$z_int], $direction);
                push @results, [$item, $x_weight * $y_weight * $z_weight] if defined $item;
            }
        }
    }

    my $total = 0;
    $total += $_->[1] for @results;
    return () unless $total;
    $_->[1] /= $total for @results;

    return @results;
}

=item Get_interpolated_position

Retrieve the occlusion angle for an xyz point in a given direction (fractional
point coordinates will be interpolated, direction is range 0 .. 15):

  $rad_inclination = $field->Get_interpolated_position ([$x,$y,$z], $int_direction);

=cut

sub Get_interpolated_position
{
    my $self = shift;
    my $coor = shift;
    my $direction = shift;

    my $result = 0;
    for my $item ($self->Get_weighted ($coor, $direction))
    {
        my ($scalar, $weight) = @{$item};
        $result += $scalar * $weight;
    }
    return $result;
}

=item Walls_intersecting

Intersect a list of walls with a 'ray' defined by xyz coordinate and plan angle
and retrieve inclination and type for any that intersect:

  $list_ref = $field->Walls_intersecting ([1.0,2.0,3.0], 1.570796);

=cut

sub Walls_intersecting
{
    my $self = shift;
    my $point = shift;
    my $angle = shift;
    my $result = [];
    for my $wall (@{$self->{_walls}})
    {
        my $d_height = $wall->[1]->[2] - $point->[2];
        next if $d_height < 0; # not interested in lower walls

        my $vec_a = subtract_2d ($wall->[0], $point);
        my $vec_b = subtract_2d ($wall->[1], $point);
        next if ($vec_a->[0] == 0 and $vec_b->[0] == 0);

        my $heading_a = atan2 ($vec_a->[1], $vec_a->[0]);
        my $heading_b = atan2 ($vec_b->[1], $vec_b->[0]);
        next unless is_angle_between ($angle, $heading_a, $heading_b);

        my $line_ray = angle_2line ($point, $angle);
        my $line_wall = points_2line ($wall->[0], $wall->[1]);
        my $intersection = line_intersection ($line_ray, $line_wall) || next;
        my $distance = distance_2d ($intersection, $point);
        next if $distance <= 0.01; # close walls may be a self reference

        my $inclination = atan2 ($d_height, $distance);
        push @{$result}, [$inclination, $wall->[2]];
    }
    return $result;
}

=item Walls_lowhigh

Get the z-range of a list of walls:

  ($min_z, $max_z) = $field->Walls_lowhigh;

=cut

sub Walls_lowhigh
{
    my $self = shift;
    my $walls = $self->{_walls};
    return unless scalar @{$walls};
    my $min = $walls->[0]->[0]->[2];
    my $max = $walls->[0]->[1]->[2];
    for my $wall (@{$walls})
    {
        $min = $wall->[0]->[2] if $wall->[0]->[2] < $min;
        $max = $wall->[1]->[2] if $wall->[1]->[2] > $max;
    }
    return ($min, $max);
}

=item CIEsky_horizontal

Calculate a CIE average sky illumination for a horizontal surface at $coor
using the current occlusion field as a horizon-line mask:

  $illumination = $field->CIEsky_horizontal ($coor);

Result is illumination beteen 0.0 and 2.44346095279206

=cut

sub CIEsky_horizontal
{
    my $self = shift;
    my $coor = shift;
    my $illumination = 0;

    for my $int_direction (0 .. 15)
    {
        my $theta;
        $theta = ($PI/2) - $self->Get_interpolated_position ($coor, $int_direction)
            unless $self->{nocache};
        $theta = ($PI/2) - $self->Get_float ($coor, $int_direction)
            if $self->{nocache};
        $illumination += Urb::Misc::CIESky::average_horizontal ($theta) /8;
    }
    return $illumination;
}

=item Sun_horizontal

Calculate number of minutes of sunlight for a point at $coor using current
occlusion field as the horizon-line mask:

  $minutes = $field->Sun_horizontal ($sun, $coor);

($sun is an Urb::Misc::Sun object for the current latitude)

=cut

sub Sun_horizontal
{
    my $self = shift;
    my $sun = shift;
    my $coor = shift;
    my $minutes = 0;

    for my $int_direction (0 .. 15)
    {
        my $theta;
        $theta = $self->Get_interpolated_position ($coor, $int_direction)
            unless $self->{nocache};
        $theta = $self->Get_float ($coor, $int_direction)
            if $self->{nocache};
        # oops, sun table is north/clockwise so convert from east/anticlockwise
        my $rad_direction = (4 -$int_direction) * $PI/8;
        $rad_direction += (2*$PI) if $rad_direction <= (0-$PI);
        $minutes += $sun->Minutes ($rad_direction, $theta);
        #$minutes += $sun->Minutes ((4 -$int_direction) * $PI/8, $theta);
    }
    return $minutes;
}

=item CIEsky_vertical

Calculate a CIE average sky illumination for a vertical surface at $coor,
facing in direction $radians using the current occlusion field as a horizon-line
mask:

  $illumination = $field->CIEsky_vertical ($coor, $radians);

Result is illumination beteen 0.0 and 0.608335238775694

=back

=cut

sub CIEsky_vertical
{
    my $self = shift;
    my $coor = shift;
    my $rad_angle = shift;
    my $illumination = 0;

    for my $int_direction (0 .. 15)
    {
        my $rad_direction = $PI * $int_direction / 8;
        my $rad_incidence = abs ($rad_angle - $rad_direction);
        $rad_incidence = abs ($rad_incidence - (2*$PI)) if $rad_incidence > $PI;
        next if $rad_incidence > $PI/2;

        my $theta;
        $theta = ($PI/2) - $self->Get_interpolated_position ($coor, $int_direction)
            unless $self->{nocache};
        $theta = ($PI/2) - $self->Get_float ($coor, $int_direction)
            if $self->{nocache};
        $illumination += cos ($rad_incidence) * Urb::Misc::CIESky::average_vertical ($theta) /8;
    }
    return $illumination;
}

1;

__END__

=head1 SEE ALSO

L<Urb::Dom>
L<Urb::Field>
L<Urb::Misc::CIESky>

=head1 AUTHOR

Bruno Postle <bruno@postle.net>

=head1 COPYRIGHT AND LICENSE

Copyright (C) 2013 by Bruno Postle

This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself, either Perl version 5.10.1 or,
at your option, any later version of Perl 5 you may have available.

=cut
