profileArea

PROFILEAREA  Compute the area eclosed by two profiles.

  Syntax:
    A = PROFILEAREA(X1, Y1, X2, Y2)

  Description:
    A = PROFILEAREA(X1, Y1, X2, Y2) returns the area A enclosed by consecutive 
    profiles with opposite directions in vectors X1 and Y1, and X2 and Y2.

  Notes:
    This function is a simpler rewording of a previous function by Tomeu Garau,
    called BUILDPOLYGON. He is the true glider man.

    The union of the two profiles may be a complex polygon (self-intersecting).
    Hence, the area is computed decomposing it in triangles with the function
    POLY2TRY, and adding the absolute value of the area of each triangular 
    component returned by POLYAREA.

    Profile points with invalid coordinates (NaN) are ignored when building the
    polygonal contour.

  Examples:
    x1 = [ 0 -1  1  0]
    y1 = [ 2  1 -1 -2]
    x2 = [ 0 -1  1  0]
    y2 = [-2 -1  1  2]
    figure
    hold on
    plot(x1, y1, 'b', x2, y2, 'r')
    a = profileArea(x1, y1, x2, y2)
    % POLYAREA would fail because of complex polygon:
    a = polyarea([x1(:); x2(:)], [y1(:); y2(:)])

  See also:
    POLYAREA

  Authors:
    Joan Pau Beltran  <joanpau.beltran@socib.cat>

0001 function a = profileArea(x1, y1, x2, y2)
0002 %PROFILEAREA  Compute the area eclosed by two profiles.
0003 %
0004 %  Syntax:
0005 %    A = PROFILEAREA(X1, Y1, X2, Y2)
0006 %
0007 %  Description:
0008 %    A = PROFILEAREA(X1, Y1, X2, Y2) returns the area A enclosed by consecutive
0009 %    profiles with opposite directions in vectors X1 and Y1, and X2 and Y2.
0010 %
0011 %  Notes:
0012 %    This function is a simpler rewording of a previous function by Tomeu Garau,
0013 %    called BUILDPOLYGON. He is the true glider man.
0014 %
0015 %    The union of the two profiles may be a complex polygon (self-intersecting).
0016 %    Hence, the area is computed decomposing it in triangles with the function
0017 %    POLY2TRY, and adding the absolute value of the area of each triangular
0018 %    component returned by POLYAREA.
0019 %
0020 %    Profile points with invalid coordinates (NaN) are ignored when building the
0021 %    polygonal contour.
0022 %
0023 %  Examples:
0024 %    x1 = [ 0 -1  1  0]
0025 %    y1 = [ 2  1 -1 -2]
0026 %    x2 = [ 0 -1  1  0]
0027 %    y2 = [-2 -1  1  2]
0028 %    figure
0029 %    hold on
0030 %    plot(x1, y1, 'b', x2, y2, 'r')
0031 %    a = profileArea(x1, y1, x2, y2)
0032 %    % POLYAREA would fail because of complex polygon:
0033 %    a = polyarea([x1(:); x2(:)], [y1(:); y2(:)])
0034 %
0035 %  See also:
0036 %    POLYAREA
0037 %
0038 %  Authors:
0039 %    Joan Pau Beltran  <joanpau.beltran@socib.cat>
0040 
0041 %  Copyright (C) 2013-2016
0042 %  ICTS SOCIB - Servei d'observacio i prediccio costaner de les Illes Balears
0043 %  <http://www.socib.es>
0044 %
0045 %  This program is free software: you can redistribute it and/or modify
0046 %  it under the terms of the GNU General Public License as published by
0047 %  the Free Software Foundation, either version 3 of the License, or
0048 %  (at your option) any later version.
0049 %
0050 %  This program is distributed in the hope that it will be useful,
0051 %  but WITHOUT ANY WARRANTY; without even the implied warranty of
0052 %  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0053 %  GNU General Public License for more details.
0054 %
0055 %  You should have received a copy of the GNU General Public License
0056 %  along with this program.  If not, see <http://www.gnu.org/licenses/>.
0057 
0058   error(nargchk(4, 4, nargin, 'struct'));
0059 
0060   % Join both profiles discarding invalid coordinates.
0061   % The resulting contour may be a complex polygon.
0062   % Decompose it in triangular components and sum up their areas.
0063   % We could use ISFINITE instead of ISNAN to discard all non-numerical values.
0064   % However, this may not be practical because the decomposition would omit
0065   % infinite triangles, and their contribution to the total area would be 0.
0066   xy = [x1(:) y1(:); x2(:) y2(:)];
0067   xy = xy(~any(isnan(xy), 2), :);
0068   [x, y] = poly2tri(xy(:,1), xy(:,2));
0069   a = sum(polyarea(x, y));
0070 
0071 end