function crs=pol2crs(fileIn,fileOut,grd,varargin)
%POL2CRS Convert tekal file to .crs file or vice versa.
%BETA VERSION. Not garanteed to work if polygon is outside grid. 
%
%Syntax:
%	pol2obs(fileIn,fileOut,grd,<keyword>,<value>,...)
%
%Input:
%	fileIn:		[string] Filepath of input .crs,.pol or .ldb file.
%	fileOut:	[string] Filepath of output .crs,.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: pol2dry, pol2obs, pol2thd, pol2bnd, delft3d_io_crs

%   --------------------------------------------------------------------
%   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: pol2crs.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/pol2crs.m $

%% Preprocessing

%keywords
OPT.Rmax=1e3; %maximum distance observation point grid point
OPT.dstep=1e2; %stepsize interpolation polygon
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

if strcmpi(fext,'.crs')
%% convert crs to polygon

	%read crs
	crs=delft3d_io_crs('read',fileIn); 
	
	Field=[];
	for k=1:length(crs.DATA)
		
		%name of cross section
		Field1.Name=sprintf('crs%-.3d',k); 
		Field1.Comments=strtrim(crs.DATA(k).name);

		%position of cross-section
		mn=crs.DATA(k).mn; 
		mn=reshape(mn,[2 2]); 
		%[m_A,m_B;n_A n_B]
		dmn=diff(mn,[],2);
		dmn=sign(dmn); 
		
		if dmn(1)==1 | dmn(2)==1
			mn(1,1)=mn(1,1)-dmn(1); 
			mn(2,1)=mn(2,1)-dmn(2); 
		elseif dmn(1)==-1 | dmn(2) == -1
			mn(1,2)=mn(1,2)+dmn(1); 
			mn(2,2)=mn(2,2)+dmn(2); 
		end
		
		dmn(dmn==0)=1; 
		xLine=grd.cor.x(mn(2,1):dmn(2):mn(2,2),mn(1,1):dmn(1):mn(1,2)); 
		yLine=grd.cor.y(mn(2,1):dmn(2):mn(2,2),mn(1,1):dmn(1):mn(1,2)); 
		if size(xLine,1)<size(xLine,2)
			xLine=xLine';
		end
		if size(yLine,1)<size(yLine,2)
			yLine=yLine';
		end
		Field1.Data=[xLine,yLine]; 
		
		Field=[Field,Field1]; 
	end %end for k
	
	%write to polygon
	INFO.Field=Field;
	tekal('write',fileOut,INFO);
	
elseif strcmpi(fext,'.pol') |strcmpi(fext,'.ldb')
%% convert polygon to one or more cross sections


    %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
	
	%load polygons
	pol=tekal('read',fileIn); 
	pol=pol.Field; 
	
	DATA=[]; 
	for k=1:length(pol)
		
		xLine=pol(k).Data(:,1); 
		yLine=pol(k).Data(:,2); 
		name=pol(k).Comments; 
		name=cell2mat(name); 
		
		%create not closed polygon
		if xLine(1)==xLine(end) & yLine(1)==yLine(end)
			xLine=xLine(1:end-1); yLine=yLine(1:end-1); 
		end
		
		%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
		if length(xLine)<2
			continue; 
		end
		
		[xLine yLine]=local_polyInterp(xLine,yLine,OPT.dstep); 
		
		DATA1=local_xy2mn(grd,xLine,yLine,OPT.Rmax); 
		
		if length(DATA1)==1
			DATA1.name=sprintf('%-20s',name(2:end)); 
			DATA=[DATA,DATA1];
		else
			for l=1:length(DATA1)
				lname=length(name(:)); 
				lname=min(lname,18); 
				DATA1(l).name=sprintf('%-17s %-.2d',name(2:lname),l); 
			end
			DATA=[DATA,DATA1]; 
		end
	end %end for k
	
	crs.DATA=DATA;
	
	%reshape
	mn=[crs.DATA(:).mn] ;
	mn=reshape(mn,4,[]); mn=mn';
	crs.m=mn(:,[1,3]); 
	crs.n=mn(:,[2,4]); 
	namst={crs.DATA(:).name}; 
	namst=cell2mat(namst); 
	namst=reshape(namst,20,[]); namst=namst';
	crs.namst=namst; 
    delft3d_io_crs('write',fileOut,crs); 



else
	error('Not a valid file format')
end %end if strcmpi(fext

end %end function pol2crs

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
		
		%Find for every pair of points the position and
		%orientation of the dam
		if length(mi)<2
			warning('Insufficient points to create cross section')
		else
			
			%differences between points
			dmn=[diff(mi),diff(ni)];  
			
			DATA1=struct('name',[],'mn',[mi(1) ni(1)],'m',[],'n',[]); 
			mn0=[mi(1),ni(1)]; mem=[0 0];  
			for k=1:size(dmn,1)
				
				if sum(abs(dmn(k,:)))==0 & k~=size(dmn,1)
					continue;
				end
				
				if abs(dmn(k,1))>1 | abs(dmn(k,2))>1
					dmn(k,:)
					error('Cannot create continuous cross section. Try using smaller dstep.'); 
				end
				
				mem=mem+dmn(k,:); 
				
				if abs(mem(1))>=abs(mem(2))
					%v cross section
					if ( abs( mem(2) )>1 & abs(mem(1))>1 )| k==size(dmn,1)
						DATA1=local_write_cross(mn0,mn0+[mem(1) 0]); 
						DATA=[DATA,DATA1]; 
						mem(1)=0; mn0=mn0+[mem(1) 0]; 
					end
					if isnan( grd.cor.x(mn0(2),mn0(1)+mem(1)) )
						DATA1=local_write_cross(mn0,mn0+[mem(1) 0]-[dmn(k,1) 0]); 
						DATA=[DATA,DATA1]; 
						mem(1)=0; mn0=mn0+[mem(1) 0]-[dmn(k,1) 0];
					end
				elseif abs(mem(2))>=abs(mem(1))
					%u cross section
					if ( abs( mem(1) )>1 & abs(mem(2))>1 )| k==size(dmn,1)
						DATA1=local_write_cross(mn0,mn0+[0 mem(2)]); 
						DATA=[DATA,DATA1]; 
						mem(2)=0; mn0=mn0+[0 mem(2)]; 
					end
					if isnan( grd.cor.x(mn0(2)+mem(2),mn0(1)) )
						DATA1=local_write_cross(mn0,mn0+[0 mem(2)]-[0 dmn(k,2)]);
						DATA=[DATA,DATA1]; 
						mem(2)=0; mn0=mn0+[0 mem(2)]-[0 dmn(k,2)]; 
					end
				end
					
				
					
			end
			
		end %end if
		
	end %end if isempty(mi)
end %end function local_xy2mn

function data=local_write_cross(mnA,mnB)
%write cross sections using dep indices mnA en mnB

mnA=reshape(mnA,[],1); 
mnB=reshape(mnB,[],1); 

dmn=diff([mnA,mnB],[],2); 

if sum(abs(dmn)~=0)>1 | sum(abs(dmn))<2
	error('Input does not define a proper cross section')
end
dmn=sign(dmn);

if sum(dmn)==1
	mnA=mnA+dmn;
else
	mnB=mnB-dmn; 
end

if dmn(1)~=0
	m=[mnA(1):dmn(1):mnB(1)]; 
	n=mnA(2)*ones(size(m)); 
end
if dmn(2)~=0
	n=[mnA(2):dmn(2):mnB(2)]; 
	m=mnA(1)*ones(size(n)); 
end
data=struct('name',[],'mn',[mnA',mnB'],'m',m,'n',n);


end %end function local_write_cross
