/* 
* 
* 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 "FlowDirGrid.h"
#include "Topology.h"
#include "esriheader.h"
#include "Erreur.h"
#include <set>
#include <map>
#include <iterator>

int verifyUpstreamDownstream ( FlowDirGrid * F , int * _src , int * _dst )
{
    int nr_src, nc_src, nr_dst, nc_dst , Dr , Dc ;

    (void) rec2RC ( &(F->Header) , _src , & nr_src , & nc_src );
    (void) rec2RC ( &(F->Header) , _dst , & nr_dst , & nc_dst );

    Dc = nr_dst - nr_src ;
    Dr = nc_dst - nc_src ;

    if ( ( Dc*Dc +Dr*Dr ) < 2 )
    {
        return ( (int)( F->Data[nr_src][nc_src] == F->Compass[Dr+1][Dc+1] ) );
    }
    else
    {
        return(0); /* Les cellules src et dst ne sont pas voisines -> forcement faux */
    }
}

int GetDownstream ( FlowDirGrid * F , int * _cell )
{
    int nr, nc, Dr, Dc, ncols , nrows , nc1 , nr1 , down ;
    int code;

    (void) rec2RC ( &(F->Header) , _cell , & nr , & nc );
     code = (int) F->Data[nr][nc];
     Dc = code % 3 ;
     Dr = ( code - Dc ) / 3 ;

     nr1 = nr+Dr-1;
     nc1 = nc+Dc-1;

     nrows = F->Header.nrows;
     ncols = F->Header.ncols;

     if ( nr1<0 | nr1>=nrows | nc1<0 | nc1>=ncols )
        down = -1;
     else
        down = nr1*ncols+nc1;
        
     return(down);
}

int GetUpstream ( FlowDirGrid * F , int * _cell , int * _nup , int ** _UPSTREAM )
{
    int Dr,Dc,nr,nc,Drmin,Drmax,Dcmin,Dcmax,irec,nup;
    
    *_UPSTREAM = NULL;
    nup = 0;
    
    (void) rec2RC ( &(F->Header) , _cell , & nr , & nc );
    (void) GetNeighborhood ( &(F->Header) , & nr , & nc , & Drmin , & Drmax , & Dcmin , & Dcmax );

    for(Dr=Drmin;Dr<=Drmax;Dr++)
    {
        for(Dc=Dcmin;Dc<=Dcmax;Dc++)
        {
			if ( Dr !=0 | Dc !=0 )
			{
				if ( F->Data[nr+Dr][nc+Dc] == F->Compass[1-Dr][1-Dc] )
				{
					irec = (nr+Dr)*(F->Header.ncols)+(nc+Dc);

					if ( !( (* _UPSTREAM) = (int *) realloc( (* _UPSTREAM) ,(nup+1)*sizeof(int) ) ) )
					{
						SetErreur("GetUpstream : could not reallocate integer pointer\n");
						return(-1);
					}

					(*_UPSTREAM)[nup] = irec;
					nup += 1;
				}	
            }
        }
    }

    *_nup = nup;

    return(0);

}



int GetCatchmentSize_ ( FlowDirGrid * F , int * _outlet, set<int> & UpNodes, int ** _BB , int * _nup , set<int> & ixUpNodes , int * _nr , int * _nc )
{
    int nup_i , i , nup1 ;
    set<int>::iterator isi;
    int * UPSTREAM;

    if ( GetUpstream ( F , _outlet , & nup1 , & UPSTREAM ) )
    {
            SetErreur("GetCatchmentSize_ : could not get upstream cells\n");
            return(-1);
    }

    *_nup = 1;
    (void) rec2RC ( &(F->Header) , _outlet , _nr , _nc );

    // Mise Ã  jour de la Bounding Box

    if ( (*_nc) < (*_BB)[0] )
        (*_BB)[0] = (*_nc);     // col min

    if ( (*_nr) < (*_BB)[1] )
        (*_BB)[1] = (*_nr);     // row min

    if ( (*_nc) > (*_BB)[2] )
        (*_BB)[2] = (*_nc);     // col max

    if ( (*_nr) > (*_BB)[3] )
        (*_BB)[3] = (*_nr);     // row max

    // Boucle sur les pixels amont

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

        isi = UpNodes.find(UPSTREAM[i]);

        if ( isi == UpNodes.end() )       // Si le pixel UPSTREAM[i] n'est pas trouvÃ© dans la liste des noeuds d'arrÃªt
        {
            if ( GetCatchmentSize_ ( F , UPSTREAM + i , UpNodes , _BB , & nup_i , ixUpNodes , _nr , _nc ) )     // Appel rÃ©cursif
            {
                SetErreur("GetCatchmentSize_ : recursive call failed\n");
                return(-1);
            }

            (*_nup) += nup_i ;
        }
        else
        {
            ixUpNodes.insert ( (int) distance ( UpNodes.begin() , isi ) );
        }
    }
 
    free(UPSTREAM);

    return (0);
}


int GetCatchmentSize ( int HandleFlowDir , int outlet, int nupnodes, int * UPNODES , int ** _BB , int * _nup , int ** _IXUPNODES , int * _nixupnodes )
{
    int i;
    FlowDirGrid * F;
    set<int> UpNodes,ixUpNodes;
    set<int>::iterator isi;
    int nr;
    int nc;

    *_nup = -1;

    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
        SetErreur("GetCatchmentSize :flow direction grid %d does not exist\n",HandleFlowDir);
	return(-1);
    }

    // Recopie des Ã©lÃ©ments du tableau UPNODES dans le set UpNodes

    UpNodes.clear();
    ixUpNodes.clear();

    for (i=0;i<nupnodes;i++)
    {
        UpNodes.insert( UPNODES[i] );
    }

    (*_BB)[0] = F->Header.ncols - 1;
    (*_BB)[1] = F->Header.nrows - 1;
    (*_BB)[2] = 0;
    (*_BB)[3] = 0;

    if ( GetCatchmentSize_ ( F , & outlet , UpNodes, _BB , _nup , ixUpNodes , & nr , & nc ) )
    {
        SetErreur("GetCatchmentSize : recursive call failed\n");
	return(-1);
    }

    // Recopie des Ã©lÃ©ments du set ixUpNodes dans le tableau IXUPNODES

    (*_IXUPNODES) = (int *) malloc ( ixUpNodes.size()*sizeof(int) );
    (*_nixupnodes) = ixUpNodes.size();

    i=0;
    for (isi=ixUpNodes.begin();isi!=ixUpNodes.end();isi++)
    {
        (*_IXUPNODES)[i] = (*isi);
        i++;
    }

    // Vidage des containers

    ixUpNodes.clear();
    UpNodes.clear();

    return(0);
}

int GetCatchmentCells_ ( FlowDirGrid * F , int * _outlet , set<int> & Inlets , set<int> & CATCHMENT , set<int> & SelectInlets )
{
    int i , nup1 ;
    int * UPSTREAM;
	set<int>::iterator isi;

    if ( GetUpstream ( F , _outlet , & nup1 , & UPSTREAM ) )
    {
            SetErreur("GetCatchmentCells_ : could not get upstream cells\n");
            return(-1);
    }

    for(i=0;i<nup1;i++)
    {
		isi = Inlets.find(UPSTREAM[i]);
		
		if (isi==Inlets.end())
		{
			/* Recuperation des cellules du sous-bv i */

			if ( GetCatchmentCells_ ( F ,UPSTREAM + i , Inlets , CATCHMENT , SelectInlets ) )
			{
				SetErreur("GetCatchmentCells_ : recursive call failed\n");
				return(-1);
			}
		}
		else
		{
			SelectInlets.insert((*isi));
		}
    }

    CATCHMENT.insert(*_outlet);  /* Initialisation avec la cellule exutoire */

    return (0);
}

int GetCatchmentCells ( int HandleFlowDir , int outlet, int * INLETS , int ninlets , int * _nup , int ** _CATCHMENT , int ** _bb , int ** _SELINLETS , int * _nselinlets )
{
    FlowDirGrid * F;
    set<int> CATCHMENT;
    set<int>::iterator isi;
    int i,rec,nc,ncols,nrmin,nrmax,ncmin,ncmax;
	set<int> Inlets , SelectInlets ;
	
    /* Initialisation du pointeur de retour en cas d'echec */

    *_CATCHMENT = NULL;

    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
        SetErreur("GetCatchmentCells : flow direction grid %d does not exist\n",HandleFlowDir);
		return(-1);
    }

    ncols = F->Header.ncols;
	
	/* Recopie du tableau INLETS dans le container Inlets */ 
	
	Inlets.clear();
	
	for(i=0;i<ninlets;i++)
		Inlets.insert(INLETS[i]);
		
	SelectInlets.clear();
	
    if ( GetCatchmentCells_ ( F , & outlet , Inlets , CATCHMENT , SelectInlets ) )
    {
        SetErreur("GetCatchmentCells : recursive call failed\n");
		return(-1);
    }

    /* (RÃ©-)allocation du tableau d'entiers */

    if ( !( (* _CATCHMENT) = (int *) realloc( (* _CATCHMENT) , CATCHMENT.size()*sizeof(int) ) ) )
    {
        SetErreur("GetCatchmentCells : could not reallocate integer pointer\n");
		return(-1);
    }

    /* Recopie des valeurs du set dans le tableau d'entiers et calcul de la "bounding-box" */

    *_nup = CATCHMENT.size();

    ncmin = ncols;
    ncmax = 0;

    isi = CATCHMENT.begin();
    rec = (*isi);
    (void) rec2RC( & (F->Header) , & rec , & nrmin , & nc );
    
    isi = CATCHMENT.end();isi--;
    rec = (*isi);
    (void) rec2RC( & (F->Header) , & rec , & nrmax , & nc );

    i=0;

    for (isi=CATCHMENT.begin();isi!=CATCHMENT.end();isi++)
    {
        (*_CATCHMENT)[i] = (*isi);
        i++;

        nc = (*isi) % ncols;

        if (nc<ncmin)
            ncmin = nc;

        if (nc>ncmax)
            ncmax = nc;
    }

    CATCHMENT.clear();

    (*_bb)[0] = ncmin ;
    (*_bb)[1] = ncmax ;
    (*_bb)[2] = nrmin ;
    (*_bb)[3] = nrmax ;
    
	/* Recopie des valeurs des injections recontrees */
	 
	if ( !( (*_SELINLETS) = (int *) realloc( (*_SELINLETS) , SelectInlets.size()*sizeof(int) ) ) )
    {
        SetErreur("GetCatchmentCells : could not reallocate integer pointer\n");
		return(-1);
    }
	
	i=0;
	*_nselinlets = SelectInlets.size();
	
	for (isi=SelectInlets.begin();isi!=SelectInlets.end();isi++)
	{
		(*_SELINLETS)[i] = (*isi) ;
		i++ ;
	}
	
	Inlets.clear();
	SelectInlets.clear();
	
    return(0);
}


int GetDownstreamSet ( const int HandleFlowDir , int ncells , int * CELLS , int ** _DOWNSTREAM )
{
    int i, nr, nc, Dr, Dc, ncols , nrows , nc1 , nr1 , down ;
    int code;
    FlowDirGrid * F;

    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
        SetErreur("GetDownstreamSet : flow direction grid %d does not exist\n",HandleFlowDir);
	return(-1);
    }

    for (i=0;i<ncells;i++)
    {
        (void) rec2RC ( &(F->Header) , CELLS + i , & nr , & nc );
        code = (int)( F->Data[nr][nc] );
        Dc = code % 3 ;
        Dr = ( code - Dc ) / 3 ;

        nr1 = nr+Dr-1;
        nc1 = nc+Dc-1;

        nrows = F->Header.nrows;
        ncols = F->Header.ncols;

        if ( nr1<0 | nr1>=nrows | nc1<0 | nc1>=ncols )
            down = -1;
        else
            down = nr1*ncols+nc1;

        (*_DOWNSTREAM)[i] = down;
    }

    return(0);

}

int MakeFlowLineBetweenPoints ( FlowDirGrid * F , int A , int B, int ** _FlowLine , int * _length )
{
    int Amont,Aval,n ;

    Amont = A;

    if ( ! ( (*_FlowLine) = (int *) malloc ( sizeof(int) ) ) )
    {
        SetErreur("MakeFlowLineBetweenPoints : could not allocate integer pointer\n",A,B);
        return(-1);
    }

    (*_FlowLine)[0] = Amont;
    n=1;

    Aval = GetDownstream ( F , & Amont );
    
    while ( Aval>=0 & Amont!=B & Aval!=Amont )
    {
        n++;
        
        if ( ! ( (*_FlowLine) = (int *) realloc( (*_FlowLine), n*sizeof(int)) ) )
        {
            SetErreur("MakeFlowLineBetweenPoints : could not reallocate integer pointer\n",A,B);
            return(-1);
        }

        (*_FlowLine)[n-1] = Aval;
        Amont = Aval;
        Aval = GetDownstream ( F , & Amont );
    }

    (*_length) = n;

    // Si on n'arrive pas au point B demandÃ©, erreur

    if ( (*_FlowLine)[n-1] != B )
    {
        SetErreur("MakeFlowLineBetweenPoints : pixel %d is not located downstream of pixel %d\n",A,B);
        free(*_FlowLine);
	return(-1);
    }

    return(0);

}

int MakeFlowLine ( int HandleFlowDir , int A , int B, int ** _FlowLine , int * _length )
{
    FlowDirGrid * F;

    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
        SetErreur("MakeFlowLine : flow direction grid %d does not exist\n",HandleFlowDir);
	return(-1);
    }

    if ( MakeFlowLineBetweenPoints ( F , A , B, _FlowLine , _length ) )
    {
        SetErreur("MakeFlowLine : could not build flowline between cells %d and %d\n",A,B);
	return(-1);
    }

    return(0);

}

int MakeArcs_ ( FlowDirGrid * F , int * _root , map<int,int> & Nodes , int ** _FLAG , int * _downid , int * & UP , int * & DOWN , int * _narcs )
{
    int i, nup1 , upid , narcs ;
    int * UPSTREAM;
    map<int,int>:: iterator imii;

    if ( GetUpstream ( F , _root , & nup1 , & UPSTREAM ) )
    {
            SetErreur("MakeArcs_ : could not get upstream cells\n");
            return(-1);
    }

    narcs = (*_narcs);

    for(i=0;i<nup1;i++)
    {
            /* Test pour savoir si UPSTREAM[i] est un noeud */

            imii = Nodes.find ( UPSTREAM[i] );

            if ( imii != Nodes.end() )
            {
                    upid = (*imii).second;

                    if ( (*_FLAG) [upid] )
                    {
                        /* Si le noeud n'a pas encore été parcouru, on poursuit récursivement */

                        if ( MakeArcs_ ( F , UPSTREAM + i , Nodes , _FLAG , & upid , UP , DOWN , & narcs ) )
                        {
                            SetErreur("MakeArcs_ : recursive call failed\n");
                            return(-1);
                        }

                        (*_FLAG) [upid] = 0; /* on marque le noeud comme parcouru */
                    }

                    /* On ajoute l'arc upid->downid */

                    narcs = narcs + 1 ;

                    if ( !(UP = (int *) realloc ( UP , narcs * sizeof (int) )) )
		    {
			SetErreur("MakeArcs_: could not reallocation pointer *_UP, size %d\n",narcs);
			return(-1);
		    }
                    UP[narcs-1] = upid ;

                    if ( !(DOWN = (int *) realloc ( DOWN , narcs * sizeof (int) )) );
		    {
			SetErreur("MakeArcs_: could not reallocate pointer *_DOWN, size %d\n",narcs);
			return(-1);
		    }
                    DOWN[narcs-1] = (*_downid) ;
            }
            else
            {
                    if ( MakeArcs_ ( F , UPSTREAM + i , Nodes , _FLAG , _downid , UP , DOWN , & narcs ) )
                    {
                        SetErreur("MakeArcs_ : recursive call failed\n");
                        return(-1);
                    }
            }
    }

    free(UPSTREAM);
    (*_narcs) = narcs;

    return(0);

}

int MakeArcs ( int HandleFlowDir , int * NODES , int nnodes , int ** _UP , int ** _DOWN , int * _narcs )
{
    FlowDirGrid * F ;
    map<int,int> Nodes ;
    int * FLAG , * UP , * DOWN ;
    int i,downid,root,narcs ;

    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
        SetErreur("MakeArcs : flow direction grid %d does not exist\n",HandleFlowDir);
	return(-1);
    }

    /* Construction du set et du tableau de travail */

    Nodes.clear() ;

    if ( !(FLAG = (int * ) malloc ( nnodes * sizeof(int) ) ) )
    {
        SetErreur ("MakeArcs : could not allocate integer pointer FLAG\n");
        return(-1) ;
    }

    for(i=0;i<nnodes;i++)
    {
        Nodes.insert ( make_pair ( NODES[i] , i ) ) ;
        FLAG[i] = 1 ;
    }

    narcs = 1 ;
    UP = (int *) malloc ( 1*sizeof(int));
    DOWN = (int *) malloc ( 1*sizeof(int));

    /* Boucle sur les noeuds */

    for(i=0;i<nnodes;i++)
    {
        downid = i ;
//        root = NODES[i];

        /* Construction de tous les arcs à l'amont du noeud i */

	if ( FLAG[i] )
	{
		if ( MakeArcs_ ( F , NODES + i , Nodes , & FLAG , & downid , UP , DOWN , & narcs ) )
		{
		    SetErreur("MakeArcs : initial recursive call failed at cell %d\n",downid);
		    return(-1);
		}

		FLAG[i] = 0;
	}
    }

    (*_narcs) = narcs;
    (*_UP) = UP ;
    (*_DOWN) = DOWN ;

    free(FLAG) ;
    Nodes.clear();

    return(0) ;

}


int GetDownstreamPath ( int HandleFlowDir , int start , int ** _FlowLine , int * _length  )
{
    FlowDirGrid * F ;
	int n,Amont,Aval ;

    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
		SetErreur("GetDownstreamPath : flow direction grid %d does not exist\n",HandleFlowDir);
		return(-1);
    }

	if ( ! ( (*_FlowLine) = (int *) malloc ( sizeof(int) ) ) )
    {
        SetErreur("GetDownstreamPath : could not allocate integer pointer\n");
        return(-1);
    }
	
	
    Amont = start ; 
	(*_FlowLine)[0] = Amont;
    n=1;

    Aval = GetDownstream ( F , & Amont );
		
	while ( Aval>=0 & Aval!=Amont )
    {
        n++;
        
        if ( ! ( (*_FlowLine) = (int *) realloc( (*_FlowLine), n*sizeof(int)) ) )
        {
            SetErreur("GetDownstreamPath : could not reallocate integer pointer\n");
            return(-1);
        }

        (*_FlowLine)[n-1] = Aval;
        Amont = Aval;
        Aval = GetDownstream ( F , & Amont );
    }
	
	(*_length) = n ;
	
    return(0) ;

}

