mask.m

MDPSAS工具箱是马里兰大学开发的

function B = mask(A,shape,params)

% mask.m Creates a mask defined as a scalarfield object
%        with value 1 inside the given vertices and zero
%        elsewhere. Called as
%
%        B = mask(A,shape,params)
%
%        INPUT parameters:
%           A      : quadgrid object
%           shape  : 'lineseg'
%           params : [x0, xfinal]
%           shape  : 'disk'
%           params : [radius, xcenter, ycenter]
%                    note that xcenter = r; ycenter = theta if
%                    the quadgrid is in polar coordinates

switch shape

    case 'lineseg'

        if length(A(1).name) ~= 1, error('Must be a 1D quadgrid'), end
        if prod(size(params)) ~= 2, error('Must specify 2 vertices'), end

        qp = A.qp{1};

        for i = 1:length(qp)

            if qp(i) < params(1) | qp(i) > params(2)
                val(i,1) = 0;
            else
                val(i,1) = 1;
            end

        end

    case 'disk'

        if length(A(1).name) ~= 2, error('Must be a 2D quadgrid'), end
        if prod(size(params)) ~= 3, error('Must specify 3 parameters'), end

        x = A.qp{1};
        y = A.qp{2};

        radius  = params(1);
        xcenter = params(2);
        ycenter = params(3);

        gflag = 0;
        if A.geomcode == 11, gflag = 1; end
        if A.geomcode == 24, gflag = 2; end

        if gflag == 0, error('geometry not supported'), end

        for i = 1:length(x)
            for j = 1:length(y)

                switch gflag
                    case 1
                        if sqrt( (x(i)-xcenter)^2 + (y(j)-ycenter)^2 ) > radius
                            val(i,j) = 0;
                        else
                            val(i,j) = 1;
                        end

                    case 2
                        if sqrt( (x(i)*cos(y(j))-xcenter*cos(ycenter))^2 ...
                                + (x(i)*sin(y(j))-xcenter*sin(ycenter))^2 ) > radius
                            val(i,j) = 0;
                        else
                            val(i,j) = 1;
                        end
                end

            end
        end

end % ends cases

B = scalarfield(A,val);