use strict;
use warnings;
use feature qw/ say /;

use POSIX qw/ round /;      # use 5.022;
use List::Util qw/ sum /;

use Math::Geometry::Delaunay qw/ TRI_CONSTRAINED /;

use constant EPSILON_NEGATIVE => -1e-6;

sub key { pack 'II', @{ $_[ 0 ]}}

my @red_in  = ( [58,48], [108,155], [186,80], [255,191], [331,48] );
my @red_out = ( [471,15], [531,141], [603,90], [682,227], [747,107] );

my %mapped  = map { key( $red_in[ $_ ]), $red_out[ $_ ]} 0 .. $#red_in;

my $tri = Math::Geometry::Delaunay-> new( TRI_CONSTRAINED );
$tri-> addPoints( \@red_in );
$tri-> triangulate;

my @blue_in = ( [125,73], [197,158], [282,94] );

BLUE: for my $blue ( @blue_in ) {
    TRI: for my $elem ( @{ $tri-> elements }) {
        my $y23 = $elem-> [ 1 ][ 1 ] -
                  $elem-> [ 2 ][ 1 ];
        my $x32 = $elem-> [ 2 ][ 0 ] -
                  $elem-> [ 1 ][ 0 ];
        my $x13 = $elem-> [ 0 ][ 0 ] -
                  $elem-> [ 2 ][ 0 ];
        my $y13 = $elem-> [ 0 ][ 1 ] -
                  $elem-> [ 2 ][ 1 ];
        
        my $denominator = $y23 * $x13 + $x32 * $y13;
        
        my $xx3 = $blue-> [ 0 ] - $elem-> [ 2 ][ 0 ];
        my $yy3 = $blue-> [ 1 ] - $elem-> [ 2 ][ 1 ];

        my @weights;
        next TRI if EPSILON_NEGATIVE > ( $weights[ 0 ] =
            ( $y23 * $xx3 + $x32 * $yy3 ) / $denominator );
        next TRI if EPSILON_NEGATIVE > ( $weights[ 1 ] =
            ( -$y13 * $xx3 + $x13 * $yy3 ) / $denominator );
        next TRI if EPSILON_NEGATIVE > ( $weights[ 2 ] =
            1 - $weights[ 0 ] - $weights[ 1 ]);

        printf "[%3d,%3d] => [%3d,%3d]\n", 
            @$blue,
            map {
                my $coord = $_;
                round sum map {
                    $weights[ $_ ] * $mapped{ key $elem-> [ $_ ]}[ $coord ]
                } 0 .. 2    # 3 vertices
            } 0, 1;         # x, y
        next BLUE
    }
    die "point @$blue outside convex hull"
}

__END__

[125, 73] => [541, 63]
[197,158] => [621,174]
[282, 94] => [701,137]