%ROIMANIP - manipulate shape and size of a region of interest mask.
%
%   [BWout, xi, yi] = roimanip(BWin,'method',par)
%
%   BWout = roimanip(BWin) modifies the mask image BWin by adding a 10
%   pixel wide border around the existing regions in BWin and returns the
%   resulting mask image in BWout.
%
%   BWout = roimanip(BWin,'method') modifies the mask image BWin according
%   to <method> and returns the result in BWout. For valid methods, see
%   below.
%
%   BWout = roimanip(BWin,'method', par) modifies the mask image BWin
%   according to <method> using the paramter <par>. For details, see below.
%
%   [BWout, xi, yi] = roimanip(...) also returns the bounding polygon
%   coordinates of the regions in BWout. If only one bounding polygon is
%   present, xi and yi are the vectors containing those coordinates. For
%   several region boundaries (segmented images), xi and yi are n-by-1 cell
%   arrays where each cell contains the coordinates for one of the n
%   bounding polygons.
%
%   Input arguments:
%   ----------------
%
%   BWin:   the mask image to be modified. You must convert the binary
%           image BWin into a label matrix before calling roimanip. 
%           There are two common ways to convert a binary image into a
%           label matrix:
%               a) Using the bwlabel function
%                   L = bwlabel(BW);
%               b) Using the double function
%                   L = double(BW);
%           Note, however, that these functions produce different but
%           equally valid label matrices from the same binary image. 
%           For example, given the following binary image, BWin,
%
%           BWin = 
%                   0  0  0  0  0  0  0  0
%                   0  1  1  0  0  0  0  0
%                   1  1  1  1  0  0  0  0
%                   1  0  1  0  0  0  0  0
%                   1  0  1  0  0  1  1  0
%                   0  0  0  0  0  1  1  0
%                   0  0  0  0  0  1  1  0
%                   0  0  0  0  0  0  0  0
%
%           the two following label matrices are obtained:
%
%           BWD = double(BWin) =                BWL = bwlabel(BWin) =
%                   0  0  0  0  0  0  0  0          0  0  0  0  0  0  0  0
%                   0  1  1  0  0  0  0  0          0  1  1  0  0  0  0  0
%                   1  1  1  1  0  0  0  0          1  1  1  1  0  0  0  0
%                   1  0  1  0  0  0  0  0          1  0  1  0  0  0  0  0
%                   1  0  1  0  0  1  1  0          1  0  1  0  0  2  2  0
%                   0  0  0  0  0  1  1  0          0  0  0  0  0  2  2  0
%                   0  0  0  0  0  1  1  0          0  0  0  0  0  2  2  0
%                   0  0  0  0  0  0  0  0          0  0  0  0  0  0  0  0
%
%           The output of double(BWin) is interpreted as one single
%           (segmented) region, while the result of bwlabel(BWin) is
%           interpreted as an image containing multiple regions which are
%           to be processed individually.
%
%   method: Can be one of the following:
%
%           'boundingBox' -     finds the tightest bounding box around each
%                               region in BWin and scales it by a factor
%                               <par> with respect to its center of mass
%                               coordinate.
%
%           'convexHull' -      finds the smallest convex enclosing polygon
%                               around each region in BWin and scales it
%                               by a factor <par> with respect to its
%                               center of mass coordinate.
%
%           'enclosingEllipse' - finds the smallest enclosing ellipse
%                               containing all the points for each region
%                               in BWin and scales it by a factor <par>
%                               with respect to its center of mass
%                               coordinate.
%
%           'border' -          moves the border of each region outwards by
%                               <par> pixels. Negative values of <par> move
%                               the border inwards. No scaling is applied.
%                               Sharp corners are rounded by this
%                               operation.
%
%           'scale' -           scales each region in BWin by a factor
%                               <par> with respect to its center of mass
%                               coordinate. <par> must be positive, values
%                               smaller than 1 shrink the regions.
%
%           'rotate' -          rotates each region in BWin anti-clockwise
%                               by an angle <par> with respect to its
%                               center of mass. <par> must be in degrees.
%
%           For the example image BWin from above, note the different
%           output of a 'convexHull' operation for the two label matrix
%           types:
%
%           roimanip(BWD,'convexHull') =    roimanip(BWL,'convexHull')=
%                   0  1  1  0  0  0  0  0          0  0  0  0  0  0  0  0
%                   1  1  1  1  0  0  0  0          1  1  1  1  0  0  0  0
%                   1  1  1  1  1  0  0  0          1  1  1  1  0  0  0  0
%                   1  1  1  1  1  1  1  0          1  1  1  1  0  0  0  0
%                   1  1  1  1  1  1  1  1          1  1  1  0  0  1  1  0
%                   0  1  1  1  1  1  1  1          0  0  0  0  0  1  1  0
%                   0  0  0  1  1  1  1  1          0  0  0  0  0  1  1  0
%                   0  0  0  0  0  0  1  0          0  0  0  0  0  0  0  0
%
%
%   par:    parameter which controls the manipulation process.
%           For <method> = 'border', <par> specifies the number of pixels
%           with which the border is moved outwards (inwards for negative
%           values of <par>). Default = 10.
%           For <method> = 'rotate', <par> specifies the anti-clockwise
%           rotation-angle in degrees. Default = 5.
%           For all other methods, <par> controls the amount of scaling
%           which is applied to the regions with respect to their center of
%           mass coordinate.
%           If <method> = 'scale', the default value of <par> is 1.2. For
%           <method> = 'boundingBox', 'convexHull', and 'enclosingEllipse',
%           no scaling is applied by default (<par> = 1).
%
%
%   Output arguments:
%   -----------------
%
%   BWout:  Binary mask image (of type logical, not a label matrix) with
%           the modified region of interest.
%
%   xi, yi: Coordinates of the bounding polygons of the regions in BWout.
%           If more than one distinct bounding polygon is present (for
%           example in segmented images, or regions containing holes), xi
%           and yi are cell arrays of size n-by-1 where each cell contains
%           the coordinates for one of the n bounding polygons. For only
%           one polygon, xi and yi are directly returned as standard
%           vectors (1-by-k arrays).
%           To obtain the i-th bounding polygon coordinates, use the
%           following syntax,
%               x = xi{i}, y = yi{i}
%           which returns the two vectors x and y of size 1-by-k, where k
%           is the number of vertices of the i-th bounding polygon.
%
%
%   See also:
%   ---------
%
%   REGIONPROPS (details about the label matrix conversion),
%   BWLABEL, IMDILATE (for the 'border' method), BWBOUNDARIES (details
%   about image objects and their children, i.e., holes and enclosed
%   objects).
%

%%
%==========================================================================
% Further info for function ROIMANIP
% ----------------------------------
%
% Function subroutines:
% --------------------
% parse_inputs        - parse the input paramters to the main function.
%
% findBoundingPolygon - find the bounding polygons in a BW-image. Used to
%                       restrict the bounding polygon coordinates to the
%                       image frame after a scaling operation.
%
% contourc2Lines      - convert contour matrix C (see MATLAB command
%                       'contourc') to cell array containing the vertices
%                       of each contour line in a separate cell.
%
% MinVolEllipse       - finds the minimum volume enclsing ellipsoid (MVEE)
%                       of a set of data points
%
%
% CVS-Info:
% ---------
% $Date: 2006/11/09 15:56:01 $
% $Author: schlepuetz $
% $Revision: 1.5 $
% $Source: /import/cvs/X/PILATUS/App/lib/X_PILATUS_Matlab/roimanip.m,v $
% $Tag: $
%
% Author(s):            C. M. Schlepuetz (CS)
% Co-author(s):         D. Martoccia (DM)
%                       S. A. Pauli (SP)
% Address:              Surface Diffraction Station
%                       Materials Science Beamline X04SA
%                       Swiss Light Source (SLS)
%                       Paul Scherrer Institut
%                       CH - 5232 Villigen PSI
% Created:              2006/06/02
%
% Change Log:
% -----------
% 2006/06/02 (CS):
% - created first version of this file.
%
% 2006/09/12 (CS):
% - changed default values of <par> to 1 for <method> = 'boundingBox',
%   'convexHull', and 'enclosingEllipse'.
%
% 2006/10/23 (CS):
% - fixed bug in 'scale' method which caused the returned BWout to have the
%   wrong size when removing the padded border.
%
% 2006/10/25 (SP,CS):
% - added method 'rotate'.
%
% 2006/11/09 (CS):
% - fixed bug in method 'border': When passing a bwlabel matrix, the output
%   was also a bwlabel matrix (values other than 0 and 1 possible).

%%
%==========================================================================
%  Main function - ROIMANIP
%                  ========

function varargout = roimanip(varargin)

% check number of output arguments
error(nargoutchk(1, 3, nargout))

% parse input arguments
[BWin, method, parms] = parse_inputs(varargin{:});

% check number of regions contained in label matrix BWin.
nregions = max(BWin(:));

% find dimensions of BWin
dim = size(BWin);
xdim = dim(1);
ydim = dim(2);

%%
%-------------------------------------------------------------
% perform region of interest manipulation according to method.
switch method

    %-----------------
    case 'boundingbox'
        BWout = zeros(dim);
        props = regionprops(BWin,'boundingBox','centroid');
        
        for i = 1:nregions
            bb = props(i,1).BoundingBox;
            cen = props(i,1).Centroid;
            xi = [bb(1) bb(1)+bb(3) bb(1)+bb(3) bb(1)];
            yi = [bb(2) bb(2) bb(2)+bb(4) bb(2)+bb(4)];

            % scale xi and yi with respect to the center of mass
            xi = cen(1) + (xi-cen(1)).*parms;
            yi = cen(2) + (yi-cen(2)).*parms;

            BWout = BWout | poly2mask(xi,yi,xdim,ydim);
        end


    %----------------
    case 'convexhull'
        BWout = zeros(dim);
        props = regionprops(BWin,'convexHull','centroid');
        
        for i = 1:nregions
            cen = props(i,1).Centroid;
            xi = props(i,1).ConvexHull(:,1);
            yi = props(i,1).ConvexHull(:,2);

            % scale xi and yi with respect to the center of mass
            xi = cen(1) + (xi-cen(1)).*parms;
            yi = cen(2) + (yi-cen(2)).*parms;

            BWout = BWout | poly2mask(xi,yi,xdim,ydim);

        end

        
    %-------------
    case 'enclosingellipse'

    % approach: determine convexHull of region. Using the corner points of
    % this hull is enough to determine the ellipse and results in a huge
    % reduction of data points.

        BWout = zeros(size(BWin));
        props = regionprops(BWin,'convexHull','Centroid');
        
        for i = 1:nregions
            xp = props(i,1).ConvexHull(:,1);
            yp = props(i,1).ConvexHull(:,2);
            points = [xp yp]';
            [xi,yi] = MinVolEllipse(points, 0.001);

            % scale xi and yi with respect to the center of mass
            cen = props(i,1).Centroid;
            xi = cen(1) + (xi-cen(1)).*parms;
            yi = cen(2) + (yi-cen(2)).*parms;

            BWout = BWout | poly2mask(xi,yi,xdim,ydim);

        end
    
    %------------
    case 'border'

        if parms > 0
            s = strel('disk',parms);
            BWout = imdilate(BWin,s);
        elseif parms < 0
            s = strel('disk',-parms);
            BWout = imerode(BWin,s);
        else
            eid = sprintf('Images:%s:InvalidParameter',mfilename);
            error(eid,'%s %s','''parms'' cannot be 0 for method',...
                '''border'' in call to ROIMANIP');
        end

        BWout = BWout > 0;

    %-----------
    case 'scale'

        % pad BWin with border of zeros to avoid that a roi touches the
        % border.
        BWin = padarray(BWin,[1 1]);
        BWout = zeros(size(BWin));

        % B=boundaries, L=label matrix, N=number of objects,
        % A=adjacency matrix
        [B,L,N,A] = bwboundaries(BWin);
        size(B);
        props = regionprops(L,'centroid');
        xi = cell(N+sum(A(:)),1);
        yi = cell(N+sum(A(:)),1);
        for k=1:N

            C = contourc(double(imfill(L==k,'holes')),1);
            [xi{k} yi{k}] = contourc2Lines(C);
            xi{k} = cell2mat(xi{k});
            yi{k} = cell2mat(yi{k});

            % scale xi and yi OUTWARDS with respect to the center of mass
            cen = props(k,1).Centroid;
            xi{k} = cen(1) + (xi{k}-cen(1)).*parms;
            yi{k} = cen(2) + (yi{k}-cen(2)).*parms;

            BWout = BWout | poly2mask(xi{k},yi{k},xdim+2,ydim+2);

            for l=find(A(:,k))' % the corresponding hole boundaries

                C = contourc(double(imfill(L==l,'holes')),1);
                [xi{l} yi{l}] = contourc2Lines(C);
                xi{l} = cell2mat(xi{l});
                yi{l} = cell2mat(yi{l});

                % scale xi and yi INWARDS by parms with respect to
                % the center of mass.
                cen = props(l,1).Centroid;
                xi{l} = cen(1) + (xi{l}-cen(1))./parms;
                yi{l} = cen(2) + (yi{l}-cen(2))./parms;

                BWout = BWout & ~poly2mask(xi{l},yi{l},xdim+2,ydim+2);

            end
        end

        % remove padded border from BWout
        BWout = BWout(2:xdim+1, 2:ydim+1);
        

    %-----------
    case 'rotate'
 
        % pad BWin with border of zeros to avoid that a roi touches the
        % border.
        BWin = padarray(BWin,[1 1]);
        BWout = zeros(size(BWin));
 
        % B=boundaries, L=label matrix, N=number of objects,
        % A=adjacency matrix
        [B,L,N,A] = bwboundaries(BWin);
        size(B);
        props = regionprops(L,'centroid');
        xi = cell(N+sum(A(:)),1);
        yi = cell(N+sum(A(:)),1);
        rotang = parms*pi/180;
        for k=1:N
 
            C = contourc(double(imfill(L==k,'holes')),1);
            [xi{k} yi{k}] = contourc2Lines(C);
            xi{k} = cell2mat(xi{k});
            yi{k} = cell2mat(yi{k});
 
            % rotate xi and yi with respect to the center of mass
            cen = props(k,1).Centroid;
            % store old values to xi_old and yi_old
            xi_old = xi{k};
            yi_old = yi{k};
            xi{k} = cen(1) + cos(rotang)*(xi_old-cen(1)) + ...
                sin(rotang)*(yi_old-cen(2));
            yi{k} = cen(2) - sin(rotang)*(xi_old-cen(1)) + ...
                cos(rotang)*(yi_old-cen(2));
 
            BWout = BWout | poly2mask(xi{k},yi{k},xdim+2,ydim+2);
 
            for l=find(A(:,k))' % the corresponding hole boundaries
 
                C = contourc(double(imfill(L==l,'holes')),1);
                [xi{l} yi{l}] = contourc2Lines(C);
                xi{l} = cell2mat(xi{l});
                yi{l} = cell2mat(yi{l});
 
                % rotate xi and yi by rotang with respect to
                % the center of mass.
                cen = props(l,1).Centroid;
                % store old values to xi_old and yi_old
                xi_old = xi{l};
                yi_old = yi{l};
                xi{l} = cen(1) + cos(rotang)*(xi_old-cen(1)) + ...
                    sin(rotang)*(yi_old-cen(2));
                yi{l} = cen(2) - sin(rotang)*(xi_old-cen(1)) + ...
                    cos(rotang)*(yi_old-cen(2));
 
                BWout = BWout & ~poly2mask(xi{l},yi{l},xdim+2,ydim+2);
 
            end
        end
 
        % remove padded border from BWout
        BWout = BWout(2:xdim+1, 2:ydim+1);
end

%%
%-------------------------
% assign output parameters

varargout{1} = BWout;

if nargout > 1

    [xi yi] = findBoundingPolygon(BWout);

    % convert xi and yi to 1-dimensional array instead of cell array if
    % only one border line is present.
    if length(xi) == 1
        xi = xi{1};
        yi = yi{1};
    end

    varargout{2} = xi;
    varargout{3} = yi;

end

%%
%==========================================================================
%  Sub-function - parse_inputs
%  ===========================
%
%  used to parse the input paramters to the main function.

function [BWin, method, parms] = parse_inputs(varargin)

error(nargchk(1, 3, nargin))

switch nargin

    case 1
        BWin = varargin{1};
        method = 'border';
        parms = 10;

    case 2
        BWin = varargin{1};
        method = varargin{2};
        parms = [];

    case 3
        BWin = varargin{1};
        method = varargin{2};
        parms = varargin{3};
end

% check whether BWin contains at least one region
if max(max(bwlabel(BWin)))<1
    eid = sprintf('Images:%s:NoDefinedRegion',mfilename);
    error(eid,'%s','BWin does not contain any regions in call to ROIMANIP');
end


% check whether method is a valid manipulation method.
% valid methods: boundingBox, boundingPolygon, convexHull, border,
% and enclosingellipse.
method = lower(method);
if ~(strcmp(method,'boundingbox') || strcmp(method,'enclosingellipse') ||...
        strcmp(method,'convexhull') || strcmp(method,'border') || ...
        strcmp(method,'scale') || strcmp(method,'rotate'))
    eid = sprintf('Images:%s:invalidMethod',mfilename);
    error(eid,'%s','Invalid string for ''method'' in ROIMANIP');
end

% check whether parms has legal value according to method
if (strcmp(method,'border'))
    if parms ~= round(parms)
        wid = sprintf('Images:%s:NonIntegerValue',mfilename);
        warning(wid,'%s %s %s','borders can only be moved by an',...
            'integer number of pixels when using method ''border''',...
            'in ROIMANIP. Rounding input to nearest integer.');
        parms = round(parms);
    end
end

% if parms is empty, use default parameter according to method:
% expand borders by 10 pixels for method border and scale with factor 1.2
% for all other methods.
if isempty(parms)
    if (strcmp(method,'border'))
        parms = 10;
    elseif (strcmp(method,'scale'))
        parms = 1.2;
    elseif (strcmp(method,'rotate'))
        parms = 5;
    else 
        parms = 1.0;
    end
end

%%
%==========================================================================
%  Sub-function - findBoundingPolygon
%  ==================================
%
% find the bounding polygons in a BW-image. Used to restrict the bounding
% polygon coordinates to the image frame after a scaling operation.

function [x,y] = findBoundingPolygon(BW)

% pad with zeros to have well defined border
BW = padarray(BW,[1 1]);

% find contour lines
C = contourc(double(BW),1);
[x y] = contourc2Lines(C);

% shift contourlines back to original position before the padding occured
for i = 1:length(x)
    x{i} = x{i}-1;
    y{i} = y{i}-1;
end


%%
%==========================================================================
%  Sub-function - contourc2Lines
%  =============================
%
% convert contour matrix C (see MATLAB command 'contourc') to cell
% array containing the vertices of each contour line in a separate cell.

function [xi,yi] = contourc2Lines(C)

dim = size(C);
cols = dim(2);

% find number of contour lines
ind = 1;
nlines = 0;
while ind < cols
    nvertices = C(2,ind);
    ind = ind+nvertices+1;
    nlines = nlines+1;
end

xi = cell(nlines,1);
yi = cell(nlines,1);

ind = 1;
i = 1;
while ind < cols
    nvertices = C(2,ind);
    xi{i} = C(1,ind+1:ind+nvertices);
    yi{i} = C(2,ind+1:ind+nvertices);
    ind = ind+nvertices+1;
    i = i+1;
end


%%
%==========================================================================
%  Sub-function - MINVOLELLIPSE
%  ============================
%
% [X , Y] = MinVolEllipse(P, tolerance)
% Finds the minimum volume enclsing ellipsoid (MVEE) of a set of data
% points stored in matrix P. The following optimization problem is solved: 
%
% minimize       log(det(A))
% subject to     (P_i - c)' * A * (P_i - c) <= 1
%                
% in variables A and c, where P_i is the i-th column of the matrix P. 
% The solver is based on Khachiyan Algorithm, and the final solution 
% is different from the optimal value by the pre-spesified amount of 'tolerance'.
%
% inputs:
%---------
% P : (d x N) dimnesional matrix containing N points in R^d.
% tolerance : error in the solution with respect to the optimal value.
%
% outputs:
%---------
% X,Y: x- and y-coordinates of the MVEE as a polygon with 720 points.
% 
% example:
% --------
%      P = rand(5,100);
%      [X,Y] = MinVolEllipse(P, .01)
%
%
% this subfunction is based upon the MinVolEllipse routine obtained from
% the MATLAB exchange file server and originally posted by:
%
% Nima Moshtagh (nima@seas.upenn.edu)
% University of Pennsylvania
% December 2005

function [X , Y] = MinVolEllipse(P, tolerance)

% data points 
[d N] = size(P);

Q = zeros(d+1,N);
Q(1:d,:) = P(1:d,1:N);
Q(d+1,:) = ones(1,N);

% initializations
count = 1;
err = 1;
u = (1/N) * ones(N,1);          % 1st iteration

% Khachiyan Algorithm
while err > tolerance,
    X = Q * diag(u) * Q';       % X = \sum_i ( u_i * q_i * q_i')  is a (d+1)x(d+1) matrix
    M = diag(Q' * inv(X) * Q);  % M the diagonal vector of an NxN matrix
    [maximum j] = max(M);
    step_size = (maximum - d -1)/((d+1)*(maximum-1));
    new_u = (1 - step_size)*u ;
    new_u(j) = new_u(j) + step_size;
    count = count + 1;
    err = norm(new_u - u);
    u = new_u;
end

% Find the ellipse equation in the 'center form': 
% (x-c)' * A * (x-c) = 1
% It computes a dxd matrix 'A' and a d dimensional vector 'c' as the center
% of the ellipse. 
U = diag(u);

% the A matrix for the ellipse
A = (1/d) * inv(P * U * P' - (P * u)*(P*u)' );

% center of the ellipse 
c = P * u;

% eigenvalues and eigenvectors of A:
[L,V] = eig(A);

% axes a and b
a = sqrt(1/V(1,1));
b = sqrt(1/V(2,2));

% angle between major axis and coordinate system.
phi = acos([-1 0]*L(:,1));

theta = 0:2*pi/720:2*pi;

% Parametric equation of the ellipse
%----------------------------------------
 x = a*cos(theta);
 y = b*sin(theta);

% Coordinate transform 
%----------------------------------------
 X = cos(phi)*x - sin(phi)*y;
 Y = sin(phi)*x + cos(phi)*y;
 X = X + c(1);
 Y = Y + c(2);
  
%%
%==========================================================================
%
%---------------------------------------------------%
% emacs setup:  force text mode to get no help with %
%               indentation and force use of spaces %
%               when tabbing.                       %
% Local Variables:                                  %
% mode:text                                         %
% indent-tabs-mode:nil                              %
% End:                                              %
%---------------------------------------------------%
%
% $Log: roimanip.m,v $
% Revision 1.5  2006/11/09 15:56:01  schlepuetz
% fixed bug in method border
%
% Revision 1.4  2006/10/25 15:50:05  schlepuetz
% added method rotate
%
% Revision 1.3  2006/10/24 08:22:44  schlepuetz
% fixed bug in 'scale' method which caused the returned BWout to have the wrong size when removing the padded border
%
% Revision 1.2  2006/10/20 07:37:43  schlepuetz
% cleanup, cross-checked by co-author
%
% Revision 1.1  2006/09/13 15:07:53  schlepuetz
% first release of this file
%
%
%============================= End of $RCSfile: roimanip.m,v $ ===