function e = area3 ( s_num )

%*****************************************************************************80
%
%% AREA3 estimates the area of the hand in the hand data file.
%
%  Discussion:
%
%    This program assumes that the file 'HAND_NODES.TXT' is available.
%
%    It uses sampling to estimate the area contained within the hand outline.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2015
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Cleve Moler,
%    Numerical Computing with MATLAB,
%    SIAM, 2004,
%    ISBN13: 978-0-898716-60-3,
%    LC: QA297.M625. 
%
%  Parameters:
%
%    Input, integer S_NUM, the number of sample points.
%
%    Output, real E, the error in the area estimate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'AREA3:\n' );
  fprintf ( 1, '  Estimate the area inside a curve that approximates the\n' );
  fprintf ( 1, '  shape of a hand, using sampling.\n' );
%
%  1: DATA
%

%
%  Read the hand data.
%  It comes in as an Nx2 array.
%
  xy = load ( 'hand_nodes.txt' );
%
%  For convenience, transpose it to a 2xN array.
%
  xy = xy';
%
%  Append a copy of the first column to the end of the array.
%  This will make the PLOT program happy when it draws the outline.
%
  xy = [ xy, [ xy(:,1) ] ];
%
%  We will want to know how many columns there are.  Use the SIZE command.
%
  [ two, xy_num ] = size ( xy );
%
%  Determine the range of the data.
%
  x_min = min ( xy(1,:) );
  x_max = max ( xy(1,:) );
  y_min = min ( xy(2,:) );
  y_max = max ( xy(2,:) );
%
%  Compute the area of the rectangle that contains the data.
%
  area_total = ( x_max - x_min ) * ( y_max - y_min );
%
%  2.  SAMPLING
%

%
%  Get an array of 2xS_NUM random values, each in [0,1].
%
  s = rand ( 2, s_num );
%
%  Shift and stretch the points to fit in the rectangle.
%
  s(1,:) = x_min + s(1,:) * ( x_max - x_min );
  s(2,:) = y_min + s(2,:) * ( y_max - y_min );
%
%  Set up a vector Z to record if each sample point is inside (1)
%  or outside (0) the hand curve.
%
  z = zeros ( s_num, 1 ); 
  for j = 1 : s_num
    p = [ s(1,j); s(2,j) ];
    z(j) = polygon_contains_point_2d ( xy_num, xy, p );
  end
%
%  Use FIND to make a list of the points inside and outside the curve.
%
  j = find ( z == 1 );
  j_num = length ( j );

  k = find ( z == 0 );
  k_num = length ( k );
%
%  Estimate the hand area using the proportion of sample points that landed
%  inside the hand.
%
  area_hand = area_total * j_num / s_num;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of sample points = %d\n', s_num );
  fprintf ( 1, '  Total area = %g\n', area_total );
  fprintf ( 1, '  Estimated hand area = %g\n', area_hand );
%
%  Since the hand is a polygon, we can actually compute the exact area.
%
  area_hand2 = polygon_area_2d ( xy_num, xy );
  fprintf ( 1, '  Exact hand area = %g\n', area_hand2 );
  e = abs ( area_hand - area_hand2 );
  fprintf ( 1, '  Error = %g\n', e );
%
%  3: GRAPHICS
%

%
%  Clear any old graphics information from the screen.
%
  clf
%
%  Prepare for multiple plot commands on one image.
%
  hold on
%
%  Plot:
%    the hand curve (cyan)
%    the hand data points (blue)
%    the sample points inside the hand (red)
%    the sample points outside the hand (green)
%
  plot ( xy(1,:), xy(2,:), 'c-', 'LineWidth', 4 );
  plot ( xy(1,:), xy(2,:), 'b.', 'MarkerSize', 15 );
  plot ( s(1,j), s(2,j), 'r.', 'Markersize', 15 );
  plot ( s(1,k), s(2,k), 'g.', 'Markersize', 15 );
%
%  Use the same scale for X and Y.
%  Draw grid lines.
%  Print a title.
%
  axis equal
  grid on
  title ( 'Sample points inside (red) and outside (green) the hand data', 'Fontsize', 16 )
%
%  The plot is done.
%
  hold off
%
%  Save a copy to a file.
%
  filename = 'hand_area.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot saved to file "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'AREA3\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function area = polygon_area_2d ( n, v )

%*****************************************************************************80
%
%% POLYGON_AREA_2D computes the area of a polygon in 2D.
%
%  Discussion:
%
%    AREA = 1/2 * abs ( sum ( 1 <= I <= N ) X(I) * ( Y(I+1) - Y(I-1) ) )
%    where Y(0) should be replaced by Y(N), and Y(N+1) by Y(1).
%
%    If the vertices are given in counterclockwise order, the area
%    will be positive.  If the vertices are given in clockwise order,
%    the area will be negative.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of vertices of the polygon.
%
%    Input, real V(2,N), the vertices.
%
%    Output, real AREA, the area of the polygon.
%
  area = 0.0;

  for i = 1 : n

    im1 = i - 1;
    ip1 = i + 1;

    if ( im1 < 1 )
      im1 = n;
    end

    if ( n < ip1 )
      ip1 = 1;
    end

    area = area + v(1,i) * ( v(2,ip1) - v(2,im1) );

  end

  area = 0.5 * area;

  return
end
function inside = polygon_contains_point_2d ( n, v, p )

%*****************************************************************************80
%
%% POLYGON_CONTAINS_POINT_2D finds if a point is inside a simple polygon in 2D.
%
%  Discussion:
%
%    A simple polygon is one whose boundary never crosses itself.
%    The polygon does not need to be convex.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    M Shimrat,
%    Position of Point Relative to Polygon,
%    ACM Algorithm 112,
%    Communications of the ACM,
%    Volume 5, Number 8, page 434, August 1962.
%
%  Parameters:
%
%    Input, integer N, the number of nodes or vertices in the polygon.
%    N must be at least 3.
%
%    Input, real V(2,N), the coordinates of the vertices of the polygon.
%
%    Input, real P(2,1), the coordinates of the point to be tested.
%
%    Output, logical INSIDE, is TRUE if the point is inside the polygon.
%
  inside = 0;

  for i = 1 : n

    x1 = v(1,i);
    y1 = v(2,i);

    if ( i < n )
      x2 = v(1,i+1);
      y2 = v(2,i+1);
    else
      x2 = v(1,1);
      y2 = v(2,1);
    end

    if ( ( y1 < p(2,1) & p(2,1) <= y2 ) | ( p(2,1) <= y1 & y2 < p(2,1) ) )
      if ( ( p(1,1) - x1 ) - ( p(2,1) - y1 ) * ( x2 - x1 ) / ( y2 - y1 ) < 0.0 )
        inside = ~inside;
      end
    end

  end

  return
end
