00001 /***************************************************************************
00002                           sumagts.c  -  description
00003                              -------------------
00004     begin                : Mon Jun 14 2004
00005     copyright            : (C) 2004 by Jakub Otwinowski
00006     email                : jakub@hurin
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017 #ifdef HAVE_CONFIG_H
00018 #include "gts/config.h"
00019 #endif
00020 
00021 #include "SUMA_suma.h"
00022 #include "SUMA_gts.h"
00023 
00024 extern SUMA_CommonFields *SUMAg_CF;
00025 extern SUMA_DO *SUMAg_DOv;
00026 extern SUMA_SurfaceViewer *SUMAg_SVv;
00027 extern int SUMAg_N_SVv; 
00028 extern int SUMAg_N_DOv;  
00029 
00030 #if 0
00031 /* Normally, functions in SUMA_gts_insert.c would be here */
00032 #endif 
00033 
00034 GtsSurface* SumaToGts( SUMA_SurfaceObject *SO)
00035 {
00036         static char FuncName[]={"SumaToGts"};
00037         GtsSurface* s = NULL;
00038         GtsVertex ** vertices = NULL;
00039         GtsEdge ** edges = NULL;
00040         int i = 0; /*counters */
00041         int n = 0;
00042         int *MapToGts = NULL; /*used to map EL vector to edge vector*/
00043    SUMA_Boolean LocalHead = NOPE;
00044    
00045    SUMA_ENTRY;
00046    
00047         if (!SO->EL) {
00048       SUMA_SL_Err("Null Edgeist!");
00049       SUMA_RETURN(s);
00050    }
00051    
00052    SUMA_LH("In with ");
00053    if (LocalHead) SUMA_Print_Surface_Object(SO, stderr);
00054    
00055    s = gts_surface_new( gts_surface_class (),
00056                                                                          gts_face_class (),
00057                                                                          gts_edge_class (),
00058                                                                          gts_vertex_class ());
00059         vertices = (GtsVertex **)g_malloc (SO->N_Node * sizeof (GtsVertex *));
00060         n = 0;
00061    for ( i=0; i< SO->N_Node*3; i+=3)
00062         {
00063                 vertices[n] = gts_vertex_new( s->vertex_class,
00064                         (gdouble)SO->NodeList[i],
00065                         (gdouble)SO->NodeList[i+1],
00066                         (gdouble)SO->NodeList[i+2]);
00067       if (LocalHead) 
00068          fprintf(SUMA_STDERR, "Added vertex (%d/%d)%f %f %f\n"
00069                                           "Stored gts (fails for some reason, but surface is OK)  %f %f %f\n"
00070                                           , n, SO->N_Node-1, SO->NodeList[i], SO->NodeList[i+1], SO->NodeList[i+2]
00071                                           , (float)(vertices[n]->p.x), (float)(vertices[n]->p.y), (float)(vertices[n]->p.z));  
00072            ++n;
00073    }
00074    edges = (GtsEdge**)g_malloc (SO->EL->N_Distinct_Edges * sizeof (GtsEdge*));
00075         n = 0;
00076         MapToGts = (int *)g_malloc ( SO->EL->N_EL * sizeof(int));
00077    for ( i=0; i< SO->EL->N_EL; i++)
00078         {
00079                 if (SO->EL->ELps[i][2] > 0)
00080                 { /* a unique edge*/
00081                         GtsVertex* v1 = vertices[SO->EL->EL[i][0]];
00082                         GtsVertex* v2 = vertices[SO->EL->EL[i][1]];
00083                         if (SO->EL->ELps[i][0] == 1) {
00084                                 edges[n++] = gts_edge_new ( s->edge_class, v2, v1 );
00085                            if (LocalHead) fprintf(SUMA_STDERR,"Added edge of vertices: %d %d \n", SO->EL->EL[i][1], SO->EL->EL[i][0]);
00086          } else {
00087             edges[n++] = gts_edge_new ( s->edge_class, v1, v2 );
00088                                 if (LocalHead) fprintf(SUMA_STDERR,"Added edge of nodes: %d %d \n", SO->EL->EL[i][0], SO->EL->EL[i][1]);
00089          }
00090                 }
00091       MapToGts[i] = n-1; /* n-1 bc n was just incremented */
00092                 if (LocalHead) fprintf(SUMA_STDERR,"SUMA edge %d is mapped to GTS edge %d\n", i, n-1);
00093         }
00094 
00095         if (n != SO->EL->N_Distinct_Edges)
00096         {
00097                 fprintf(SUMA_STDERR, "distinct edges didn't equal N_Distinct_Edges");
00098                 exit(1);
00099         }
00100         /* this loop isn't working, stupid tri_limb
00101         for ( i=0; i< SO->N_FaceSet; i++)
00102         {
00103                 gts_surface_add_face(s,
00104                    gts_face_new ( s->face_class,
00105                         edges[MapToGts[SO->EL->Tri_limb[i][0]]],
00106                         edges[MapToGts[SO->EL->Tri_limb[i][1]]],
00107                         edges[MapToGts[SO->EL->Tri_limb[i][2]]]));
00108         }
00109         */
00110         for ( i=0; i< SO->N_FaceSet*3; i+=3)
00111         {
00112                 int n1 = SO->FaceSetList[i];
00113                 int n2 = SO->FaceSetList[i+1];
00114                 int n3 = SO->FaceSetList[i+2];
00115                 GtsEdge* e1 = edges[MapToGts[SUMA_FindEdge( SO->EL, n1, n2)]];
00116                 GtsEdge* e2 = edges[MapToGts[SUMA_FindEdge( SO->EL, n2, n3)]];
00117                 GtsEdge* e3 = edges[MapToGts[SUMA_FindEdge( SO->EL, n3, n1)]];
00118                 gts_surface_add_face(s,
00119                    gts_face_new ( s->face_class, e1, e2, e3));
00120       if (LocalHead) fprintf(SUMA_STDERR, "Added face of SUMA edges: %d %d %d\n"
00121                                           "              MapToGTS  : %d %d %d\n", 
00122                                           SUMA_FindEdge( SO->EL, n1, n2), SUMA_FindEdge( SO->EL, n2, n3), SUMA_FindEdge( SO->EL, n3, n1),
00123                                           MapToGts[SUMA_FindEdge( SO->EL, n1, n2)], MapToGts[SUMA_FindEdge( SO->EL, n2, n3)], MapToGts[SUMA_FindEdge( SO->EL, n3, n1)]);
00124         }
00125    
00126         g_free (vertices);
00127         g_free (edges);
00128         g_free (MapToGts);
00129         SUMA_RETURN( s );
00130 }
00131 
00132 
00133 #if 0
00134 /* These functions fail on OSX, see function gts_surface_suma and comments preceding it */ 
00135 SUMA_SurfaceObject* GtsToSuma( GtsSurface *gs)
00136 {
00137         static char FuncName[]={"GtsToSuma"};
00138    GHashTable *hash = g_hash_table_new(NULL,NULL);
00139         SUMA_SurfaceObject* SO = NULL; 
00140         int i = 0;
00141         gpointer data1[3];
00142    
00143    SUMA_ENTRY;
00144    
00145    GTS_OUT("shit.surf", gs);
00146    
00147    SO = (SUMA_SurfaceObject *)SUMA_malloc(sizeof(SUMA_SurfaceObject));
00148         SO->N_Node = gts_surface_vertex_number(gs);
00149         SO->NodeList = (float*) SUMA_calloc(SO->N_Node * 3, sizeof(float));
00150         SO->N_FaceSet = gts_surface_face_number(gs);
00151         SO->FaceSetList = (int*) SUMA_calloc(SO->N_FaceSet * 3, sizeof(int));
00152 
00153         /*
00154         foreach vertex
00155                 get xyz,
00156                 make hash key=pvertex value=integer
00157                 make NodeList from integer and xyz
00158         foreach face
00159                 get 3 vertexes
00160                 find their integers from hash
00161                 make FaceSetList from integers
00162         */
00163         data1[0] = hash;
00164         data1[1] = SO;
00165         data1[2] = &i;
00166         i = 0;
00167         gts_surface_foreach_vertex( gs, (GtsFunc) MakeNodeList_foreach_vertex, data1);
00168         i = 0;
00169         gts_surface_foreach_face( gs, (GtsFunc) MakeFaceList_foreach_face, data1);
00170         g_hash_table_destroy(hash);
00171         SUMA_RETURN(SO);
00172 }
00173 
00174 
00175 void MakeNodeList_foreach_vertex ( GtsPoint *p, gpointer* data)
00176 {
00177         GHashTable *hash = data[0];
00178         SUMA_SurfaceObject* SO = data[1];
00179         int* i = data[2];
00180         fprintf(SUMA_STDERR,"this %f %f %f\n", p->x, p->y, p->z);
00181    SO->NodeList[(*i)*3] = p->x;
00182         SO->NodeList[(*i)*3+1] = p->y;
00183         SO->NodeList[(*i)*3+2] = p->z;
00184         if (*i >= SO->N_Node || *i < 0)
00185         {
00186                 fprintf(SUMA_STDERR, "node %i out of range", *i);
00187                 exit(1);
00188         }
00189         
00190         if (g_hash_table_lookup(hash, p) == NULL)
00191                 g_hash_table_insert(hash, p, GINT_TO_POINTER(*i));
00192         else
00193         {
00194                 fprintf(SUMA_STDERR, "something already in hash table??");
00195                 exit(1);
00196         }
00197         (*i)++;
00198 }
00199 
00200 void MakeFaceList_foreach_face ( GtsFace* f, gpointer* data)
00201 {
00202         GHashTable *hash = data[0];
00203         SUMA_SurfaceObject* SO = data[1];
00204         int* i = data[2];
00205         GtsVertex *v1,*v2,*v3;
00206         gpointer presult = NULL;
00207         int iresult = 0;
00208         gpointer temp = NULL;
00209         gts_triangle_vertices(&(f->triangle), &v1, &v2, &v3);
00210         if ((*i) >= SO->N_FaceSet)
00211         {
00212                 fprintf(SUMA_STDERR, "counter %i passed N_FaceSet %i",(*i),SO->N_FaceSet);
00213                 exit(1);
00214         }
00215 
00216         if (!g_hash_table_lookup_extended(hash, v1, temp, &presult))
00217         {
00218                 fprintf(SUMA_STDERR, "hash table lookup failed");
00219                 exit(1);
00220         }
00221         iresult = GPOINTER_TO_INT (presult);
00222         if (iresult >= SO->N_Node || iresult < 0)
00223         {
00224                 fprintf(SUMA_STDERR, "node %i out of range", iresult);
00225                 exit(1);
00226         }
00227         SO->FaceSetList[(*i)*3] = iresult;
00228         
00229         if (!g_hash_table_lookup_extended(hash, v2, temp, &presult))
00230         {
00231                 fprintf(SUMA_STDERR, "hash table lookup failed");
00232                 exit(1);
00233         }
00234         iresult = GPOINTER_TO_INT (presult);
00235         if (iresult >= SO->N_Node || iresult < 0)
00236         {
00237                 fprintf(SUMA_STDERR, "node %i out of range", iresult);
00238                 exit(1);
00239         }
00240         SO->FaceSetList[(*i)*3+1] = iresult;
00241         
00242         if (!g_hash_table_lookup_extended(hash, v3, temp, &presult))
00243         {
00244                 fprintf(SUMA_STDERR, "hash table lookup failed");
00245                 exit(1);
00246         }
00247         iresult = GPOINTER_TO_INT (presult);
00248         if (iresult >= SO->N_Node || iresult < 0)
00249         {
00250                 fprintf(SUMA_STDERR, "node %i out of range", iresult);
00251                 exit(1);
00252         }
00253         SO->FaceSetList[(*i)*3+2] = iresult;
00254         (*i)++;
00255 }
00256 #endif
00257 
00258 void coarsen( GtsSurface* s, int stop)
00259 {
00260    GtsVolumeOptimizedParams params = { 0.5, 0.5, 0. };
00261    gdouble fold = PI/180.;
00262    gts_surface_coarsen (s,
00263       (GtsKeyFunc) gts_volume_optimized_cost,
00264       &params,
00265       (GtsCoarsenFunc) gts_volume_optimized_vertex,
00266       &params,
00267       (GtsStopFunc) gts_coarsen_stop_number,
00268       &stop,
00269       fold);
00270 }
00271 
00272 gboolean stop_number (gdouble cost, guint number, guint * max)
00273 {
00274   if (number > *max)
00275     return TRUE;
00276   return FALSE;
00277 }
00278 
00279 void refine( GtsSurface* s, int stop)
00280 {
00281    gts_surface_refine(s, NULL, NULL, NULL, NULL,
00282       (GtsStopFunc)stop_number, &stop);
00283 }
00284 
00285 /*!
00286    \brief resample a mesh so that the resultant surface has edge_factor as many edges as the original one
00287 */
00288 SUMA_SurfaceObject *SUMA_Mesh_Resample (SUMA_SurfaceObject *SO, float edge_factor)
00289 {
00290    static char FuncName[]={"SUMA_Mesh_Resample"};
00291    SUMA_SurfaceObject *S2=NULL;
00292    GtsSurface* s  = NULL;
00293    SUMA_Boolean LocalHead = NOPE;
00294    
00295    SUMA_ENTRY;
00296 
00297    if (!SO) {
00298       SUMA_SL_Err("NULL SO");
00299       SUMA_RETURN(S2);
00300    }
00301    if (!SO->EL) {
00302       SUMA_S_Warn("NULL Edge List, computing it");
00303       if (!SUMA_SurfaceMetrics(SO, "EdgeList", NULL)) {
00304          SUMA_SL_Err("Failed to create EdgeList");
00305          SUMA_RETURN(S2);
00306       }
00307    }
00308    /* create a GTS surface */
00309    s  = SumaToGts(SO);
00310    if (!s) {
00311       SUMA_S_Err("Failed to change SO to GTS surface");
00312       SUMA_RETURN(S2);
00313    }
00314    if (LocalHead) { /* see if surface was created well in GTS format */
00315       SUMA_S_Note("Writing initial surface in GTS form");
00316       GTS_OUT("gtsout.surf",s);
00317       GTS_VTK_OUT("vtkout.surf",s);
00318    }
00319 
00320    /* resample */
00321    if (1) {
00322       SUMA_S_Note("Changing mesh density\n");
00323       if (edge_factor < 1)
00324          coarsen(s, SO->EL->N_Distinct_Edges * edge_factor);
00325       else
00326          refine(s, SO->EL->N_Distinct_Edges * edge_factor);
00327    } else {
00328       SUMA_S_Note("Leaving surface untouched\n");
00329    }
00330    if (!s) {
00331       SUMA_SL_Err("Failed to refine");
00332       SUMA_RETURN(S2);
00333    }
00334   
00335    /* change the surface back to SUMA */
00336    S2 = SUMA_Alloc_SurfObject_Struct(1);
00337    gts_surface_suma (s, 
00338                      &(S2->NodeList), &(S2->N_Node), &(S2->NodeDim), 
00339                            &(S2->FaceSetList), &(S2->N_FaceSet), &(S2->FaceSetDim));
00340    
00341    gts_object_destroy((GtsObject*)s); s = NULL;
00342 
00343    SUMA_RETURN(S2);
00344 }