1 /***************************************************************************
2 sumagts.c - description
3 -------------------
4 begin : Mon Jun 14 2004
5 copyright : (C) 2004 by Jakub Otwinowski
6 email : jakub@hurin
7 ***************************************************************************/
8
9 /***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17 #ifdef HAVE_CONFIG_H
18 #include "gts/config.h"
19 #endif
20
21 #include "SUMA_suma.h"
22 #include "SUMA_gts.h"
23
24 extern SUMA_CommonFields *SUMAg_CF;
25 extern SUMA_DO *SUMAg_DOv;
26 extern SUMA_SurfaceViewer *SUMAg_SVv;
27 extern int SUMAg_N_SVv;
28 extern int SUMAg_N_DOv;
29
30 #if 0
31 /* Normally, functions in SUMA_gts_insert.c would be here */
32 #endif
33
SumaToGts(SUMA_SurfaceObject * SO)34 GtsSurface* SumaToGts( SUMA_SurfaceObject *SO)
35 {
36 static char FuncName[]={"SumaToGts"};
37 GtsSurface* s = NULL;
38 GtsVertex ** vertices = NULL;
39 GtsEdge ** edges = NULL;
40 int i = 0; /*counters */
41 int n = 0;
42 int *MapToGts = NULL; /*used to map EL vector to edge vector*/
43 SUMA_Boolean LocalHead = NOPE;
44
45 SUMA_ENTRY;
46
47 if (!SO->EL) {
48 SUMA_SL_Err("Null Edgeist!");
49 SUMA_RETURN(s);
50 }
51
52 SUMA_LH("In with ");
53 if (LocalHead) SUMA_Print_Surface_Object(SO, stderr);
54
55 s = gts_surface_new( gts_surface_class (),
56 gts_face_class (),
57 gts_edge_class (),
58 gts_vertex_class ());
59 vertices = (GtsVertex **)g_malloc (SO->N_Node * sizeof (GtsVertex *));
60 n = 0;
61 for ( i=0; i< SO->N_Node*3; i+=3)
62 {
63 vertices[n] = gts_vertex_new( s->vertex_class,
64 (gdouble)SO->NodeList[i],
65 (gdouble)SO->NodeList[i+1],
66 (gdouble)SO->NodeList[i+2]);
67 if (LocalHead)
68 fprintf(SUMA_STDERR, "Added vertex (%d/%d)%f %f %f\n"
69 "Stored gts (fails for some reason, but surface is OK) %f %f %f\n"
70 , n, SO->N_Node-1, SO->NodeList[i], SO->NodeList[i+1], SO->NodeList[i+2]
71 , (float)(vertices[n]->p.x), (float)(vertices[n]->p.y), (float)(vertices[n]->p.z));
72 ++n;
73 }
74 edges = (GtsEdge**)g_malloc (SO->EL->N_Distinct_Edges * sizeof (GtsEdge*));
75 n = 0;
76 MapToGts = (int *)g_malloc ( SO->EL->N_EL * sizeof(int));
77 for ( i=0; i< SO->EL->N_EL; i++)
78 {
79 if (SO->EL->ELps[i][2] > 0)
80 { /* a unique edge*/
81 GtsVertex* v1 = vertices[SO->EL->EL[i][0]];
82 GtsVertex* v2 = vertices[SO->EL->EL[i][1]];
83 if (SO->EL->ELps[i][0] == 1) {
84 edges[n++] = gts_edge_new ( s->edge_class, v2, v1 );
85 if (LocalHead) fprintf(SUMA_STDERR,"Added edge of vertices: %d %d \n", SO->EL->EL[i][1], SO->EL->EL[i][0]);
86 } else {
87 edges[n++] = gts_edge_new ( s->edge_class, v1, v2 );
88 if (LocalHead) fprintf(SUMA_STDERR,"Added edge of nodes: %d %d \n", SO->EL->EL[i][0], SO->EL->EL[i][1]);
89 }
90 }
91 MapToGts[i] = n-1; /* n-1 bc n was just incremented */
92 if (LocalHead) fprintf(SUMA_STDERR,"SUMA edge %d is mapped to GTS edge %d\n", i, n-1);
93 }
94
95 if (n != SO->EL->N_Distinct_Edges)
96 {
97 fprintf(SUMA_STDERR, "distinct edges didn't equal N_Distinct_Edges");
98 exit(1);
99 }
100 /* this loop isn't working, stupid tri_limb
101 for ( i=0; i< SO->N_FaceSet; i++)
102 {
103 gts_surface_add_face(s,
104 gts_face_new ( s->face_class,
105 edges[MapToGts[SO->EL->Tri_limb[i][0]]],
106 edges[MapToGts[SO->EL->Tri_limb[i][1]]],
107 edges[MapToGts[SO->EL->Tri_limb[i][2]]]));
108 }
109 */
110 for ( i=0; i< SO->N_FaceSet*3; i+=3)
111 {
112 int n1 = SO->FaceSetList[i];
113 int n2 = SO->FaceSetList[i+1];
114 int n3 = SO->FaceSetList[i+2];
115 GtsEdge* e1 = edges[MapToGts[SUMA_FindEdge( SO->EL, n1, n2)]];
116 GtsEdge* e2 = edges[MapToGts[SUMA_FindEdge( SO->EL, n2, n3)]];
117 GtsEdge* e3 = edges[MapToGts[SUMA_FindEdge( SO->EL, n3, n1)]];
118 gts_surface_add_face(s,
119 gts_face_new ( s->face_class, e1, e2, e3));
120 if (LocalHead) fprintf(SUMA_STDERR, "Added face of SUMA edges: %d %d %d\n"
121 " MapToGTS : %d %d %d\n",
122 SUMA_FindEdge( SO->EL, n1, n2), SUMA_FindEdge( SO->EL, n2, n3), SUMA_FindEdge( SO->EL, n3, n1),
123 MapToGts[SUMA_FindEdge( SO->EL, n1, n2)], MapToGts[SUMA_FindEdge( SO->EL, n2, n3)], MapToGts[SUMA_FindEdge( SO->EL, n3, n1)]);
124 }
125
126 g_free (vertices);
127 g_free (edges);
128 g_free (MapToGts);
129 SUMA_RETURN( s );
130 }
131
132
133 #if 0
134 /* These functions fail on OSX, see function gts_surface_suma and comments preceding it */
135 SUMA_SurfaceObject* GtsToSuma( GtsSurface *gs)
136 {
137 static char FuncName[]={"GtsToSuma"};
138 GHashTable *hash = g_hash_table_new(NULL,NULL);
139 SUMA_SurfaceObject* SO = NULL;
140 int i = 0;
141 gpointer data1[3];
142
143 SUMA_ENTRY;
144
145 GTS_OUT("fails.surf", gs);
146
147 SO = (SUMA_SurfaceObject *)SUMA_malloc(sizeof(SUMA_SurfaceObject));
148 SO->N_Node = gts_surface_vertex_number(gs);
149 SO->NodeList = (float*) SUMA_calloc(SO->N_Node * 3, sizeof(float));
150 SO->N_FaceSet = gts_surface_face_number(gs);
151 SO->FaceSetList = (int*) SUMA_calloc(SO->N_FaceSet * 3, sizeof(int));
152
153 /*
154 foreach vertex
155 get xyz,
156 make hash key=pvertex value=integer
157 make NodeList from integer and xyz
158 foreach face
159 get 3 vertexes
160 find their integers from hash
161 make FaceSetList from integers
162 */
163 data1[0] = hash;
164 data1[1] = SO;
165 data1[2] = &i;
166 i = 0;
167 gts_surface_foreach_vertex( gs, (GtsFunc) MakeNodeList_foreach_vertex, data1);
168 i = 0;
169 gts_surface_foreach_face( gs, (GtsFunc) MakeFaceList_foreach_face, data1);
170 g_hash_table_destroy(hash);
171 SUMA_RETURN(SO);
172 }
173
174
175 void MakeNodeList_foreach_vertex ( GtsPoint *p, gpointer* data)
176 {
177 GHashTable *hash = data[0];
178 SUMA_SurfaceObject* SO = data[1];
179 int* i = data[2];
180 fprintf(SUMA_STDERR,"this %f %f %f\n", p->x, p->y, p->z);
181 SO->NodeList[(*i)*3] = p->x;
182 SO->NodeList[(*i)*3+1] = p->y;
183 SO->NodeList[(*i)*3+2] = p->z;
184 if (*i >= SO->N_Node || *i < 0)
185 {
186 fprintf(SUMA_STDERR, "node %i out of range", *i);
187 exit(1);
188 }
189
190 if (g_hash_table_lookup(hash, p) == NULL)
191 g_hash_table_insert(hash, p, GINT_TO_POINTER(*i));
192 else
193 {
194 fprintf(SUMA_STDERR, "something already in hash table??");
195 exit(1);
196 }
197 (*i)++;
198 }
199
200 void MakeFaceList_foreach_face ( GtsFace* f, gpointer* data)
201 {
202 GHashTable *hash = data[0];
203 SUMA_SurfaceObject* SO = data[1];
204 int* i = data[2];
205 GtsVertex *v1,*v2,*v3;
206 gpointer presult = NULL;
207 int iresult = 0;
208 gpointer temp = NULL;
209 gts_triangle_vertices(&(f->triangle), &v1, &v2, &v3);
210 if ((*i) >= SO->N_FaceSet)
211 {
212 fprintf(SUMA_STDERR, "counter %i passed N_FaceSet %i",(*i),SO->N_FaceSet);
213 exit(1);
214 }
215
216 if (!g_hash_table_lookup_extended(hash, v1, temp, &presult))
217 {
218 fprintf(SUMA_STDERR, "hash table lookup failed");
219 exit(1);
220 }
221 iresult = GPOINTER_TO_INT (presult);
222 if (iresult >= SO->N_Node || iresult < 0)
223 {
224 fprintf(SUMA_STDERR, "node %i out of range", iresult);
225 exit(1);
226 }
227 SO->FaceSetList[(*i)*3] = iresult;
228
229 if (!g_hash_table_lookup_extended(hash, v2, temp, &presult))
230 {
231 fprintf(SUMA_STDERR, "hash table lookup failed");
232 exit(1);
233 }
234 iresult = GPOINTER_TO_INT (presult);
235 if (iresult >= SO->N_Node || iresult < 0)
236 {
237 fprintf(SUMA_STDERR, "node %i out of range", iresult);
238 exit(1);
239 }
240 SO->FaceSetList[(*i)*3+1] = iresult;
241
242 if (!g_hash_table_lookup_extended(hash, v3, temp, &presult))
243 {
244 fprintf(SUMA_STDERR, "hash table lookup failed");
245 exit(1);
246 }
247 iresult = GPOINTER_TO_INT (presult);
248 if (iresult >= SO->N_Node || iresult < 0)
249 {
250 fprintf(SUMA_STDERR, "node %i out of range", iresult);
251 exit(1);
252 }
253 SO->FaceSetList[(*i)*3+2] = iresult;
254 (*i)++;
255 }
256 #endif
257
coarsen(GtsSurface * s,int stop)258 void coarsen( GtsSurface* s, int stop)
259 {
260 GtsVolumeOptimizedParams params = { 0.5, 0.5, 0. };
261 gdouble fold = PI/180.;
262 gts_surface_coarsen (s,
263 (GtsKeyFunc) gts_volume_optimized_cost,
264 ¶ms,
265 (GtsCoarsenFunc) gts_volume_optimized_vertex,
266 ¶ms,
267 (GtsStopFunc) gts_coarsen_stop_number,
268 &stop,
269 fold);
270 }
271
stop_number(gdouble cost,guint number,guint * max)272 gboolean stop_number (gdouble cost, guint number, guint * max)
273 {
274 if (number > *max)
275 return TRUE;
276 return FALSE;
277 }
278
refine(GtsSurface * s,int stop)279 void refine( GtsSurface* s, int stop)
280 {
281 gts_surface_refine(s, NULL, NULL, NULL, NULL,
282 (GtsStopFunc)stop_number, &stop);
283 }
284
285 /* A convenience function for SUMA_Mesh_Resample.
286 if new_N_Nodes is between 0 and 1 then it is taken
287 to be the fraction of nodes in the original surface.*/
SUMA_Mesh_Resample_nodes(SUMA_SurfaceObject * SO,float new_N_Nodes)288 SUMA_SurfaceObject *SUMA_Mesh_Resample_nodes(SUMA_SurfaceObject *SO,
289 float new_N_Nodes)
290 {
291 static char FuncName[]={"SUMA_Mesh_Resample_nodes"};
292 int new_N_Edges = 0;
293 float edge_factor = 0;
294
295 if (new_N_Nodes > 0.0 && new_N_Nodes < 1.0) {
296 new_N_Nodes = (int)(new_N_Nodes*SO->N_Node);
297 } else new_N_Nodes = (int)new_N_Nodes;
298
299 /* approximate number of edges */
300 new_N_Edges = 3*new_N_Nodes-6;
301 edge_factor = (float)new_N_Edges/(float)SO->EL->N_Distinct_Edges;
302
303 return(SUMA_Mesh_Resample (SO, edge_factor));
304 }
305
306 /*!
307 \brief resample a mesh so that the resultant surface has edge_factor as many edges as the original one
308 */
SUMA_Mesh_Resample(SUMA_SurfaceObject * SO,float edge_factor)309 SUMA_SurfaceObject *SUMA_Mesh_Resample (SUMA_SurfaceObject *SO,
310 float edge_factor)
311 {
312 static char FuncName[]={"SUMA_Mesh_Resample"};
313 SUMA_SurfaceObject *S2=NULL;
314 GtsSurface* s = NULL;
315 SUMA_Boolean LocalHead = NOPE;
316
317 SUMA_ENTRY;
318
319 if (!SO) {
320 SUMA_SL_Err("NULL SO");
321 SUMA_RETURN(S2);
322 }
323 if (!SO->EL) {
324 SUMA_S_Warn("NULL Edge List, computing it");
325 if (!SUMA_SurfaceMetrics(SO, "EdgeList", NULL)) {
326 SUMA_SL_Err("Failed to create EdgeList");
327 SUMA_RETURN(S2);
328 }
329 }
330 /* create a GTS surface */
331 s = SumaToGts(SO);
332 if (!s) {
333 SUMA_S_Err("Failed to change SO to GTS surface");
334 SUMA_RETURN(S2);
335 }
336 if (LocalHead) { /* see if surface was created well in GTS format */
337 SUMA_S_Note("Writing initial surface in GTS form");
338 GTS_OUT("gtsout.surf",s);
339 GTS_VTK_OUT("vtkout.surf",s);
340 }
341
342 /* resample */
343 if (1) {
344 SUMA_LH("Changing mesh density\n");
345 if (edge_factor < 1)
346 coarsen(s, SO->EL->N_Distinct_Edges * edge_factor);
347 else
348 refine(s, SO->EL->N_Distinct_Edges * edge_factor);
349 } else {
350 SUMA_LH("Leaving surface untouched\n");
351 }
352 if (!s) {
353 SUMA_SL_Err("Failed to refine");
354 SUMA_RETURN(S2);
355 }
356
357 /* change the surface back to SUMA */
358 S2 = SUMA_Alloc_SurfObject_Struct(1);
359 gts_surface_suma (s,
360 &(S2->NodeList), &(S2->N_Node), &(S2->NodeDim),
361 &(S2->FaceSetList), &(S2->N_FaceSet), &(S2->FaceSetDim));
362
363 gts_object_destroy((GtsObject*)s); s = NULL;
364
365 SUMA_RETURN(S2);
366 }
367