function pol2thd(fileIn,fileOut,grd,varargin)
%POL2THD Convert .pol to .thd file and vice versa
%
%Syntax:
%	pol2thd(fileIn,fileOut,grd,<keyword>,<value>)
%
%Input:
%	fileIn:		[string] Filepath of input .thd,.pol or .ldb file.
%	fileOut:	[string] Filepath of output .thd,.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, pol2bnd, delft3d_io_thd

%   --------------------------------------------------------------------
%   Copyright (C) 2013 Arcadis
%       Ivo Pasmans
%
%       ivo.pasmans@arcadis.nl	
%
%		Arcadis Zwolle
%
%   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: pol2thd.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/pol2thd.m $

%% Preprocessing

%keywords
OPT.dstep=1e2; %numbers of points per dam for conversion .pol->.thd 
OPT.Rmax=5e3; %maximum distance between points on polygon and grid points
OPT=setproperty(OPT,varargin);

%Check existence
if ~exist(fileIn,'file')
    error(sprintf('Matlab cannot find %s',fileIn)); 
else
    [fpath fname fext]=fileparts(fileIn); 
end %end if exist



%% Conversion .thd to .pol

if strcmpi(fext,'.thd')
    
	%open .thd file
    thd=delft3d_io_thd('read',fileIn,grd); 
    
	%read all dams in .thd file
    icount=1; 
    Field=[]; 
    for k=1:size(thd.m,2)
       
	    %determine direction of dam
        if thd.DATA(k).direction=='V'
            shift=[0 1]; 
        else
            shift=[1 0]; 
        end
        
		%find cartesian coordinates dam
        x1=[grd.cor.x(thd.n(1,k)-shift(1),thd.m(1,k)-shift(2)),grd.cor.x(thd.n(2,k),thd.m(2,k))];
        y1=[grd.cor.y(thd.n(1,k)-shift(1),thd.m(1,k)-shift(2)),grd.cor.y(thd.n(2,k),thd.m(2,k))];
		
        
		%write cartesian coordinates to file. Only break dams if end point not equal to start point. 
        if icount==1 
            Field1=struct('Name',sprintf('thd%-.3d',icount),'Data',[x1',y1']); 
            icount=icount+1; 
        elseif (Field1.Data(end,1)==x1(1) && Field1.Data(end,2)==y1(1))
            Field1.Data=[Field1.Data; [x1(2),y1(2)]]; 
        else
            Field=[Field,Field1]; 
            Field1=struct('Name',sprintf('thd%-.3d',icount),'Data',[x1',y1']); 
            icount=icount+1; 
        end
            
        
    end
    
	%Write cartesian coordinates to output
    if ~isempty(Field)
        INFO.Field=Field;
        tekal('write',fileOut,INFO);
	else
		error('No data for output'); 
    end %end if isempty(Field)
    

%% Conversion .pol to .thd

elseif strcmpi(fext,'.pol') | strcmpi(fext,'.ldb')

	%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

	
	%open input
	pol=tekal('read',fileIn); 
	pol=pol.Field; 
	mpol=[]; npol=[]; 
	
	DATA=[]; hold on
	display('Finding grid indices.');
	for k=1:length(pol)
		
		xpol=pol(k).Data(:,1);
		ypol=pol(k).Data(:,2);
		
		%removing points outside grid
		for l=1:length(enc)
			if isempty(xpol)
				break; 
			end
			inpoly=inpolygon(xpol,ypol,enc(l).x,enc(l).y); 
			if l==1
				%external enclosure
				xpol=xpol(inpoly); ypol=ypol(inpoly); 
			else
				%internal enclosure
				xpol=xpol(~inpoly); ypol=ypol(~inpoly); 
			end
		end
		
		if length(xpol)<2
			continue; 
		end
			
		
		if xpol(1)==xpol(end) & ypol(1)==ypol(end)
			xpol=xpol(1:end-1); 
			ypol=ypol(1:end-1); 
		end
		
		%split polygon at npoint positions
		[xpol ypol]=local_polyInterp(xpol,ypol,OPT.dstep);
		
		%Get m,n coordinates
		data1=local_xy2mn(grd,xpol,ypol,OPT.Rmax); 
		
		DATA=[DATA,data1]; 
		
	end %end for k
	
	%Remove superfluous m,n-coordinates
	%TO DO
	
	%Write to file
	thd=struct('DATA',DATA); 
	delft3d_io_thd('write',fileOut,thd); 
	
else
    error('Matlab cannot recognise file format'); 
end %end if strcmpi(fext


end %end function pol2thd

function [xout yout]=local_polyInterp(x,y,dstep)
%Splits a polygon with stepsize dstep

	%check input
	if ( sum(isnan(x))+sum(isnan(y))+sum(isinf(x))+sum(isinf(y)) )>0
		error('Polygon contains NaN values')
	end

	%calculate tangentiable coordinates input
	d= [0;cumsum( hypot(diff(x),diff(y)) )]; 

	%calculate tangentiable coordinates output
	di=[0:dstep:d(end)]; di=di';
	if di(end)~=d(end)
		di=[di; d(end)]; 
	end %end if

	%calculate cartesian coordinates output
	xout=interp1(d,x,di); 
	yout=interp1(d,y,di); 

end %end function local_polyInterp

function DATA=local_xy2mn(grd,x,y,Rmax)
%Convert cartesian coordinates to grid indices

	[ni mi]=xy2mn(grd.cor.x,grd.cor.y,x,y,'Rmax',Rmax);
	DATA=[]; 
	if isempty(mi)
		warning('Polygon lies outside grid.'); 
	else
		%Remove double points
		%[u1 u2 u3]=unique([mi,ni],'rows');
		%[u4 u5]=sort(u2);
		%mi=mi(u5,1); ni=ni(u5,1); 
		
		%Find for every pair of points the position and
		%orientation of the dam
		if length(mi)<2
			warning('Insufficient points to create thin dam.')
		else
			
			%differences between points
			dmn=[diff(mi,[],1),diff(ni,[],1)];  
			dmn=dmn'; dmn=dmn(:); 
			
			m0=mi(1); n0=ni(1); 
			for k=1:length(dmn)
				if abs(dmn(k))>1
					error('Cannot create continuous thin dam. Try using smaller dstep.'); 
				end
				if mod(k,2)==1
					%m changes
					n1=n0; m1=m0+dmn(k); 
					data1.m=max(m0,m1); 
					data1.n=n0;
					data1.direction='V';
				else
					%n changes
					n1=n0+dmn(k); m1=m0;
					data1.m=m1;
					data1.n=max(n0,n1); 
					data1.direction='U';
				end
				
				%remove double points
				if n0~=n1 | m0~=m1
					data1.mn=[data1.m,data1.n,data1.m,data1.n]; 
					DATA=[DATA,data1]; 
					n0=n1; m0=m1; 
				end
			end
			
		end %end if
		
	end %end if isempty(mi)

end %end function local_xy2mn