function pol2bnd(fileIn,fileOut,grd,varargin)
%POL2BND Convert tekal file to bnd and vice versa.
%BETA VERSION
%
%Syntax:
%	pol2bnd(fileIn,fileOut,grd,<keyword>,<value>)
%
%Input:
%	fileIn:		[string] Filepath of input .bnd,.pol or .ldb file.
%	fileOut:	[string] Filepath of output .bnd,.pol or .ldb file.
%	grd:		[struct] grid struct as generated by delft3d_io_grd
%
%Keyword:	
%	Rmax		[double] parts of the polygon which are further than Rmax
%				from a grid point will not be taken into account. Rmax should
%				normally be larger than maximum cell size.
%	dstep		[double] Polygon is matched to grid by dividing it in sections
%				of length dstep. Normally dstep should be 4x (or more) smaller 
%				than minimal cell size
%
%See also: pol2crs, pol2dry, pol2obs, pol2thd, delft3d_io_bnd

%   --------------------------------------------------------------------
%   Copyright (C) 2013 Arcadis
%       Ivo Pasmans
%
%       ivo.pasmans@arcadis.nl	
%
%		Arcadis, Zwolle, The Netherlands
%
%   This library is free software; you can redistribute it and/or
%   modify it under the terms of the GNU Lesser General Public
%   License as published by the Free Software Foundation; either
%   version 2.1 of the License, or (at your option) any later version.
%
%   This library is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%   Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this library; if not, write to the Free Software
%   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
%   USA
%   or http://www.gnu.org/licenses/licenses.html, http://www.gnu.org/, http://www.fsf.org/
%   --------------------------------------------------------------------

% $Id: pol2bnd.m 9598 2013-11-07 13:50:02Z boer_g $
% $Date: 2013-11-07 14:50:02 +0100 (Thu, 07 Nov 2013) $
% $Author: boer_g $
% $Revision: 9598 $
% $HeadURL: https://svn.oss.deltares.nl/repos/openearthtools/trunk/matlab/applications/delft3d/flow/pol2bnd.m $

%find settings
OPT.location=0;
OPT.d_step=50;
OPT.Rmax=10e3;
[OPT OPTset]=setproperty(OPT,varargin);

%check if input exist
if ~exist(fileIn,'file')
    error(sprintf('Cannot find file %s',fileIn)); 
else
    [fpath fname fext]=fileparts(fileIn); 
end

if strcmpi(fext,'.bnd')
    %% Transform .bnd into .pol 

    %Read data in bnd file
	bnd=delft3d_io_bnd('read',fileIn);
    bnd=bnd.DATA; 
	
    %initiate output
	Field=[]; 
    
    %initiate matrix containing points used in boundary
    u_mn=[]; v_mn=[];
    
     %Enclosure. Needed to remove points outside enclosure
     enc=polyinterp(grd.Enclosure(:,1),grd.Enclosure(:,2),'d_step',NaN,'seperator','equal');
     for l=1:length(enc)
         enc(l).m=enc(l).x; enc(l).n=enc(l).y;
         for k=1:length(enc(l).x)
             enc(l).x(k)=grd.cend.x(enc(l).n(k),enc(l).m(k));
             enc(l).y(k)=grd.cend.y(enc(l).n(k),enc(l).m(k));
         end
     end
    
    %expand grid
    corx=nan(size(grd.cor.x')+[2 2]); corx(2:end-1,2:end-1)=grd.cor.x';
    cory=nan(size(grd.cor.y')+[2 2]); cory(2:end-1,2:end-1)=grd.cor.y';
    cendx=grd.cend.x'; cendy=grd.cend.y';
    
    for l=1:length(enc)
        inpoly=inpolygon(corx,cory,enc(l).x,enc(l).y); 
        if l==1
            %external enclosure
            corx(~inpoly)=NaN; cory(~inpoly)=NaN; 
        else
            %internal enclosure
            corx(inpoly)=NaN; cory(inpoly)=NaN;
        end
    end
    
    %read boundary by boundary section
	for k=1:length(bnd)
		
        %mn=[m_A, m_B; n_A n_B]
        mn=reshape(bnd(k).mn,[2 2]); 
        %direction boundary
		dmn=diff(mn,[],2); dmn=sign(dmn);         
		
        if dmn(1)==0
            %u boundary

           

            %indices WL-points section
            clear mn_section; 
            mn_section(:,2)=[mn(2,1):dmn(2):mn(2,2)];
            mn_section(:,1)=mn(1,1);
            %remove indices outside grid
            mn_section=mn_section(mn_section(:,1)>=1 & mn_section(:,1)<=size(cendx,1)...
                & mn_section(:,2)>=1 & mn_section(:,2)<=size(cendx,2),:);
            
            %remove points that are already in other sections
            u_mn=[mn_section; u_mn;]; 
            [u_mn i_unique]=unique(u_mn,'rows');
            i_unique=i_unique(i_unique<=size(mn_section,1)); i_unique=sort(i_unique); 
            %no uniques test: mn_section=mn_section(i_unique,:);
            if size(mn_section,1)<2
                display(sprintf('Section %s is empty',bnd(k).name));
                continue; 
            end
            
            %convert n-coordinates from WL to depth point
            n_section=[mn_section(1,2)-max(0,dmn(2)):dmn(2):mn_section(end,2)+min(0,dmn(2))];
            
            %convert m-coordinates from depth point
            m_section1=corx(mn(1,1)+1,n_section+1); 
            m_section2=corx(mn(1,1),n_section+1);
            if sum(isnan(m_section1))<=sum(isnan(m_section2))
                m_section=mn(1,1); 
            else
                m_section=(mn(1,1)-1);
            end

            
            
            %convert m,n-coordinates to x,y-coordinates
            xLine=0.5*corx(m_section+1,n_section(1:end-1)+1)+...
                0.5*corx(m_section+1,n_section(2:end)+1); 
            yLine=0.5*cory(m_section+1,n_section(1:end-1)+1)+...
                0.5*cory(m_section+1,n_section(2:end)+1); 
            inan=~isnan(xLine&yLine); 
            xLine=xLine(inan); yLine=yLine(inan); 
            if length(xLine)<2
                %too few points to make boundary
                display(sprintf('Section %s is empty',bnd(k).name));
                continue; 
            else
                %convert output to column vectors
                xLine=reshape(xLine,[],1); yLine=reshape(yLine,[],1); 
            end %end if line           
            
		elseif dmn(2)==0
        %v boundary
            
            %indices WL-points section
            clear mn_section; 
            mn_section(:,1)=[mn(1,1):dmn(1):mn(1,2)];
            mn_section(:,2)=mn(2,1);
            %remove indices outside grid
            mn_section=mn_section(mn_section(:,1)>=1 & mn_section(:,1)<=size(cendx,1)...
                & mn_section(:,2)>=1 & mn_section(:,2)<=size(cendx,2),:);
            
            %remove points that are already in other sections
            v_mn=[mn_section; v_mn]; 
            [v_mn i_unique]=unique(v_mn,'rows'); 
            i_unique=i_unique(i_unique<=size(mn_section,1)); 
            %no uniqueness test: mn_section=mn_section(i_unique,:); 
            i_unique=sort(i_unique);
            if size(mn_section,1)<2
                display(sprintf('Section %s is empty',bnd(k).name));
                continue; 
            end
            
            %convert m-coordinates from WL to depth point
            m_section=[mn_section(1,1)-max(0,dmn(1)):dmn(1):mn_section(end,1)+min(0,dmn(1))];
            
            %convert n-coordinates from depth point
            n_section1=corx(m_section+1,mn(2,1)+1); 
            n_section2=corx(m_section+1,mn(2,1)); 
            if sum(isnan(n_section1))<=sum(isnan(n_section2))
                n_section=mn(2,1); 
            else
                n_section=(mn(2,1)-1);
            end
            
            %convert m,n-coordinates to x,y-coordinates
            xLine=0.5*corx(m_section(1:end-1)+1,n_section+1)+...
                0.5*corx(m_section(2:end)+1,n_section+1); 
            yLine=0.5*cory(m_section(1:end-1)+1,n_section+1)+...
               0.5*cory(m_section(2:end)+1,n_section+1);  
            inan=~isnan(xLine&yLine); 
            xLine=xLine(inan); yLine=yLine(inan); 
            if length(xLine)<2
                %too few points to make boundary
                display(sprintf('Section %s is empty',bnd(k).name));
                continue; 
            else
                %convert output to column vectors
                xLine=reshape(xLine,[],1); yLine=reshape(yLine,[],1); 
            end
            
        end %end if dmn

		Field1.Data=[xLine,yLine]; 
		Field1.Name=sprintf('%s; %s; %s; %f; %s; %s; %d',bnd(k).name,bnd(k).bndtype,...
            bnd(k).datatype,bnd(k).alfa,bnd(k).vert_profile,bnd(k).pli_name,bnd(k).pli_nr); 
		Field=[Field,Field1];
		clear Field1;
	end %end for k
	
    %write data to tekal file
	INFO.Field=Field; 
	tekal('write',fileOut,INFO); 
	
	elseif strcmpi(fext,'.ldb') | strcmpi(fext,'.pol')
        %% convert .pol or .ldb to .bnd
        
        %initiate output
        bnd=struct('DATA',[],'m',[],'n',[],'mn',[]); 
		
        %read tekal file
		pol=tekal('read',fileIn); 
		pol=pol.Field;
        
        %Enclosure. Needed to remove points outside enclosure and find find
        %points where interpolation takes place
        enc=polyinterp(grd.Enclosure(:,1),grd.Enclosure(:,2),'d_step',NaN,'seperator','equal');
        encx=[]; ency=[]; encm=[]; encn=[]; encl=[]; 
        for l=1:length(enc)
            enc(l).m=enc(l).x; enc(l).n=enc(l).y; 
            for k=1:length(enc(l).x)
                enc(l).x(k)=grd.cend.x(enc(l).n(k),enc(l).m(k)); 
                enc(l).y(k)=grd.cend.y(enc(l).n(k),enc(l).m(k)); 
                
            end
            encx=[encx;enc(l).x];
            ency=[ency;enc(l).y];
            encm=[encm;enc(l).m];
            encn=[encn;enc(l).n];
            encl=[encl;l*ones(size(enc(l).x))];
        end	
        		
        for k=1:length(pol)
            
            %get x,y-coordinates
            xLine=pol(k).Data(:,1);
            yLine=pol(k).Data(:,2);            

            %removing points outside grid
            for l=1:length(enc)
                if isempty(xLine)
                    break;
                end
                inpoly=inpolygon(xLine,yLine,enc(l).x,enc(l).y);
                if l==1
                    %external enclosure
                    xLine=xLine(inpoly); yLine=yLine(inpoly);
                else
                    %internal enclosure
                    xLine=xLine(~inpoly); yLine=yLine(~inpoly);
                end
            end %enf for l

            %check if there are sufficient points to create boundary
            if length(xLine)<2
                display(sprintf('Polygon %d does not contain sufficient points for a boundary',k)); 
                continue;
            end

            %convert x,y to points in enclosure
            [xi yi]=local_lineinterp(xLine,yLine,OPT.d_step);
            [enci_m enci_n]=xy2mn(encx,ency,xi,yi);
            enci=unique(enci_m);
            n=encn(enci); m=encm(enci);
            
            if n(1)~=n(end) & m(1)~=m(end)
                display(sprintf('Cannot write boundary %d, because its begin and end point are not on 1 line.',k));
                [m'; n']
                m=[0 0]; n=[0 0]; 
            end

            %read meta data
            meta=regexp(pol(k).Name,';','split');
            meta=strtrim(meta);
            data1.mn=[m(1),n(1),m(end),n(end)];
            data1.m=m;
            data1.n=n;
            data1.name=meta{1};
            data1.bndtype=meta{2};
            data1.datatype=meta{3};
            data1.alfa=str2double(meta{4});
            data1.vert_profile=meta{5};
            data1.pli_name=meta{6};
            data1.pli_nr=str2double(meta{7});
            bnd.m=[bnd.m; m(1), m(end)];
            bnd.n=[bnd.n; n(1), n(end)];
            bnd.DATA=[bnd.DATA,data1]; 
            

        end %end for k
		
		delft3d_io_bnd('write',fileOut,bnd); 
	end
end


function [xi yi]=local_lineinterp(x,y,dstep)
%LOCAL_LINEINTERP

%tangential coordinate
dist=distance(x,y); 
di=[dist(1):dstep:dist(end)]; 
di=unique([di,dist]); 

xi=interp1(dist,x,di); 
yi=interp1(dist,y,di); 

end %end local_lineinterp

