/***************************************************************************
sumagts.c - description
-------------------
begin : Mon Jun 14 2004
copyright : (C) 2004 by Jakub Otwinowski
email : jakub@hurin
***************************************************************************/
/***************************************************************************
* *
* This program 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 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
#ifdef HAVE_CONFIG_H
#include "gts/config.h"
#endif
#include "SUMA_suma.h"
#include "SUMA_gts.h"
extern SUMA_CommonFields *SUMAg_CF;
extern SUMA_DO *SUMAg_DOv;
extern SUMA_SurfaceViewer *SUMAg_SVv;
extern int SUMAg_N_SVv;
extern int SUMAg_N_DOv;
#if 0
/* Normally, functions in SUMA_gts_insert.c would be here */
#endif
GtsSurface* SumaToGts( SUMA_SurfaceObject *SO)
{
static char FuncName[]={"SumaToGts"};
GtsSurface* s = NULL;
GtsVertex ** vertices = NULL;
GtsEdge ** edges = NULL;
int i = 0; /*counters */
int n = 0;
int *MapToGts = NULL; /*used to map EL vector to edge vector*/
SUMA_Boolean LocalHead = NOPE;
SUMA_ENTRY;
if (!SO->EL) {
SUMA_SL_Err("Null Edgeist!");
SUMA_RETURN(s);
}
SUMA_LH("In with ");
if (LocalHead) SUMA_Print_Surface_Object(SO, stderr);
s = gts_surface_new( gts_surface_class (),
gts_face_class (),
gts_edge_class (),
gts_vertex_class ());
vertices = (GtsVertex **)g_malloc (SO->N_Node * sizeof (GtsVertex *));
n = 0;
for ( i=0; i< SO->N_Node*3; i+=3)
{
vertices[n] = gts_vertex_new( s->vertex_class,
(gdouble)SO->NodeList[i],
(gdouble)SO->NodeList[i+1],
(gdouble)SO->NodeList[i+2]);
if (LocalHead)
fprintf(SUMA_STDERR, "Added vertex (%d/%d)%f %f %f\n"
"Stored gts (fails for some reason, but surface is OK) %f %f %f\n"
, n, SO->N_Node-1, SO->NodeList[i], SO->NodeList[i+1], SO->NodeList[i+2]
, (float)(vertices[n]->p.x), (float)(vertices[n]->p.y), (float)(vertices[n]->p.z));
++n;
}
edges = (GtsEdge**)g_malloc (SO->EL->N_Distinct_Edges * sizeof (GtsEdge*));
n = 0;
MapToGts = (int *)g_malloc ( SO->EL->N_EL * sizeof(int));
for ( i=0; i< SO->EL->N_EL; i++)
{
if (SO->EL->ELps[i][2] > 0)
{ /* a unique edge*/
GtsVertex* v1 = vertices[SO->EL->EL[i][0]];
GtsVertex* v2 = vertices[SO->EL->EL[i][1]];
if (SO->EL->ELps[i][0] == 1) {
edges[n++] = gts_edge_new ( s->edge_class, v2, v1 );
if (LocalHead) fprintf(SUMA_STDERR,"Added edge of vertices: %d %d \n", SO->EL->EL[i][1], SO->EL->EL[i][0]);
} else {
edges[n++] = gts_edge_new ( s->edge_class, v1, v2 );
if (LocalHead) fprintf(SUMA_STDERR,"Added edge of nodes: %d %d \n", SO->EL->EL[i][0], SO->EL->EL[i][1]);
}
}
MapToGts[i] = n-1; /* n-1 bc n was just incremented */
if (LocalHead) fprintf(SUMA_STDERR,"SUMA edge %d is mapped to GTS edge %d\n", i, n-1);
}
if (n != SO->EL->N_Distinct_Edges)
{
fprintf(SUMA_STDERR, "distinct edges didn't equal N_Distinct_Edges");
exit(1);
}
/* this loop isn't working, stupid tri_limb
for ( i=0; i< SO->N_FaceSet; i++)
{
gts_surface_add_face(s,
gts_face_new ( s->face_class,
edges[MapToGts[SO->EL->Tri_limb[i][0]]],
edges[MapToGts[SO->EL->Tri_limb[i][1]]],
edges[MapToGts[SO->EL->Tri_limb[i][2]]]));
}
*/
for ( i=0; i< SO->N_FaceSet*3; i+=3)
{
int n1 = SO->FaceSetList[i];
int n2 = SO->FaceSetList[i+1];
int n3 = SO->FaceSetList[i+2];
GtsEdge* e1 = edges[MapToGts[SUMA_FindEdge( SO->EL, n1, n2)]];
GtsEdge* e2 = edges[MapToGts[SUMA_FindEdge( SO->EL, n2, n3)]];
GtsEdge* e3 = edges[MapToGts[SUMA_FindEdge( SO->EL, n3, n1)]];
gts_surface_add_face(s,
gts_face_new ( s->face_class, e1, e2, e3));
if (LocalHead) fprintf(SUMA_STDERR, "Added face of SUMA edges: %d %d %d\n"
" MapToGTS : %d %d %d\n",
SUMA_FindEdge( SO->EL, n1, n2), SUMA_FindEdge( SO->EL, n2, n3), SUMA_FindEdge( SO->EL, n3, n1),
MapToGts[SUMA_FindEdge( SO->EL, n1, n2)], MapToGts[SUMA_FindEdge( SO->EL, n2, n3)], MapToGts[SUMA_FindEdge( SO->EL, n3, n1)]);
}
g_free (vertices);
g_free (edges);
g_free (MapToGts);
SUMA_RETURN( s );
}
#if 0
/* These functions fail on OSX, see function gts_surface_suma and comments preceding it */
SUMA_SurfaceObject* GtsToSuma( GtsSurface *gs)
{
static char FuncName[]={"GtsToSuma"};
GHashTable *hash = g_hash_table_new(NULL,NULL);
SUMA_SurfaceObject* SO = NULL;
int i = 0;
gpointer data1[3];
SUMA_ENTRY;
GTS_OUT("shit.surf", gs);
SO = (SUMA_SurfaceObject *)SUMA_malloc(sizeof(SUMA_SurfaceObject));
SO->N_Node = gts_surface_vertex_number(gs);
SO->NodeList = (float*) SUMA_calloc(SO->N_Node * 3, sizeof(float));
SO->N_FaceSet = gts_surface_face_number(gs);
SO->FaceSetList = (int*) SUMA_calloc(SO->N_FaceSet * 3, sizeof(int));
/*
foreach vertex
get xyz,
make hash key=pvertex value=integer
make NodeList from integer and xyz
foreach face
get 3 vertexes
find their integers from hash
make FaceSetList from integers
*/
data1[0] = hash;
data1[1] = SO;
data1[2] = &i;
i = 0;
gts_surface_foreach_vertex( gs, (GtsFunc) MakeNodeList_foreach_vertex, data1);
i = 0;
gts_surface_foreach_face( gs, (GtsFunc) MakeFaceList_foreach_face, data1);
g_hash_table_destroy(hash);
SUMA_RETURN(SO);
}
void MakeNodeList_foreach_vertex ( GtsPoint *p, gpointer* data)
{
GHashTable *hash = data[0];
SUMA_SurfaceObject* SO = data[1];
int* i = data[2];
fprintf(SUMA_STDERR,"this %f %f %f\n", p->x, p->y, p->z);
SO->NodeList[(*i)*3] = p->x;
SO->NodeList[(*i)*3+1] = p->y;
SO->NodeList[(*i)*3+2] = p->z;
if (*i >= SO->N_Node || *i < 0)
{
fprintf(SUMA_STDERR, "node %i out of range", *i);
exit(1);
}
if (g_hash_table_lookup(hash, p) == NULL)
g_hash_table_insert(hash, p, GINT_TO_POINTER(*i));
else
{
fprintf(SUMA_STDERR, "something already in hash table??");
exit(1);
}
(*i)++;
}
void MakeFaceList_foreach_face ( GtsFace* f, gpointer* data)
{
GHashTable *hash = data[0];
SUMA_SurfaceObject* SO = data[1];
int* i = data[2];
GtsVertex *v1,*v2,*v3;
gpointer presult = NULL;
int iresult = 0;
gpointer temp = NULL;
gts_triangle_vertices(&(f->triangle), &v1, &v2, &v3);
if ((*i) >= SO->N_FaceSet)
{
fprintf(SUMA_STDERR, "counter %i passed N_FaceSet %i",(*i),SO->N_FaceSet);
exit(1);
}
if (!g_hash_table_lookup_extended(hash, v1, temp, &presult))
{
fprintf(SUMA_STDERR, "hash table lookup failed");
exit(1);
}
iresult = GPOINTER_TO_INT (presult);
if (iresult >= SO->N_Node || iresult < 0)
{
fprintf(SUMA_STDERR, "node %i out of range", iresult);
exit(1);
}
SO->FaceSetList[(*i)*3] = iresult;
if (!g_hash_table_lookup_extended(hash, v2, temp, &presult))
{
fprintf(SUMA_STDERR, "hash table lookup failed");
exit(1);
}
iresult = GPOINTER_TO_INT (presult);
if (iresult >= SO->N_Node || iresult < 0)
{
fprintf(SUMA_STDERR, "node %i out of range", iresult);
exit(1);
}
SO->FaceSetList[(*i)*3+1] = iresult;
if (!g_hash_table_lookup_extended(hash, v3, temp, &presult))
{
fprintf(SUMA_STDERR, "hash table lookup failed");
exit(1);
}
iresult = GPOINTER_TO_INT (presult);
if (iresult >= SO->N_Node || iresult < 0)
{
fprintf(SUMA_STDERR, "node %i out of range", iresult);
exit(1);
}
SO->FaceSetList[(*i)*3+2] = iresult;
(*i)++;
}
#endif
void coarsen( GtsSurface* s, int stop)
{
GtsVolumeOptimizedParams params = { 0.5, 0.5, 0. };
gdouble fold = PI/180.;
gts_surface_coarsen (s,
(GtsKeyFunc) gts_volume_optimized_cost,
¶ms,
(GtsCoarsenFunc) gts_volume_optimized_vertex,
¶ms,
(GtsStopFunc) gts_coarsen_stop_number,
&stop,
fold);
}
gboolean stop_number (gdouble cost, guint number, guint * max)
{
if (number > *max)
return TRUE;
return FALSE;
}
void refine( GtsSurface* s, int stop)
{
gts_surface_refine(s, NULL, NULL, NULL, NULL,
(GtsStopFunc)stop_number, &stop);
}
/*!
\brief resample a mesh so that the resultant surface has edge_factor as many edges as the original one
*/
SUMA_SurfaceObject *SUMA_Mesh_Resample (SUMA_SurfaceObject *SO, float edge_factor)
{
static char FuncName[]={"SUMA_Mesh_Resample"};
SUMA_SurfaceObject *S2=NULL;
GtsSurface* s = NULL;
SUMA_Boolean LocalHead = NOPE;
SUMA_ENTRY;
if (!SO) {
SUMA_SL_Err("NULL SO");
SUMA_RETURN(S2);
}
if (!SO->EL) {
SUMA_S_Warn("NULL Edge List, computing it");
if (!SUMA_SurfaceMetrics(SO, "EdgeList", NULL)) {
SUMA_SL_Err("Failed to create EdgeList");
SUMA_RETURN(S2);
}
}
/* create a GTS surface */
s = SumaToGts(SO);
if (!s) {
SUMA_S_Err("Failed to change SO to GTS surface");
SUMA_RETURN(S2);
}
if (LocalHead) { /* see if surface was created well in GTS format */
SUMA_S_Note("Writing initial surface in GTS form");
GTS_OUT("gtsout.surf",s);
GTS_VTK_OUT("vtkout.surf",s);
}
/* resample */
if (1) {
SUMA_S_Note("Changing mesh density\n");
if (edge_factor < 1)
coarsen(s, SO->EL->N_Distinct_Edges * edge_factor);
else
refine(s, SO->EL->N_Distinct_Edges * edge_factor);
} else {
SUMA_S_Note("Leaving surface untouched\n");
}
if (!s) {
SUMA_SL_Err("Failed to refine");
SUMA_RETURN(S2);
}
/* change the surface back to SUMA */
S2 = SUMA_Alloc_SurfObject_Struct(1);
gts_surface_suma (s,
&(S2->NodeList), &(S2->N_Node), &(S2->NodeDim),
&(S2->FaceSetList), &(S2->N_FaceSet), &(S2->FaceSetDim));
gts_object_destroy((GtsObject*)s); s = NULL;
SUMA_RETURN(S2);
}
syntax highlighted by Code2HTML, v. 0.9.1