poi_area.m

toolbox of BVQX, This is the access between BV and matlab. It will help you to analysis data from BV

function [pas, pa] = poi_Area(hfile, srf, poin)
% POI::Area  - calcuate area of POI
%
% FORMAT:       [pas, pa] = poi.Area(srf [, poin])
%
% Input fields:
%
%       srf         surface file object (for coordinates)
%       poin        number of POI (if not given set to 1)
%
% Output fields:
%
%       pas         sum of POI area
%       pa          area of within-POI triangles

% Version:  v0.7b
% Build:    7090213
% Date:     Sep-02 2007, 1:05 PM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'poi') || ...
    numel(srf) ~= 1 || ...
   ~isBVQXfile(srf, 'srf')
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if numel(bc.POI) < 1
    error( ...
        'BVQXfile:InvalidObject', ...
        'POI object must contain at least one POI definition.' ...
    );
end
if nargin < 3 || ...
   ~isa(poin, 'double') || ...
    numel(poin) ~= 1 || ...
    isinf(poin) || ...
    isnan(poin) || ...
    poin < 1 || ...
    poin ~= fix(poin) || ...
    poin > numel(bc.POI)
    poin = 1;
end
srfc = bvqxfile_getcont(srf.L);

% get POI vertices
poiv = bc.POI(poin).Vertices(:);

% find triangles that match
tv = srfc.TriangleVertex;
for tvc = 1:numel(poiv)
    tv(tv == poiv(tvc)) = Inf;
end
ti = find(sum(isinf(tv),2) == 3);
if numel(ti) == 0
    error( ...
        'BVQXfile:InternalError', ...
        'No triangles in POI.' ...
    );
end

% get vertices of triangles and replace with numbers from 1 to numel(poiv)
tv = srfc.TriangleVertex(ti, :);
for tvc = 1:numel(poiv), tv(tv == poiv(tvc)) = tvc; end

% get coordinates
tc = srfc.VertexCoordinate(poiv, :);

% calculate sides for Heron / Archimedes formula
% (see http://www.mste.uiuc.edu/dildine/heron/triarea.html for details)
ts = [ ...
    sqrt((tc(tv(:,1),1) - tc(tv(:,2),1)) .^ 2 + ...
         (tc(tv(:,1),2) - tc(tv(:,2),2)) .^ 2 + ...
         (tc(tv(:,1),3) - tc(tv(:,2),3)) .^ 2), ...
    sqrt((tc(tv(:,2),1) - tc(tv(:,3),1)) .^ 2 + ...
         (tc(tv(:,2),2) - tc(tv(:,3),2)) .^ 2 + ...
         (tc(tv(:,2),3) - tc(tv(:,3),3)) .^ 2), ...
    sqrt((tc(tv(:,3),1) - tc(tv(:,1),1)) .^ 2 + ...
         (tc(tv(:,3),2) - tc(tv(:,1),2)) .^ 2 + ...
         (tc(tv(:,3),3) - tc(tv(:,1),3)) .^ 2)];

% use Heron / Archimedes formula
tcs = sum(ts, 2) / 2;
pa = sqrt(tcs .* (tcs-ts(:,1)) .* (tcs-ts(:,2)) .* (tcs-ts(:,3)));

% global area
pas = sum(pa);