/* 
* 
* 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 "CellTree.h"
#include "Erreur.h"
#include "FlowDirGrid.h"
#include "FloatRaster.h"
#include "esriheader.h"
#include "Strahler.h"

using namespace std;

int StrahlerOrdering_ ( CellTree * T , float * _threshold , FloatRaster * R , int * _nsrc , int ** _SRC )
{
     int nr,nc,s1,s2,s,rec;
     s1=0;
     s2=0;

    if ( T->FirstChild )    // Si l'arbre T a un fils, on continue à construire la numérotation
    {

         if ( StrahlerOrdering_ ( T->FirstChild , _threshold , R , _nsrc , _SRC ) )
         {
             SetErreur("StrahlerOrdering_ : recursive call failed for first-child tree\n");
             return(-1);
         }

         rec = (T->FirstChild)->RootCell;
         (void) rec2RC ( &(R->Header) , & rec  , & nr , & nc );
         s1 = (int) (R->Data[nr][nc]);

         if ( (T->FirstChild)->NextSibling )    // Si l'arbre T a un deuxième fils
         {
            rec = ((T->FirstChild)->NextSibling)->RootCell;
            (void) rec2RC ( &(R->Header) , & rec , & nr , & nc );
            s2 = (int) (R->Data[nr][nc]);
         }

    }

    /* Calcul de l'ordre de Strahler de l'arbre T à partir des ordres de ses 2 plus grands fils */

    if (s1==0 & s2==0)
    {   
        if ( (T->Size) >= (*_threshold) )   // L'arbre T est le début d'un ordre 1
        {
            s = 1;
            // On ajoute la cellule à la liste des sources :
            (*_nsrc) += 1;
            (*_SRC) = (int *) realloc( (*_SRC), (*_nsrc)*sizeof(int));
            (*_SRC)[(*_nsrc)-1] = T->RootCell ;
        }
        else
            s = 0;
    }
    else
    {
        if (s1==s2)     // Confluence de 2 ordres de Strahler identiques => incrément de 1
            s = s1+1;
        else
        {
            if (s1>s2)
                s = s1;     // Confluence de 2 ordres de Strahler différents => on garde le max
            else
                s = s2;
        }
    }

    (void) rec2RC ( &(R->Header) , &(T->RootCell) , & nr , & nc );
    R->Data[nr][nc] = (float)s;

    /* Si l'arbre T a un arbre frère, appel sur l'arbre frère */

    if ( T->NextSibling )
    {
        if ( StrahlerOrdering_ ( T->NextSibling , _threshold , R , _nsrc , _SRC ) )
        {
            SetErreur("StrahlerOrdering_ : recursive call failed for next-sibling tree\n");
            return(-1);
        }
    }

    return(0);
}

int StrahlerOrdering ( int HandleTree , float threshold , int HandleRaster  , int * _nsrc , int ** _SRC )
{
    CellTree * T;
    FloatRaster * R;
//   float * _threshold;

    if( !(T = GetCellTree(HandleTree) ) )
    {
        SetErreur("StrahlerOrdering : tree # %d does not exist\n",HandleTree);
	return(-1);
    }

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

//    (*_threshold) = threshold;

    (*_nsrc) = 0;
    (*_SRC) = (int *) malloc( 0 * sizeof(int) );
    
    if ( StrahlerOrdering_ ( T , & threshold , R , _nsrc , _SRC ) )
    {
        SetErreur("StrahlerOrdering : recursive call failed at root cell\n");
        return(-1);
    }

    return(0);
}

int MakeStrahlerStreams_ ( CellTree * T , float * _threshold , CellTree * Ancestor , int * _s , int * _ssib , int * _nstreams , int ** _StreamLength , int ** _StreamOrder , int ** _StreamConnect , int *** _StreamCells , int * _ncurr , int * _nsib )
{
    int s1,s2,n1,n2,len;

    (*_s) = 0;
    (*_ssib) = 0;

    if ( T->Size < (*_threshold) )
        return(0);

    if ( T->FirstChild )    // Si l'arbre T a un fils, on continue à construire le drain courant
    {

         if ( MakeStrahlerStreams_ ( T->FirstChild , _threshold , T , & s1 , & s2 , _nstreams , _StreamLength , _StreamOrder , _StreamConnect , _StreamCells , & n1 , & n2 )  )
         {
             SetErreur("MakeStrahlerStreams_ : recursive call failed for first-child tree\n");
             return(-1);
         }
    }

    /* Calcul de l'ordre de Strahler de l'arbre T à partir des ordres de ses 2 plus grands fils */

    if (s1==s2) // Confluence de 2 ordres de Strahler identiques => incrément de 1, création d'un nouveau drain
    {
        // L'arbre T est le début d'un ordre s1+1

        (*_s) = s1+1;

        // On crée un nouveau drain

        (*_nstreams) += 1;

        (*_StreamLength) = (int *) realloc ( *_StreamLength , (*_nstreams) * sizeof (int) );
        (*_StreamLength)[(*_nstreams)-1] = 1;

        (*_StreamOrder) = (int *) realloc ( *_StreamOrder , (*_nstreams) * sizeof (int) );
        (*_StreamOrder)[(*_nstreams)-1] = s1+1;

        (*_StreamConnect) = (int *) realloc ( *_StreamConnect , (*_nstreams) * sizeof (int) );

        (*_StreamCells) = (int **) realloc ( *_StreamCells , (*_nstreams) * sizeof (int *) );
        (*_StreamCells)[(*_nstreams)-1] = (int *) malloc ( sizeof(int) ) ;
        (*_StreamCells)[(*_nstreams)-1][0] = T->RootCell ;

        (*_ncurr) = (*_nstreams) ;
    }
    else
    {
        if (s1>s2)
        {
            (*_s) = s1;     // Confluence de 2 ordres de Strahler différents => on garde le max
            (*_ncurr) = n1;
        }
        else
        {
            (*_s) = s2;
            (*_ncurr) = n2;
        }

        len = (*_StreamLength)[(*_ncurr)-1] ;
        len += 1;
        (*_StreamLength)[(*_ncurr)-1] = len;
        (*_StreamCells)[(*_ncurr)-1] = (int *) realloc ( (*_StreamCells)[(*_ncurr)-1] , len * sizeof(int) ) ;
        (*_StreamCells)[(*_ncurr)-1][len-1] = T->RootCell ;

    }
      
    (*_StreamConnect)[(*_ncurr)-1] = Ancestor->RootCell ;

    /* Si l'arbre T a un arbre cadet, appel sur cet arbre frère, de même ancêtre */

    if ( T->NextSibling )
    {
        if ( MakeStrahlerStreams_ ( T->NextSibling , _threshold , Ancestor , _ssib , & s2 , _nstreams , _StreamLength , _StreamOrder , _StreamConnect , _StreamCells , _nsib , & n2 )  )
        {
            SetErreur("MakeStrahlerStreams_ : recursive call failed for next-sibling tree\n");
            return(-1);
        }
    }

    return(0);
}

int MakeStrahlerStreams ( int HandleTree , float * _threshold , int * _nstreams , int ** _StreamLength , int ** _StreamOrder , int ** _StreamConnect , int *** _StreamCells )
{
    CellTree * T;
    int s,ssib,ncurr,nsib ;

    if( !(T = GetCellTree(HandleTree) ) )
    {
        SetErreur("MakeStrahlerStreams : tree # %d does not exist\n",HandleTree);
	return(-1);
    }

    ncurr = 0;
    (*_nstreams) = 0;
    (*_StreamLength) = NULL ;
    (*_StreamOrder) = NULL ;
    (*_StreamConnect) = NULL ;
    (*_StreamCells) = NULL ;

    if ( MakeStrahlerStreams_ ( T , _threshold , T , & s , & ssib , _nstreams , _StreamLength , _StreamOrder , _StreamConnect , _StreamCells , & ncurr , & nsib ) )
    {
        SetErreur("MakeStrahlerStreams : recursive call failed at root cell\n");
        return(-1);
    }

    return(0);
}
