/*   mktransform    Create by  Juan Velasquez S                            */
/*                  email:  jvelasqu@dcc.uchile.cl                         */
/*                  University of  Chile                                   */
/*                                                                         */
/*                                                                         */
/* This program finds the transformation matrix between two frames.        */
/* Mktransform finds the object configurations in the frames that have     */
/* the best match and derives the transformation matrix.                   */
/* For this purpose, mktransform uses Stetson's method "Matching stars     */
/* between frames" based on triangle similarity of the brightest objects   */
/* ======================================================================= */
#include "focas1.h"

#define TOL		0.002   /* Default matching tolerance */
#define	NOBJS		20	/* Default number of objects */
#define MIN_MATCH	6	/* Min # of matched objects for transform */
#define MAX_MATCH	120	/* Max # of matches to use (~MAX_OBJS) */
#define MIN_OBJS	10	/* Min # of objects to use */
#define MAX_OBJS	100	/* Max # of objects to use */
#define	CLIP		3	/* Sigma clipping factor */
#define PM		57.2958 /* Radian to degree conversion */

static char *help[] = {
  " Finds the transformation matrix between two catalogs.",
  " Usage:    mktransform [-t tol] [-n nobj] target reference [filter]",
  " Argument: target    - catalog to find transformation and set coords",
  "           reference - catalog giving reference coordinates",
  "           -t  tol  - matching tolerance in triangle space (default 0.002)",
  "           -n  nobj - number of objects to use (default 20)",
  "           filter   - filter options",
  " Output:   The transformation matrix and reference coordinates are set",
  "           in the target catalog.  The standard output gives the",
  "           transformation coefficients and estimates of the offsets,",
  "           scales, and rotations.",
  0
};


/* Structure definitions */
typedef struct object {
	     float x, y; /* Object's position */
	     float m;    /* Object's magnitude */
	} OBJECT;
typedef struct triangulo { 
	    short  s1;	/* Object index at vertex between a and b */
	    short  s2;	/* Object index at vertex between a and c */
	    short  s3;	/* Object index at vertex between b and c */
	    float x, y;	/* Triangle space coordinate */
	} TRIANG;


/* Function Definitions */
void	Read_Sort(); 
int	Generate_List_Triang();
void	Found_Match();
void	Tol_Check();
void	Find_Transform();
void	Set_Catalog();
void	Chcoords();
int	Sides_Pos();
void	Quicksort();
void	Bin_Search();

int	DEBUG = 0;


main (argc, argv)
int	argc;
char	*argv[];
{
	int	nobjs;
	float	tolerance;
	float	xfrm[2][3];
	int	i, ntriang, argcat1, argcat2;
        int     ntriangA, ntriangB;
	TRIANG	*List_triangA, *List_triangB;
	OBJECT	*List1, *List2;

	/* Print usage help */
	if ((argc < 3) || (argv[1][0] == '^')) {
	    for (i = 0; help[i] != 0; i++)
		fprintf (stderr, "%s\n", help[i]);
	    exit (0);
	}

	/* Initialize */
	nobjs = NOBJS;
	tolerance = TOL;

	/* Get switch parameters */
	for (i = 1; (i < argc - 1) && (argv[i][0] == '-'); i++) {
	    switch (argv[i][1]) {
	    case 'd': /* Debug */
		DEBUG = 1;
		break;
	    case 't': /* Get tolerance */
		sscanf (argv[++i], "%f", &tolerance);              
		break;
	    case 'n': /* Get number of objects */
		sscanf (argv[++i], "%d", &nobjs);
		if (nobjs > MAX_OBJS) {
		    fprintf (stderr, "nobj too large: maximum is %d\n",
			MAX_OBJS);
		    exit (1);
		}
		if (nobjs < MIN_OBJS) {
		    fprintf (stderr, "nobj too small: minimum is %d\n",
			MIN_OBJS);
		    exit (2);
		}
		break;
	    default:
		for (i = 0; help[i] != 0; i++)
		    fprintf (stderr, "%s\n", help[i]);
		exit (0);
	    }
	}

	/* Set catalog indices */
	argcat1 = i++;
	argcat2 = i++;

	/* Set filter options */
	if (argc - i > 0)
	    setflt (argc-i, &argv[i]);

	/* Setup triangle lists */
	ntriang = (nobjs - 2) * (nobjs - 1) * nobjs / 6;
	List1 = (OBJECT *) malloc (nobjs * sizeof (OBJECT));
	List2 = (OBJECT *) malloc (nobjs * sizeof (OBJECT));
	List_triangA = (TRIANG *)malloc (ntriang * sizeof (TRIANG));
	List_triangB = (TRIANG *)malloc (ntriang * sizeof (TRIANG));
	for (i = 0; i < nobjs; i++)
	    List1[i].m = List2[i].m = 1000;

	/* Read the catalogs for the brightest objects */
	Read_Sort (argv[argcat1], argv[argcat2], nobjs, List1, List2);

	/* Make the triangle lists */
        ntriangA = Generate_List_Triang (nobjs, List1, List_triangA);
	ntriangB = Generate_List_Triang (nobjs, List2, List_triangB);
	if (DEBUG)
	    printf(" ntriangA = %d, ntriangB = %d\n",ntriangA, ntriangB);

	/* Match the triangles and compute the transformation */
	Found_Match (nobjs, ntriangA, ntriangB, tolerance,
	    List_triangA, List_triangB, List1, List2, xfrm);

	/* Set the catalog and print output information */
	Set_Catalog (argv[argcat1], xfrm, argc, argv);
}


/* ==================================================================== */
/*  Read_Sort: Read the catalogs and sort by magnitude			*/
/* ==================================================================== */

void Read_Sort (target, reference, nobjs, List1, List2)
char	*target, *reference;
int	nobjs;
OBJECT	*List1, *List2; 
{
	float	x, y, m, tempm, tempx, tempy;
	int	i;
	struct objrec ob;

	/* Read and sort the target catalog. */
	/* Only non-multiple objects matching the filter are selected. */

	if ((cfd = catopen (target, 0)) == -1) {
	    fprintf (stderr, "mktransform: can't open catalog %s\n", target);
	    exit (0);
	}
	     
	for (;;) {
	    if (rdcatob (cfd, 0L, &ob))
		break;
	    if (ffilter (&ob) !=1)
		continue;
	    if (ob.class[0] == 'm')
	    continue;
	      
	    x = ob.xc;
	    y = ob.yc;
	    m = ob.mag;

	    for (i=0; i < nobjs; i++)
		if  (m < List1[i].m) {
		    tempm = List1[i].m;
		    tempx = List1[i].x;
		    tempy = List1[i].y;
		    List1[i].m = m;
		    List1[i].x = x;
		    List1[i].y = y;
		    m = tempm;
		    x = tempx;
		    y = tempy;
		    break;
		}

	    for (i=i+1; i < nobjs; i++) {
		tempm = List1[i].m;
		tempx = List1[i].x;
		tempy = List1[i].y;
		List1[i].m = m;
		List1[i].x = x;
		List1[i].y = y;
		m = tempm;
		x = tempx;
		y = tempy;
	    }
	}

	fclose (cfd);

	/* Read and sort the reference catalog. */
	/* Only non-multiple objects are selected. */

	if ((cfd = catopen (reference, 0)) == -1) {
	    fprintf (stderr, "mktransform: can't open catalog %s\n", reference);
	    exit (0);
	}

	for (;;) {
	    if (rdcatob (cfd, 0L, &ob))
		break;
	    if (ob.class[0] == 'm')
		continue;
	      
	    x =  ob.ra;
	    y =  ob.dec;
	    m =  ob.mag;

	    for (i=0; i < nobjs; i++)
		if  (m < List2[i].m) {
		    tempm = List2[i].m;
		    tempx = List2[i].x;
		    tempy = List2[i].y;
		    List2[i].m  = m;
		    List2[i].x  = x;
		    List2[i].y  = y;
		    m = tempm;
		    x = tempx;
		    y = tempy;
		    break;
		}

	    for (i=i+1; i < nobjs; i++) {
		tempm = List2[i].m;
		tempx = List2[i].x;
		tempy = List2[i].y;
		List2[i].m  = m;
		List2[i].x  = x;
		List2[i].y  = y;
		m = tempm;
		x = tempx;
		y = tempy;
	    }
	}
	fclose (cfd);
}


/* ========================================================================== */
/* Rutina Generate_list_triang : Se encarga de generar las listas de          */
/* triangulos  a partir de la cantidad de estrellas que se consideraran de    */
/* ambas imagenes.                                                            */
/* ========================================================================== */

int Generate_List_Triang (nobjs, List, List_triang)
int nobjs;
OBJECT  *List;
TRIANG  *List_triang;
{
	int	i, j, k, p=0;
	float	di, dj, dk, a, b, c, h1, h2; 
        float   *sides;

	/* Minimize the number of triangle sides to compute */
	sides = (float *) calloc (nobjs * (nobjs - 1)/2 + 1, sizeof(float));
	for (i = 0; i < nobjs; i++) {
	    for (j = i + 1; j < nobjs; j++) {
		h1 = (List[i].x - List[j].x);
		h2 = (List[i].y - List[j].y);
		k = Sides_Pos (i, j, nobjs);
		sides[k] = sqrt(h1*h1 + h2*h2);
	    }
	}

	/* Order the triangle sides and compute the triangle space coords */
	for (i = 0; i < (nobjs-2); i++) {
	    for (j = i+1; j < (nobjs-1); j++) {
		for (k = j+1; k < nobjs; k++) {
		    di = sides [Sides_Pos(i, j, nobjs)] ;
		    dj = sides [Sides_Pos(j, k, nobjs)] ;
		    dk = sides [Sides_Pos(k, i, nobjs)] ;

		    if (dk > dj) {
			if (dk > di) {
			    if (dj > di) {	/* kji */
				a = dk; b = dj; c = di;
				List_triang[p].s1 = k;
				List_triang[p].s2 = i;
				List_triang[p].s3 = j;
			    } else {		/* kij */
				a = dk; b = di; c = dj;
				List_triang[p].s1 = i;
				List_triang[p].s2 = k;
				List_triang[p].s3 = j;
			    }
			} else {		/* ikj */
			    a = di; b = dk; c = dj;
			    List_triang[p].s1 = i;
			    List_triang[p].s2 = j;
			    List_triang[p].s3 = k;
			}
		    } else if (dj > di) {
			if (di > dk) {		/* jik */
			    a = dj; b = di; c = dk;
			    List_triang[p].s1 = j;
			    List_triang[p].s2 = k;
			    List_triang[p].s3 = i;
			} else {		/* jki */
			    a = dj; b = dk; c = di;
			    List_triang[p].s1 = k;
			    List_triang[p].s2 = j;
			    List_triang[p].s3 = i;
			}
		    } else {			/* ijk */
			a = di; b = dj; c = dk;
			List_triang[p].s1 = j;
			List_triang[p].s2 = i;
			List_triang[p].s3 = k;
		    }

		    /* Include only triangles with b/a < 0.9 */
		    if ( b/a < 0.9 ) {
			List_triang[p].x = b/a;
			List_triang[p].y = c/a;
			p++;
		    }
		}
	    }
	}

	free(sides);
	return (p);
}


/* ======================================================================== */
/*  Found_Match: Find the objects with the best match between two catalogs  */
/*  Compute the transformation matrix and set the catalog header and        */
/*  reference coordinates.  Output the results.                             */
/* =======================================================================  */

void Found_Match (nobjs, ntriangA, ntriangB, tolerance, List_triangA,
	List_triangB, List1, List2, xfrm)
int	nobjs, ntriangA, ntriangB;
float	tolerance;
TRIANG	*List_triangA, *List_triangB;
OBJECT	*List1, *List2;
float	xfrm[2][3];
{
	short	matches[3][MAX_MATCH];
	int	i, j, k=0, l, n, first, last, *Table_match;
        float	tolerance2 = tolerance * tolerance;

	/* Allocate memory for the match table. */
      	Table_match = (int *) calloc (nobjs * nobjs, sizeof(int));
	if( Table_match == NULL) {
	    printf("  Error: Can't create Table_match array.\n");
	    exit(0);
	}
    
	/* Sort List_triangB by x coordinate. */
	Quicksort (List_triangB, 0, ntriangB - 1);
		   
	/* Find objects within tolerance distance in triangle space. */
	for (i = 0; i < ntriangA; i++) {
            Bin_Search (List_triangA[i].x, List_triangB, &first, &last, 
                tolerance, ntriangB);
	    Tol_Check (nobjs, tolerance, tolerance2, List_triangA[i],
		List_triangB, first, last, Table_match);
	}

	if (DEBUG) {
	    n = 0;
	    for (i=0; i<nobjs*nobjs; i++)
		n += Table_match[i];
	    printf ("Number of trianges matched = %d\n", n/3);
	}

	/* Find the nobjs points with the most matches. */
	k = 0;
	for (i = 0; i < nobjs; i++) {
	    for (j = 0; j < nobjs; j++) {
		n = Table_match[i*nobjs+j];
		if (n > 0) {
		    if (k < nobjs) {
			for (l=k; l>0 && n>matches[2][l-1]; l--) {
			    matches[0][l] = matches[0][l-1];
			    matches[1][l] = matches[1][l-1];
			    matches[2][l] = matches[2][l-1];
			}
			matches[0][l] = i;
			matches[1][l] = j;
			matches[2][l] = n;
			k++;
		    } else if (n >= matches[2][nobjs-1]) {
			for (l=k; l>0 && n>matches[2][l-1]; l--) {
			    matches[0][l] = matches[0][l-1];
			    matches[1][l] = matches[1][l-1];
			    matches[2][l] = matches[2][l-1];
			}
			matches[0][l] = i;
			matches[1][l] = j;
			matches[2][l] = n;
			l = k < MAX_MATCH ? k + 1: k;
			n = matches[2][nobjs-1];
			for (k=nobjs; k<l && n==matches[2][k]; k++);
		    }
		}

		/*
		if (DEBUG) {
		    printf ("%3d : %6.2f %6.2f : %6.2f %6.2f\n",
			Table_match[i*nobjs+j],
			List1[i].x, List1[i].y, List2[j].x, List2[j].y);
		}
		*/
		if (DEBUG)
		    printf ("%2d ", Table_match[i*nobjs+j]);
	    }
	    if (DEBUG)
		printf ("\n");
	}

	if (DEBUG) {
	    for (i = 0; i < k; i++)
		printf ("%3d %3d %4d\n", matches[1][i]+1, matches[0][i]+1,
		    matches[2][i]);
	}

	/* Compute transformation matrix from the nobjs best matches. */
	Find_Transform (matches, k, List1, List2, xfrm);
}


/* ============================================================	*/
/* Find triangles within specified tolerance distance.		*/
/* ============================================================	*/

void Tol_Check (nobjs, tolerance, tolerance2, List_triangA,
	List_triangB, first, last, Table_match)
int	nobjs;
float	tolerance, tolerance2;
TRIANG	List_triangA, *List_triangB;
int	first, last;
int	*Table_match;
{
	float	temp1, temp2, distance;
	int	i, h1, h2;

	for (i = first ; i <= last; i++) {
	    temp2 = (List_triangA.y - List_triangB[i].y);
	    if (temp2 < tolerance) {
		temp1 = (List_triangA.x - List_triangB[i].x);
		distance = temp1 * temp1 + temp2 * temp2;
		if (distance < tolerance2) {
		    h1 = List_triangA.s1;
		    h2 = List_triangB[i].s1;
		    Table_match[h1*nobjs+h2]++;

		    h1 = List_triangA.s2;
		    h2 = List_triangB[i].s2;
		    Table_match[h1*nobjs+h2]++;

		    h1 = List_triangA.s3;
		    h2 = List_triangB[i].s3;
		    Table_match[h1*nobjs+h2]++;
		}
	    }
	}
}


/* ============================================================ */
/* Compute transformation matrix with iterative rejection	*/
/* ============================================================ */

void Find_Transform (matches, m, List1, List2, xfrm)
short	matches[3][MAX_MATCH];
int	m;
OBJECT	*List1, *List2;
float	xfrm[2][3];
{
	int	i, j, i1, i2;
	int	mda = MAX_MATCH, mdb = MAX_MATCH, n = 3, nb = 2;
	int	krank, ip[3];
	float	tau = 0.1, rnorm[2], h[3], g[3];
	float	a[3][MAX_MATCH], b[2][MAX_MATCH];
	float	x, y, r2, sum, rms;

	/* Require a minimum number of points. */
	if (m < MIN_MATCH) {
	    printf ("Match not found: use more objects or larger tolerance\n");
	    exit (0);
	}

	/* Compute the initial transformation with the 12 best matches. */
	j = (m < 12) ? m : 12;
	for (i=0; i<j; i++) {
	    a[0][i] = List1[matches[0][i]].x;
	    a[1][i] = List1[matches[0][i]].y;
	    a[2][i] = 1.;
	    b[0][i] = List2[matches[1][i]].x;
	    b[1][i] = List2[matches[1][i]].y;
	}
#ifndef NOUS
	hfti_ (a, &mda, &j, &n, b, &mdb, &nb, &tau, &krank, rnorm, h, g, ip);
#else
	hfti (a, &mda, &j, &n, b, &mdb, &nb, &tau, &krank, rnorm, h, g, ip);
#endif
	for (i=0; i<2; i++)
	    for (j=0; j<3; j++)
		xfrm[i][j] = b[i][j];

	/* Start with all matches compute RMS and reject outliers.	      */
	/* The outliers are found from the 60% point in the sorted residuals. */
	for (;;) {
	    sum = 0.;
	    for (i=0; i<m; i++) {
		i1 = matches[0][i];
		i2 = matches[1][i];
		r2 = List1[i1].x;
		y = List1[i1].y;
		x = xfrm[0][0] * r2 + xfrm[0][1] * y + xfrm[0][2];
		y = xfrm[1][0] * r2 + xfrm[1][1] * y + xfrm[1][2];
		x -= List2[i2].x;
		y -= List2[i2].y;
		r2 = x * x + y * y;
		for (j=i; j>0 && r2<a[1][j-1]; j--)
		    a[1][j] = a[1][j-1];
		a[0][i] = r2;
		a[1][j] = r2;
		sum += r2;
	    }

	    /* Set clipping limit and quit when no points are clipped. */
	    i = 0.6 * m;
	    r2 = CLIP * a[1][i];
	    if (r2 >= a[1][m-1]) {
		rms = sqrt (sum / m);
		break;
	    }

	    /* Clip outliers and redo the fit. */
	    j = 0;
	    for (i=0; i<m; i++) {
		if (a[0][i] < r2) {
		    i1 = matches[0][i];
		    i2 = matches[1][i];
		    matches[0][j] = i1;
		    matches[1][j] = i2;
		    a[0][j] = List1[i1].x;
		    a[1][j] = List1[i1].y;
		    a[2][j] = 1.;
		    b[0][j] = List2[i2].x;
		    b[1][j] = List2[i2].y;
		    j++;
		}
	    }
	    m = j;

	    if (m < MIN_MATCH) {
		printf (
		    "Match not found: use more objects or larger tolerance\n");
		exit (0);
	    }

#ifndef NOUS
	    hfti_ (a, &mda, &m, &n, b, &mdb, &nb, &tau, &krank, rnorm, h, g,ip);
#else
	    hfti (a, &mda, &m, &n, b, &mdb, &nb, &tau, &krank, rnorm, h, g,ip);
#endif
	    for (i=0; i<2; i++)
		for (j=0; j<3; j++)
		    xfrm[i][j] = b[i][j];
	}

	printf ("Number of matches = %d, RMS of fit = %8.2f\n", m, rms);
}


/* ========================================================================= */
/* Set_Catalog - Set the catalog transformation and reference coordinates.   */
/* Also print diagnostic information to the standard output.                 */
/* ========================================================================= */

void Set_Catalog (target, xfrm, argc, argv)
char	*target;
float	xfrm[2][3];
int	argc;
char	*argv[];
{
	int	i, j;
	float	rotation, scale;
	float	A, B, C, D, E, F, angle1, angle2, scalex, scaley;

	/* Open catalog */
	cfd = catopen (target, 2);

	/* Set transformation matrix and reference coordinates */
	for (i = 0; i < 2; i++)
	    for (j = 0; j < 3; j++)
		sp.xfrm[i][j] = xfrm[i][j];
	for (j = 0; j < 3; j++)
	    sp.xfrm[2][j] = 0.;

	Chcoords ();

	/* Update catalog header and close */
	cmmnt (argc, argv);
	wtcathdr (cfd, 0);
	fclose (cfd);

	/* Output results */
	A = xfrm[0][0];
	B = xfrm[0][1];
	C = xfrm[0][2];
	D = xfrm[1][0];
	E = xfrm[1][1];
	F = xfrm[1][2];

	printf ("Transformation coefficients: \n"); 
	printf ("          %10f  %10f %10f\n", A, B, C);
	printf ("          %10f  %10f %10f\n", D, E, F);

	scalex = sqrt (A * A + B * B);
	scaley = sqrt (D * D + E * E);
	scale = (scalex + scaley) / 2;
	angle1 = atan2 (B, A)  * PM;
	angle2 = atan2 (-D, E) * PM;

	if (angle1 - angle2 > 90.)
	    rotation = (angle1 + angle2 + 180.)/2;
	else if (angle1 - angle2 < -90.)
	    rotation = (angle1 + angle2 - 180.)/2;
	else
	    rotation = (angle1 + angle2)/2;

printf ("\n");
printf ("                         Xtransform   Ytransform      Average\n");
printf ("X Offset (pixels):                                 %10.2f\n", C);
printf ("Y Offset (pixels):                                 %10.2f\n", F);
printf ("Scale:                   %10.4f   %10.4f   %10.4f\n", scalex, scaley, scale);
printf ("Rotation (degrees):      %10.2f   %10.2f   %10.2f\n",
    angle1, angle2,rotation);
if (fabs (angle1 - angle2) > 90.)
    printf ("Flip:                                                     YES\n");
else
    printf ("Flip:                                                      NO\n");
}


/* ============================================================	*/
/* CHCOORDS - Set the reference coordinates			*/
/* ============================================================	*/

void Chcoords () {
	long    pos, lseek();
	struct	objrec ob;

	for (;;) {
	    pos = lseek (cfd, 0L, 1);
	    if (rdcatob (cfd, 0L, &ob))
		break;
	    coords (&ob.xc, &ob.ra);
	    lseek (cfd, pos, 0);
	    wtcat (&ob, cfd);
	}
}


/* ============================================================ */
/* Routine to select index in array of side lengths		*/
/* ============================================================ */

int Sides_Pos (i, j, n)
int	i, j, n;
{
	if (i < j)
	    return(i*(2*n-i-3)/2 + j);
	else
	    return(j*(2*n-j-3)/2 + i);
}


/* ============================================================	*/
/*  Bin_Search : Binary Search in a array order by x coordinate	*/
/* ============================================================	*/

void Bin_Search (key, List_triang, first, last, tolerance, ntriang)
float   key;
TRIANG	*List_triang;
int	*first, *last;
float	tolerance;
int	ntriang;
{
	int	min = 0, max = ntriang - 1, middle ;
	int	found = 0, i;  

	while ((!found) && (max - min > 1)) {
	    middle  = (min + max ) / 2;
	    if (fabs(List_triang[middle].x - key) < tolerance)
		found = 1;
	    else if ( key < ( List_triang[middle ].x - tolerance))
		max = middle ;
	    else if ( key > (List_triang[middle ].x + tolerance))
		min = middle;
	}

	/* Not found */
	if (!found) {
	    *first = 2;
	    *last  = 1;
	    return;
	}

	for (i = middle; i > 0; i--) {
	    if (fabs (List_triang[i].x - key) > tolerance)
		break;
	} 
	*first = i;

	for (i = middle; i < ntriang-1; i++) {
	    if (fabs (List_triang[i].x - key) > tolerance)
		break;
	} 
	*last = i;

}


/* ============================================================== */
/*  Quicksort : Sort by x coordinate                              */
/* ============================================================== */

void Quicksort (List_Triang, l, r)
TRIANG	*List_Triang;
int	l, r;
{
	TRIANG	v, t;
	int	i, j;

	if( r > l ) {
	    v.x = List_Triang[r].x;
	    i = l-1; 
	    j = r;

	    for(;;) {
		while(List_Triang[++i].x < v.x );
		while((j > 1) && (List_Triang[--j].x > v.x));
		if(i >= j)
		    break;

		t.x = List_Triang[i].x; 
		List_Triang[i].x = List_Triang[j].x;
		List_Triang[j].x = t.x;

		t.y = List_Triang[i].y; 
		List_Triang[i].y = List_Triang[j].y;
		List_Triang[j].y = t.y;

		t.s1 = List_Triang[i].s1; 
		List_Triang[i].s1 = List_Triang[j].s1;
		List_Triang[j].s1 = t.s1;

		t.s2 = List_Triang[i].s2; 
		List_Triang[i].s2 = List_Triang[j].s2;
		List_Triang[j].s2 = t.s2;

		t.s3 = List_Triang[i].s3; 
		List_Triang[i].s3 = List_Triang[j].s3;
		List_Triang[j].s3 = t.s3;
	    }

	    t.x = List_Triang[i].x; 
	    List_Triang[i].x = List_Triang[r].x;
	    List_Triang[r].x = t.x;

	    t.y = List_Triang[i].y; 
	    List_Triang[i].y = List_Triang[r].y;
	    List_Triang[r].y = t.y;

	    t.s1= List_Triang[i].s1; 
	    List_Triang[i].s1 = List_Triang[r].s1;
	    List_Triang[r].s1 = t.s1;	

	    t.s2 = List_Triang[i].s2; 
	    List_Triang[i].s2 = List_Triang[r].s2;
	    List_Triang[r].s2 = t.s2;

	    t.s3 = List_Triang[i].s3; 
	    List_Triang[i].s3 = List_Triang[r].s3;
	    List_Triang[r].s3 = t.s3;

	    Quicksort (List_Triang, l, i-1);
	    Quicksort (List_Triang, i+1, r);
	}
}