/* 
* 
* Copyright (C) 2013 UMPC
* 
* This file is part of CLaRiNet.
*  
* CLaRiNet is free software: you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* 
* CLaRiNet 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 General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with CLaRiNet.  If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <stddef.h>
#include "GeoMorph.h"
#include "Correction.h"
#include "Erreur.h"
#include "esriheader.h"
#include "FlowDirGrid.h"
#include "FloatRaster.h"
#include "Topology.h"
#include<vector>
#include<map>
#include<fstream>
#include<algorithm>

using namespace std;

/* int MakeElevConsistentWithFlowDir ( const int HandleRaster ,
                                    const int HandleFlowDir ,
                                    int root )
{
    // Met en coherence les altitudes des cellules du MNT avec un plan de drainage prescris.
    // L'analyse est effectuee a partir de la cellule "root"
    
    FloatRaster * R;
    FlowDirGrid * F;

    int recurs;

    if ( !(R = GetFloatRaster(HandleRaster)) )
    {
        SetErreur("RasterFindTransition : le raster # %d n'existe pas\n",HandleRaster);
        return(-1);
    }

    if ( !(F = GetFlowDirGrid(HandleFlowDir)) )
    {
        SetErreur("RasterFindTransition : la grille d'ecoulement # %d n'existe pas\n",HandleFlowDir);
        return(-1);
    }

    if ( CheckSameHeader( &(F->Header) , &(R->Header) ) )
    {
        SetErreur("RasterFindTransition : le raster # %d et la grille d'ecoulement # %d n'ont pas la meme reference spatiale\n",HandleRaster, HandleFlowDir);
        return(-1);
    }

    // Initialisation des valeurs de retour

    recurs = 1;

    if ( MakeElevConsistentWithFlowDir_ ( R, F , & root , & recurs ) )
    {
        SetErreur("MakeElevConsistentWithFlowDir : echec du 1er appel de la recursion\n");
        return(-1);
    }

    return(0);
}


int MakeElevConsistentWithFlowDir_ ( FloatRaster * R, FlowDirGrid * F , int * _root , int * _recurs )
{

    int i,nup1,recmin;
    int nr,nc;
    float Zi , Zmin , Zroot ;
    int zero = 0;
    int * UPSTREAM ;

    // Lecture de la valeur de la cellule courante

    (void) rec2RC ( &(R->Header) , _root , & nr, & nc );
    Zroot = R->Data[nr][nc];
    Zmin = Zroot;
    recmin = (*_root);

    // Recuperation des indices des cellules amont

    if ( GetUpstream ( F , _root , & nup1 , & UPSTREAM ) )
    {
            SetErreur("MakeElevConsistentWithFlowDir_ : echec de recuperation des cellules amont\n");
            return(-1);
    }

    if(!nup1)
        return(0);

    for(i=0;i<nup1;i++)
    {

        if ( *_recurs )
        {
            if ( MakeElevConsistentWithFlowDir_ ( R, F , UPSTREAM + i , _recurs ) )
            {
                SetErreur("MakeElevConsistentWithFlowDir_ : echec de l''appel recursif\n");
                return(-1);
            }
        }

        (void) rec2RC ( &(R->Header) , UPSTREAM + i , & nr, & nc );
        Zi = R->Data[nr][nc];

        if ( Zi < Zmin )
        {
            recmin = UPSTREAM[i];
            Zmin = Zi;
        }

    }

    if ( recmin!=(*_root) )       // Si le point bas du bv n'est pas la racine de l'appel
    {
        (void) rec2RC ( &(R->Header) , _root , & nr, & nc );    // On echange les altitude des cellules root et recmin
        R->Data[nr][nc] = Zmin;
        (void) rec2RC ( &(R->Header) , & recmin , & nr, & nc );
        R->Data[nr][nc] = Zroot;

        // On rappelle la fonction en recmin (puisqu'on a remonte son altitude, cette cellule n'est peut-etre plus le point bas de son bv)

        if ( MakeElevConsistentWithFlowDir_ ( R, F , & recmin , & zero ) )
        {
                SetErreur("MakeElevConsistentWithFlowDir_ : echec de l''appel correctif\n");
                return(-1);
        }

    }

    return(0);

}
*/

//==============================================================================

int MakeLinkTree ( LinkTree *** _ptrLTree , char strFichGeom[120], char strFichLinks[120], char strFichNodes[120] , double bb[4] , int * _nlinks )
{
        double x,y,scal,dxgeom,dxtopo,dygeom,dytopo,xstart,xend,ystart,yend;
        float f,LAm ;
        complex<double> z;
        complex<double> I(0,1);
        bool inbb;
        int i , id , hylcov , hylcov0 , id_Am , id_Av , nlinks , uplink , nchildren , npoints , nodesinbb ;
        map<int, pair<complex<double>,bool> > Nodes ;
        map<int, pair<complex<double>,bool> > ::iterator itnodes ;
        multimap<int,int> EndsAtNode ;
        multimap<int,int>::iterator immii;
        pair<multimap<int,int>::iterator,multimap<int,int>::iterator> ret;
        map<int,int> MapLinkId ;
        map<int,int>::iterator imii;
        LinkTree ** _LTree ;

        // Verification de l'existence des 3 fichiers : noeuds, topologie, et geometrie des troncons

        ifstream FichNodes ( strFichNodes , ios::in);

        if (FichNodes.fail())
        {
            SetErreur("MakeLinkTree : could not open file containing nodes (%s)\n",strFichNodes);
            return(-1);
        }

        ifstream FichLinks (strFichLinks,ios::in);

        if (FichLinks.fail())
        {
            SetErreur("MakeLinkTree : could not open file containing link topology (%s)\n",strFichLinks);
            return(-1);
        }

        ifstream FichGeom (strFichGeom,ios::in);

        if (FichGeom.fail())
        {
            SetErreur("MakeLinkTree : could not open file containing link geometry (%s)\n",strFichGeom);
            return(-1);
        }

        // Chargement des noeuds -----------------------------------------------

        Nodes.clear();

	while (FichNodes >> id >> x >> y)
	{
                z = x + I*y;
                inbb = ( x>=bb[0] ) & ( x<=bb[1] ) & ( y>=bb[2] ) & ( y<=bb[3] ) ;
		Nodes.insert ( make_pair (id , make_pair(z,inbb) ) );
	}

	FichNodes.close();

	// Chargement des troncons ---------------------------------------------

        _LTree = NULL ;
        nlinks = 0;

	while (FichLinks >> hylcov >> id_Am >> id_Av >> LAm)
	{
            nodesinbb = 0;

            // Identification du noeud amont

            itnodes = Nodes.find(id_Am);

            if ( itnodes == Nodes.end() )
            {
                SetErreur("MakeLinkTree : node %d, which is the start node of link %d, is not listed in file %s\n",id_Am,hylcov,strFichLinks);
                return(-1);
            }
            xstart = (*itnodes).second.first.real();
            ystart = (*itnodes).second.first.imag();

            if ((*itnodes).second.second)
                nodesinbb++ ;

            // Identification du noeud aval

            itnodes = Nodes.find(id_Av);

            if ( itnodes == Nodes.end() )
            {
                SetErreur("MakeLinkTree : node %d, which is the end node of link %d, is not listed in file %s\n",id_Av,hylcov,strFichLinks);
                return(-1);
            }
            xend = (*itnodes).second.first.real();
            yend = (*itnodes).second.first.imag();

            if ((*itnodes).second.second)
                nodesinbb++;

            if ( nodesinbb ) // Si au moins l'une des extrÃ©mitÃ©s du troncon est dans le domaine
            {
                    nlinks++;

                    if ( ! (_LTree = (LinkTree **) realloc ( _LTree , nlinks * sizeof ( LinkTree *) ) ) )
                    {
                        SetErreur("MakeLinkTree : could not reallocate pointer _LTree (size %d)\n",nlinks * sizeof ( LinkTree *) );
                        return(-1);
                    }

                    if ( ! (_LTree[nlinks-1] = (LinkTree * ) malloc ( sizeof (LinkTree) ) ) )
                    {
                        SetErreur("MakeLinkTree : could not reallocate pointer _LTree[%d] (size %d)\n",nlinks,sizeof ( LinkTree ) );
                        return(-1);
                    }

                    _LTree[nlinks-1]->id = hylcov ;

                    _LTree[nlinks-1]->startnode = id_Am ;
                    _LTree[nlinks-1]->xstart = xstart ;
                    _LTree[nlinks-1]->ystart = ystart ;

//                    _LTree[nlinks-1]->endnode = id_Av ;
                    _LTree[nlinks-1]->xend = xend ;
                    _LTree[nlinks-1]->yend = yend ;

                    (_LTree[nlinks-1]->Points)[0] = NULL ;
                    (_LTree[nlinks-1]->Points)[1] = NULL ;
                    _LTree[nlinks-1]->Child = NULL ;
                    _LTree[nlinks-1]->Parent = NULL ;

                    EndsAtNode.insert(make_pair(id_Av,nlinks-1));
                    MapLinkId.insert(make_pair(hylcov,nlinks-1));
            }
	}

	FichLinks.close();
        Nodes.clear();
        
        // Lecture de la geometrie ---------------------------------------------

        hylcov0=-1;

	while (FichGeom >> x >> y >> f)
	{
		hylcov=(int)f;

                if(hylcov!=hylcov0)
		{
                    imii = MapLinkId.find(hylcov);
                    
                    if ( imii!=MapLinkId.end() )
                    {
                        i = (*imii).second;
                    }
                    else
                    {
                        i=-1;
                    }
                    
                    hylcov0=hylcov;
                    npoints = 0;
		}

                if (i>=0)
                {
                    npoints++;
                    _LTree[i]->npoints = npoints ;

                    if ( !( _LTree[i]->Points[0] = (double *) realloc ( _LTree[i]->Points[0] , npoints * sizeof ( double ) ) ) )
                    {
                        SetErreur("MakeLinkTree : could not reallocate pointer _LTree[%d]->Points[0] (size %d)\n",i,npoints);
                        return(-1);                        
                    }

                    if ( !( _LTree[i]->Points[1] = (double *) realloc ( _LTree[i]->Points[1] , npoints * sizeof ( double ) ) ) )
                    {
                        SetErreur("MakeLinkTree : could not reallocate pointer _LTree[%d]->Points[1] (size %d)\n",i,npoints);
                        return(-1);
                    }

                    (_LTree[i]->Points)[0][npoints-1] = x;
                    (_LTree[i]->Points)[1][npoints-1] = y;
                }
	}

	FichGeom.close();

        // Creation de la topologie --------------------------------------------

        for (i=0;i<nlinks;i++)
        {
            
            ret = EndsAtNode.equal_range( _LTree[i]->startnode );
            
            nchildren = 0;
            _LTree[i]->Child = NULL ;

            for (immii=ret.first; immii!=ret.second; immii++)
            {
                nchildren++;
                uplink = (*immii).second ;

                if ( ! (_LTree[i]->Child = (LinkTree **) realloc ( _LTree[i]->Child , nchildren * sizeof(LinkTree *) ) ) )
                {
                    SetErreur("MakeLinkTree : could not reallocate pointer _LTree[%d]->Child (size %d)\n",i,nchildren * sizeof ( LinkTree *) );
                    return(-1);
                }
                (_LTree[i]->Child)[nchildren-1] = _LTree[uplink] ;

                _LTree[uplink]->Parent = _LTree[i] ;
            }
            
            _LTree[i]->nchildren = nchildren ;

            // Si la geometrie est manquante pour le troncon, ligne droite entre startnode et endnode

            if ( ! ((_LTree[i]->Points)[0]) | ! ((_LTree[i]->Points)[1]) )
            {
                _LTree[i]->npoints = 2;

                if ( ! ( (_LTree[i]->Points)[0] = (double *) realloc ( (_LTree[i]->Points)[0] , 2 * sizeof ( double ) ) ) )
                {
                    SetErreur("MakeLinkTree : could not reallocate pointer _LTree[%d]->Points[0] (size 2)\n" );
                    return(-1);
                }

                if ( ! ( (_LTree[i]->Points)[1] = (double *) realloc ( (_LTree[i]->Points)[1] , 2 * sizeof ( double ) ) ) )
                {
                    SetErreur("MakeLinkTree : could not reallocate pointer _LTree[%d]->Points[1] (size 2)\n" );
                    return(-1);
                }
                
                (_LTree[i]->Points)[0][0] = _LTree[i]->xstart ;
                (_LTree[i]->Points)[1][0] = _LTree[i]->ystart ;

                (_LTree[i]->Points)[0][1] = _LTree[i]->xend ;
                (_LTree[i]->Points)[1][1] = _LTree[i]->yend ;
            }

            // Verification de l'orientation du troncon
            // (coherence entre les points de depart / arrivee topologique et geometrique)

            dxtopo =  _LTree[i]->xend -  _LTree[i]->xstart ;
            dytopo =  _LTree[i]->yend -  _LTree[i]->ystart ;

            npoints = _LTree[i]->npoints ;
            dxgeom =  (_LTree[i]->Points)[0][npoints-1] -  (_LTree[i]->Points)[0][0] ;
            dygeom =  (_LTree[i]->Points)[1][npoints-1] -  (_LTree[i]->Points)[1][0] ;

/*
                             ---------------->   ---------------->
            Produit scalaire DepartTopoFinTopo . DepartGeomFinGeom = orientation du troncon
*/
            scal = dxgeom * dxtopo + dygeom * dytopo ;

            if (scal<0)
            {
                _LTree[i]->flip = 1;
            }
            else
            {
                _LTree[i]->flip = 0;
            }

            // Calcul de la longueur du troncon

            if ( PolyLineLength (  & (_LTree[i]->Points) , _LTree[i]->npoints , & (_LTree[i]->Length) ) )
            {
                SetErreur("CreateLinkTree : could not compute length of segment %d\n",i);
                return(-1);
            }

        }

        EndsAtNode.clear();
        MapLinkId.clear();
        
        // Mise a jour des PK --------------------------------------------------

        if ( UpdatePK ( _LTree , nlinks ) )
        {
            SetErreur("CreateLinkTree : could not update PK\n");
            return(-1);
        }
              
        //----------------------------------------------------------------------

        (*_nlinks) = nlinks ;
        (*_ptrLTree) = _LTree ;
        return(0);   
}

//====================================================================================================================================

int UpdatePK ( LinkTree ** _LTree , int nlinks )
{
    // Parcours de l'arbre pour mettre a jour les PK aval de chaque troncon

    int i;

    for(i=0;i<nlinks;i++)
    {
        if ( !(_LTree[i]->Parent) ) // Le troncon i est n'a pas de troncon aval => PK = 0 puis construction recursive
        {
            _LTree[i]->PKend = 0;

            if ( UpdatePK_ ( _LTree[i] ) )
            {
                SetErreur("UpdatePK : recursice call of function UpdatePK_ failed at root node %d\n",i);
                return(-1);
            }

        }
    }

    return(0);
}

//====================================================================================================================================

int UpdatePK_ ( LinkTree * LTree )
{
    int i;

    for(i=0;i<LTree->nchildren;i++)
    {
/*        if ( !( (LTree->Child)[i] ) )
        {
            SetErreur("UpdatePK_ : %d fils declares pour le lien, mais (LTree->Child)[%d] est NULL\n",LTree->nchildren,i);
            return(-1);
        }
*/
        (LTree->Child)[i]->PKend = LTree->PKend + LTree->Length ;

        if ( UpdatePK_ ( (LTree->Child)[i] ) )
        {
            SetErreur("UpdatePK_ : recursive call failed\n");
            return(-1);
        }
    }

    return(0);
}

//====================================================================================================================================

int RaiseLinkTree ( LinkTree ** _LTree , int nlinks , double slope , FloatRaster * R )
{
    double z,zbase,Dz,N,x;
    int * CELLS ;
    int i,j,n,nr,nc,irec ;

    for(i=0;i<nlinks;i++)
    {
        zbase = (_LTree[i]->PKend) * slope ;
        Dz = (_LTree[i]->Length) * slope ;
        CELLS = NULL;

        if ( RasterizePolyline ( & (_LTree[i]->Points) , _LTree[i]->npoints , & (R->Header) , & CELLS , & n ) )
        {
            SetErreur("RaiseLinkTree : could not rasterize link %d\n",i);
            return(-1);
        }

        N = (double)(n-1);

        for(j=0;j<n;j++ )
        {
            if ( CELLS[j] >= 0 )
            {
                x = ((double) j ) / N;

                if ( !(_LTree[i]->flip) )
                    x = (1.0-x);

                z = zbase + x*Dz ;

                (void) rec2RC ( &(R->Header) , CELLS + j , & nr, & nc );
                R->Data[nr][nc] = (float)z;
            }
        }

        free(CELLS);
    }

    return(0);
}

//====================================================================================================================================

int PolyLineLength ( double * (*_PolyLine)[2] , int n , double * _length )
{
        int i;
    	complex<double> z0,z1;
        complex<double> I(0,1);

        z0 = (*_PolyLine)[0][0] + I * (*_PolyLine)[1][0];

        (*_length) = 0;

	for(i=1;i<n;i++)
	{
		z1 = (*_PolyLine)[0][i] + I * (*_PolyLine)[1][i];
		(*_length) += abs(z1-z0);
                z0 = z1;
        }

        return(0);
}

//====================================================================================================================================

int RasterizePolyline ( double * (*_PolyLine)[2] , int npoints , esriheader * Header, int ** _CELLS , int * _ncells )
{
	complex<double>z0,z1,Vt,z;
        complex<double> I(0,1);
        vector<int> PolyLineCells;
	vector<int>::iterator im;
	int rec,i,j;
        double pas ;
	double L,nt;

	pas = Header->cellsize ;
	z0 = (*_PolyLine)[0][0] + I * (*_PolyLine)[1][0];
	(void) c2rec ( Header , z0 , & rec );
        PolyLineCells.clear();

	if(rec>=0)
		PolyLineCells.push_back(rec);

	for(i=1;i<npoints;i++)
	{
		z1=(*_PolyLine)[0][i] + I * (*_PolyLine)[1][i];
		L=abs(z1-z0);
		nt=(double)ceil(L/pas);
		Vt=(z1-z0)/nt;
		for(j=0;j<(int)nt;j++)
		{
			z=z0+(double)(j+1)*Vt;
			(void) c2rec ( Header , z , & rec );
			im = find(PolyLineCells.begin(),PolyLineCells.end(),rec);

			// Si on a une boucle non-descriptible a la resolution de la grille,
			// on supprime toutes les cellules comprises entre la premiere et la deuxieme occurrence
			// de la cellule "i"
			if(im!=PolyLineCells.end() & rec>=0)
			{
				PolyLineCells.erase(im,PolyLineCells.end());
			}

			if(rec>=0)
			PolyLineCells.push_back(rec);

		}
		z0=z1;
	}

        (*_ncells) = PolyLineCells.size();
        
        if ( ! ((*_CELLS) = (int *) malloc( PolyLineCells.size() * sizeof(int) ) ) )
        {
            SetErreur("RasterizePolyLine : could not allocate pointer *_CELLS (size %d)\n",PolyLineCells.size() * sizeof(int) );
            return(-1);
        }

        for (i=0;i<PolyLineCells.size();i++)
        {
            (*_CELLS)[i] = PolyLineCells[i];
        }
        PolyLineCells.clear();

        return(0);
}

//====================================================================================================================================

int RaiseStreamNetwork ( char strFichGeom[120] , char strFichLinks[120] , char strFichNodes[120] , int HandleRaster , double slope )
{
    LinkTree ** _LTree ;
    FloatRaster * R ;
    int nlinks;
    double bb[4];

    _LTree = NULL ;

    if ( !(R = GetFloatRaster(HandleRaster)) )
    {
        SetErreur("RaiseStreamNetwork : raster # %d does not exist\n",HandleRaster);
        return(-1);
    }

    bb[0] = R->Header.xll ;
    bb[1] = bb[0] + (R->Header.ncols) * (R->Header.cellsize) ;
    bb[2] = R->Header.yll ;
    bb[3] = bb[2] + (R->Header.nrows) * (R->Header.cellsize) ;

    if ( MakeLinkTree ( & _LTree , strFichGeom , strFichLinks , strFichNodes , bb , & nlinks ) )
    {
        SetErreur("RaiseStreamNetwork : could not create link tree\n");
        return(-1);
    }

    if ( RaiseLinkTree ( _LTree , nlinks , slope , R ) )
    {
        SetErreur("RaiseStreamNetwork : could not compute elevation of links\n");
        return(-1);
    }

    if ( DestroyLinkTree ( & _LTree , nlinks ) )
    {
        SetErreur("RaiseStreamNetwork : could not destroy link tree\n");
        return(-1);
    }

    return(0);

}

//====================================================================================================================================

int DestroyLinkTree ( LinkTree *** _ptrLTree , int nlinks )
{
    int i;

    for (i=0;i<nlinks;i++)
    {
        free ( ((*_ptrLTree)[i]->Points)[0] );
        free ( ((*_ptrLTree)[i]->Points)[1] );
        free ( (*_ptrLTree)[i]->Child );
        free ( (*_ptrLTree)[i] );
    }

    // Liberation du tableau complet

    free ( (*_ptrLTree) );

    return(0);
}

