00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
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
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;
00041 int n = 0;
00042 int *MapToGts = NULL;
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 {
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;
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
00101
00102
00103
00104
00105
00106
00107
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
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
00155
00156
00157
00158
00159
00160
00161
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 ¶ms,
00265 (GtsCoarsenFunc) gts_volume_optimized_vertex,
00266 ¶ms,
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
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
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) {
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
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
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 }