//==============================================================================
/* Nouvelle-Fonction */

// code By Reza Ansari


NTuple *  MkAssocNTuple (StarList* stl1, StarList* stl2)
{
    // Association des listes d'objets (sources)
    stl1->Complete();
    stl2->Complete();
    StarList::SetPrtLevel(2);
    stl1->AlGTransfo().SetIdentite(1);
    int rc = stl1->AlignGeom(stl2, 3, 50, AlG_Transl);
    cout << " stl1->AlignGeom() -> " << rc << endl;
    if (rc) { stl1->resTr1_2 = 1;  stl1->AlGTransfo().SetIdentite(1); }
    int nassoc = stl1->AssocStars(stl2, 3);
    cout << " stl->AssocStars() NAssoc= " << nassoc << endl;
    
    char* nomnt[6] = { "x", "y", "ysz", "yszdet", "dismin", "dism2" } ;
    NTuple* nt = new NTuple(6, nomnt);
    float xnt[10];
    BStar *s1, *s2;
    int i,j;
    for(i=0; i<stl1->NbStars(); i++) 
    {
	for(j=0; j<10; j++) xnt[j] = -0.;
	s1 = stl1->Star(i);
	if (!s1->Nice())  continue;
	xnt[0] = s1->PosX();
	xnt[1] = s1->PosY();
	xnt[2] = s1->Flux();
	j = stl1->XRef(i);        // probleme ici a cause de aligngeom ou assocstar 
	if (j >= 0) 
	{
	    
	    s2 = stl2->Star(j);
	    if (s2->Nice()) 
	    {

		xnt[3] = s2->Flux();
		xnt[4] = stl1->DAss(i);
		xnt[5] = stl1->DAss2(i);
	    }
        }
	nt->Fill(xnt); 
    }
    return(nt);
}

//==============================================================================
/* Nouvelle-Fonction */ 
// la fonction repond a la question : le point de coordonnees (x,y) est-
// il un maximum local ( sur un pave 9*9 centre sur (x,y) )
// ceci quelques soient x et y

void MaxImg(ImageR4* im,double x,double y,bool& max)
{
  
   	
    if ( ( x>0 ) && ( x<im->XSize()-1 ) && ( y>0 ) && ( y<im->YSize()-1 ) )
    // pour le centre de l'image (8 points a verifier)
    {
		
	if (( (*im)(x,y) > (*im)(x-1,y-1) )&& ( (*im)(x,y) > (*im)(x,y-1) ))
	{
	    if (( (*im)(x,y) > (*im)(x+1,y-1) ) && ( (*im)(x,y) > (*im)(x-1,y+1) ))
	    {
		if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x+1,y+1) ))
		{
		    if (( (*im)(x,y) > (*im)(x-1,y) ) && ( (*im)(x,y) > (*im)(x+1,y) ))
		    {
			max= true;
		    }
		}
	    }
	}
    }
    if ( (x==0) &&  ( y>0 ) && ( y<im->YSize()-1 ) )
    // bord en x a "gauche" (5 points a verifier)
    {
	if (( (*im)(x,y) > (*im)(x+1,y-1) ) && ( (*im)(x,y) > (*im)(x,y-1) ))
	{
	    if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x+1,y+1) ))
	    {	
		if ( (*im)(x,y) > (*im)(x+1,y) )
		{
		    max= true;
		}
	    }
	}
    }
    if ( (y==0) &&  ( x>0 ) && ( x<im->XSize()-1 ) )
    // bord en y en haut (5 points )
    {
	if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x+1,y+1) ))
	{
	    if (( (*im)(x,y) > (*im)(x-1,y) ) && ( (*im)(x,y) > (*im)(x+1,y) ))
	    {
		if ( (*im)(x,y) > (*im)(x-1,y+1) )
		{
		    max= true;
		}
	    }
	}
    }
    if ( (x==im->XSize()-1 ) &&  ( y>0 ) && ( y<im->YSize()-1 ) )
    // bord en x a "droite" (5 points a verifier)
    {
	if (( (*im)(x,y) > (*im)(x-1,y-1) ) && ( (*im)(x,y) > (*im)(x,y-1) ))
	{
	    if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x-1,y+1) ))
	    {	
		if ( (*im)(x,y) > (*im)(x-1,y) )
	       	{
		    max= true;
		}
	    }
	}
    }
    if ( (y==im->YSize()-1 ) &&  ( x>0 ) && ( x<im->XSize()-1 ) )
    // bord en y en bas (5 points )
    {
	if (( (*im)(x,y) > (*im)(x,y-1) ) && ( (*im)(x,y) > (*im)(x+1,y-1) ))
	{
	    if (( (*im)(x,y) > (*im)(x-1,y) ) && ( (*im)(x,y) > (*im)(x+1,y) ))
	    {
		if ( (*im)(x,y) > (*im)(x-1,y-1) )
		{
		    max= true;
		}
	    }
	}
    }
    
    if ( (x==0) &&  ( y==0 ) )
    // coin en haut a "gauche" (3 points a verifier)
    {
	if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x+1,y+1) ))
	{	
	    if ( (*im)(x,y) > (*im)(x+1,y) )
	    {
		max= true;
	    }
	}
    }
    

    if ( (x==im->XSize()-1 ) &&  ( y==0 ) )
    // coin en haut a "droite" (3 points a verifier)
    {
	if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x-1,y+1) ))
	{	
	    if ( (*im)(x,y) > (*im)(x-1,y) )
	    {
		max= true;
	    }
	}
    }
		
    if ( (x==im->XSize()-1 ) &&  ( y==im->YSize()-1 ) )
    // coin en bas a "droite" (3 points a verifier)
    {
	if (( (*im)(x,y) > (*im)(x,y-1) ) && ( (*im)(x,y) > (*im)(x-1,y-1) ))
	{	
	    if ( (*im)(x,y) > (*im)(x-1,y) )
	    {
		max= true;
	    }
	}
    }
    if ( (x==0) &&  ( y==im->YSize()-1 ) )
    // coin en bas a "gauche" (3 points a verifier)
    {
	if (( (*im)(x,y) > (*im)(x,y-1) ) && ( (*im)(x,y) > (*im)(x+1,y-1) ))
	{	
	    if ( (*im)(x,y) > (*im)(x+1,y) )
	    {
		max= true;
	    }
	}
    }
}


//==============================================================================
/* Nouvelle-Fonction */ 
// calcule la somme des pixels dans un carre 9*9 centre sur (x,y)
// dans l'image im, ceci quelques soient x et y

void SomPix(ImageR4& im,double x,double y,float& som)
{
  
   	
    if ( ( x>0 ) && ( x<im.XSize()-1 ) && ( y>0 ) && ( y<im.YSize()-1 ))
    // pour le centre de l'image (9 points)
    {
		
	som=  im(x,y) +  im(x,y+1) + 
		im(x,y-1) + im(x+1,y) +
		im(x-1,y) + im(x-1,y-1) +
		im(x-1,y+1) + im(x+1,y+1) +
		im(x+1,y-1);
    }
    if ( (x==0) &&  ( y>0 ) && ( y<im.YSize()-1 ) )
    // bord en x a "gauche" (6 points )
    {
	som= im(x,y) +  im(x,y+1) + im(x,y-1) + im(x+1,y) + im(x+1,y+1) +
		im(x+1,y-1);
    }

    if ( (y==0) &&  ( x>0 ) && ( x<im.XSize()-1 ) )
    // bord en y en haut (6 points )
    {
	som=  im(x,y) +  im(x,y+1) + im(x+1,y) +
		im(x-1,y) +	im(x-1,y+1) + im(x+1,y+1) ;
    }

    if ( (x==im.XSize()-1 ) &&  ( y>0 ) && ( y<im.YSize()-1 ) )
    // bord en x a "droite" (6 points)
    {
	som= im(x,y) +  im(x,y+1) + 
	       im(x,y-1) +	im(x-1,y) + im(x-1,y-1) +
	       im(x-1,y+1) ;

    }
    if ( (y==im.YSize()-1 ) &&  ( x>0 ) && ( x<im.XSize()-1 ) )
    // bord en y en bas (6 points )
    {
	som=  im(x,y) +	im(x,y-1) + im(x+1,y) +
		im(x-1,y) + im(x-1,y-1) +	im(x+1,y-1);
    }
    
    if ( (x==0) &&  ( y==0 ) )
    // coin en haut a "gauche" (4 points )
    {
	som= im(x,y) +im(x+1,y) + im(x,y+1) + im(x+1,y+1) ;
    }
    

    if ( (x==im.XSize()-1 ) &&  ( y==0 ) )
    // coin en haut a "droite" (4 points )
    {
	som= im(x,y) + im(x-1,y) + im(x-1,y+1) + im(x,y+1) ;
    }
		
    if ( (x==im.XSize()-1 ) &&  ( y==im.YSize()-1 ) )
    // coin en bas a "droite" (4 points )
    {
	som= im(x,y) + im(x-1,y) + im(x-1,y-1) + im(x,y-1) ;
    }
    if ( (x==0) &&  ( y==im.YSize()-1 ) )
    // coin en bas a "gauche" (4 points)
    {
	som= im(x,y) + im(x+1,y) + im(x+1,y-1) + im(x,y-1) ;
    }
    

}

/* ------Nouvelle fonction -----------*/
// detection de l'effet SZ
// la carte des sources detectees est mise dans nb_det (d'ou la tramsmission par reference)
// le nombre de sources detectees est mis dans nb (d'ou la tramsmission par reference)
// la fonction retourne une StarList

StarList * detec_SZ(ImageR4 & corin, ImageR4& nb_det,double& nb,int taille_f)
{
    ImageR4 * cor=NULL; 
    // On cree une liste d'etoiles pour stocker les sources detectees
    StarList* stl = new StarList(NULL, 1000, 250); 
    BStar* star;
    
    if (taille_f <= 1 ) cor = &corin;             // si l'on ne veut pas de filtrage il suffit d'appeler
                                                  // la fonction detec_SZ(..,.,.,1)
    else 
    { 
	ImageR4 filt_c(taille_f,taille_f);        // declaration et ...
	filt_c=1.;                                // initialisation du filtre 
	cor = FilterImage(&corin,cor,&filt_c);    // filtrage
    }

    nb=0;
    nb_det=0;
    
    
    double seuil=9/taille_f;                        // on garde 3 sigma pour chaque dimension x et y
    //    cout << "in seuil!!!" << seuil << endl;
    double taille=cor->XSize()*cor->YSize();

    bool tmp=false;
    bool& fgdet=tmp;
    float posx, posy, flux, fond, sompix, pixmax;   // parametre d'un etoile (source SZ)
    


    for (int i=0;i<cor->XSize();i++)    
    {
	for (int j=0;j<cor->YSize();j++)
	{
            fgdet = false;
	    
	    if ( (*cor)(i,j)>seuil )        // image filtree ou non filtree  
	    { 
		MaxImg(cor,i,j,fgdet);
		if ( fgdet )
		{
		    
		    nb=nb+1;          // calcul du nombre de sources
		    nb_det(i,j)=10;   // valeur arbitraire qui pourra etre remplacee par le y mesure par exemple
		    //fgdet = true;
		 
		}
	    }

	   
            if ( !fgdet ) continue;                     // creation de la liste d'etoiles
            posx = i+0.5;  posy = j+0.5;                      // on utilise les caracteristiques de la carte
	                                              // non filtree
            pixmax = corin(i,j);
	    SomPix(corin,i,j,sompix);
	    flux = sompix;
            fond = 0.;
	    star = stl->Next();
	    star->SetX(posx, posy, flux, fond, pixmax, sompix, 
                       posx, posy, kCourRef);
	    star->ClearFlag(BStar::flagAll, kCourRef);
	    star->SetFlag(BStar::flagOK, kCourRef);
	}
    }
    return(stl);
}