00001 #include "SUMA_suma.h"
00002 #include "SUMA_Macros.h"
00003 #if 0
00004
00005 #include "malloc.h"
00006 #endif
00007
00008 #undef STAND_ALONE
00009
00010 #if defined SUMA_CreateIcosahedron_STAND_ALONE
00011 #define STAND_ALONE
00012 #elif defined SUMA_MapIcosahedron_STAND_ALONE
00013 #define STAND_ALONE
00014 #elif defined SUMA_Map_SurfacetoSurface_STAND_ALONE
00015 #define STAND_ALONE
00016 #elif defined SUMA_AverageMaps_STAND_ALONE
00017 #define STAND_ALONE
00018 #elif defined SUMA_GiveColor_STAND_ALONE
00019 #define STAND_ALONE
00020 #elif defined SUMA_Test_STAND_ALONE
00021 #define STAND_ALONE
00022 #endif
00023
00024
00025 #ifdef STAND_ALONE
00026
00027 SUMA_SurfaceViewer *SUMAg_cSV;
00028 SUMA_SurfaceViewer *SUMAg_SVv;
00029 int SUMAg_N_SVv = 0;
00030 SUMA_DO *SUMAg_DOv;
00031 int SUMAg_N_DOv = 0;
00032 SUMA_CommonFields *SUMAg_CF;
00033 #else
00034 extern SUMA_CommonFields *SUMAg_CF;
00035 extern int SUMAg_N_DOv;
00036 extern SUMA_DO *SUMAg_DOv;
00037 #endif
00038
00039 float ep = 1e-4;
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066 SUMA_Boolean SUMA_SphereQuality(SUMA_SurfaceObject *SO, char *Froot, char *shist)
00067 {
00068 static char FuncName[]={"SUMA_SphereQuality"};
00069 float *dist = NULL, mdist, *dot=NULL, nr, r[3], *bad_dot = NULL;
00070 float dmin, dmax, dminloc, dmaxloc;
00071 int i, i3, *isortdist = NULL;
00072 int *bad_ind = NULL, ibad =-1;
00073 FILE *fid;
00074 char *fname;
00075 SUMA_COLOR_MAP *CM;
00076 SUMA_SCALE_TO_MAP_OPT * OptScl;
00077 SUMA_COLOR_SCALED_VECT * SV;
00078 SUMA_Boolean LocalHead = NOPE;
00079
00080 SUMA_ENTRY;
00081
00082 if (!SO) {
00083 SUMA_SL_Err("NULL SO");
00084 SUMA_RETURN(NOPE);
00085 }
00086
00087
00088 OptScl = SUMA_ScaleToMapOptInit();
00089 if (!OptScl) {
00090 fprintf (SUMA_STDERR,"Error %s: Could not get scaling option structure.\n", FuncName);
00091 exit (1);
00092 }
00093
00094
00095 CM = SUMA_GetStandardMap (SUMA_CMAP_MATLAB_DEF_BYR64);
00096 if (CM == NULL) {
00097 fprintf (SUMA_STDERR,"Error %s: Could not get standard colormap.\n", FuncName);
00098 if (OptScl) SUMA_free(OptScl);
00099 exit (1);
00100 }
00101
00102
00103 dist = (float *)SUMA_calloc(SO->N_Node, sizeof(float));
00104 mdist = 0.0;
00105 for (i=0; i<SO->N_Node; ++i) {
00106 i3 = 3*i;
00107 dist[i] = sqrt ( pow((double)(SO->NodeList[i3] - SO->Center[0]), 2.0) +
00108 pow((double)(SO->NodeList[i3+1] - SO->Center[1]), 2.0) +
00109 pow((double)(SO->NodeList[i3+2] - SO->Center[2]), 2.0) );
00110 mdist += dist[i];
00111 }
00112 mdist /= (float)SO->N_Node;
00113
00114
00115 for (i=0; i<SO->N_Node; ++i) dist[i] = fabs(dist[i] - mdist);
00116
00117
00118
00119 SV = SUMA_Create_ColorScaledVect(SO->N_Node);
00120 if (!SV) {
00121 fprintf (SUMA_STDERR,"Error %s: Could not allocate for SV.\n", FuncName);
00122 if (dist) SUMA_free(dist);
00123 if (CM) SUMA_Free_ColorMap (CM);
00124 if (OptScl) SUMA_free(OptScl);
00125 exit(1);
00126 }
00127 SUMA_MIN_MAX_VEC(dist, SO->N_Node, dmin, dmax, dminloc, dmaxloc);
00128 if (!SUMA_ScaleToMap (dist, SO->N_Node, dmin, dmax, CM, OptScl, SV)) {
00129 fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_ScaleToMap.\n", FuncName);
00130 if (dist) SUMA_free(dist);
00131 if (CM) SUMA_Free_ColorMap (CM);
00132 if (OptScl) SUMA_free(OptScl);
00133 exit(1);
00134 }
00135
00136
00137 fname = SUMA_append_string(Froot, "_Dist.1D.dset");
00138 if (LocalHead) fprintf (SUMA_STDERR,"%s:\nWriting %s...\n", FuncName, fname);
00139 fid = fopen(fname, "w");
00140 fprintf(fid,"#Node distance from center of mass.\n"
00141 "#col 0: Node Index\n"
00142 "#col 1: distance\n");
00143 if (shist) fprintf(fid,"#History:%s\n", shist);
00144 for (i=0; i<SO->N_Node; ++i) fprintf(fid,"%d\t%f\n", i, dist[i]);
00145 fclose(fid);
00146 SUMA_free(fname); fname = NULL;
00147
00148
00149 fname = SUMA_append_string(Froot, "_Dist.1D.col");
00150 if (LocalHead) fprintf (SUMA_STDERR,"%s:\nWriting %s...\n", FuncName, fname);
00151 fid = fopen(fname, "w");
00152 fprintf(fid,"#Color file of node distance from center of mass.\n"
00153 "#col 0: Node Index\n"
00154 "#col 1: R\n"
00155 "#col 2: G\n"
00156 "#col 3: B\n");
00157 if (shist) fprintf(fid,"#History:%s\n", shist);
00158 for (i=0; i<SO->N_Node; ++i) fprintf(fid,"%d\t%f\t%f\t%f\n", i, SV->cV[3*i ], SV->cV[3*i+1], SV->cV[3*i+2]);
00159 fclose(fid);
00160 SUMA_free(fname); fname = NULL;
00161 if (SV) SUMA_Free_ColorScaledVect (SV);
00162
00163
00164 isortdist = SUMA_z_qsort ( dist , SO->N_Node );
00165
00166
00167 fprintf (SUMA_STDERR,"\n");
00168 fprintf (SUMA_STDERR,"%s: \n"
00169 "Reporting on Spheriosity of %s\n", FuncName, SO->Label);
00170 fprintf (SUMA_STDERR," Mean distance from center (estimated radius): %f\n", mdist);
00171 fprintf (SUMA_STDERR," Largest 10 absolute departures from estimated radius:\n"
00172 " See output files for more detail.\n");
00173 for (i=SO->N_Node-1; i > SO->N_Node - 10; --i) {
00174 fprintf (SUMA_STDERR,"dist @ %d: %f\n", isortdist[i], dist[i]);
00175 }
00176
00177
00178
00179
00180
00181
00182
00183 dot = (float *)SUMA_calloc(SO->N_Node, sizeof(float));
00184 bad_ind = (int *) SUMA_calloc(SO->N_Node, sizeof(int) );
00185 bad_dot = (float *)SUMA_calloc(SO->N_Node, sizeof(float));
00186 ibad = 0;
00187 for (i=0; i<SO->N_Node; ++i) {
00188 i3 = 3*i;
00189 r[0] = SO->NodeList[i3] - SO->Center[0];
00190 r[1] = SO->NodeList[i3+1] - SO->Center[1];
00191 r[2] = SO->NodeList[i3+2] - SO->Center[2];
00192 nr = sqrt ( r[0] * r[0] + r[1] * r[1] + r[2] * r[2] );
00193 r[0] /= nr; r[1] /= nr; r[2] /= nr;
00194
00195 dot[i] = r[0]*SO->NodeNormList[i3] +
00196 r[1]*SO->NodeNormList[i3+1] +
00197 r[2]*SO->NodeNormList[i3+2] ;
00198
00199 if (fabs(dot[i]) < 0.9) {
00200 bad_ind[ibad] = i;
00201 bad_dot[ibad] = dot[i];
00202 ++ibad;
00203 }
00204 }
00205
00206 bad_ind = (int *) SUMA_realloc(bad_ind, ibad * sizeof(int));
00207 bad_dot = (float *)SUMA_realloc(bad_dot, ibad * sizeof(float));
00208
00209 fname = SUMA_append_string(Froot, "_dotprod.1D.dset");
00210 if (LocalHead) fprintf (SUMA_STDERR,"%s:\nWriting %s...\n", FuncName, fname);
00211 fid = fopen(fname, "w");
00212 fprintf(fid,"#Cosine of node normal angles with radial direction\n"
00213 "#col 0: Node Index\n"
00214 "#col 1: cos(angle)\n"
00215 );
00216 if (shist) fprintf(fid,"#History:%s\n", shist);
00217 for (i=0; i<SO->N_Node; ++i) fprintf(fid,"%d\t%f\n", i, dot[i]);
00218 fclose(fid);
00219 SUMA_free(fname); fname = NULL;
00220
00221
00222 SV = SUMA_Create_ColorScaledVect(SO->N_Node);
00223 if (!SV) {
00224 fprintf (SUMA_STDERR,"Error %s: Could not allocate for SV.\n", FuncName);
00225 if (dot) SUMA_free(dot);
00226 if (bad_dot) SUMA_free(bad_dot);
00227 if (bad_ind) SUMA_free(bad_ind);
00228 if (isortdist) SUMA_free(isortdist);
00229 if (dist) SUMA_free(dist);
00230 if (CM) SUMA_Free_ColorMap (CM);
00231 if (OptScl) SUMA_free(OptScl);
00232 exit(1);
00233 }
00234
00235 if (!SUMA_ScaleToMap (dot, SO->N_Node, -1.0, 1.0, CM, OptScl, SV)) {
00236 fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_ScaleToMap.\n", FuncName);
00237 exit(1);
00238 }
00239 fname = SUMA_append_string(Froot, "_dotprod.1D.col");
00240 if (LocalHead) fprintf (SUMA_STDERR,"%s:\nWriting %s...\n", FuncName, fname);
00241 fid = fopen(fname, "w");
00242 fprintf(fid,"#Color file of cosine of node normal angles with radial direction\n"
00243 "#col 0: Node Index\n"
00244 "#col 1: R\n"
00245 "#col 2: G\n"
00246 "#col 3: B\n"
00247 );
00248 if (shist) fprintf(fid,"#History:%s\n", shist);
00249 for (i=0; i<SO->N_Node; ++i) fprintf(fid,"%d\t%f\t%f\t%f\n", i, SV->cV[3*i ], SV->cV[3*i+1], SV->cV[3*i+2]);
00250 fclose(fid);
00251 SUMA_free(fname); fname = NULL;
00252 if (SV) SUMA_Free_ColorScaledVect (SV);
00253
00254 fname = SUMA_append_string(Froot, "_BadNodes.1D.dset");
00255 if (LocalHead) fprintf (SUMA_STDERR,"%s:\nWriting %s...\n", FuncName, fname);
00256 fid = fopen(fname, "w");
00257 fprintf(fid,"#Nodes with normals at angle with radial direction: abs(dot product < 0.9)\n"
00258 "#col 0: Node Index\n"
00259 "#col 1: cos(angle)\n"
00260 );
00261 if (shist) fprintf(fid,"#History:%s\n", shist);
00262 for (i=0; i<ibad; ++i) fprintf(fid,"%d\t%f\n", bad_ind[i], bad_dot[i]);
00263 fclose(fid);
00264 SUMA_free(fname); fname = NULL;
00265
00266
00267 SV = SUMA_Create_ColorScaledVect(ibad);
00268 if (!SV) {
00269 fprintf (SUMA_STDERR,"Error %s: Could not allocate for SV.\n", FuncName);
00270 if (dot) SUMA_free(dot);
00271 if (bad_dot) SUMA_free(bad_dot);
00272 if (bad_ind) SUMA_free(bad_ind);
00273 if (isortdist) SUMA_free(isortdist);
00274 if (dist) SUMA_free(dist);
00275 if (CM) SUMA_Free_ColorMap (CM);
00276 if (OptScl) SUMA_free(OptScl);
00277 exit(1);
00278 }
00279
00280 if (!SUMA_ScaleToMap (bad_dot, ibad, -1.0, 1.0, CM, OptScl, SV)) {
00281 fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_ScaleToMap.\n", FuncName);
00282 if (dot) SUMA_free(dot);
00283 if (bad_dot) SUMA_free(bad_dot);
00284 if (bad_ind) SUMA_free(bad_ind);
00285 if (isortdist) SUMA_free(isortdist);
00286 if (dist) SUMA_free(dist);
00287 if (CM) SUMA_Free_ColorMap (CM);
00288 if (OptScl) SUMA_free(OptScl);
00289 exit(1);
00290 }
00291 fname = SUMA_append_string(Froot, "_BadNodes.1D.col");
00292 if (LocalHead) fprintf (SUMA_STDERR,"%s:\nWriting %s...\n", FuncName, fname);
00293 fid = fopen(fname, "w");
00294 fprintf(fid,"#Color file of nodes with normals at angle with radial direction: abs(dot product < 0.9)\n"
00295 "#col 0: Node Index\n"
00296 "#col 1: R\n"
00297 "#col 2: G\n"
00298 "#col 3: B\n" );
00299 if (shist) fprintf(fid,"#History:%s\n", shist);
00300 for (i=0; i<ibad; ++i) fprintf(fid,"%d\t%f\t%f\t%f\n", bad_ind[i], SV->cV[3*i ], SV->cV[3*i+1], SV->cV[3*i+2]);
00301 fclose(fid);
00302 SUMA_free(fname); fname = NULL;
00303 if (SV) SUMA_Free_ColorScaledVect (SV);
00304
00305
00306
00307 if (ibad > 10) ibad = 10;
00308 fprintf (SUMA_STDERR,"%d of the nodes with normals at angle with radial direction\n"
00309 " i.e. abs(dot product < 0.9)\n"
00310 " See output files for full list\n", ibad);
00311 for (i=0; i < ibad; ++i) {
00312 fprintf (SUMA_STDERR,"cos(ang) @ node %d: %f\n", bad_ind[i], bad_dot[i]);
00313 }
00314
00315 if (dot) SUMA_free(dot);
00316 if (bad_dot) SUMA_free(bad_dot);
00317 if (bad_ind) SUMA_free(bad_ind);
00318 if (isortdist) SUMA_free(isortdist);
00319 if (dist) SUMA_free(dist);
00320 if (CM) SUMA_Free_ColorMap (CM);
00321 if (OptScl) SUMA_free(OptScl);
00322
00323
00324 SUMA_RETURN(YUP);
00325 }
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344 void SUMA_binTesselate(float *nodeList, int *triList, int *nCtr, int *tCtr, int recDepth, int depth, int n1, int n2, int n3)
00345 {
00346 double x1=0,y1=0,z1=0, x2=0,y2=0,z2=0, x3=0,y3=0,z3=0;
00347 double x12=0, y12=0, z12=0;
00348 double x23=0, y23=0, z23=0;
00349 double x31=0, y31=0, z31=0;
00350 int currIndex, index1, index2, index3;
00351 int i=0, j=0, m=0, k=0;
00352 static char FuncName[]={"SUMA_binTesselate"};
00353
00354 SUMA_ENTRY;
00355
00356 currIndex = (nCtr[0]-2)/3;
00357
00358 x1=(double)nodeList[3*n1]; y1=(double)nodeList[3*n1+1]; z1=(double)nodeList[3*n1+2];
00359 x2=(double)nodeList[3*n2]; y2=(double)nodeList[3*n2+1]; z2=(double)nodeList[3*n2+2];
00360 x3=(double)nodeList[3*n3]; y3=(double)nodeList[3*n3+1]; z3=(double)nodeList[3*n3+2];
00361
00362 x12=(x1+x2)/2.0; y12=(y1+y2)/2.0; z12=(z1+z2)/2.0;
00363 x23=(x2+x3)/2.0; y23=(y2+y3)/2.0; z23=(z2+z3)/2.0;
00364 x31=(x3+x1)/2.0; y31=(y3+y1)/2.0; z31=(z3+z1)/2.0;
00365
00366
00367 index1 = -1; index2 = -1; index3 = -1;
00368 i=0; j=0;
00369 for (i=0; i<=currIndex; ++i) {
00370 j = 3*i;
00371 if ( fabs(nodeList[j]-x12)<ep && fabs(nodeList[j+1]-y12)<ep && fabs(nodeList[j+2]-z12)<ep ) {
00372 index1 = i;
00373 }
00374 if ( fabs(nodeList[j]-x23)<ep && fabs(nodeList[j+1]-y23)<ep && fabs(nodeList[j+2]-z23)<ep ) {
00375 index2 = i;
00376 }
00377 if ( fabs(nodeList[j]-x31)<ep && fabs(nodeList[j+1]-y31)<ep && fabs(nodeList[j+2]-z31)<ep ) {
00378 index3 = i;
00379 }
00380 }
00381
00382 if (index1==-1) {
00383 ++currIndex;
00384 index1 = currIndex;
00385 SUMA_addNode( nodeList, nCtr, (float)x12, (float)y12, (float)z12);
00386 }
00387 if (index2==-1) {
00388 ++currIndex;
00389 index2 = currIndex;
00390 SUMA_addNode( nodeList, nCtr, (float)x23, (float)y23, (float)z23);
00391 }
00392 if (index3==-1) {
00393 ++currIndex;
00394 index3 = currIndex;
00395 SUMA_addNode( nodeList, nCtr, (float)x31, (float)y31, (float)z31);
00396 }
00397
00398
00399 if (depth>=recDepth) {
00400 SUMA_addTri( triList, tCtr, n1, index1, index3);
00401 SUMA_addTri( triList, tCtr, index1, n2, index2);
00402 SUMA_addTri( triList, tCtr, index3, index2, n3);
00403 SUMA_addTri( triList, tCtr, index3, index2, index1);
00404 }
00405
00406
00407 else {
00408 ++depth;
00409 SUMA_binTesselate( nodeList, triList, nCtr, tCtr, recDepth, depth, n1, index1, index3 );
00410 SUMA_binTesselate( nodeList, triList, nCtr, tCtr, recDepth, depth, index1, n2, index2 );
00411 SUMA_binTesselate( nodeList, triList, nCtr, tCtr, recDepth, depth, index3, index2, n3 );
00412 SUMA_binTesselate( nodeList, triList, nCtr, tCtr, recDepth, depth, index3, index2, index1 );
00413 }
00414
00415 SUMA_RETURNe;
00416 }
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433 void SUMA_tesselate( float *nodeList, int *triList, int *nCtr, int *tCtr, int N_Div, int n0, int n1, int n2) {
00434
00435 int i=0, j=0;
00436 int *edge01=NULL, *edge12=NULL, *edge20=NULL, *currFloor=NULL;
00437 static char FuncName[]={"SUMA_tesselate"};
00438
00439 SUMA_ENTRY;
00440
00441 edge01 = SUMA_divEdge( nodeList, nCtr, n0, n1, N_Div);
00442 edge12 = SUMA_divEdge( nodeList, nCtr, n2, n1, N_Div);
00443 edge20 = SUMA_divEdge( nodeList, nCtr, n0, n2, N_Div);
00444 if (!edge01 || !edge12 || !edge20) {
00445 fprintf (SUMA_STDERR, "Error %s: Failed in SUMA_divEdge.\n", FuncName);
00446 SUMA_RETURNe;
00447 }
00448
00449 currFloor = edge20;
00450
00451 for (i=1; i<N_Div; ++i) {
00452 SUMA_triangulateRow( nodeList, triList, nCtr, tCtr, N_Div-i, currFloor, edge01[i], edge12[i]);
00453 }
00454
00455 SUMA_addTri( triList, tCtr, currFloor[1], n1, currFloor[0]);
00456
00457 if (edge01) SUMA_free(edge01);
00458 if (edge12) SUMA_free(edge12);
00459 if (edge20) SUMA_free(edge20);
00460
00461 SUMA_RETURNe;
00462 }
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476 int * SUMA_divEdge( float *nodeList, int *nCtr, int node1, int node2, int N_Div) {
00477
00478 float *newNodes = NULL;
00479 float n1[3], n2[3];
00480 int *edge = NULL;
00481 int i=0, j=0, k=0, m=0;
00482 int currIndex = (nCtr[0]-2)/3;
00483 static char FuncName[]={"SUMA_divEdge"};
00484
00485 SUMA_ENTRY;
00486
00487
00488 edge = (int *) SUMA_calloc(N_Div+1, sizeof(int));
00489 newNodes = (float *)SUMA_calloc (3*(N_Div-1), sizeof(float));
00490
00491 if (!edge || !newNodes) {
00492 fprintf (SUMA_STDERR, "Error %s: Failed to allocate.\n", FuncName);
00493 SUMA_RETURN (edge);
00494 }
00495
00496 for(i=0; i<N_Div+1; ++i) {
00497 edge[i] = -1;
00498 }
00499
00500 edge[0] = node1; edge[N_Div] = node2;
00501
00502 n1[0] = nodeList[3*node1]; n1[1] = nodeList[3*node1+1]; n1[2] = nodeList[3*node1+2];
00503 n2[0] = nodeList[3*node2]; n2[1] = nodeList[3*node2+1]; n2[2] = nodeList[3*node2+2];
00504
00505
00506 for(i=0; i<N_Div-1; ++i) {
00507 j = 3*i;
00508 newNodes[j] = ((i+1.0)/(float)N_Div)*(n2[0]-n1[0]) + n1[0];
00509 newNodes[j+1] = ((i+1.0)/(float)N_Div)*(n2[1]-n1[1]) + n1[1];
00510 newNodes[j+2] = ((i+1.0)/(float)N_Div)*(n2[2]-n1[2]) + n1[2];
00511 }
00512
00513
00514 for (i=0; i<=currIndex; ++i) {
00515 j = 3*i;
00516 for (m=0; m<N_Div-1; ++m) {
00517 k = 3*m;
00518 if ( fabs(nodeList[j]-newNodes[k])<ep && fabs(nodeList[j+1]-newNodes[k+1])<ep &&
00519 fabs(nodeList[j+2]-newNodes[k+2])<ep ) {
00520 edge[m+1] = i;
00521 }
00522 }
00523 }
00524
00525 for (i=1; i<N_Div; ++i) {
00526 if (edge[i]==-1) {
00527 SUMA_addNode( nodeList, nCtr, newNodes[3*(i-1)], newNodes[3*(i-1)+1], newNodes[3*(i-1)+2]);
00528 edge[i] = (nCtr[0]-2)/3;
00529 }
00530 }
00531
00532 if (newNodes) SUMA_free(newNodes);
00533
00534 SUMA_RETURN (edge);
00535 }
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553 void SUMA_triangulateRow( float *nodeList, int *triList, int *nCtr, int *tCtr, int N_Div, int *currFloor, int node1, int node2) {
00554
00555 int i=0, j=0;
00556 float n1[3], n2[3], newNode[3];
00557 int *newArray = NULL;
00558 static char FuncName[]={"SUMA_triangulateRow"};
00559
00560 SUMA_ENTRY;
00561
00562 newArray = (int *)SUMA_calloc(N_Div+1, sizeof(int));
00563 if (!newArray) {
00564 fprintf (SUMA_STDERR, "Error %s: Failed to allocate.\n", FuncName);
00565 SUMA_RETURNe;
00566 }
00567
00568 n1[0] = nodeList[3*node1]; n1[1] = nodeList[3*node1+1]; n1[2] = nodeList[3*node1+2];
00569 n2[0] = nodeList[3*node2]; n2[1] = nodeList[3*node2+1]; n2[2] = nodeList[3*node2+2];
00570 newArray[0] = node1; newArray[N_Div] = node2;
00571
00572 SUMA_addTri( triList, tCtr, currFloor[1], currFloor[0], newArray[0]);
00573
00574 for (i=1; i<N_Div; ++i) {
00575 newNode[0] = ((float)i/(float)N_Div)*(n2[0]-n1[0]) + n1[0];
00576 newNode[1] = ((float)i/(float)N_Div)*(n2[1]-n1[1]) + n1[1];
00577 newNode[2] = ((float)i/(float)N_Div)*(n2[2]-n1[2]) + n1[2];
00578
00579 SUMA_addNode( nodeList, nCtr, newNode[0], newNode[1], newNode[2]);
00580 newArray[i] = (nCtr[0]-2)/3;
00581 SUMA_addTri( triList, tCtr, newArray[i-1], currFloor[i], newArray[i]);
00582 SUMA_addTri( triList, tCtr, currFloor[i+1], newArray[i], currFloor[i]);
00583 }
00584 SUMA_addTri( triList, tCtr, newArray[N_Div-1], currFloor[N_Div], newArray[N_Div]);
00585 SUMA_addTri( triList, tCtr, newArray[N_Div], currFloor[N_Div+1], currFloor[N_Div]);
00586
00587 for (i=0; i<N_Div+1; ++i) {
00588 currFloor[i] = newArray[i];
00589 }
00590
00591 if (newArray) SUMA_free(newArray);
00592
00593 SUMA_RETURNe;
00594 }
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606 void SUMA_addNode(float *nodeList, int *ctr, float x, float y, float z) {
00607
00608 static char FuncName[]={"SUMA_addNode"};
00609
00610 SUMA_ENTRY;
00611
00612 ++*ctr;
00613 nodeList[*ctr] = x;
00614 ++*ctr;
00615 nodeList[*ctr] = y;
00616 ++*ctr;
00617 nodeList[*ctr] = z;
00618
00619 SUMA_RETURNe;
00620 }
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630 void SUMA_addTri(int *triList, int *ctr, int n1, int n2, int n3) {
00631
00632 static char FuncName[]={"SUMA_addTri"};
00633
00634 SUMA_ENTRY;
00635
00636 ++*ctr;
00637 triList[*ctr] = n1;
00638 ++*ctr;
00639 triList[*ctr] = n2;
00640 ++*ctr;
00641 triList[*ctr] = n3;
00642
00643 SUMA_RETURNe;
00644 }
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661 SUMA_SurfaceObject * SUMA_CreateIcosahedron (float r, int depth, float ctr[3], char bin[], int ToSphere)
00662 {
00663 static char FuncName[]={"SUMA_CreateIcosahedron"};
00664 SUMA_SurfaceObject *SO = NULL;
00665 int i, numNodes=0, numTri=0, j, i3;
00666 float a,b, lgth;
00667 int nodePtCt, triPtCt, *icosaTri=NULL;
00668 float *icosaNode=NULL;
00669 SUMA_SURF_NORM SN;
00670 SUMA_NODE_FIRST_NEIGHB *firstNeighb=NULL;
00671 SUMA_Boolean DoWind = YUP;
00672 int n=0, m=0, in=0, trouble;
00673 SUMA_Boolean LocalHead = NOPE;
00674
00675 SUMA_ENTRY;
00676
00677 SO = SUMA_Alloc_SurfObject_Struct(1);
00678 if (SO == NULL) {
00679 fprintf (SUMA_STDERR,"Error %s: Failed to allocate for Surface Object.", FuncName);
00680 SUMA_RETURN (NULL);
00681 }
00682
00683
00684 if (strcmp(bin, "y") == 0) { numTri = 20*pow(2,2*depth); }
00685 else {
00686 if (depth !=0) { numTri = 20*pow(depth, 2); }
00687 else numTri = 20;
00688 }
00689 if (depth != 0) { numNodes = 3*numTri; }
00690 else numNodes = 12;
00691
00692
00693 if (LocalHead) fprintf(SUMA_STDERR,"%s: Allocated for %d Nodes, %d numTri\n", FuncName, numNodes, numTri);
00694
00695
00696 SUMA_ICOSAHEDRON_DIMENSIONS(r, a, b, lgth);
00697
00698 if (LocalHead) {
00699 fprintf(SUMA_STDERR,"%s: a = %f, b=%f, rad = %f, lgth = %f\n", FuncName, a, b, r, lgth);
00700 }
00701
00702
00703 if (strcmp(bin, "y") == 0) {
00704 ep = lgth / pow(2, depth+1);
00705 }
00706 else ep = lgth / (2*depth);
00707
00708
00709 nodePtCt = -1;
00710 icosaNode = (float *) SUMA_calloc(3*numNodes, sizeof(float));
00711 icosaTri = (int *) SUMA_calloc(3*numTri, sizeof(int));
00712
00713 if (!icosaNode || !icosaTri) {
00714 fprintf (SUMA_STDERR,"Error %s: Could not allocate for icosaNode and/or icosaTri.\n",FuncName);
00715 SUMA_Free_Surface_Object (SO);
00716 SUMA_RETURN (NULL);
00717 }
00718
00719 SUMA_addNode( icosaNode, &nodePtCt, 0+ctr[0], b+ctr[1], -a+ctr[2] );
00720 SUMA_addNode( icosaNode, &nodePtCt, 0+ctr[0], b+ctr[1], a+ctr[2] );
00721 SUMA_addNode( icosaNode, &nodePtCt, 0+ctr[0], -b+ctr[1], a+ctr[2] );
00722 SUMA_addNode( icosaNode, &nodePtCt, 0+ctr[0], -b+ctr[1], -a+ctr[2] );
00723 SUMA_addNode( icosaNode, &nodePtCt, -b+ctr[0], a+ctr[1], 0+ctr[2] );
00724 SUMA_addNode( icosaNode, &nodePtCt, -b+ctr[0], -a+ctr[1], 0+ctr[2] );
00725 SUMA_addNode( icosaNode, &nodePtCt, b+ctr[0], a+ctr[1], 0+ctr[2] );
00726 SUMA_addNode( icosaNode, &nodePtCt, b+ctr[0], -a+ctr[1], 0+ctr[2] );
00727 SUMA_addNode( icosaNode, &nodePtCt, a+ctr[0], 0+ctr[1], b+ctr[2] );
00728 SUMA_addNode( icosaNode, &nodePtCt, -a+ctr[0], 0+ctr[1], -b+ctr[2] );
00729 SUMA_addNode( icosaNode, &nodePtCt, -a+ctr[0], 0+ctr[1], b+ctr[2] );
00730 SUMA_addNode( icosaNode, &nodePtCt, a+ctr[0], 0+ctr[1], -b+ctr[2] );
00731
00732
00733
00734 triPtCt = -1;
00735
00736
00737 if (depth==0) {
00738
00739 SUMA_addTri( icosaTri, &triPtCt, 0, 4, 6 ); SUMA_addTri( icosaTri, &triPtCt, 1, 6, 4 );
00740 SUMA_addTri( icosaTri, &triPtCt, 0, 9, 4 ); SUMA_addTri( icosaTri, &triPtCt, 1, 8, 6 );
00741 SUMA_addTri( icosaTri, &triPtCt, 0, 3, 9 ); SUMA_addTri( icosaTri, &triPtCt, 1, 2, 8 );
00742 SUMA_addTri( icosaTri, &triPtCt, 0, 11, 3 ); SUMA_addTri( icosaTri, &triPtCt, 1, 10, 2 );
00743 SUMA_addTri( icosaTri, &triPtCt, 0, 6, 11 ); SUMA_addTri( icosaTri, &triPtCt, 1, 4, 10 );
00744 SUMA_addTri( icosaTri, &triPtCt, 2, 7, 8 ); SUMA_addTri( icosaTri, &triPtCt, 3, 11, 7 );
00745 SUMA_addTri( icosaTri, &triPtCt, 2, 5, 7 ); SUMA_addTri( icosaTri, &triPtCt, 3, 7, 5 );
00746 SUMA_addTri( icosaTri, &triPtCt, 2, 10, 5 ); SUMA_addTri( icosaTri, &triPtCt, 3, 5, 9 );
00747 SUMA_addTri( icosaTri, &triPtCt, 4, 9, 10 ); SUMA_addTri( icosaTri, &triPtCt, 6, 8, 11 );
00748 SUMA_addTri( icosaTri, &triPtCt, 5, 10, 9 ); SUMA_addTri( icosaTri, &triPtCt, 7, 11, 8 );
00749 }
00750
00751 else {
00752 if (strcmp(bin, "y") == 0) {
00753
00754 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 0, 4, 6);
00755 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 0, 9, 4 );
00756 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 0, 3, 9 );
00757 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 0, 11, 3 );
00758 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 0, 6, 11 );
00759
00760 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 1, 6, 4 );
00761 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 1, 8, 6 );
00762 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 1, 2, 8 );
00763 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 1, 10, 2 );
00764 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 1, 4, 10 );
00765
00766 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 2, 7, 8 );
00767 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 2, 5, 7 );
00768 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 2, 10, 5 );
00769 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 4, 9, 10 );
00770 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 5, 10, 9 );
00771
00772 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 3, 11, 7 );
00773 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 3, 7, 5 );
00774 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 3, 5, 9 );
00775 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 6, 8, 11 );
00776 SUMA_binTesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 7, 11, 8 );
00777 }
00778
00779 else {
00780
00781 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 0, 4, 6);
00782 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 0, 9, 4 );
00783 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 0, 3, 9 );
00784 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 0, 11, 3 );
00785 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 0, 6, 11 );
00786
00787 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 6, 4 );
00788 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 8, 6 );
00789 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 2, 8 );
00790 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 10, 2 );
00791 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 1, 4, 10 );
00792
00793 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 2, 7, 8 );
00794 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 2, 5, 7 );
00795 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 2, 10, 5 );
00796 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 4, 9, 10 );
00797 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 5, 10, 9 );
00798
00799 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 3, 11, 7 );
00800 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 3, 7, 5 );
00801 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 3, 5, 9 );
00802 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 6, 8, 11 );
00803 SUMA_tesselate(icosaNode, icosaTri, &nodePtCt, &triPtCt, depth, 7, 11, 8 );
00804 }
00805 }
00806
00807 numNodes = (nodePtCt+1)/3;
00808 numTri = (triPtCt+1)/3;
00809
00810 if (LocalHead) fprintf(SUMA_STDERR,"%s: There are %d nodes, %d triangles in the icosahedron.\n", FuncName, numNodes, numTri);
00811
00812
00813 SO->NodeList = icosaNode;
00814 SO->FaceSetList = icosaTri;
00815 SO->N_Node = numNodes;
00816 SO->N_FaceSet = numTri;
00817 SO->NodeDim = 3;
00818 SO->FaceSetDim = 3;
00819 SO->idcode_str = (char *)SUMA_calloc (SUMA_IDCODE_LENGTH, sizeof(char));
00820 UNIQ_idcode_fill (SO->idcode_str);
00821
00822
00823 if (DoWind) {
00824 if (LocalHead) fprintf(SUMA_STDOUT, "%s: Making Edge list ....\n", FuncName);
00825 SO->EL = SUMA_Make_Edge_List_eng (SO->FaceSetList, SO->N_FaceSet, SO->N_Node, SO->NodeList, 0, SO->idcode_str);
00826 if (SO->EL == NULL) {
00827 fprintf(SUMA_STDERR, "Error %s: Failed in SUMA_Make_Edge_List. Neighbor list will not be created\n", FuncName);
00828 SUMA_Free_Surface_Object (SO);
00829 SUMA_RETURN (NULL);
00830 } else {
00831 }
00832
00833 if (!SUMA_MakeConsistent (SO->FaceSetList, SO->N_FaceSet, SO->EL, 0, &trouble)) {
00834 fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_MakeConsistent.\n", FuncName);
00835 SUMA_Free_Surface_Object (SO);
00836 SUMA_RETURN (NULL);
00837 }
00838 else {
00839 if (LocalHead) fprintf(SUMA_STDERR,"%s: Eeeexcellent. All triangles consistent.\n", FuncName);
00840 }
00841
00842 if (LocalHead) fprintf(SUMA_STDOUT, "%s: Determining MemberFaceSets ...\n", FuncName);
00843 SO->MF = SUMA_MemberFaceSets(SO->N_Node, SO->FaceSetList, SO->N_FaceSet, SO->FaceSetDim, SO->idcode_str);
00844 if (SO->MF->NodeMemberOfFaceSet == NULL) {
00845 fprintf(SUMA_STDERR,"Error %s: Error in SUMA_MemberFaceSets\n", FuncName);
00846 SUMA_Free_Surface_Object (SO);
00847 SUMA_RETURN (NULL);
00848 }else {
00849 }
00850
00851
00852 }
00853
00854
00855 if (ToSphere) {
00856 float dv, uv[3], U[2][3], *p1;
00857 for (i=0; i<SO->N_Node; ++i) {
00858 i3 = 3*i;
00859 p1 = &(SO->NodeList[i3]);
00860
00861 uv[0] = p1[0] - ctr[0]; uv[1] = p1[1] - ctr[1]; uv[2] = p1[2] - ctr[2];
00862 SUMA_POINT_AT_DISTANCE(uv, ctr, r, U);
00863 SO->NodeList[i3 ] = U[0][0]; SO->NodeList[i3+1] = U[0][1]; SO->NodeList[i3+2] = U[0][2];
00864 }
00865 }
00866
00867
00868 SN = SUMA_SurfNorm( SO->NodeList, SO->N_Node, SO->FaceSetList, SO->N_FaceSet);
00869 SO->NodeNormList = SN.NodeNormList;
00870 SO->FaceNormList = SN.FaceNormList;
00871
00872
00873 if ( SO->EL && SO->N_Node )
00874 firstNeighb = SUMA_Build_FirstNeighb( SO->EL, SO->N_Node, SO->idcode_str);
00875
00876 if (!DoWind) SO->EL = SUMA_Make_Edge_List (SO->FaceSetList, SO->N_FaceSet, SO->N_Node, SO->NodeList, SO->idcode_str);
00877 SO->FN = SUMA_Build_FirstNeighb( SO->EL, SO->N_Node, SO->idcode_str);
00878 if(SO->FN==NULL) {
00879 fprintf(SUMA_STDERR, "Error %s: Failed in creating neighb list.\n", FuncName);
00880 } else {
00881 }
00882
00883
00884 SUMA_RETURN (SO);
00885 }
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904 SUMA_Boolean SUMA_inNodeNeighb( SUMA_SurfaceObject *surf, float *nodeList, int *node, float *P0, float *P1) {
00905
00906 int i=0, j=0, k=0, examinedNum=0;
00907 SUMA_Boolean found=NOPE;
00908 float hitOnSurf[3];
00909 int incidentTri[100], N_incident = 0, itry;
00910 int examinedTri[100], ifound, i_node0 = -1, i_node1 = -1, i_node2 = -1;
00911 SUMA_Boolean LocalHead = NOPE;
00912 static char FuncName[]={"SUMA_inNodeNeighb"};
00913
00914 SUMA_ENTRY;
00915
00916 if (nodeList==NULL) {
00917 fprintf (SUMA_STDERR, "Warning %s: Assigning surf->NodeList to nodeList.\n", FuncName);
00918 nodeList = surf->NodeList;
00919 }
00920
00921 if (LocalHead) fprintf(SUMA_STDERR, "%s: P0-P1 [%f, %f, %f] - [%f, %f, %f]\n",
00922 FuncName, P0[0], P0[1], P0[2], P1[0], P1[1], P1[2]);
00923
00924 found = NOPE;
00925 itry = 0;
00926 examinedNum = 0;
00927 while (itry < 3 && node[itry] >= 0 && !found) {
00928 if (LocalHead) fprintf(SUMA_STDERR, "%s: Trying neighbors of node %d.\n", FuncName, node[itry]);
00929 i = 0;
00930 while ((i < surf->FN->N_Neighb[node[itry]] ) && !found) {
00931
00932 if (!SUMA_Get_Incident( node[itry], surf->FN->FirstNeighb[node[itry]][i], surf->EL, incidentTri, &N_incident, 1)) {
00933 fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_Get_Incident.\n", FuncName);
00934 SUMA_RETURN (NOPE);
00935 }
00936
00937
00938 j = 0;
00939 while ((j < N_incident) && !found) {
00940
00941
00942 SUMA_IS_IN_VEC(examinedTri, examinedNum, incidentTri[j], ifound);
00943
00944
00945 if (ifound < 0) {
00946 examinedTri[examinedNum] = incidentTri[j];
00947 ++examinedNum;
00948
00949 i_node0 = surf->FaceSetList[ 3*incidentTri[j] ];
00950 i_node1 = surf->FaceSetList[ 3*incidentTri[j]+1 ];
00951 i_node2 = surf->FaceSetList[ 3*incidentTri[j]+2 ];
00952
00953 if (SUMA_MT_isIntersect_Triangle (P0, P1, &(nodeList[3*i_node0]), &(nodeList[3*i_node1]),
00954 &(nodeList[3*i_node2]), hitOnSurf, NULL, NULL)) {
00955 found = YUP;
00956 node[0] = i_node0;
00957 node[1] = i_node1;
00958 node[2] = i_node2;
00959 if (LocalHead) {
00960 fprintf(SUMA_STDERR, "%s: Triangle %d [%d, %d, %d] is intersected at (%f, %f, %f)\n",
00961 FuncName, incidentTri[j], node[0], node[1], node[2], hitOnSurf[0], hitOnSurf[1], hitOnSurf[2]);
00962 fprintf(SUMA_STDERR, "%s: Coordinates of nodes forming triangle are:\n", FuncName);
00963 fprintf(SUMA_STDERR, "%f, %f, %f\n", nodeList[3*i_node0], nodeList[3*i_node0+1], nodeList[3*i_node0+2]);
00964 fprintf(SUMA_STDERR, "%f, %f, %f\n", nodeList[3*i_node1], nodeList[3*i_node1+1], nodeList[3*i_node1+2]);
00965 fprintf(SUMA_STDERR, "%f, %f, %f\n", nodeList[3*i_node2], nodeList[3*i_node2+1], nodeList[3*i_node2+2]);
00966 }
00967 #if 0
00968 {
00969
00970 SUMA_MT_INTERSECT_TRIANGLE *MTI;
00971 MTI = SUMA_MT_intersect_triangle (P1, P0, nodeList, surf->N_Node, surf->FaceSetList, surf->N_FaceSet, NULL);
00972 if (MTI) {
00973 if (LocalHead)fprintf(SUMA_STDERR, "%s: Meth2-Triangle %d [%d, %d, %d] is intersected at (%f, %f, %f)\n",
00974 FuncName, MTI->ifacemin, surf->FaceSetList[3*MTI->ifacemin], surf->FaceSetList[3*MTI->ifacemin+1],
00975 surf->FaceSetList[3*MTI->ifacemin+2], MTI->P[0], MTI->P[1], MTI->P[2]);
00976
00977 if (MTI->N_hits) {
00978
00979 if (MTI->ifacemin != incidentTri[j]) {
00980 fprintf (SUMA_STDERR,"Error %s: Warning, mismatch in results of triangle intersection. This should not be\n", FuncName);
00981 exit(1);
00982 }
00983 }
00984
00985 MTI = SUMA_Free_MT_intersect_triangle(MTI);
00986 }
00987
00988 }
00989 #endif
00990
00991 P1[0] = hitOnSurf[0]; P1[1] = hitOnSurf[1]; P1[2] = hitOnSurf[2];
00992 }else {
00993 if (LocalHead)fprintf(SUMA_STDERR, "%s: Triangle %d [%d, %d, %d] is not intersected.\n",
00994 FuncName, incidentTri[j], i_node0, i_node1, i_node2);
00995 }
00996 }
00997 ++j;
00998 }
00999 ++i;
01000 }
01001 ++itry;
01002 }
01003
01004 SUMA_RETURN (found);
01005 }
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020 float * SUMA_detWeight (float node0[3], float node1[3], float node2[3], float ptHit[3]) {
01021
01022 int i=0;
01023 float triNode0[3], triNode1[3], triNode2[3];
01024 float p00[3], p01[3], p02[3];
01025 float p10[3], p11[3], p12[3];
01026 float p20[3], p21[3], p22[3];
01027 float tri0[3], tri1[3], tri2[3], triOrig[3];
01028 float s0=0, s1=0, s2=0, sOrig=0, A0=0, A1=0, A2=0, Aorig=0;
01029 float wsum=0, *weight=NULL;
01030 static char FuncName[]={"SUMA_detWeight"};
01031
01032 SUMA_ENTRY;
01033
01034
01035
01036
01037 p00[0] = node0[0]; p00[1] = node0[1]; p00[2] = node0[2];
01038 p11[0] = node1[0]; p11[1] = node1[1]; p11[2] = node1[2];
01039 p22[0] = node2[0]; p22[1] = node2[1]; p22[2] = node2[2];
01040
01041
01042
01043
01044
01045 for (i=0; i<3; ++i) {
01046
01047 if (p00[i]==p22[i]) { p01[i] = intersection_map( p11[i], p22[i], p00[i], p11[i], ptHit[i] ); }
01048 else { p01[i] = intersection_map( p11[i], p22[i], p11[i], p00[i], ptHit[i] ); }
01049
01050 if (p11[i]==p00[i]) { p02[i] = intersection_map( p11[i], p22[i], p22[i], p00[i], ptHit[i] ); }
01051 else { p02[i] = intersection_map( p11[i], p22[i], p00[i], p22[i], ptHit[i] ); }
01052
01053 if (p22[i]==p11[i]) { p10[i] = intersection_map( p22[i], p00[i], p00[i], p11[i], ptHit[i] ); }
01054 else { p10[i] = intersection_map( p22[i], p00[i], p11[i], p00[i], ptHit[i] ); }
01055
01056 if (p11[i]==p00[i]) { p12[i] = intersection_map( p22[i], p00[i], p11[i], p22[i], ptHit[i] ); }
01057 else { p12[i] = intersection_map( p22[i], p00[i], p22[i], p11[i], ptHit[i] ); }
01058
01059 if (p22[i]==p11[i]) { p20[i] = intersection_map( p00[i], p11[i], p22[i], p00[i], ptHit[i] ); }
01060 else { p20[i] = intersection_map( p00[i], p11[i], p00[i], p22[i], ptHit[i] ); }
01061
01062 if (p00[i]==p22[i]) { p21[i] = intersection_map( p00[i], p11[i], p11[i], p22[i], ptHit[i] ); }
01063 else { p21[i] = intersection_map( p00[i], p11[i], p22[i], p11[i], ptHit[i] ); }
01064 }
01065
01066
01067
01068 tri0[0] = sqrt( pow(p01[0]-p00[0],2) + pow(p01[1]-p00[1],2) + pow(p01[2]-p00[2],2) );
01069 tri0[1] = sqrt( pow(p02[0]-p01[0],2) + pow(p02[1]-p01[1],2) + pow(p02[2]-p01[2],2) );
01070 tri0[2] = sqrt( pow(p00[0]-p02[0],2) + pow(p00[1]-p02[1],2) + pow(p00[2]-p02[2],2) );
01071
01072 tri1[0] = sqrt( pow(p11[0]-p10[0],2) + pow(p11[1]-p10[1],2) + pow(p11[2]-p10[2],2) );
01073 tri1[1] = sqrt( pow(p12[0]-p11[0],2) + pow(p12[1]-p11[1],2) + pow(p12[2]-p11[2],2) );
01074 tri1[2] = sqrt( pow(p10[0]-p12[0],2) + pow(p10[1]-p12[1],2) + pow(p10[2]-p12[2],2) );
01075
01076 tri2[0] = sqrt( pow(p21[0]-p20[0],2) + pow(p21[1]-p20[1],2) + pow(p21[2]-p20[2],2) );
01077 tri2[1] = sqrt( pow(p22[0]-p21[0],2) + pow(p22[1]-p21[1],2) + pow(p22[2]-p21[2],2) );
01078 tri2[2] = sqrt( pow(p20[0]-p22[0],2) + pow(p20[1]-p22[1],2) + pow(p20[2]-p22[2],2) );
01079
01080
01081
01082 s0 = .5*(tri0[0] + tri0[1] + tri0[2]);
01083 s1 = .5*(tri1[0] + tri1[1] + tri1[2]);
01084 s2 = .5*(tri2[0] + tri2[1] + tri2[2]);
01085
01086 A0 = sqrt( s0*(s0-tri0[0])*(s0-tri0[1])*(s0-tri0[2]) );
01087 A1 = sqrt( s1*(s1-tri1[0])*(s1-tri1[1])*(s1-tri1[2]) );
01088 A2 = sqrt( s2*(s2-tri2[0])*(s2-tri2[1])*(s2-tri2[2]) );
01089
01090
01091
01092 triOrig[0] = sqrt( pow(p11[0]-p00[0],2) + pow(p11[1]-p00[1],2) + pow(p11[2]-p00[2],2) );
01093 triOrig[1] = sqrt( pow(p22[0]-p11[0],2) + pow(p22[1]-p11[1],2) + pow(p22[2]-p11[2],2) );
01094 triOrig[2] = sqrt( pow(p00[0]-p22[0],2) + pow(p00[1]-p22[1],2) + pow(p00[2]-p22[2],2) );
01095
01096 sOrig = .5*(triOrig[0] + triOrig[1] + triOrig[2]);
01097 Aorig = sqrt( sOrig*(sOrig-triOrig[0])*(sOrig-triOrig[1])*(sOrig-triOrig[2]) );
01098
01099
01100 weight = (float *)SUMA_calloc( 3, sizeof(float) );
01101 weight[0] = (Aorig-A0)/Aorig; weight[1] = (Aorig-A1)/Aorig; weight[2] = (Aorig-A2)/Aorig;
01102 wsum = weight[0] + weight[1] + weight[2];
01103 weight[0] = weight[0]/wsum; weight[1] = weight[1]/wsum; weight[2] = weight[2]/wsum;
01104
01105
01106
01107 SUMA_RETURN (weight);
01108
01109 }
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122 SUMA_Boolean SUMA_binSearch( float *nodeList, float target, int *seg) {
01123
01124 int mid=0;
01125 int beg = seg[0], end = seg[1];
01126 SUMA_Boolean found=YUP;
01127 static char FuncName[]={"SUMA_binSearch"};
01128
01129 SUMA_ENTRY;
01130
01131 if ( end<beg) {
01132 fprintf(SUMA_STDERR, "Error %s: Segment must be passed with seg[0] being of lower index of seg[1].\n\n", FuncName);
01133 SUMA_RETURN (found = NOPE);
01134 }
01135 if ( nodeList[end]<nodeList[beg] ) {
01136 fprintf(SUMA_STDERR, "Error %s: Nodelist must be passed sorted and in ascending order.\n\n", FuncName);
01137 SUMA_RETURN (found = NOPE);
01138 }
01139 if ( (nodeList[beg]>target) || (nodeList[end]<target) ) {
01140 fprintf(SUMA_STDERR, "Error %s: Target does not lie within segment!\n\n", FuncName);
01141 SUMA_RETURN (found = NOPE);
01142 }
01143
01144 if (beg!=end) {
01145 mid =(end-beg)/2 + beg;
01146
01147 if (beg+1==end) {
01148 seg[0] = beg;
01149 seg[1] = end;
01150 }
01151 else if (target==nodeList[mid]) {
01152 seg[0] = mid;
01153 seg[1] = mid;
01154 }
01155
01156 else if ( target < nodeList[mid]) {
01157 seg[0] = beg; seg[1] = mid;
01158 found = SUMA_binSearch( nodeList, target, seg);
01159 }
01160 else if ( target > nodeList[mid]) {
01161 seg[0] = mid; seg[1] = end;
01162 found = SUMA_binSearch( nodeList, target, seg);
01163 }
01164 }
01165
01166 else {
01167 seg[0] = mid;
01168 seg[1] = mid;
01169 }
01170
01171 SUMA_RETURN(found);
01172 }
01173
01174
01175 float intersection_map(float a, float b, float c, float d, float val) {
01176
01177 float sol = (val*(c-d) - d*(a-b)) / (c+b-a-d);
01178
01179 return sol;
01180 }
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197 SUMA_MorphInfo * SUMA_MapSurface (SUMA_SurfaceObject *surf1, SUMA_SurfaceObject *surf2)
01198 {
01199 static char FuncName[]={"SUMA_MapSurface"};
01200
01201
01202 int numNodes_1=0, numFace_1=0;
01203 float *nodeList_1=NULL, *ctrNodeList_1=NULL;
01204 int *faceList_1=NULL;
01205
01206
01207 int numNodes_2=0, numFace_2=0;
01208 float *nodeList_2=NULL, *ctrNodeList_2=NULL;
01209 int *faceList_2=NULL;
01210
01211 int i=0, j=0, k=0, m=0, j_srtd;
01212 float *weight=NULL;
01213 int *clsNodes=NULL;
01214 SUMA_MorphInfo *MI;
01215 float ctr1[3], ctr2[3], zero[3], r2, dist_tmp;
01216 float *justX_2=NULL, *justX_1=NULL, *srtdX_ctrNodeList_2=NULL;
01217 int *i_SrtdX_2=NULL;
01218 int N_outliers;
01219 float currNode[3], ptHit[3], currDist=0, avgDist=0.0, pi=3.14159265359;
01220 int seg[2], i_node[3];
01221 float min_dist[3], curr_restr;
01222
01223 SUMA_Boolean found=NOPE;
01224 float *triNode0, *triNode1, *triNode2, weight_tot;
01225 SUMA_SO_map *SO=NULL;
01226 SUMA_Boolean LocalHead = NOPE;
01227
01228 SUMA_ENTRY;
01229
01230 MI = SUMA_Create_MorphInfo();
01231 if (MI == NULL) {
01232 fprintf (SUMA_STDERR,"Error %s: Failed to allocate for MorphInfo.\n", FuncName);
01233 SUMA_RETURN (NULL);
01234 }
01235
01236
01237 nodeList_1 = surf1->NodeList;
01238 faceList_1 = surf1->FaceSetList;
01239 numNodes_1 = surf1->N_Node;
01240 numFace_1 = surf1->N_FaceSet;
01241
01242
01243 nodeList_2 = surf2->NodeList;
01244 faceList_2 = surf2->FaceSetList;
01245 numNodes_2 = surf2->N_Node;
01246 numFace_2 = surf2->N_FaceSet;
01247
01248 clsNodes = (int *)SUMA_calloc( 3*numNodes_1, sizeof(int) );
01249 weight = (float *)SUMA_calloc( 3*numNodes_1, sizeof(float) );
01250 if (!clsNodes || !weight) {
01251 if (clsNodes) SUMA_free(clsNodes);
01252 if (weight) SUMA_free(weight);
01253 fprintf (SUMA_STDERR,"Error %s: Failed to allocate for clsNodes || weight.\n", FuncName);
01254 SUMA_RETURN (NULL);
01255 }
01256
01257
01258
01259
01260
01261 ctr1[0]=0; ctr1[1]=0; ctr1[2]=0;
01262 ctr2[0]=0; ctr2[1]=0; ctr2[2]=0;
01263 zero[0]=0; zero[1]=0; zero[2]=0;
01264
01265 for (i=0; i<numNodes_1; ++i) {
01266 j = 3*i;
01267 ctr1[0] = ctr1[0] + nodeList_1[j];
01268 ctr1[1] = ctr1[1] + nodeList_1[j+1];
01269 ctr1[2] = ctr1[2] + nodeList_1[j+2];
01270 }
01271 ctr1[0] = ctr1[0]/numNodes_1;
01272 ctr1[1] = ctr1[1]/numNodes_1;
01273 ctr1[2] = ctr1[2]/numNodes_1;
01274
01275 for (i=0; i<numNodes_2; ++i) {
01276 j = 3*i;
01277 ctr2[0] = ctr2[0] + nodeList_2[j];
01278 ctr2[1] = ctr2[1] + nodeList_2[j+1];
01279 ctr2[2] = ctr2[2] + nodeList_2[j+2];
01280 }
01281 ctr2[0] = ctr2[0]/numNodes_2;
01282 ctr2[1] = ctr2[1]/numNodes_2;
01283 ctr2[2] = ctr2[2]/numNodes_2;
01284
01285
01286 zero[0] = ctr2[0];
01287 zero[1] = ctr2[1];
01288 zero[2] = ctr2[2];
01289
01290 ctrNodeList_1 = (float *) SUMA_calloc( 3*numNodes_1, sizeof(float) );
01291 ctrNodeList_2 = (float *) SUMA_calloc( 3*numNodes_2, sizeof(float) );
01292 if (!ctrNodeList_1 || !ctrNodeList_2) {
01293 if (ctrNodeList_1) SUMA_free(ctrNodeList_1);
01294 if (ctrNodeList_2) SUMA_free(ctrNodeList_2);
01295 if (clsNodes) SUMA_free(clsNodes);
01296 if (weight) SUMA_free(weight);
01297 if (i_SrtdX_2) SUMA_free(i_SrtdX_2);
01298 if (justX_2) SUMA_free(justX_2);
01299 fprintf (SUMA_STDERR,"Error %s: Failed to allocate for ctrNodeList_1 || ctrNodeList_2.\n", FuncName);
01300 SUMA_RETURN (NULL);
01301 }
01302
01303
01304 for (i=0; i<numNodes_1; ++i) {
01305 j = 3*i;
01306 ctrNodeList_1[j] = nodeList_1[j] - ctr1[0] + zero[0];
01307 ctrNodeList_1[j+1] = nodeList_1[j+1] - ctr1[1] + zero[1];
01308 ctrNodeList_1[j+2] = nodeList_1[j+2] - ctr1[2] + zero[2];
01309 }
01310 for (i=0; i<numNodes_2; ++i) {
01311 j = 3*i;
01312 ctrNodeList_2[j] = nodeList_2[j] - ctr2[0] + zero[0];
01313 ctrNodeList_2[j+1] = nodeList_2[j+1] - ctr2[1] + zero[1];
01314 ctrNodeList_2[j+2] = nodeList_2[j+2] - ctr2[2] + zero[2];
01315 }
01316
01317
01318
01319 r2 = 0.0;
01320 for (i=0; i<numNodes_2; ++i) {
01321 j = 3*i;
01322 r2 = r2 +
01323 sqrt( pow( ctrNodeList_2[j]-zero[0], 2) + pow( ctrNodeList_2[j+1]-zero[1], 2) + pow( ctrNodeList_2[j+2]-zero[2], 2) );
01324 }
01325 r2 /= numNodes_2;
01326
01327 avgDist = (4*pi*pow(r2,2))/numNodes_2;
01328
01329
01330
01331 N_outliers = 0;
01332 for (i=0; i<numNodes_2; ++i) {
01333 j = 3*i;
01334 dist_tmp = sqrt( pow( ctrNodeList_2[j]-zero[0], 2) + pow( ctrNodeList_2[j+1]-zero[1], 2)
01335 + pow( ctrNodeList_2[j+2]-zero[2], 2) );
01336 if ( abs(dist_tmp-r2)>r2/10) {
01337
01338 if ( N_outliers>(numNodes_2/1000)) {
01339
01340 fprintf(SUMA_STDERR, "\nError %s: Too many outliers. Surface considered to be non-spherical.\n\n", FuncName);
01341 SUMA_RETURN(NULL);
01342 }
01343 fprintf(SUMA_STDERR, "Warning %s: Outlier detected! Resetting to lie on sphere...\n", FuncName);
01344 N_outliers = N_outliers+1;
01345 ctrNodeList_2[j] = (r2/dist_tmp)*ctrNodeList_2[j];
01346 ctrNodeList_2[j+1] = (r2/dist_tmp)*ctrNodeList_2[j+1];
01347 ctrNodeList_2[j+2] = (r2/dist_tmp)*ctrNodeList_2[j+2];
01348 }
01349 }
01350
01351
01352
01353
01354
01355 justX_2 = (float *) SUMA_calloc( numNodes_2, sizeof(float) );
01356 if (!justX_2 ) {
01357 fprintf (SUMA_STDERR,"Error %s: Failed to allocate for justX_2.\n", FuncName);
01358 if (ctrNodeList_1) SUMA_free(ctrNodeList_1);
01359 if (ctrNodeList_2) SUMA_free(ctrNodeList_2);
01360 if (clsNodes) SUMA_free(clsNodes);
01361 if (weight) SUMA_free(weight);
01362 SUMA_RETURN (NULL);
01363 }
01364
01365 for (i=0; i<numNodes_2; ++i) {
01366 j = 3*i;
01367 justX_2[i] = ctrNodeList_2[j];
01368 }
01369
01370
01371 i_SrtdX_2 = SUMA_z_qsort( justX_2, numNodes_2 );
01372
01373 if (!i_SrtdX_2) {
01374 fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_z_qsort.\n", FuncName);
01375 if (ctrNodeList_1) SUMA_free(ctrNodeList_1);
01376 if (ctrNodeList_2) SUMA_free(ctrNodeList_2);
01377 if (clsNodes) SUMA_free(clsNodes);
01378 if (weight) SUMA_free(weight);
01379 if (justX_2) SUMA_free(justX_2);
01380
01381 SUMA_RETURN (NULL);
01382 }
01383
01384
01385 srtdX_ctrNodeList_2 = SUMA_calloc( 3*numNodes_2, sizeof(float));
01386 for (i=0; i<numNodes_2; ++i) {
01387 j = 3*i;
01388 j_srtd = 3*i_SrtdX_2[i];
01389 srtdX_ctrNodeList_2[j] = ctrNodeList_2[j_srtd];
01390 srtdX_ctrNodeList_2[j+1] = ctrNodeList_2[j_srtd+1];
01391 srtdX_ctrNodeList_2[j+2] = ctrNodeList_2[j_srtd+2];
01392 }
01393
01394
01395
01396
01397
01398 fprintf(SUMA_STDERR,"\nComputing intersections...\n\n");
01399 ptHit[0]=0; ptHit[1]=0; ptHit[2]=0;
01400 triNode0=0; triNode1=0; triNode2=0;
01401
01402 for (i=0; i<numNodes_1; ++i) {
01403
01404 j=3*i;
01405 currNode[0]=ctrNodeList_1[j];
01406 currNode[1]=ctrNodeList_1[j+1];
01407 currNode[2]=ctrNodeList_1[j+2];
01408 currDist = sqrt( pow( currNode[0]-zero[0], 2) + pow( currNode[1]-zero[1], 2) + pow( currNode[2]-zero[2], 2) );
01409
01410
01411
01412 ptHit[0] = (r2/currDist)*currNode[0];
01413 ptHit[1] = (r2/currDist)*currNode[1];
01414 ptHit[2] = (r2/currDist)*currNode[2];
01415
01416
01417
01418
01419
01420 found = NOPE;
01421 for (k=0; k<3; ++k) {
01422 min_dist[k] = 2*r2;
01423 i_node[k] = -1;
01424 }
01425 curr_restr = (float)12.0*avgDist;
01426
01427
01428
01429 seg[0] = 0;
01430 seg[1] = numNodes_2-1;
01431
01432 if ( ptHit[0] < justX_2[seg[0]] )
01433 seg[1] = seg[0];
01434 else if ( ptHit[0] > justX_2[seg[1]] )
01435 seg[0] = seg[1];
01436 else {
01437 if ( !SUMA_binSearch( justX_2, ptHit[0], seg )) {
01438 fprintf(SUMA_STDERR, "Error %s: Failed in binary search !(%f < %f < %f).\n\n", FuncName, justX_2[seg[0]], ptHit[0], justX_2[seg[1]]);
01439 if (ctrNodeList_1) SUMA_free(ctrNodeList_1);
01440 if (ctrNodeList_2) SUMA_free(ctrNodeList_2);
01441 if (clsNodes) SUMA_free(clsNodes);
01442 if (weight) SUMA_free(weight);
01443 if (i_SrtdX_2) SUMA_free(i_SrtdX_2);
01444 if (justX_2) SUMA_free(justX_2);
01445 if (srtdX_ctrNodeList_2) SUMA_free(srtdX_ctrNodeList_2);
01446 SUMA_RETURN (NULL);
01447 }
01448 }
01449
01450
01451 while ( (ptHit[0] - srtdX_ctrNodeList_2[3*seg[0]]) < curr_restr && seg[0]>0) {
01452 if ( seg[0]>10 ) seg[0] = seg[0]-10;
01453 else --seg[0];
01454 }
01455 while ( (srtdX_ctrNodeList_2[3*seg[1]] - ptHit[0]) < curr_restr && seg[1]<(numNodes_2-1) ) {
01456 if ( seg[1]<(numNodes_2-11) ) seg[1] = seg[1]+10;
01457 else ++seg[1];
01458 }
01459
01460
01461 while ( !found && seg[1]-seg[0]<numNodes_2 && curr_restr<3*r2 ) {
01462
01463
01464 SUMA_Search_Min_Dist( ptHit, srtdX_ctrNodeList_2, seg, curr_restr, min_dist, i_node );
01465
01466 if ( i_node[0]==-1 || i_node[1]==-1 || i_node[2]==-1 ) {
01467
01468 curr_restr = (float) 1.5*curr_restr;
01469 found = NOPE;
01470 while ( ptHit[0] - srtdX_ctrNodeList_2[3*seg[0]] < curr_restr && seg[0]>0) {
01471 if (seg[0]>10) seg[0] = seg[0]-10;
01472 else --seg[0];
01473 }
01474 while ( srtdX_ctrNodeList_2[3*seg[1]] - ptHit[0] < curr_restr && seg[1]<numNodes_2-1) {
01475 if (k<numNodes_2-11) seg[1] = seg[1]+10;
01476 else ++seg[1];
01477 }
01478 }
01479 else found = YUP;
01480 }
01481
01482
01483 if ( i_node[0]==-1 || i_node[1]==-1 || i_node[2]==-1 ) {
01484
01485 fprintf(SUMA_STDERR, "Error %s: Unable to acquire 3 closest nodes ?!?\n\n", FuncName);
01486 if (ctrNodeList_1) SUMA_free(ctrNodeList_1);
01487 if (ctrNodeList_2) SUMA_free(ctrNodeList_2);
01488 if (clsNodes) SUMA_free(clsNodes);
01489 if (weight) SUMA_free(weight);
01490 if (i_SrtdX_2) SUMA_free(i_SrtdX_2);
01491 if (justX_2) SUMA_free(justX_2);
01492 if (srtdX_ctrNodeList_2) SUMA_free(srtdX_ctrNodeList_2);
01493 SUMA_RETURN (NULL);
01494 }
01495
01496
01497 i_node[0] = i_SrtdX_2[i_node[0]];
01498 i_node[1] = i_SrtdX_2[i_node[1]];
01499 i_node[2] = i_SrtdX_2[i_node[2]];
01500
01501 if (LocalHead) {
01502 fprintf(SUMA_STDERR,"----------------------------------------\n");
01503 fprintf(SUMA_STDERR, "%s: PtHit: [%f, %f, %f].\n", FuncName, ptHit[0], ptHit[1], ptHit[2]);
01504 fprintf(SUMA_STDERR, "%s: Node %d [%f, %f, %f], distances %f.\n",
01505 FuncName, i_node[0], ctrNodeList_2[3*i_node[0]], ctrNodeList_2[3*i_node[0]+1], ctrNodeList_2[3*i_node[0]+2], min_dist[0]);
01506 fprintf(SUMA_STDERR, "%s: Node %d [%f, %f, %f], distances %f.\n",
01507 FuncName, i_node[1], ctrNodeList_2[3*i_node[1]], ctrNodeList_2[3*i_node[1]+1], ctrNodeList_2[3*i_node[1]+2], min_dist[1]);
01508 fprintf(SUMA_STDERR, "%s: Node %d [%f, %f, %f], distances %f.\n",
01509 FuncName, i_node[2], ctrNodeList_2[3*i_node[2]], ctrNodeList_2[3*i_node[2]+1], ctrNodeList_2[3*i_node[2]+2], min_dist[2]);
01510 fprintf(SUMA_STDERR, "%s: orig ptHit (%f, %f, %f)\n", FuncName, ptHit[0], ptHit[1], ptHit[2]);
01511 fprintf(SUMA_STDERR, "%s: Trying 1- node %d\n", FuncName, i_node[0]);
01512 }
01513
01514
01515
01516
01517 if (surf2->FN == NULL) {
01518 fprintf(SUMA_STDERR, "%s: Surf2->FN is NULL.\n", FuncName);
01519 if (ctrNodeList_1) SUMA_free(ctrNodeList_1);
01520 if (ctrNodeList_2) SUMA_free(ctrNodeList_2);
01521 if (clsNodes) SUMA_free(clsNodes);
01522 if (weight) SUMA_free(weight);
01523 if (i_SrtdX_2) SUMA_free(i_SrtdX_2);
01524 if (justX_2) SUMA_free(justX_2);
01525 if (srtdX_ctrNodeList_2) SUMA_free(srtdX_ctrNodeList_2);
01526 SUMA_RETURN (NULL);
01527 }
01528
01529
01530 found = SUMA_inNodeNeighb( surf2, ctrNodeList_2, i_node, zero, ptHit);
01531
01532 if (!found) {
01533
01534 if (LocalHead) fprintf(SUMA_STDERR, "%s: Trying Brute force. (%d)\n", FuncName, i);
01535 {
01536 SUMA_MT_INTERSECT_TRIANGLE *MTI;
01537
01538 MTI = SUMA_MT_intersect_triangle(ptHit, zero, ctrNodeList_2, numNodes_2, faceList_2, numFace_2, NULL);
01539 if (MTI) {
01540 if (MTI->N_hits) {
01541 if (LocalHead) fprintf(SUMA_STDERR, "%s: Brute force-Triangle %d [%d, %d, %d] is intersected at (%f, %f, %f)\n",
01542 FuncName, MTI->ifacemin, surf2->FaceSetList[3*MTI->ifacemin], surf2->FaceSetList[3*MTI->ifacemin+1],
01543 surf2->FaceSetList[3*MTI->ifacemin+2], MTI->P[0], MTI->P[1], MTI->P[2]);
01544 found = YUP;
01545 ptHit[0] = MTI->P[0];
01546 ptHit[1] = MTI->P[1];
01547 ptHit[2] = MTI->P[2];
01548 i_node[0] = surf2->FaceSetList[3*MTI->ifacemin];
01549 i_node[1] = surf2->FaceSetList[3*MTI->ifacemin+1];
01550 i_node[2] = surf2->FaceSetList[3*MTI->ifacemin+2];
01551 }
01552 MTI = SUMA_Free_MT_intersect_triangle(MTI);
01553 }
01554 }
01555 }
01556
01557 if (!found) {
01558 fprintf(SUMA_STDERR, "Error %s: !!!!!!!!!! intersected triangle not found.\n", FuncName);
01559 if (ctrNodeList_1) SUMA_free(ctrNodeList_1);
01560 if (ctrNodeList_2) SUMA_free(ctrNodeList_2);
01561 if (clsNodes) SUMA_free(clsNodes);
01562 if (weight) SUMA_free(weight);
01563 if (i_SrtdX_2) SUMA_free(i_SrtdX_2);
01564 if (justX_2) SUMA_free(justX_2);
01565 if (srtdX_ctrNodeList_2) SUMA_free(srtdX_ctrNodeList_2);
01566 SUMA_RETURN (NULL);
01567 }
01568
01569 if (LocalHead) fprintf (SUMA_STDERR, "%s: (%d : %d : %d)\n ptHit(%f, %f, %f)\n", FuncName, i_node[0], i_node[1], i_node[2], ptHit[0], ptHit[1], ptHit[2]);
01570
01571
01572 clsNodes[j] = i_node[0]; clsNodes[j+1] = i_node[1]; clsNodes[j+2] = i_node[2];
01573
01574
01575 triNode0 = &(ctrNodeList_2[ 3*i_node[0] ]);
01576 triNode1 = &(ctrNodeList_2[ 3*i_node[1] ]);
01577 triNode2 = &(ctrNodeList_2[ 3*i_node[2] ]);
01578
01579
01580 SUMA_TRI_AREA( ptHit, triNode1, triNode2, weight[j]);
01581 SUMA_TRI_AREA( ptHit, triNode0, triNode2, weight[j+1]);
01582 SUMA_TRI_AREA( ptHit, triNode0, triNode1, weight[j+2]);
01583
01584
01585
01586 weight_tot = weight[j] + weight[j+1] + weight[j+2];
01587 if (weight_tot) {
01588 weight[j] /= weight_tot;
01589 weight[j+1] /= weight_tot;
01590 weight[j+2] /= weight_tot;
01591 }else {
01592 weight[j] = weight[j+1] = weight[j+2] = 1.0/3.0;
01593 }
01594
01595 }
01596
01597 MI->N_Node = numNodes_1;
01598 MI->N_FaceSet = numFace_1;
01599 MI->Weight = weight;
01600 MI->ClsNodes = clsNodes;
01601 MI->FaceSetList = (int *) SUMA_calloc( 3*numFace_1, sizeof(int));
01602 if (!MI->FaceSetList) {
01603 fprintf(SUMA_STDERR, "Error %s: Failed to allocate for MI->FaceSetList.\n", FuncName);
01604 if (ctrNodeList_1) SUMA_free(ctrNodeList_1);
01605 if (ctrNodeList_2) SUMA_free(ctrNodeList_2);
01606 if (clsNodes) SUMA_free(clsNodes);
01607 if (weight) SUMA_free(weight);
01608 if (i_SrtdX_2) SUMA_free(i_SrtdX_2);
01609 if (justX_2) SUMA_free(justX_2);
01610 if (srtdX_ctrNodeList_2) SUMA_free(srtdX_ctrNodeList_2);
01611 SUMA_RETURN (NULL);
01612 }
01613 for (i=0; i<numFace_1; ++i) {
01614 j = 3*i;
01615 MI->FaceSetList[j] = faceList_1[j];
01616 MI->FaceSetList[j+1] = faceList_1[j+1];
01617 MI->FaceSetList[j+2] = faceList_1[j+2];
01618 }
01619
01620 if (ctrNodeList_1) SUMA_free(ctrNodeList_1);
01621 if (ctrNodeList_2) SUMA_free(ctrNodeList_2);
01622 if (i_SrtdX_2) SUMA_free(i_SrtdX_2);
01623 if (justX_2) SUMA_free(justX_2);
01624 if (srtdX_ctrNodeList_2) SUMA_free(srtdX_ctrNodeList_2);
01625
01626 SUMA_RETURN (MI);
01627 }
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643 void SUMA_Search_Min_Dist( float* pt, float* nodeList, int* seg, float restr, float *dist, int *i_dist ) {
01644
01645 static char FuncName[]={"SUMA_Search_Min_Dist"};
01646 float tempD;
01647 int j, k;
01648
01649 if ( !dist[0] || !dist[1] || !dist[2] ) {
01650 tempD = 3*pow(restr,2);
01651 dist[0] = tempD; dist[1] = tempD; dist[2] = tempD;
01652 i_dist[0] = -1; i_dist[1] = -1; i_dist[2] = -1;
01653 }
01654 else tempD = dist[2]+1;
01655
01656 for (k=seg[0]; k<=seg[1]; ++k) {
01657 j = 3*k;
01658 if (pt[0]-nodeList[j] < restr) {
01659 if (pt[0]-nodeList[j] > -restr) {
01660 if (pt[1]-nodeList[j+1] < restr) {
01661 if (pt[1]-nodeList[j+1] > -restr) {
01662 if (pt[2]-nodeList[j+2] < restr) {
01663 if (pt[2]-nodeList[j+2] > -restr) {
01664
01665 tempD = sqrt( pow(pt[0]-nodeList[j],2) + pow(pt[1]-nodeList[j+1],2) +
01666 pow(pt[2]-nodeList[j+2],2) );
01667
01668 if (tempD < dist[2]) {
01669 if (tempD < dist[1]) {
01670 if (tempD < dist[0]) {
01671 dist[2] = dist[1]; i_dist[2] = i_dist[1];
01672 dist[1] = dist[0]; i_dist[1] = i_dist[0];
01673 dist[0] = tempD; i_dist[0] = k;
01674 }
01675 else {
01676 dist[2] = dist[1]; i_dist[2] = i_dist[1];
01677 dist[1] = tempD; i_dist[1] = k;
01678 }
01679 }
01680 else {
01681 dist[2] = tempD; i_dist[2] = k;
01682 }
01683 }
01684 }
01685 }
01686 }
01687 }
01688 }
01689 }
01690 }
01691
01692 SUMA_RETURNe;
01693 }
01694
01695
01696
01697
01698
01699
01700 SUMA_SO_map *SUMA_Create_SO_map (void)
01701 {
01702 static char FuncName[]={"SUMA_Create_SO_map"};
01703 SUMA_SO_map *SOM = NULL;
01704
01705 SUMA_ENTRY;
01706
01707 SOM = (SUMA_SO_map *) SUMA_malloc (sizeof(SUMA_SO_map));
01708 if (!SOM) {
01709 fprintf (SUMA_STDERR, "Error %s: Failed to allocate for SOM.\n", FuncName);
01710 SUMA_RETURN (NULL);
01711 }
01712
01713 SOM->N_Node = 0;
01714 SOM->NewNodeList = NULL;
01715 SOM->NodeVal = NULL;
01716 SOM->NodeDisp = NULL;
01717 SOM->NodeCol = NULL;
01718
01719 SUMA_RETURN (SOM);
01720 }
01721
01722
01723
01724
01725 SUMA_Boolean SUMA_Free_SO_map (SUMA_SO_map *SOM)
01726 {
01727 static char FuncName[]={"SUMA_Free_SO_map"};
01728
01729 SUMA_ENTRY;
01730
01731 if (!SOM) {
01732 SUMA_RETURN (YUP);
01733 }
01734
01735 if (SOM->NewNodeList) SUMA_free (SOM->NewNodeList);
01736 if (SOM->NodeVal) SUMA_free (SOM->NodeVal);
01737 if (SOM->NodeDisp) SUMA_free (SOM->NodeDisp);
01738 if (SOM->NodeCol) SUMA_free(SOM->NodeCol);
01739
01740 SUMA_free (SOM);
01741
01742 SUMA_RETURN (YUP);
01743 }
01744
01745
01746
01747
01748 SUMA_Boolean SUMA_Show_SO_map (SUMA_SO_map *SOM, FILE *out)
01749 {
01750 static char FuncName[]={"SUMA_Show_SO_map"};
01751 int i=0, imax;
01752
01753 SUMA_ENTRY;
01754
01755 if (!out) out = SUMA_STDERR;
01756
01757 fprintf (out, "\n%s: Showing contents of SUMA_SO_map structure:\n", FuncName);
01758 if (!SOM) {
01759 fprintf (out, "\tpointer is NULL.\n");
01760 SUMA_RETURN (YUP);
01761 }
01762
01763 if (SOM->N_Node > 5) imax = 5;
01764 else imax = SOM->N_Node;
01765
01766 fprintf (SUMA_STDERR, "NodeList, (1st %d elements):\n", imax);
01767 for (i=0; i<imax; ++i) {
01768 fprintf (SUMA_STDERR, "\t%f, %f, %f\n", SOM->NewNodeList[3*i], SOM->NewNodeList[3*i+1], SOM->NewNodeList[3*i+2]);
01769 }
01770
01771 SUMA_RETURN (YUP);
01772 }
01773
01774
01775
01776
01777
01778 SUMA_MorphInfo *SUMA_Create_MorphInfo (void)
01779 {
01780 static char FuncName[]={"SUMA_Create_MorphInfo"};
01781 SUMA_MorphInfo *MI = NULL;
01782
01783 SUMA_ENTRY;
01784
01785 MI = (SUMA_MorphInfo *) SUMA_malloc (sizeof(SUMA_MorphInfo));
01786 if (!MI) {
01787 fprintf (SUMA_STDERR, "Error %s: Failed to allocate for MI.\n", FuncName);
01788 SUMA_RETURN (NULL);
01789 }
01790
01791 MI->IDcode = NULL;
01792 MI->N_Node = 0;
01793 MI->N_FaceSet = 0;
01794 MI->Weight = NULL;
01795 MI->ClsNodes = NULL;
01796 MI->FaceSetList = NULL;
01797
01798 SUMA_RETURN (MI);
01799 }
01800
01801
01802
01803
01804 SUMA_Boolean SUMA_Free_MorphInfo (SUMA_MorphInfo *MI)
01805 {
01806 static char FuncName[]={"SUMA_Free_MorphInfo"};
01807
01808 SUMA_ENTRY;
01809
01810 if (!MI) {
01811 SUMA_RETURN (YUP);
01812 }
01813
01814 if (MI->IDcode) SUMA_free (MI->IDcode);
01815 if (MI->Weight) SUMA_free (MI->Weight);
01816 if (MI->ClsNodes) SUMA_free (MI->ClsNodes);
01817 if (MI->FaceSetList) SUMA_free (MI->FaceSetList);
01818
01819 SUMA_free (MI);
01820
01821 SUMA_RETURN (YUP);
01822 }
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835 SUMA_SurfaceObject* SUMA_morphToStd (SUMA_SurfaceObject *SO, SUMA_MorphInfo *MI, SUMA_Boolean nodeChk) {
01836
01837 float *newNodeList = NULL;
01838 int *tmp_newFaceSetList = NULL, *newFaceSetList = NULL, *inclNodes=NULL;
01839 int i, j, N_FaceSet;
01840 SUMA_SurfaceObject *SO_new=NULL;
01841 static char FuncName[] = {"SUMA_morphToStd"};
01842
01843 SUMA_ENTRY;
01844
01845 SO_new = SUMA_Alloc_SurfObject_Struct(1);
01846 if (SO_new == NULL) {
01847 fprintf (SUMA_STDERR,"Error %s: Failed to allocate for Surface Object.", FuncName);
01848 SUMA_RETURN (NULL);
01849 }
01850
01851 newNodeList = (float *) SUMA_calloc( 3*MI->N_Node, sizeof(float));
01852 if (!newNodeList) {
01853 fprintf (SUMA_STDERR, "Error %s: Failed to allocate. \n", FuncName);
01854 SUMA_RETURN (NULL);
01855 }
01856 N_FaceSet = 0;
01857
01858 if ( !nodeChk ) {
01859
01860 fprintf(SUMA_STDERR, "Warning %s: Assuming face sets of surface %s to contain all nodes indicated in morphing to standard mesh.\n\n",
01861 FuncName, SO->State);
01862
01863 for (i=0; i<(MI->N_Node); ++i){
01864 j = 3*i;
01865
01866 newNodeList[j] = (MI->Weight[j])*SO->NodeList[3*(MI->ClsNodes[j])] +
01867 (MI->Weight[j+1])*SO->NodeList[3*(MI->ClsNodes[j+1])] +
01868 (MI->Weight[j+2])*SO->NodeList[3*(MI->ClsNodes[j+2])];
01869 newNodeList[j+1] = (MI->Weight[j])*SO->NodeList[3*(MI->ClsNodes[j])+1] +
01870 (MI->Weight[j+1])*SO->NodeList[3*(MI->ClsNodes[j+1])+1] +
01871 (MI->Weight[j+2])*SO->NodeList[3*(MI->ClsNodes[j+2])+1];
01872 newNodeList[j+2] = (MI->Weight[j])*SO->NodeList[3*(MI->ClsNodes[j])+2] +
01873 (MI->Weight[j+1])*SO->NodeList[3*(MI->ClsNodes[j+1])+2] +
01874 (MI->Weight[j+2])*SO->NodeList[3*(MI->ClsNodes[j+2])+2];
01875 }
01876
01877 newFaceSetList = MI->FaceSetList;
01878 N_FaceSet = MI->N_FaceSet;
01879 }
01880
01881 else {
01882
01883
01884 if ( !SO->FN ) {
01885 fprintf(SUMA_STDERR, "Error %s: No First Neighbor information passed.\n", FuncName);
01886 return (NULL);
01887 }
01888
01889
01890 inclNodes = SUMA_calloc( MI->N_Node, sizeof(int));
01891 for (i=0; i<MI->N_Node; ++i) {
01892 inclNodes[i] = 0;
01893 }
01894
01895 for (i=0; i<(MI->N_Node); ++i) {
01896
01897 j = 3*i;
01898 if ( (MI->ClsNodes[j])<=(SO->N_Node) && (MI->ClsNodes[j+1])<=(SO->N_Node) && (MI->ClsNodes[j+2])<=(SO->N_Node) ) {
01899
01900 if ( SO->FN->N_Neighb[MI->ClsNodes[j]]>0 && SO->FN->N_Neighb[MI->ClsNodes[j+1]]>0 && SO->FN->N_Neighb[MI->ClsNodes[j+2]]>0 ) {
01901
01902
01903 inclNodes[i] = 1;
01904 newNodeList[j] = (MI->Weight[j])*SO->NodeList[3*(MI->ClsNodes[j])] +
01905 (MI->Weight[j+1])*SO->NodeList[3*(MI->ClsNodes[j+1])] +
01906 (MI->Weight[j+2])*SO->NodeList[3*(MI->ClsNodes[j+2])];
01907 newNodeList[j+1] = (MI->Weight[j])*SO->NodeList[3*(MI->ClsNodes[j])+1] +
01908 (MI->Weight[j+1])*SO->NodeList[3*(MI->ClsNodes[j+1])+1] +
01909 (MI->Weight[j+2])*SO->NodeList[3*(MI->ClsNodes[j+2])+1];
01910 newNodeList[j+2] = (MI->Weight[j])*SO->NodeList[3*(MI->ClsNodes[j])+2] +
01911 (MI->Weight[j+1])*SO->NodeList[3*(MI->ClsNodes[j+1])+2] +
01912 (MI->Weight[j+2])*SO->NodeList[3*(MI->ClsNodes[j+2])+2];
01913 }
01914 }
01915
01916 }
01917
01918
01919 tmp_newFaceSetList = SUMA_calloc( 3*MI->N_FaceSet, sizeof(int));
01920 if (!tmp_newFaceSetList) {
01921 fprintf (SUMA_STDERR, "Error %s: Failed to allocate. \n", FuncName);
01922 SUMA_RETURN (NULL);
01923 }
01924
01925 for (i=0; i<MI->N_FaceSet; ++i) {
01926 j = 3*i;
01927 if ( inclNodes[MI->FaceSetList[j]]==1 && inclNodes[MI->FaceSetList[j+1]]==1 &&
01928 inclNodes[MI->FaceSetList[j+2]]==1) {
01929
01930 tmp_newFaceSetList[3*N_FaceSet] = MI->FaceSetList[j];
01931 tmp_newFaceSetList[3*N_FaceSet+1] = MI->FaceSetList[j+1];
01932 tmp_newFaceSetList[3*N_FaceSet+2] = MI->FaceSetList[j+2];
01933 N_FaceSet++;
01934 }
01935 }
01936
01937
01938 if ( N_FaceSet == MI->N_FaceSet ) {
01939
01940 newFaceSetList = tmp_newFaceSetList;
01941 }
01942 else {
01943
01944 newFaceSetList = SUMA_calloc( 3*N_FaceSet, sizeof(int));
01945 if (!newFaceSetList) {
01946 fprintf (SUMA_STDERR, "Error %s: Failed to allocate. \n", FuncName);
01947 SUMA_RETURN (NULL);
01948 }
01949 for (i=0; i<3*N_FaceSet; ++i) {
01950 newFaceSetList[i] = tmp_newFaceSetList[i];
01951 }
01952 SUMA_free (tmp_newFaceSetList);
01953 }
01954 SUMA_free (inclNodes);
01955 }
01956
01957
01958 SO_new->NodeList = newNodeList;
01959 SO_new->FaceSetList = newFaceSetList;
01960 SO_new->N_Node = MI->N_Node;
01961 SO_new->N_FaceSet = N_FaceSet;
01962 SO_new->NodeDim = 3;
01963 SO_new->FaceSetDim = 3;
01964 SO_new->idcode_str = (char *)SUMA_calloc (SUMA_IDCODE_LENGTH, sizeof(char));
01965 UNIQ_idcode_fill (SO_new->idcode_str);
01966
01967 SUMA_RETURN( SO_new );
01968 }
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980 float* SUMA_readColor (int numNodes, char* colFileNm) {
01981
01982 float *colArray=NULL;
01983 FILE *colFile=NULL;
01984 char *line=NULL, *temp=NULL;
01985 int i=0, j=0, k=0, index=0;
01986 static char FuncName[]={"SUMA_readColor"};
01987
01988 SUMA_ENTRY;
01989
01990 colArray = (float *) SUMA_calloc( 3*numNodes, sizeof(float) );
01991 line = (char *) SUMA_calloc( 10000, sizeof(char));
01992 temp = (char *) SUMA_calloc( 10000, sizeof(char));
01993
01994 if( (colFile = fopen(colFileNm, "r"))==NULL) {
01995 fprintf (SUMA_STDERR, "Failed in opening %s for reading.\n", colFileNm);
01996 if (colArray) SUMA_free (colArray);
01997 if (line) SUMA_free (line);
01998 if (temp) SUMA_free (temp);
01999 exit(1);
02000 }
02001 else {
02002 fgets( line, 1000, colFile);
02003 while( !feof(colFile) ) {
02004
02005 j = 3*index;
02006 i = 0;
02007 while ( isdigit(line[i]) ) ++i;
02008
02009 ++i; k=0;
02010 while ( !isspace(line[i])) {
02011 temp[k] = line[i];
02012 ++i; ++k;
02013 }
02014 colArray[j] = atof(temp);
02015 SUMA_free(temp);
02016 temp = SUMA_calloc(10000, sizeof(char));
02017
02018 ++i; k=0;
02019 while ( !isspace(line[i])) {
02020 temp[k] = line[i];
02021 ++i; ++k;
02022 }
02023 colArray[j+1] = atof(temp);
02024 SUMA_free(temp);
02025 temp = SUMA_calloc( 10000, sizeof(char));
02026
02027 ++i; k=0;
02028 while ( !isspace(line[i])) {
02029 temp[k] = line[i];
02030 ++i; ++k;
02031 }
02032 colArray[j+2] = atof(temp);
02033 SUMA_free(temp);
02034 temp = SUMA_calloc( 10000, sizeof(char));
02035
02036 fgets( line, 10000, colFile );
02037 ++index;
02038 }
02039 }
02040 SUMA_free(line);
02041 SUMA_free(temp);
02042
02043 SUMA_RETURN( colArray);
02044 }
02045
02046
02047
02048
02049
02050
02051
02052
02053
02054
02055
02056
02057 void SUMA_writeColorFile (float *array, int numNode, int *index, char fileNm[]) {
02058
02059 FILE *outFile=NULL;
02060 int i=0, j=0;
02061 static char FuncName[] = {"SUMA_writeColorFile"};
02062
02063 SUMA_ENTRY;
02064
02065 for (i=0; i<numNode; ++i) {
02066 j = 3*i;
02067 }
02068
02069 for (i=0; i<numNode; ++i) {
02070 j = 3*i;
02071 }
02072
02073 if((outFile = fopen(fileNm, "w"))==NULL) {
02074 fprintf(SUMA_STDERR, "Could not open file %s.\n", fileNm);
02075 exit(1);
02076 }
02077 else {
02078 if (index!=NULL) {
02079
02080 for (i=0; i<numNode; ++i) {
02081 j = 3*i;
02082 fprintf (outFile, "%d\t%f\t%f\t%f\n", index[i], array[j], array[j+1], array[j+2]);
02083 }
02084 }
02085 else {
02086
02087 for (i=0; i < numNode; ++i) {
02088 j = i*3;
02089 fprintf (outFile, "%d\t%f\t%f\t%f\n", i, array[j], array[j+1], array[j+2]);
02090 }
02091 }
02092 fclose (outFile);
02093 }
02094 SUMA_RETURNe;
02095 }
02096
02097
02098
02099
02100
02101
02102
02103
02104
02105
02106
02107
02108
02109
02110
02111 void SUMA_writeFSfile (SUMA_SurfaceObject *SO, char firstLine[], char fileNm[]) {
02112
02113 FILE *outFile=NULL;
02114 int i=0, j=0;
02115 static char FuncName[]={"SUMA_writeFSfile"};
02116
02117 SUMA_ENTRY;
02118
02119 outFile = fopen(fileNm, "w");
02120 if (!outFile) {
02121 fprintf (SUMA_STDERR, "Error %s: Failed in opening %s for writing.\n",FuncName, fileNm);
02122 exit(1);
02123 }
02124 else {
02125 if ( firstLine!=NULL )
02126 fprintf (outFile,"# %s\n", firstLine);
02127 else fprintf (outFile, "#\n");
02128 fprintf (outFile, "%d %d\n", SO->N_Node, SO->N_FaceSet);
02129
02130 j=0;
02131 for (i=0; i<SO->N_Node; ++i) {
02132 j=3*i;
02133 fprintf (outFile, "%f %f %f 0\n", SO->NodeList[j], SO->NodeList[j+1], SO->NodeList[j+2]);
02134 }
02135
02136 j=0;
02137 for (i=0; i<SO->N_FaceSet; ++i) {
02138 j = 3*i;
02139 fprintf (outFile, "%d %d %d 0\n", SO->FaceSetList[j], SO->FaceSetList[j+1], SO->FaceSetList[j+2]);
02140 }
02141
02142 fclose(outFile);
02143 }
02144
02145 SUMA_RETURNe;
02146 }
02147
02148
02149
02150
02151
02152
02153
02154
02155
02156
02157
02158
02159
02160
02161 void SUMA_writeSpecFile (SUMA_SpecSurfInfo *surfaces, int numSurf, char program[], char group[], char specFileNm[]) {
02162
02163 FILE *outFile=NULL;
02164 int i=0, k=0, tag=0, ifSmwm=0, p=0;
02165 static char FuncName[]={"SUMA_writeSpecFile"};
02166
02167 SUMA_ENTRY;
02168
02169 outFile = fopen(specFileNm, "w");
02170 if (!outFile) {
02171 fprintf (SUMA_STDERR, "Failed in opening %s for writing.\n", specFileNm);
02172 exit (1);
02173 }
02174 else {
02175 fprintf (outFile, "# %s spec file for %s\n\n", program, group);
02176 fprintf (outFile, "#define the group\n\tGroup = %s\n\n", group);
02177 fprintf (outFile, "#define various States\n");
02178 for (i=0; i<numSurf; ++i) {
02179 tag = 0;
02180 for (k=0; k<i; ++k) {
02181 if ( strcmp( surfaces[k].state, surfaces[i].state ) == 0) tag = -1;
02182 }
02183 if (tag==0) {
02184 fprintf( outFile, "\tStateDef = %s\n", surfaces[i].state);
02185 }
02186 }
02187
02188 for (i=0; i<numSurf; ++i) {
02189 fprintf (outFile, "\nNewSurface\n\tSurfaceFormat = %s\n\tSurfaceType = %s\n", surfaces[i].format, surfaces[i].type);
02190 fprintf (outFile, "\tFreeSurferSurface = %s\n\tLocalDomainParent = %s\n", surfaces[i].fileToRead, surfaces[i].mapRef );
02191 fprintf (outFile, "\tSurfaceState = %s\n\tEmbedDimension = %s\n", surfaces[i].state, surfaces[i].dim);
02192 }
02193
02194 fclose(outFile);
02195 }
02196 SUMA_RETURNe;
02197 }
02198
02199 #ifdef SUMA_CreateIcosahedron_STAND_ALONE
02200
02201 void SUMA_CreateIcosahedron_usage ()
02202
02203 {
02204 static char FuncName[]={"SUMA_CreateIcosahedron_usage"};
02205 char * s = NULL;
02206
02207 printf ( "\n"
02208 "Usage: CreateIcosahedron [-rad r] [-rd recDepth] [-ld linDepth] \n"
02209 " [-ctr ctr] [-prefix fout] [-help]\n"
02210 "\n"
02211 " -rad r: size of icosahedron. (optional, default 100)\n"
02212 "\n"
02213 " -rd recDepth: recursive (binary) tesselation depth for icosahedron \n"
02214 " (optional, default:3) \n"
02215 " (recommended to approximate number of nodes in brain: 6\n"
02216 " let rd2 = 2 * recDepth\n"
02217 " Nvert = 2 + 10 * 2^rd2\n"
02218 " Ntri = 20 * 2^rd2\n"
02219 " Nedge = 30 * 2^rd2\n"
02220 "\n"
02221 " -ld linDepth: number of edge divides for linear icosahedron tesselation\n"
02222 " (optional, default uses binary tesselation).\n"
02223 " Nvert = 2 + 10 * linDepth^2\n"
02224 " Ntri = 20 * linDepth^2\n"
02225 " Nedge = 30 * linDepth^2\n"
02226 "\n"
02227 " -nums: output the number of nodes (vertices), triangles, edges, total volume and total area then quit\n"
02228 "\n"
02229 " -nums_quiet: same as -nums but less verbose. For the machine in you.\n"
02230 "\n"
02231 " -ctr ctr: coordinates of center of icosahedron. \n"
02232 " (optional, default 0,0,0)\n"
02233 "\n"
02234 " -tosphere: project nodes to sphere.\n"
02235 "\n"
02236 " -prefix fout: prefix for output files. \n"
02237 " (optional, default CreateIco)\n"
02238 "\n"
02239 " -help: help message\n"
02240 "\n");
02241 s = SUMA_New_Additions(0, 1); printf("%s\n", s);SUMA_free(s); s = NULL;
02242 printf ("\n"
02243 " Brenna D. Argall LBC/NIMH/NIH bargall@codon.nih.gov \n"
02244 " Ziad S. Saad SSC/NIMH/NIH ziad@nih.gov\n");
02245 exit (0);
02246 }
02247
02248
02249
02250
02251 int main (int argc, char *argv[])
02252 {
02253
02254 static char FuncName[]={"SUMA_CreateIcosahedron-main"};
02255 int kar, depth, i, j;
02256 float r, ctr[3], a, b, lgth, A = 0.0, V = 0.0;
02257 SUMA_SurfaceObject *SO=NULL;
02258 SUMA_Boolean brk;
02259 int NumOnly, ToSphere;
02260 SUMA_Boolean LocalHead = NOPE;
02261 char fout[SUMA_MAX_DIR_LENGTH+SUMA_MAX_NAME_LENGTH];
02262 char bin[SUMA_MAX_DIR_LENGTH+SUMA_MAX_NAME_LENGTH];
02263 char surfFileNm[1000], outSpecFileNm[1000];
02264 SUMA_SpecSurfInfo *surfaces;
02265
02266
02267 if (LocalHead) fprintf (SUMA_STDERR,"%s: Calling SUMA_Create_CommonFields ...\n", FuncName);
02268
02269 SUMAg_CF = SUMA_Create_CommonFields ();
02270 if (SUMAg_CF == NULL) {
02271 fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_Create_CommonFields\n", FuncName);
02272 exit(1);
02273 }
02274 if (LocalHead) fprintf (SUMA_STDERR,"%s: SUMA_Create_CommonFields Done.\n", FuncName);
02275
02276
02277
02278 r = 100;
02279 depth = 3;
02280 ctr[0] = 0; ctr[1] = 0; ctr[2] = 0;
02281 sprintf (fout, "%s", "CreateIco");
02282 sprintf (bin, "%s", "y");
02283 NumOnly = 0;
02284 ToSphere = 0;
02285 kar = 1;
02286 brk = NOPE;
02287 while (kar < argc) {
02288 if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
02289 SUMA_CreateIcosahedron_usage ();
02290 exit (1);
02291 }
02292
02293 if (!brk && (strcmp(argv[kar], "-rad") == 0 ))
02294 {
02295 kar ++;
02296 if (kar >= argc) {
02297 fprintf (SUMA_STDERR, "need argument after -r ");
02298 exit (1);
02299 }
02300 r = atof(argv[kar]);
02301 brk = YUP;
02302 }
02303
02304 if (!brk && (strcmp(argv[kar], "-tosphere") == 0 ))
02305 {
02306 ToSphere = 1;
02307 brk = YUP;
02308 }
02309
02310 if (!brk && (strcmp(argv[kar], "-rd") == 0 ))
02311 {
02312 kar ++;
02313 if (kar >= argc) {
02314 fprintf (SUMA_STDERR, "need argument after -rd ");
02315 exit (1);
02316 }
02317 depth = atoi(argv[kar]);
02318 sprintf (bin, "y");
02319 brk = YUP;
02320
02321 }
02322 if (!brk && (strcmp(argv[kar], "-ld") == 0 ))
02323 {
02324 kar ++;
02325 if (kar >= argc) {
02326 fprintf (SUMA_STDERR, "need argument after -ld ");
02327 exit (1);
02328 }
02329 depth = atoi(argv[kar]);
02330 sprintf (bin, "n");
02331 brk = YUP;
02332 }
02333
02334 if (!brk && (strcmp(argv[kar], "-nums") == 0 ))
02335 {
02336 NumOnly = 1;
02337 brk = YUP;
02338 }
02339 if (!brk && (strcmp(argv[kar], "-nums_quiet") == 0 ))
02340 {
02341 NumOnly = 2;
02342 brk = YUP;
02343 }
02344 if (!brk && strcmp(argv[kar], "-ctr") == 0)
02345 {
02346 kar ++;
02347 if (kar >= argc) {
02348 fprintf (SUMA_STDERR, "need argument after -ctr ");
02349 exit (1);
02350 }
02351 ctr[0] = atof(argv[kar]); kar ++;
02352 ctr[1] = atof(argv[kar]); kar ++;
02353 ctr[2] = atof(argv[kar]);
02354
02355 brk = YUP;
02356 }
02357
02358 if (!brk && strcmp(argv[kar], "-prefix") == 0)
02359 {
02360 kar ++;
02361 if (kar >= argc) {
02362 fprintf (SUMA_STDERR, "need argument after -so ");
02363 exit (1);
02364 }
02365 sprintf (fout, "%s", argv[kar]);
02366
02367 brk = YUP;
02368 }
02369
02370 if (!brk) {
02371 fprintf (SUMA_STDERR,"Error %s: Option %s not understood. Try -help for usage\n", FuncName, argv[kar]);
02372 exit (1);
02373 } else {
02374 brk = NOPE;
02375 kar ++;
02376 }
02377
02378 }
02379
02380 if (LocalHead) fprintf (SUMA_STDERR, "%s: Recursion depth %d, Size %f.\n", FuncName, depth, r);
02381
02382 if (NumOnly) {
02383
02384 int Ntri, Nedge, Nvert;
02385 if (strcmp(bin, "y") == 0) {
02386 Nvert = (int)(pow(2, (2*depth)))*10 + 2;
02387 Ntri = (int)(pow(2, (2*depth)))*20;
02388 Nedge = (int)(pow(2, (2*depth)))*30;
02389 } else {
02390 Nvert = 2 + (10 * depth * depth);
02391 Ntri = 20 * depth * depth;
02392 Nedge = 30 * depth * depth;
02393 }
02394
02395 SUMA_ICOSAHEDRON_DIMENSIONS(r, a, b, lgth);
02396 A = 1/4.0 * lgth * lgth * sqrt(3.0);
02397 V = 5.0 / 12.0 * ( 3 + sqrt(5.0) ) * lgth * lgth * lgth;
02398 if (NumOnly == 1) fprintf (SUMA_STDOUT,"#Nvert\t\tNtri\t\tNedge\t\tArea\t\t\tVolume\n %d\t\t%d\t\t%d\t\t%f\t\t%f\n", Nvert, Ntri, Nedge, A, V);
02399 else fprintf (SUMA_STDOUT," %d\t\t%d\t\t%d\t\t%f\t\t%f\n", Nvert, Ntri, Nedge, A, V);
02400
02401 exit(0);
02402 }
02403
02404 sprintf (surfFileNm, "%s_surf.asc", fout);
02405 sprintf (outSpecFileNm, "%s.spec", fout);
02406
02407 if ( SUMA_filexists(surfFileNm) || SUMA_filexists(outSpecFileNm)) {
02408 fprintf (SUMA_STDERR,"Error %s: At least one of output files %s, %s exists.\nWill not overwrite.\n", \
02409 FuncName, surfFileNm, outSpecFileNm);
02410 exit(1);
02411 }
02412
02413
02414
02415 SO = SUMA_CreateIcosahedron (r, depth, ctr, bin, ToSphere);
02416 if (!SO) {
02417 fprintf (SUMA_STDERR, "Error %s: Failed in SUMA_CreateIcosahedron.\n", FuncName);
02418 exit (1);
02419 }
02420
02421 if (LocalHead) fprintf (SUMA_STDERR, "%s: Now writing surface %s to disk ...\n", FuncName, surfFileNm);
02422
02423
02424
02425 SUMA_writeFSfile (SO, "#tesselated icosahedron for SUMA_CreateIcosahedron (SUMA_SphericalMapping.c)", surfFileNm);
02426
02427
02428 surfaces = (SUMA_SpecSurfInfo *) SUMA_calloc(1, sizeof(SUMA_SpecSurfInfo));
02429
02430 strcpy (surfaces[0].format, "ASCII"); strcpy (surfaces[0].type, "FreeSurfer");
02431 sprintf (surfaces[0].fileToRead, "%s", surfFileNm); strcpy( surfaces[0].mapRef, "SAME");
02432 strcpy (surfaces[0].state, "icosahedron"); strcpy (surfaces[0].dim, "3");
02433
02434 SUMA_writeSpecFile ( surfaces, 1, FuncName, fout, outSpecFileNm );
02435 fprintf (SUMA_STDERR, "\n* To view in SUMA, run:\n suma -spec %s \n\n", outSpecFileNm);
02436
02437
02438 if (LocalHead) fprintf(SUMA_STDERR, "\n... before free surf in createIco\n\n");
02439 SUMA_Free_Surface_Object (SO);
02440 if (LocalHead) fprintf(SUMA_STDERR, "\n... after free surf in createIco\n\n");
02441 SUMA_free(surfaces);
02442
02443 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
02444
02445 SUMA_RETURN(0);
02446
02447 }
02448 #endif
02449
02450
02451
02452 #ifdef SUMA_Map_SurfacetoSurface_STAND_ALONE
02453
02454 void SUMA_Map_StoS_usage ()
02455
02456 {
02457 printf ("\nUsage: SUMA_Map_SurfacetoSurface <-spec spec surf1 surf2> [-prefix fout]\n");
02458 printf ("\n\n\tspec: spec file containing surfaces.\n");
02459 printf ("\n\tsurf1: surface state whose topology (connectivity) will be used.\n");
02460 printf ("\n\tsurf2: surface state whose geometry (shape) will be used.\n\t\t(Spherical input recommended)\n");
02461 printf ("\n\t The topology of surf1 is mapped onto the geometry of surf2.\n\n");
02462 printf ("\n\tfout: prefix for output files. (optional, default StoS)\n\n");
02463 printf ("\n\t Brenna D. Argall LBC/NIMH/NIH brenna.argall@nih.gov \n\t\t\t Fri Sept 20 14:23:42 EST 2002\n\n");
02464 exit (0);
02465 }
02466
02467
02468
02469
02470 int main (int argc, char *argv[])
02471 {
02472
02473 static char FuncName[]={"SUMA_Map_SurfacetoSurface-main"};
02474
02475 char *input=NULL;
02476 char fout[SUMA_MAX_DIR_LENGTH+SUMA_MAX_NAME_LENGTH];
02477 char surfState_1[SUMA_MAX_DIR_LENGTH+SUMA_MAX_NAME_LENGTH];
02478 char surfState_2[SUMA_MAX_DIR_LENGTH+SUMA_MAX_NAME_LENGTH];
02479 SUMA_SurfSpecFile spec;
02480 char surfFileNm[1000], outSpecFileNm[1000];
02481
02482 int kar, i, j;
02483 SUMA_SurfaceObject **surfaces_orig=NULL, *currSurf=NULL;
02484 char *specFile=NULL;
02485 SUMA_MorphInfo *MI=NULL;
02486 float r_temp, ctr[3];
02487 SUMA_SpecSurfInfo *spec_info=NULL;
02488 SUMA_SurfaceObject *morph_SO=NULL;
02489 SUMA_Boolean brk, LocalHead = NOPE, writeFile;
02490
02491
02492
02493 if (LocalHead) fprintf (SUMA_STDERR,"%s: Calling SUMA_Create_CommonFields ...\n", FuncName);
02494
02495 SUMAg_CF = SUMA_Create_CommonFields ();
02496 if (SUMAg_CF == NULL) {
02497 fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_Create_CommonFields\n", FuncName);
02498 exit(1);
02499 }
02500 SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct(SUMA_MAX_DISPLAYABLE_OBJECTS);
02501 if (LocalHead) fprintf (SUMA_STDERR,"%s: SUMA_Create_CommonFields Done.\n", FuncName);
02502
02503
02504 kar = 1;
02505 sprintf( fout, "%s", "StoS");
02506 brk = NOPE;
02507 if (argc < 4) {
02508 SUMA_Map_StoS_usage ();
02509 exit (1);
02510 }
02511 while (kar < argc) {
02512 if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
02513 SUMA_Map_StoS_usage ();
02514 exit (1);
02515 }
02516
02517 if (!brk && (strcmp(argv[kar], "-spec") == 0 ))
02518 {
02519 kar ++;
02520 if (kar >= argc) {
02521 fprintf (SUMA_STDERR, "need argument after -spec ");
02522 exit (1);
02523 }
02524
02525
02526 specFile = argv[kar];
02527
02528
02529 ++kar;
02530 input = argv[kar];
02531 if ( strcmp(input, "-spec")==0 || strcmp(input, "-prefix")==0) {
02532 fprintf(SUMA_STDERR, "\nError %s: Improper format for -spec option. Exiting.\n", FuncName);
02533 exit(1);
02534 }
02535 else strcpy( surfState_1, input);
02536
02537
02538 ++kar;
02539 input = argv[kar];
02540 if ( strcmp(input, "-spec")==0 || strcmp(input, "-prefix")==0) {
02541 fprintf(SUMA_STDERR, "\nError %s: Improper format for -spec option. Exiting.\n", FuncName);
02542 exit(1);
02543 }
02544 else strcpy( surfState_2, input);
02545
02546 brk = YUP;
02547 }
02548
02549 if (!brk && strcmp(argv[kar], "-prefix") == 0)
02550 {
02551 kar ++;
02552 if (kar >= argc) {
02553 fprintf (SUMA_STDERR, "need argument after -prefix ");
02554 exit (1);
02555 }
02556 sprintf (fout, "%s", argv[kar]);
02557
02558 brk = YUP;
02559 }
02560
02561 if (!brk) {
02562 fprintf (SUMA_STDERR,"Error %s: Option %s not understood. Try -help for usage\n", FuncName, argv[kar]);
02563 exit (1);
02564 } else {
02565 brk = NOPE;
02566 kar ++;
02567 }
02568
02569 }
02570
02571
02572 if (specFile == NULL) {
02573 fprintf (SUMA_STDERR,"Error %s: No spec file specified.\n", FuncName);
02574 exit(1);
02575 }
02576
02577
02578 sprintf (surfFileNm, "%s_mappedSurf.asc", fout);
02579 sprintf (outSpecFileNm, "%s.spec", fout);
02580
02581 if ( SUMA_filexists(surfFileNm) || SUMA_filexists(outSpecFileNm) ) {
02582 fprintf (SUMA_STDERR,"Error %s: At least one of output files %s, %s exists.\nWill not overwrite.\n", \
02583 FuncName, surfFileNm, outSpecFileNm);
02584 exit(1);
02585 }
02586
02587
02588
02589 if ( !SUMA_Read_SpecFile (specFile, &spec) ) {
02590 fprintf(SUMA_STDERR,"Error %s: Error in SUMA_Read_SpecFile\n", FuncName);
02591 exit(1);
02592 }
02593
02594
02595
02596 surfaces_orig = (SUMA_SurfaceObject **) SUMA_calloc( 2, sizeof(SUMA_SurfaceObject));
02597 surfaces_orig[0] = NULL; surfaces_orig[1] = NULL;
02598
02599 if ( !SUMA_LoadSpec_eng( &spec, SUMAg_DOv, &SUMAg_N_DOv, NULL , 0, SUMAg_CF->DsetList)) {
02600 fprintf(SUMA_STDERR, "Error %s: Error in SUMA_LoadSpec\n", FuncName);
02601 exit(1);
02602 }
02603
02604 for (i=0; i < SUMAg_N_DOv; ++i) {
02605 if (SUMA_isSO(SUMAg_DOv[i])) {
02606 currSurf = (SUMA_SurfaceObject *)(SUMAg_DOv[i].OP);
02607
02608 if ( SUMA_iswordin(currSurf->State, surfState_1) ==1 ) {
02609 if ( SUMA_iswordin(currSurf->State, "sphere.reg") ==1 ) {
02610 if ( SUMA_iswordin(surfState_1, "sphere.reg") ==1 ) {
02611
02612 surfaces_orig[0] = (SUMA_SurfaceObject *)(SUMAg_DOv[i].OP);
02613 }
02614 }
02615 else surfaces_orig[0] = (SUMA_SurfaceObject *)(SUMAg_DOv[i].OP);
02616 }
02617
02618 if ( SUMA_iswordin(currSurf->State, surfState_2) ==1 ) {
02619 if ( SUMA_iswordin(currSurf->State, "sphere.reg") ==1 ) {
02620 if ( SUMA_iswordin(surfState_2, "sphere.reg") ==1 ) {
02621
02622 surfaces_orig[1] = (SUMA_SurfaceObject *)(SUMAg_DOv[i].OP);
02623 }
02624 }
02625 else surfaces_orig[1] = (SUMA_SurfaceObject *)(SUMAg_DOv[i].OP);
02626 }
02627 }
02628 }
02629
02630 if ( surfaces_orig[0]==NULL || surfaces_orig[1]==NULL) {
02631 fprintf(SUMA_STDERR, "\nError %s: Unable to aquire SO. Exiting.\n (Perhaps your indicated surface state does not exist in the given spec file?)\n\n", FuncName);
02632 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
02633 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
02634 exit(1);
02635 }
02636
02637
02638
02639 if ( !(SUMA_iswordin(surfState_2, "sphere") ==1) )
02640 fprintf(SUMA_STDERR, "\n***\n Warning %s:\n\tIt is recommended that Surface 2 be spherical.\n\tMapping not likely to occur properly or cleanly.\n Proceed at your own risk...\n***\n\n", FuncName);
02641 if (SUMA_iswordin(surfState_1, "smoothwm") ==1 || SUMA_iswordin(surfState_1, "white") ==1 ||
02642 SUMA_iswordin(surfState_2, "smoothwm") ==1 || SUMA_iswordin(surfState_2, "white") ==1 )
02643 fprintf(SUMA_STDERR, "\n***\n Warning %s:\n\tAt least one surface is of the type smoothwm or white.\n\tSuch a surface will likely not map properly or cleanly.\n Assuming you to know what you are doing...\n***\n\n", FuncName);
02644
02645
02646 for (i=0; i<2; ++i) {
02647 if ( surfaces_orig[i]->FileType!=SUMA_FREE_SURFER &&
02648 surfaces_orig[i]->FileType!=SUMA_PLY && surfaces_orig[i]->FileType!=SUMA_VEC ) {
02649 fprintf(SUMA_STDERR, "\n***\n The SurfaceType (of surface %d) is not currently handled\n by this program due to lack of data.\n If you would like this option to be added, please contact\n ziad@nih.gov or brenna.argall@nih.gov.\n***\n\n", i);
02650 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
02651 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
02652 exit(1);
02653 }
02654 }
02655
02656
02657
02658
02659 spec_info = SUMA_calloc(3, sizeof(SUMA_SpecSurfInfo));
02660 if ( spec_info==NULL ) {
02661 fprintf(SUMA_STDERR, "Error %s: Unable to allocate spec_info. Exiting.\n", FuncName);
02662 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
02663 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
02664 exit(1);
02665 }
02666
02667 for (i=0; i<2; ++i) {
02668
02669 if ( surfaces_orig[i]->FileType==SUMA_PLY )
02670 sprintf( spec_info[2*i].type, "Ply");
02671 else if (surfaces_orig[i]->FileType==SUMA_FREE_SURFER)
02672 sprintf( spec_info[2*i].type, "FreeSurfer");
02673 else if (surfaces_orig[i]->FileType==SUMA_VEC)
02674 sprintf( spec_info[2*i].type, "Vec");
02675
02676 if ( surfaces_orig[i]->FileFormat==SUMA_ASCII )
02677 sprintf( spec_info[2*i].format, "ASCII");
02678 else if (surfaces_orig[i]->FileType==SUMA_BINARY ||
02679 surfaces_orig[i]->FileType==SUMA_BINARY_BE ||
02680 surfaces_orig[i]->FileType==SUMA_BINARY_LE)
02681 sprintf( spec_info[2*i].format, "BINARY");
02682 strcpy (spec_info[2*i].dim, "3");
02683 strcpy (spec_info[2*i].mapRef, "SAME");
02684 strcpy (spec_info[2*i].state, surfaces_orig[i]->State);
02685 sprintf (spec_info[2*i].fileToRead, surfaces_orig[i]->Name.FileName);
02686 }
02687
02688
02689 strcpy (spec_info[1].type, spec_info[0].type);
02690 strcpy (spec_info[1].format, spec_info[0].format);
02691 sprintf (spec_info[1].state, "%s_map", spec_info[0].state);
02692 strcpy (spec_info[1].dim, "3");
02693 strcpy (spec_info[1].mapRef, "SAME");
02694 strcpy (spec_info[1].fileToRead, surfFileNm);
02695
02696
02697
02698
02699
02700
02701 if ( !(SUMA_iswordin(surfaces_orig[1]->State, "sphere") ==1)) {
02702
02703
02704 fprintf(SUMA_STDERR, "\n***\n Warning %s:\n\tSurface 2 inflated to a sphere before morphing.\n\n (In many surfaces this will result in skewed connectivity between nodes.)\n (In such cases, expect unsucessful mapping and a possible program exit.)\n\tInput a spherical surface instead!\n***\n\n", FuncName);
02705 r_temp=0;
02706 ctr[0]=0; ctr[1]=0; ctr[2]=0;
02707
02708
02709 for (i=0; i<surfaces_orig[1]->N_Node; ++i) {
02710 j = 3*i;
02711 ctr[0] = ctr[0] + surfaces_orig[1]->NodeList[j];
02712 ctr[1] = ctr[1] + surfaces_orig[1]->NodeList[j+1];
02713 ctr[2] = ctr[2] + surfaces_orig[1]->NodeList[j+2];
02714 }
02715 ctr[0] = ctr[0]/surfaces_orig[1]->N_Node;
02716 ctr[1] = ctr[1]/surfaces_orig[1]->N_Node;
02717 ctr[2] = ctr[2]/surfaces_orig[1]->N_Node;
02718
02719
02720 for (i=0; i<surfaces_orig[1]->N_Node; ++i) {
02721 j = 3*i;
02722 r_temp = sqrt( pow(surfaces_orig[1]->NodeList[j]-ctr[0],2) + pow(surfaces_orig[1]->NodeList[j+1]-ctr[1],2)
02723 + pow(surfaces_orig[1]->NodeList[j+2]-ctr[2],2) );
02724 surfaces_orig[1]->NodeList[j] = (surfaces_orig[1]->NodeList[j] - ctr[0]) / r_temp * 100;
02725 surfaces_orig[1]->NodeList[j+1] = (surfaces_orig[1]->NodeList[j+1] - ctr[1]) / r_temp * 100;
02726 surfaces_orig[1]->NodeList[j+2] = (surfaces_orig[1]->NodeList[j+2] - ctr[2]) / r_temp * 100;
02727 }
02728 }
02729
02730
02731
02732
02733 if ( surfaces_orig[1]->EL==NULL)
02734 SUMA_SurfaceMetrics(surfaces_orig[1], "EdgeList", NULL);
02735 if ( surfaces_orig[1]->EL && surfaces_orig[1]->N_Node)
02736 surfaces_orig[1]->FN = SUMA_Build_FirstNeighb( surfaces_orig[1]->EL, surfaces_orig[1]->N_Node, surfaces_orig[1]->idcode_str);
02737 if ( surfaces_orig[1]->FN==NULL || surfaces_orig[1]->EL==NULL ) {
02738 fprintf(SUMA_STDERR, "Error %s: Failed in acquired Surface Metrics.\n", FuncName);
02739 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
02740 if (surfaces_orig) SUMA_free(surfaces_orig);
02741 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
02742 exit (1);
02743 }
02744
02745 MI = SUMA_MapSurface( surfaces_orig[0], surfaces_orig[1]);
02746 if (!MI) {
02747 fprintf (SUMA_STDERR, "Error %s: Failed in SUMA_MapSurface.\n", FuncName);
02748 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
02749 if (surfaces_orig) SUMA_free(surfaces_orig);
02750 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
02751 exit (1);
02752 }
02753
02754 morph_SO = SUMA_morphToStd( surfaces_orig[1], MI, YUP);
02755 if (!morph_SO) {
02756 fprintf(SUMA_STDERR, "Error %s: Fail in SUMA_morphToStd.\n", FuncName);
02757 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
02758 if (surfaces_orig) SUMA_free(surfaces_orig);
02759 if (MI) SUMA_Free_MorphInfo(MI);
02760 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
02761 exit (1);
02762 }
02763
02764
02765
02766 writeFile = NOPE;
02767 fprintf (SUMA_STDERR, "%s: Now writing surface %s to disk ...\n", FuncName, surfFileNm);
02768 if ( SUMA_iswordin(spec_info[1].type, "FreeSurfer") ==1)
02769 writeFile = SUMA_Save_Surface_Object (surfFileNm, morph_SO, SUMA_FREE_SURFER, SUMA_ASCII, NULL);
02770 else if ( SUMA_iswordin(spec_info[1].type, "Ply") ==1)
02771 writeFile = SUMA_Save_Surface_Object (surfFileNm, morph_SO, SUMA_PLY, SUMA_FF_NOT_SPECIFIED, NULL);
02772 else if ( SUMA_iswordin(spec_info[1].type, "Vec") ==1)
02773 writeFile = SUMA_Save_Surface_Object (surfFileNm, morph_SO, SUMA_VEC, SUMA_ASCII, NULL);
02774 else {
02775 fprintf(SUMA_STDERR, "\n** Surface format (%s) is not currently handled by this program due to lack of data.\n If you would like this option to be added, please contact\n ziad@nih.gov or brenna.argall@nih.gov.\n\n", spec_info[1].type);
02776 exit (0);
02777 }
02778
02779 if ( !writeFile ) {
02780 fprintf (SUMA_STDERR,"Error %s: Failed to write surface object.\n", FuncName);
02781 if (MI) SUMA_Free_MorphInfo (MI);
02782 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
02783 if (morph_SO) SUMA_Free_Surface_Object (morph_SO);
02784 SUMA_free(surfaces_orig);
02785 if (spec_info) SUMA_free(spec_info);
02786 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
02787 exit(1);
02788 }
02789
02790
02791 SUMA_writeSpecFile ( spec_info, 3, FuncName, fout, outSpecFileNm );
02792 fprintf (SUMA_STDERR, "\n**\t\t\t\t\t**\n\t To view in SUMA, run:\n\tsuma -spec %s \n**\t\t\t\t\t**\n\n", outSpecFileNm);
02793
02794
02795
02796 if (MI) SUMA_Free_MorphInfo (MI);
02797 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
02798 if (surfaces_orig) SUMA_free(surfaces_orig);
02799 if (morph_SO) SUMA_Free_Surface_Object (morph_SO);
02800 if (spec_info) SUMA_free(spec_info);
02801
02802 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
02803
02804 SUMA_RETURN(0);
02805
02806 }
02807 #endif
02808
02809
02810
02811
02812 #ifdef SUMA_MapIcosahedron_STAND_ALONE
02813
02814 void SUMA_MapIcosahedron_usage ()
02815
02816 {
02817 static char FuncName[]={"SUMA_MapIcosahedron_usage"};
02818 char * s = NULL;
02819 printf ( "\n"
02820 "Usage: MapIcosahedron <-spec specFile> \n"
02821 " [-rd recDepth] [-ld linDepth] \n"
02822 " [-morph morphSurf] \n"
02823 " [-it numIt] [-prefix fout] \n"
02824 " [-verb] [-help]\n"
02825 "\n"
02826 "Creates new versions of the original-mesh surfaces using the mesh\n"
02827 "of an icosahedron. \n"
02828 "\n"
02829 " -spec specFile: spec file containing original-mesh surfaces\n"
02830 " including the spherical and warped spherical surfaces.\n"
02831 "\n"
02832 " -rd recDepth: recursive (binary) tesselation depth for icosahedron.\n"
02833 " (optional, default:3) See CreateIcosahedron for more info.\n"
02834 "\n"
02835 " -ld linDepth: number of edge divides for linear icosahedron tesselation \n"
02836 " (optional, default uses binary tesselation).\n"
02837 " See CreateIcosahedron -help for more info.\n"
02838 "\n"
02839 " *Note: Enter -1 for recDepth or linDepth to let program \n"
02840 " choose a depth that best approximates the number of nodes in\n"
02841 " original-mesh surfaces.\n"
02842 "\n"
02843 " -morph morphSurf: surface state to which icosahedron is inflated \n"
02844 " accectable inputs are 'sphere.reg' and 'sphere'\n"
02845 " (optional, default uses sphere.reg over sphere).\n"
02846 "\n"
02847 " -it numIt: number of smoothing interations \n"
02848 " (optional, default none).\n"
02849 "\n"
02850 " -prefix fout: prefix for output files.\n"
02851 " (optional, default MapIco)\n"
02852 "\n"
02853 " NOTE: See program SurfQual -help for more info on the following 2 options.\n"
02854 " [-sph_check]: Run tests for checking the spherical surface (sphere.asc)\n"
02855 " The program exits after the checks.\n"
02856 " This option is for debugging FreeSurfer surfaces only.\n"
02857 "\n"
02858 " [-sphreg_check]: Run tests for checking the spherical surface (sphere.reg.asc)\n"
02859 " The program exits after the checks.\n"
02860 " This option is for debugging FreeSurfer surfaces only.\n"
02861 "\n"
02862 " -sph_check and -sphreg_check are mutually exclusive.\n"
02863 "\n"
02864 " -verb: When specified, includes original-mesh surfaces \n"
02865 " and icosahedron in output spec file.\n"
02866 " (optional, default does not include original-mesh surfaces)\n"
02867 "\n"
02868 "NOTE 1: The algorithm used by this program is applicable\n"
02869 " to any surfaces warped to a spherical coordinate\n"
02870 " system. However for the moment, the interface for\n"
02871 " this algorithm only deals with FreeSurfer surfaces.\n"
02872 " This is only due to user demand and available test\n"
02873 " data. If you want to apply this algorithm using surfaces\n"
02874 " created by other programs such as SureFit and Caret, \n"
02875 " Send ziad@nih.gov a note and some test data.\n"
02876 "\n"
02877 "NOTE 2: At times, the standard-mesh surfaces are visibly\n"
02878 " distorted in some locations from the original surfaces.\n"
02879 " So far, this has only occurred when original spherical \n"
02880 " surfaces had topological errors in them. \n"
02881 " See SurfQual -help and SUMA's online documentation \n"
02882 " for more detail.\n"
02883 "\n" );
02884
02885 s = SUMA_New_Additions(0, 1); printf("%s\n", s);SUMA_free(s); s = NULL;
02886
02887 printf ( "\n"
02888 " Brenna D. Argall LBC/NIMH/NIH brenna.argall@nih.gov \n"
02889 " Ziad S. Saad SSC/NIMH/NIH ziad@nih.gov\n"
02890 " Fri Sept 20 2002\n"
02891 "\n");
02892 exit (0);
02893 }
02894
02895
02896
02897
02898 int main (int argc, char *argv[])
02899 {
02900
02901 static char FuncName[]={"MapIcosahedron"};
02902 SUMA_Boolean brk, smooth=NOPE, verb=NOPE;
02903 char fout[SUMA_MAX_DIR_LENGTH+SUMA_MAX_NAME_LENGTH];
02904 char icoFileNm[10000], outSpecFileNm[10000];
02905 char bin[SUMA_MAX_DIR_LENGTH+SUMA_MAX_NAME_LENGTH];
02906 int numTriBin=0, numTriLin=0, numIt=0;
02907
02908 int kar, i, j, k, p, it, id = -1, depth;
02909 char *brainSpecFile=NULL, *OutName = NULL, *morph_surf = NULL;
02910 SUMA_SurfSpecFile brainSpec;
02911
02912 int i_surf, i_morph, mx_N_surf, N_inSpec, N_skip;
02913 float r, ctrX, ctrY, ctrZ, ctr[3];
02914 int *spec_order=NULL, *spec_mapRef=NULL;
02915 SUMA_SpecSurfInfo *spec_info=NULL;
02916 SUMA_SurfaceObject **surfaces_orig=NULL, *icoSurf=NULL, *currSurf=NULL, *currMapRef=NULL;
02917 SUMA_MorphInfo *MI=NULL;
02918 float *smNodeList=NULL, lambda, mu, bpf, *Cx = NULL;
02919 SUMA_INDEXING_ORDER d_order;
02920 SUMA_COMM_STRUCT *cs = NULL;
02921 struct timeval start_time;
02922 float etime_MapSurface;
02923 SUMA_Boolean CheckSphereReg, CheckSphere, skip, writeFile;
02924 SUMA_Boolean LocalHead = NOPE;
02925
02926 FILE *tmpFile=NULL;
02927
02928 SUMA_mainENTRY;
02929
02930
02931 if (LocalHead) fprintf (SUMA_STDERR,"%s: Calling SUMA_Create_CommonFields ...\n", FuncName);
02932
02933 SUMA_STANDALONE_INIT;
02934 #if 0
02935 SUMAg_CF = SUMA_Create_CommonFields ();
02936 if (SUMAg_CF == NULL) {
02937 fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_Create_CommonFields\n", FuncName);
02938 exit(1);
02939 }
02940 #endif
02941
02942 SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct(SUMA_MAX_DISPLAYABLE_OBJECTS);
02943 if (LocalHead) fprintf (SUMA_STDERR,"%s: SUMA_Create_CommonFields Done.\n", FuncName);
02944
02945 cs = SUMA_Create_CommSrtuct();
02946 if (!cs) exit(1);
02947
02948
02949 if (argc < 2) {
02950 SUMA_MapIcosahedron_usage ();
02951 exit (1);
02952 }
02953
02954
02955 depth = 3;
02956 morph_surf = NULL;
02957 sprintf( fout, "%s", "MapIco");
02958 sprintf( bin, "%s", "y");
02959 smooth = NOPE; numIt=0;
02960 verb = NOPE;
02961 kar = 1;
02962 brk = NOPE;
02963 CheckSphere = NOPE;
02964 CheckSphereReg = NOPE;
02965
02966 while (kar < argc) {
02967 if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
02968 SUMA_MapIcosahedron_usage ();
02969 exit (1);
02970 }
02971
02972 SUMA_SKIP_COMMON_OPTIONS(brk, kar);
02973
02974 if (!brk && (strcmp(argv[kar], "-iodbg") == 0)) {
02975 fprintf(SUMA_STDOUT,"Warning %s: SUMA running in in/out debug mode.\n", FuncName);
02976 SUMA_INOUT_NOTIFY_ON;
02977 brk = YUP;
02978 }
02979 if (!brk && (strcmp(argv[kar], "-memdbg") == 0)) {
02980 fprintf(SUMA_STDOUT,"Warning %s: SUMA running in memory trace mode.\n", FuncName);
02981 SUMAg_CF->MemTrace = YUP;
02982 brk = YUP;
02983 }
02984 if (!brk && (strcmp(argv[kar], "-spec") == 0 ))
02985 {
02986 kar ++;
02987 if (kar >= argc) {
02988 fprintf (SUMA_STDERR, "need argument after -spec \n");
02989 exit (1);
02990 }
02991 brainSpecFile = argv[kar];
02992 brk = YUP;
02993 }
02994 if (!brk && (strcmp(argv[kar], "-rd") == 0 ))
02995 {
02996 kar ++;
02997 if (kar >= argc) {
02998 fprintf (SUMA_STDERR, "need argument after -rd \n");
02999 exit (1);
03000 }
03001 depth = atoi(argv[kar]);
03002 sprintf (bin, "y");
03003 brk = YUP;
03004
03005 }
03006 if (!brk && (strcmp(argv[kar], "-ld") == 0 ))
03007 {
03008 kar ++;
03009 if (kar >= argc) {
03010 fprintf (SUMA_STDERR, "need argument after -ld \n");
03011 exit (1);
03012 }
03013 depth = atoi(argv[kar]);
03014 sprintf (bin, "n");
03015 brk = YUP;
03016 }
03017 if (!brk && strcmp(argv[kar], "-morph") == 0)
03018 {
03019 kar ++;
03020 if (kar >= argc) {
03021 fprintf (SUMA_STDERR, "need argument after -morph ");
03022 exit (1);
03023 }
03024 morph_surf = argv[kar];
03025 brk = YUP;
03026 }
03027 if (!brk && (strcmp(argv[kar], "-it") == 0 ))
03028 {
03029 kar ++;
03030 if (kar >= argc) {
03031 fprintf (SUMA_STDERR, "need argument after -it \n");
03032 exit (1);
03033 }
03034 smooth = YUP;
03035 numIt = atoi(argv[kar]);
03036 brk = YUP;
03037 }
03038 if (!brk && (strcmp(argv[kar], "-verb") == 0 ))
03039 {
03040 verb = YUP;
03041 brk = YUP;
03042
03043 }
03044 if (!brk && (strcmp(argv[kar], "-sphreg_check") == 0 ))
03045 {
03046 if (CheckSphere) {
03047 fprintf (SUMA_STDERR, "-sphreg_check & -sph_check are mutually exclusive.\n");
03048 exit (1);
03049 }
03050 CheckSphereReg = YUP;
03051 brk = YUP;
03052
03053 }
03054 if (!brk && (strcmp(argv[kar], "-sph_check") == 0 ))
03055 {
03056 if (CheckSphereReg) {
03057 fprintf (SUMA_STDERR, "-sphreg_check & -sph_check are mutually exclusive.\n");
03058 exit (1);
03059 }
03060 CheckSphere = YUP;
03061 brk = YUP;
03062
03063 }
03064 if (!brk && strcmp(argv[kar], "-prefix") == 0)
03065 {
03066 kar ++;
03067 if (kar >= argc) {
03068 fprintf (SUMA_STDERR, "need argument after -prefix ");
03069 exit (1);
03070 }
03071 sprintf (fout, "%s", argv[kar]);
03072 brk = YUP;
03073 }
03074
03075 if (!brk) {
03076 fprintf (SUMA_STDERR,"Error %s: Option %s not understood. Try -help for usage\n", FuncName, argv[kar]);
03077 exit (1);
03078 }
03079 else {
03080 brk = NOPE;
03081 kar ++;
03082 }
03083
03084 }
03085
03086
03087 if (bin[0] == 'y' && depth > 10) {
03088 fprintf (SUMA_STDERR, "%s: You cannot use a recursive depth > 10.\n", FuncName);
03089 exit(1);
03090 }
03091 if (LocalHead) fprintf (SUMA_STDERR, "%s: %s contains surfaces, tesselation depth is %d.\n", FuncName, brainSpecFile, depth);
03092 if (brainSpecFile == NULL) {
03093 fprintf (SUMA_STDERR,"Error %s: No spec file specified.\n", FuncName);
03094 exit(1);
03095 }
03096
03097
03098 if ( !SUMA_Read_SpecFile (brainSpecFile, &brainSpec)) {
03099 fprintf(SUMA_STDERR,"Error %s: Error in %s SUMA_Read_SpecFile\n", FuncName, brainSpecFile);
03100 exit(1);
03101 }
03102
03103 if ( !SUMA_LoadSpec_eng( &brainSpec, SUMAg_DOv, &SUMAg_N_DOv, NULL, 0 , SUMAg_CF->DsetList) ) {
03104 fprintf(SUMA_STDERR, "Error %s: Error in SUMA_LoadSpec\n", FuncName);
03105 exit(1);
03106 }
03107
03108
03109
03110
03111 if (CheckSphere) {
03112 fprintf(SUMA_STDERR,"%s:\n:Checking sphere surface only.\n", FuncName);
03113 }else if (CheckSphereReg) {
03114 fprintf(SUMA_STDERR,"%s:\n:Checking sphere.reg surface only.\n", FuncName);
03115 }
03116
03117
03118
03119 mx_N_surf = 10;
03120 spec_order = SUMA_calloc( mx_N_surf, sizeof(int));
03121 spec_mapRef = SUMA_calloc( mx_N_surf, sizeof(int));
03122 for (i=0; i<mx_N_surf; ++i) {
03123 spec_order[i] = -1;
03124 spec_mapRef[i] = -1;
03125 }
03126
03127
03128 if ( verb ) N_inSpec = 2*brainSpec.N_Surfs+1;
03129 else N_inSpec = brainSpec.N_Surfs;
03130
03131 spec_info = (SUMA_SpecSurfInfo *)SUMA_calloc( N_inSpec, sizeof(SUMA_SpecSurfInfo));
03132 surfaces_orig = (SUMA_SurfaceObject **) SUMA_calloc( mx_N_surf, sizeof(SUMA_SurfaceObject));
03133 N_skip = 0;
03134 id = -1;
03135 for (i=0; i<brainSpec.N_Surfs; ++i) {
03136
03137 skip = NOPE;
03138
03139 if (verb) i_surf = 2*(i-N_skip);
03140 else i_surf = i-N_skip;
03141
03142 if (SUMA_isSO(SUMAg_DOv[i]))
03143 currSurf = (SUMA_SurfaceObject *)(SUMAg_DOv[i].OP);
03144
03145
03146
03147 if (SUMA_iswordin( currSurf->State, "sphere.reg") ==1 )
03148 id = 4;
03149
03150 else if ( SUMA_iswordin( currSurf->State, "sphere") == 1 &&
03151 SUMA_iswordin( currSurf->State, "sphere.reg") == 0 )
03152 id = 3;
03153
03154 else if ((SUMA_iswordin( currSurf->State, "inflated") ==1) )
03155 id = 2;
03156
03157 else if ((SUMA_iswordin( currSurf->State, "pial") ==1) )
03158 id = 1;
03159
03160 else if ((SUMA_iswordin( currSurf->State, "smoothwm") ==1) )
03161 id = 0;
03162
03163 else if ((SUMA_iswordin( currSurf->State, "white") ==1) )
03164 id = 5;
03165
03166 else if ((SUMA_iswordin( currSurf->State, "occip.patch.3d") ==1) )
03167 id = 6;
03168
03169 else if ((SUMA_iswordin( currSurf->State, "occip.patch.flat") ==1) )
03170 id = 7;
03171
03172 else if ((SUMA_iswordin( currSurf->State, "full.patch.3d") ==1) )
03173 id = 8;
03174
03175 else if ((SUMA_iswordin( currSurf->State, "full.patch.flat") ==1) )
03176 id = 9;
03177 else {
03178
03179 fprintf(SUMA_STDERR, "\nWarning %s: Surface State %s not recognized. Skipping...\n",
03180 FuncName, currSurf->State);
03181 if ( verb ) N_inSpec = N_inSpec-2;
03182 else N_inSpec = N_inSpec-1;
03183 N_skip = N_skip+1;
03184 skip = YUP;
03185 }
03186
03187 if ( ( CheckSphere || CheckSphereReg) && (id != 3 && id !=4) ) skip = YUP;
03188
03189 if ( !skip ) {
03190
03191 if (id < 0) {
03192 SUMA_SL_Err("This cannot be.\n"
03193 "id < 0 !!!\n");
03194 exit(1);
03195 }
03196
03197 spec_order[id] = i-N_skip;
03198 sprintf(spec_info[i_surf].state, "std.%s", currSurf->State );
03199
03200
03201 surfaces_orig[id] = (SUMA_SurfaceObject *)(SUMAg_DOv[i].OP);
03202
03203
03204
03205 if ( (id==4 && CheckSphereReg) || (id==3 && CheckSphere) ){
03206 if (surfaces_orig[id]->EL==NULL)
03207 SUMA_SurfaceMetrics(surfaces_orig[id], "EdgeList", NULL);
03208 if (surfaces_orig[id]->MF==NULL)
03209 SUMA_SurfaceMetrics(surfaces_orig[id], "MemberFace", NULL);
03210 surfaces_orig[id]->Label = SUMA_SurfaceFileName(surfaces_orig[id], NOPE);
03211 OutName = SUMA_append_string (surfaces_orig[id]->Label, "_Conv_detail.1D.dset");
03212 Cx = SUMA_Convexity_Engine ( surfaces_orig[id]->NodeList, surfaces_orig[id]->N_Node,
03213 surfaces_orig[id]->NodeNormList, surfaces_orig[id]->FN, OutName);
03214 if (Cx) SUMA_free(Cx); Cx = NULL;
03215 if (surfaces_orig[id]) {
03216 if (id == 4) SUMA_SphereQuality (surfaces_orig[id], "SphereRegSurf", NULL);
03217 else if (id == 3) SUMA_SphereQuality (surfaces_orig[id], "SphereSurf", NULL);
03218 else {
03219 SUMA_SL_Err("Logic flow error.");
03220 exit(1);
03221 }
03222 }
03223 fprintf(SUMA_STDERR, "%s:\nExiting after SUMA_SphereQuality\n", FuncName);
03224 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
03225 if (surfaces_orig) SUMA_free (surfaces_orig);
03226 if (spec_order) SUMA_free(spec_order);
03227 if (spec_mapRef) SUMA_free(spec_mapRef);
03228 if (spec_info) SUMA_free(spec_info);
03229 if (OutName) SUMA_free(OutName); OutName = NULL;
03230 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
03231 exit(0);
03232 }
03233
03234
03235
03236 if ( surfaces_orig[id]->FileType!=SUMA_FREE_SURFER &&
03237 surfaces_orig[id]->FileType!=SUMA_PLY && surfaces_orig[id]->FileType!=SUMA_VEC ) {
03238 fprintf(SUMA_STDERR, "\n***\n The Surface Type is not currently handled by this program\n due to lack of data.\n If you would like this option to be added, please contact\n ziad@nih.gov or brenna.argall@nih.gov.\n***\n\n");
03239 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
03240 if (surfaces_orig) SUMA_free (surfaces_orig);
03241 if (spec_order) SUMA_free(spec_order);
03242 if (spec_mapRef) SUMA_free(spec_mapRef);
03243 if (spec_info) SUMA_free(spec_info);
03244 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
03245 exit (0);
03246 }
03247 else {
03248 if ( surfaces_orig[id]->FileType==SUMA_PLY )
03249 sprintf(spec_info[i_surf].fileToRead, "%s_%s.ply", fout, spec_info[i_surf].state);
03250 else
03251 sprintf(spec_info[i_surf].fileToRead, "%s_%s.asc", fout, spec_info[i_surf].state);
03252 if (verb) strcpy(spec_info[i_surf+1].fileToRead, surfaces_orig[id]->Name.FileName);
03253 }
03254
03255 if ( SUMA_filexists(spec_info[i_surf].fileToRead) ) {
03256 fprintf (SUMA_STDERR,"Error %s: %s exists. Will not overwrite.\n", FuncName, spec_info[i_surf].fileToRead);
03257 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
03258 if (surfaces_orig) SUMA_free (surfaces_orig);
03259 if (spec_order) SUMA_free(spec_order);
03260 if (spec_mapRef) SUMA_free(spec_mapRef);
03261 if (spec_info) SUMA_free(spec_info);
03262 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
03263 exit(1);
03264 }
03265
03266
03267
03268
03269
03270 currMapRef = (SUMA_SurfaceObject *)SUMAg_DOv[ SUMA_whichDO(surfaces_orig[id]->LocalDomainParentID, SUMAg_DOv,SUMAg_N_DOv) ].OP;
03271 if (SUMA_iswordin( currMapRef->State, "sphere.reg") ==1 )
03272 spec_mapRef[id] = 4;
03273 else if ( SUMA_iswordin( currMapRef->State, "sphere") == 1 &&
03274 SUMA_iswordin( currMapRef->State, "sphere.reg") == 0 )
03275 spec_mapRef[id] = 3;
03276 else if ( SUMA_iswordin( currMapRef->State, "inflated") ==1 )
03277 spec_mapRef[id] = 2;
03278 else if ( SUMA_iswordin( currMapRef->State, "pial") ==1 )
03279 spec_mapRef[id] = 1;
03280 else if ( SUMA_iswordin( currMapRef->State, "smoothwm") ==1 )
03281 spec_mapRef[id] = 0;
03282 else if ( SUMA_iswordin( currMapRef->State, "white") ==1 )
03283 spec_mapRef[id] = 5;
03284 else if ( SUMA_iswordin( currMapRef->State, "occip.patch.3d") ==1 )
03285 spec_mapRef[id] = 6;
03286 else if ( SUMA_iswordin( currMapRef->State, "occip.patch.flat") ==1 )
03287 spec_mapRef[id] = 7;
03288 else if ( SUMA_iswordin( currMapRef->State, "full.patch.3d") ==1 )
03289 spec_mapRef[id] = 8;
03290 else if ( SUMA_iswordin( currMapRef->State, "full.patch.flat") ==1 )
03291 spec_mapRef[id] = 9;
03292 else {
03293
03294 fprintf(SUMA_STDERR, "\nWarning %s: Mapping Reference %s has no recognized surface state in its name.\n\tSetting to default smoothwm.\n\n", FuncName, currMapRef->State);
03295 spec_mapRef[id] = 0;
03296 }
03297
03298
03299 if ( !verb ) k=1;
03300 else k=2;
03301 for ( j=0; j<k; ++j) {
03302 if ( surfaces_orig[id]->FileType==SUMA_PLY )
03303 sprintf( spec_info[i_surf+j].type, "Ply");
03304 else if (surfaces_orig[id]->FileType==SUMA_FREE_SURFER)
03305 sprintf( spec_info[i_surf+j].type, "FreeSurfer");
03306 else if (surfaces_orig[id]->FileType==SUMA_VEC)
03307 sprintf( spec_info[i_surf+j].type, "Vec");
03308 if ( surfaces_orig[id]->FileFormat==SUMA_ASCII )
03309 sprintf( spec_info[i_surf+j].format, "ASCII");
03310 else if (surfaces_orig[id]->FileType==SUMA_BINARY ||
03311 surfaces_orig[id]->FileType==SUMA_BINARY_BE ||
03312 surfaces_orig[id]->FileType==SUMA_BINARY_LE)
03313 sprintf( spec_info[i_surf+j].format, "BINARY");
03314 strcpy (spec_info[i_surf+j].dim, "3");
03315 if (j>0) strcpy(spec_info[i_surf+j].state, currSurf->State);
03316 }
03317 }
03318 }
03319
03320
03321 for (id=0; id<mx_N_surf; ++id) {
03322 if (spec_order[id]>=0) {
03323
03324
03325 if ( verb ) i_surf = 2*spec_order[id];
03326 else i_surf = spec_order[id];
03327
03328 if ( spec_order[spec_mapRef[id]] < 0 ) {
03329
03330 fprintf(SUMA_STDERR, "Warning %s: Mapping Reference for surface %d is not included in the spec file.\n\tSetting to 'SAME'.\n", FuncName, spec_order[id]);
03331 strcpy( spec_info[i_surf].mapRef, "SAME" );
03332 if (verb) strcpy( spec_info[i_surf+1].mapRef, "SAME");
03333 }
03334 else if ( spec_mapRef[id] == id ) {
03335
03336 strcpy( spec_info[i_surf].mapRef, "SAME" );
03337 if (verb) strcpy( spec_info[i_surf+1].mapRef, "SAME");
03338 }
03339 else {
03340
03341 if (verb) {
03342 strcpy( spec_info[i_surf].mapRef, spec_info[2*spec_order[spec_mapRef[id]]].fileToRead );
03343 strcpy( spec_info[i_surf+1].mapRef, spec_info[2*spec_order[spec_mapRef[id]]+1].fileToRead );
03344 }
03345 else strcpy( spec_info[i_surf].mapRef, spec_info[spec_order[spec_mapRef[id]]].fileToRead );
03346 }
03347 }
03348 }
03349
03350
03351 i_morph = -1;
03352 if ( morph_surf!=NULL ) {
03353
03354 if (SUMA_iswordin( morph_surf, "sphere.reg") ==1 )
03355 i_morph = 4;
03356 else if ( SUMA_iswordin( morph_surf, "sphere") == 1 &&
03357 SUMA_iswordin( morph_surf, "sphere.reg") == 0 )
03358 i_morph = 3;
03359 else {
03360 fprintf(SUMA_STDERR, "\nWarning %s: Indicated morphSurf (%s) is not sphere or sphere.reg.\n\tDefault path to determine morphing sphere will be taken.\n", FuncName, morph_surf);
03361 morph_surf = NULL;
03362 }
03363 if ( i_morph!=-1 ) {
03364 if ( spec_order[i_morph]==-1 ) {
03365
03366 fprintf(SUMA_STDERR, "\nWarning %s: Indicated morphSurf (%s) does not exist.\n\tDefault path to determine morphing sphere will be taken.\n", FuncName, morph_surf);
03367 morph_surf = NULL;
03368 }
03369 }
03370 }
03371 if ( morph_surf==NULL) {
03372
03373 if ( spec_order[4]!=-1 ) {
03374
03375 i_morph = 4;
03376 }
03377 else if ( spec_order[3]!=-1 ) {
03378
03379 i_morph = 3;
03380 }
03381 else {
03382
03383 fprintf(SUMA_STDERR, "\nError %s: Neither sphere.reg nor sphere brain states present in Spec file.\nWill not contintue.\n", FuncName);
03384 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
03385 if (surfaces_orig) SUMA_free (surfaces_orig);
03386 if (spec_order) SUMA_free(spec_order);
03387 if (spec_mapRef) SUMA_free(spec_mapRef);
03388 if (spec_info) SUMA_free(spec_info);
03389 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
03390 exit(1);
03391 }
03392 }
03393
03394
03395 if (surfaces_orig[i_morph]->EL==NULL)
03396 SUMA_SurfaceMetrics(surfaces_orig[i_morph], "EdgeList", NULL);
03397 if (surfaces_orig[i_morph]->MF==NULL)
03398 SUMA_SurfaceMetrics(surfaces_orig[i_morph], "MemberFace", NULL);
03399
03400
03401 for (i=0; i<mx_N_surf; ++i) {
03402 if ( spec_order[i]!=-1 && i!=6 && i!=7 && i!=8 && i!=9 &&!(surfaces_orig[i_morph]->N_Node == surfaces_orig[i]->N_Node) ) {
03403 fprintf(SUMA_STDERR, "Error %s: Surfaces (ref [%d], %d!=%d) differ in node number. Exiting.\n", FuncName, i, surfaces_orig[i_morph]->N_Node, surfaces_orig[i]->N_Node);
03404 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
03405 if (surfaces_orig) SUMA_free (surfaces_orig);
03406 if (spec_order) SUMA_free(spec_order);
03407 if (spec_mapRef) SUMA_free(spec_mapRef);
03408 if (spec_info) SUMA_free(spec_info);
03409 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
03410 exit(1);
03411 }
03412 }
03413
03414
03415
03416 if ( depth<0 ) {
03417
03418
03419 i = 0; numTriBin = 20;
03420 while ( numTriBin < surfaces_orig[i_morph]->N_FaceSet ) {
03421 ++i;
03422 numTriBin = 20*( pow(2,2*i) );
03423 }
03424
03425
03426 j = 1; numTriLin = 20;
03427 while ( numTriLin < surfaces_orig[i_morph]->N_FaceSet ) {
03428 ++j;
03429 numTriLin = 20*( pow(j,2) );
03430 }
03431
03432 if ( fabs(numTriLin-surfaces_orig[i_morph]->N_FaceSet) < fabs(numTriBin-surfaces_orig[i_morph]->N_FaceSet) ) {
03433 depth = j;
03434 sprintf (bin, "n");
03435 }
03436 else {
03437 depth = i;
03438 sprintf (bin, "y");
03439 }
03440 }
03441
03442
03443 ctrX=0; ctrY=0; ctrZ=0; j=0;
03444 for (i=0; i<surfaces_orig[i_morph]->N_Node; ++i) {
03445 j = 3*i;
03446 ctrX = ctrX + surfaces_orig[i_morph]->NodeList[j];
03447 ctrY = ctrY + surfaces_orig[i_morph]->NodeList[j+1];
03448 ctrZ = ctrZ + surfaces_orig[i_morph]->NodeList[j+2];
03449 }
03450 ctrX = ctrX/(surfaces_orig[i_morph]->N_Node);
03451 ctrY = ctrY/(surfaces_orig[i_morph]->N_Node);
03452 ctrZ = ctrZ/(surfaces_orig[i_morph]->N_Node);
03453
03454 ctr[0] = 0; ctr[1] = 0; ctr[2] = 0;
03455 r = sqrt( pow( (surfaces_orig[i_morph]->NodeList[0]-ctrX), 2) + pow( (surfaces_orig[i_morph]->NodeList[1]-ctrY), 2)
03456 + pow( (surfaces_orig[i_morph]->NodeList[2]-ctrZ), 2) );
03457
03458
03459
03460 icoSurf = SUMA_CreateIcosahedron (r, depth, ctr, bin, 0);
03461 if (!icoSurf) {
03462 fprintf (SUMA_STDERR, "Error %s: Failed in SUMA_MapIcosahedron.\n", FuncName);
03463 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
03464 if (surfaces_orig) SUMA_free (surfaces_orig);
03465 if (spec_order) SUMA_free(spec_order);
03466 if (spec_mapRef) SUMA_free(spec_mapRef);
03467 if (spec_info) SUMA_free(spec_info);
03468 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
03469 exit (1);
03470 }
03471
03472
03473 if ( verb ) {
03474 sprintf (icoFileNm, "%s_icoSurf.asc", fout);
03475 fprintf (SUMA_STDERR, "\n%s: Now writing surface %s to disk ...\n", FuncName, icoFileNm);
03476 SUMA_writeFSfile (icoSurf, "#icosahedron for SUMA_MapIcosahedron (SUMA_SphericalMapping.c)", icoFileNm);
03477
03478 strcpy (spec_info[ N_inSpec-1 ].type, "FreeSurfer");
03479 strcpy (spec_info[ N_inSpec-1 ].format, "ASCII");
03480 strcpy (spec_info[ N_inSpec-1 ].mapRef, "SAME");
03481 strcpy (spec_info[ N_inSpec-1 ].dim, "3");
03482 strcpy (spec_info[ N_inSpec-1 ].state, "icosahedron");
03483 strcpy (spec_info[ N_inSpec-1 ].fileToRead, icoFileNm);
03484 }
03485
03486
03487
03488
03489
03490 SUMA_etime(&start_time,0);
03491
03492 MI = SUMA_MapSurface( icoSurf, surfaces_orig[i_morph] ) ;
03493 if (!MI) {
03494 fprintf (SUMA_STDERR, "Error %s: Failed in SUMA_MapIcosahedron.\n", FuncName);
03495 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
03496 if (surfaces_orig) SUMA_free (surfaces_orig);
03497 if (icoSurf) SUMA_Free_Surface_Object(icoSurf);
03498 if (spec_order) SUMA_free(spec_order);
03499 if (spec_mapRef) SUMA_free(spec_mapRef);
03500 if (spec_info) SUMA_free(spec_info);
03501 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
03502 exit (1);
03503 }
03504
03505 etime_MapSurface = SUMA_etime(&start_time,1);
03506
03507
03508
03509
03510
03511 for (id=0; id<mx_N_surf; ++id) {
03512
03513 if ( spec_order[id] != -1 ) {
03514
03515
03516
03517 if ( surfaces_orig[id]->EL==NULL)
03518 SUMA_SurfaceMetrics(surfaces_orig[id], "EdgeList", NULL);
03519 if ( surfaces_orig[id]->EL && surfaces_orig[id]->N_Node)
03520 surfaces_orig[id]->FN = SUMA_Build_FirstNeighb( surfaces_orig[id]->EL, surfaces_orig[id]->N_Node, surfaces_orig[id]->idcode_str);
03521 if ( surfaces_orig[id]->FN==NULL || surfaces_orig[id]->EL==NULL ) {
03522 fprintf(SUMA_STDERR, "Error %s: Failed in acquired Surface Metrics.\n", FuncName);
03523 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
03524 if (surfaces_orig) SUMA_free (surfaces_orig);
03525 if (icoSurf) SUMA_Free_Surface_Object(icoSurf);
03526 if (currSurf) SUMA_free (currSurf);
03527 if (spec_order) SUMA_free(spec_order);
03528 if (spec_mapRef) SUMA_free(spec_mapRef);
03529 if (spec_info) SUMA_free(spec_info);
03530 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
03531 exit (1);
03532 }
03533
03534 currSurf = SUMA_morphToStd( surfaces_orig[id], MI, YUP);
03535 if ( !currSurf ) {
03536 fprintf(SUMA_STDERR, "Error %s: Failed in morphing surface object.\n", FuncName);
03537 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
03538 if (surfaces_orig) SUMA_free (surfaces_orig);
03539 if (icoSurf) SUMA_Free_Surface_Object(icoSurf);
03540 if (currSurf) SUMA_free (currSurf);
03541 if (spec_order) SUMA_free(spec_order);
03542 if (spec_mapRef) SUMA_free(spec_mapRef);
03543 if (spec_info) SUMA_free(spec_info);
03544 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
03545 exit (1);
03546 }
03547 currSurf->FileType = surfaces_orig[id]->FileType;
03548
03549
03550
03551 if ( smooth && ( id==0 || id==1 || id==5 ) ) {
03552
03553 bpf = 0.1;
03554 if ( !SUMA_Taubin_Smooth_Coef (bpf, &lambda, &mu) )
03555 fprintf(SUMA_STDERR, "Error %s: Failed in acquiring Taubin Coefficients. Surface will not be smoothed.\n\n", FuncName);
03556
03557 else {
03558 d_order = SUMA_ROW_MAJOR;
03559 currSurf->FN = icoSurf->FN;
03560
03561 smNodeList = SUMA_Taubin_Smooth (currSurf, NULL, lambda, mu, currSurf->NodeList, 2*numIt, 3, d_order, NULL, cs, NULL);
03562 if ( !smNodeList )
03563 fprintf(SUMA_STDERR, "Error %s: Failed in Taubin Smoothing. Surface will not be smoothed.\n\n", FuncName);
03564 else {
03565 SUMA_free( currSurf->NodeList);
03566 currSurf->NodeList = smNodeList;
03567 }
03568 }
03569 }
03570
03571
03572 if ( verb ) i_surf = 2*spec_order[id];
03573 else i_surf = spec_order[id];
03574
03575 fprintf (SUMA_STDERR, "%s: Now writing surface %s to disk ...\n", FuncName, spec_info[i_surf].fileToRead);
03576 writeFile = NOPE;
03577 if ( SUMA_iswordin(spec_info[i_surf].type, "FreeSurfer") ==1)
03578 writeFile = SUMA_Save_Surface_Object (spec_info[i_surf].fileToRead, currSurf, SUMA_FREE_SURFER, SUMA_ASCII, NULL);
03579 else if ( SUMA_iswordin(spec_info[i_surf].type, "Ply") ==1)
03580 writeFile = SUMA_Save_Surface_Object (spec_info[i_surf].fileToRead, currSurf, SUMA_PLY, SUMA_FF_NOT_SPECIFIED, NULL);
03581 else if ( SUMA_iswordin(spec_info[i_surf].type, "Vec") ==1)
03582 writeFile = SUMA_Save_Surface_Object (spec_info[i_surf].fileToRead, currSurf, SUMA_VEC, SUMA_ASCII, NULL);
03583 else {
03584 fprintf(SUMA_STDERR, "\n** Surface format (%s) is not currently handled by this program due to lack of data.\n If you would like this option to be added, please contact\n ziad@nih.gov or brenna.argall@nih.gov.\n\n", spec_info[i_surf].type);
03585 exit (0);
03586 }
03587
03588 if ( !writeFile ) {
03589 fprintf (SUMA_STDERR,"Error %s: Failed to write surface object.\n", FuncName);
03590 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
03591 if (surfaces_orig) SUMA_free (surfaces_orig);
03592 if (icoSurf) SUMA_Free_Surface_Object(icoSurf);
03593 if (currSurf) SUMA_free (currSurf);
03594 if (spec_order) SUMA_free(spec_order);
03595 if (spec_mapRef) SUMA_free(spec_mapRef);
03596 if (spec_info) SUMA_free(spec_info);
03597 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
03598 exit (1);
03599 }
03600
03601 if ( currSurf->FN ) currSurf->FN=NULL;
03602 if ( currSurf ) SUMA_Free_Surface_Object( currSurf );
03603 currSurf = NULL;
03604 }
03605 }
03606
03607
03608
03609 sprintf (outSpecFileNm, "%s_std.spec", fout);
03610 SUMA_writeSpecFile ( spec_info, N_inSpec, FuncName, fout, outSpecFileNm );
03611
03612
03613 fprintf (SUMA_STDERR, "\nSUMA_MapSurface took %f seconds to execute.\n", etime_MapSurface);
03614 fprintf (SUMA_STDERR, "\n**\t\t\t\t\t**\n\t To view in SUMA, run:\n\tsuma -spec %s \n**\t\t\t\t\t**\n\n", outSpecFileNm);
03615
03616
03617
03618 if (spec_order) SUMA_free(spec_order);
03619 if (spec_mapRef) SUMA_free(spec_mapRef);
03620 if (spec_info) SUMA_free(spec_info);
03621 if (MI) SUMA_Free_MorphInfo (MI);
03622
03623
03624 if (icoSurf) SUMA_Free_Surface_Object (icoSurf);
03625 if (currSurf) { SUMA_Free_Surface_Object(currSurf);}
03626 if (SUMAg_DOv) SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv);
03627 if (surfaces_orig) SUMA_free (surfaces_orig);
03628
03629 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
03630
03631 SUMA_RETURN(0);
03632
03633 }
03634 #endif
03635
03636
03637
03638
03639
03640
03641
03642
03643
03644
03645
03646
03647
03648
03649
03650 void SUMA_read1D (char* fileNm, int* i_colm, int* i_locInfo, SUMA_1dData* data) {
03651
03652 FILE *file=NULL;
03653 char *line=NULL, *frag=NULL;
03654 char scan[100];
03655 int i=0, j=0, k=0, num_node=0, lgth=0, i_curr=0, i_last=0;
03656 int num_loc=0, num_val=0, num_tot=0, valCnt=0, tempInt=0;
03657 int *i_colm_ndx=NULL, *i_colmSrtd=NULL, *i_cat=NULL;
03658 float *valArray=NULL, tempFlt=0;
03659 int *ndx_list=NULL, *vxl_list=NULL, *ijk_list=NULL, *nvox_list=NULL;
03660 SUMA_Boolean nd_given=NOPE;
03661 static char FuncName[]={"SUMA_read1D"};
03662
03663 SUMA_ENTRY;
03664
03665
03666 lgth = 500000;
03667
03668
03669 if (i_colm[0] == 0) {
03670 fprintf(SUMA_STDERR, "\nError %s: No column indices given! Exiting.\n", FuncName);
03671 exit(1);
03672 }
03673 else num_tot = i_colm[0];
03674 num_loc=0;
03675
03676 for (i=0; i<6; ++i) {
03677 for (j=0; j<num_tot-1; ++j) {
03678 if ( i_locInfo[i]==i_colm[j+1] ) {
03679
03680 ++num_loc;
03681 }
03682 }
03683 }
03684
03685 num_val = num_tot-num_loc;
03686
03687
03688 i_colmSrtd = SUMA_calloc( i_colm[0]-1, sizeof(int));
03689 for (i=0; i<i_colm[0]; ++i) {
03690
03691 i_colmSrtd[i] = i_colm[i+1];
03692 }
03693 i_colm_ndx = SUMA_z_dqsort( i_colmSrtd, num_tot );
03694
03695
03696 i_cat = SUMA_calloc( num_tot, sizeof(int));
03697 for ( i=0; i<num_tot; ++i) {
03698 if ( i_colmSrtd[i]==i_locInfo[0] ) {
03699 i_cat[i] = 0;
03700 nd_given = YUP;
03701 }
03702 else if ( i_colmSrtd[i]==i_locInfo[1] ) i_cat[i] = 1;
03703 else if ( i_colmSrtd[i]==i_locInfo[2] ) i_cat[i] = 2;
03704 else if ( i_colmSrtd[i]==i_locInfo[3] ) i_cat[i] = 3;
03705 else if ( i_colmSrtd[i]==i_locInfo[4] ) i_cat[i] = 4;
03706 else if ( i_colmSrtd[i]==i_locInfo[5] ) i_cat[i] = 5;
03707 else i_cat[i] = -1;
03708 }
03709
03710 valArray = SUMA_calloc( num_val*lgth, sizeof(float) );
03711 ndx_list = SUMA_calloc( lgth, sizeof(int) );
03712 vxl_list = SUMA_calloc( lgth, sizeof(int) );
03713 ijk_list = SUMA_calloc( 3*lgth, sizeof(int) );
03714 nvox_list = SUMA_calloc( lgth, sizeof(int) );
03715
03716 line = SUMA_calloc( 10000, sizeof(char));
03717
03718 if ( !valArray || !ndx_list || !vxl_list || !ijk_list || !nvox_list || !line) {
03719 fprintf(SUMA_STDERR, "Error %s: Failed in allocation.\n", FuncName);
03720 if (valArray) SUMA_free(valArray);
03721 if (ndx_list) SUMA_free(ndx_list);
03722 if (vxl_list) SUMA_free(vxl_list);
03723 if (ijk_list) SUMA_free(ijk_list);
03724 if (nvox_list) SUMA_free(nvox_list);
03725 if (line) SUMA_free(line);
03726 exit(1);
03727 }
03728
03729 if( (file = fopen(fileNm, "r"))==NULL) {
03730 fprintf (SUMA_STDERR, "Failed in opening %s for reading.\n", fileNm);
03731 if (valArray) SUMA_free(valArray);
03732 if (ndx_list) SUMA_free(ndx_list);
03733 if (vxl_list) SUMA_free(vxl_list);
03734 if (ijk_list) SUMA_free(ijk_list);
03735 if (nvox_list) SUMA_free(nvox_list);
03736 if (line) SUMA_free(line);
03737 exit(1);
03738 }
03739
03740 else {
03741
03742
03743 fgets( line, 1000, file);
03744 while( line[0]=='#' ) {
03745 fgets( line, 10000, file);
03746 }
03747
03748
03749 num_node = 0;
03750 while( !feof(file) && line[0]!='#' && num_node<lgth) {
03751 valCnt = 0;
03752 i_last = 0;
03753 frag = strtok(line, " \t\n\r");
03754 if (frag==NULL) {
03755 fprintf(SUMA_STDERR, "\nError %s: Indicated column for file not found. Exiting.\n", FuncName);
03756 exit(1);
03757 }
03758
03759 for ( k=0; k<num_tot; ++k ) {
03760 for ( i=0; i<i_colmSrtd[k]-i_last; ++i) {
03761 frag = strtok(NULL, " \t\n\r");
03762 if (frag==NULL) {
03763 fprintf(SUMA_STDERR, "\nError %s: Indicated column for file not found. Exiting.\n", FuncName);
03764 exit(1);
03765 }
03766 }
03767
03768 if (frag==NULL) {
03769 fprintf(SUMA_STDERR, "\nError %s: Indicated column for file not found. Exiting.\n", FuncName);
03770 exit(1);
03771 }
03772
03773 if ( i_cat[k]!=-1 ) {
03774
03775 sscanf(frag, "%d", &tempInt);
03776 }
03777 else {
03778
03779 sscanf(frag, "%f", &tempFlt);
03780 }
03781
03782 if ( i_cat[k]==0 ) ndx_list[num_node] = tempInt;
03783 else if ( i_cat[k]==1 ) vxl_list[num_node] = tempInt;
03784 else if ( i_cat[k]==2 ) ijk_list[3*num_node] = tempInt;
03785 else if ( i_cat[k]==3 ) ijk_list[3*num_node+1] = tempInt;
03786 else if ( i_cat[k]==4 ) ijk_list[3*num_node+2] = tempInt;
03787 else if ( i_cat[k]==5 ) nvox_list[num_node] = tempInt;
03788 else valArray[ (valCnt++)*lgth + num_node ] = tempFlt;
03789 i_last = i_colmSrtd[k];
03790 }
03791 fgets( line, 10000, file);
03792 ++num_node;
03793 }
03794 fclose(file);
03795 }
03796
03797
03798
03799 data->N_elem = num_node;
03800 data->nd_list = SUMA_calloc( num_node, sizeof(int));
03801 data->vxl_list = SUMA_calloc( num_node, sizeof(int));
03802 data->ijk_list = SUMA_calloc( 3*num_node, sizeof(int));
03803 data->nvox_list = SUMA_calloc( num_node, sizeof(int));
03804 data->valArray = SUMA_calloc( num_val*num_node, sizeof(float));
03805
03806 for (i=0; i<num_node; ++i) {
03807 if (nd_given) data->nd_list[i] = ndx_list[i];
03808 else data->nd_list[i] = i;
03809 data->vxl_list[i] = vxl_list[i];
03810 data->ijk_list[i] = ijk_list[i];
03811 data->nvox_list[i] = nvox_list[i];
03812 for (k=0; k<num_val; ++k) {
03813 data->valArray[ k*num_node +i ] = valArray[ k*lgth +i ];
03814 }
03815 }
03816
03817 SUMA_free(line);
03818 SUMA_free(i_colm_ndx);
03819 SUMA_free(i_colmSrtd);
03820 SUMA_free(i_cat);
03821
03822 SUMA_free(valArray);
03823 SUMA_free(ndx_list);
03824 SUMA_free(vxl_list);
03825 SUMA_free(ijk_list);
03826 SUMA_free(nvox_list);
03827
03828 SUMA_RETURNe;
03829 }
03830
03831
03832
03833
03834
03835
03836
03837
03838
03839
03840
03841
03842
03843
03844
03845 void SUMA_write1D ( int *num, float *vals, int *index, char firstline[], char outFileNm[]) {
03846
03847 FILE *outFile=NULL;
03848 int i=0, j=0, k=0;
03849 static char FuncName[]={"SUMA_write1D"};
03850
03851 SUMA_ENTRY;
03852
03853 outFile = fopen(outFileNm, "w");
03854 if (!outFile) {
03855 fprintf (SUMA_STDERR, "Failed in opening %s for writing.\n", outFileNm);
03856 exit (1);
03857 }
03858 else {
03859 if (firstline!=NULL) fprintf (outFile, "%s\n", firstline);
03860 for (i=0; i<num[0]; ++i) {
03861 if ( index!=NULL ) {
03862
03863 j = index[i] * num[1];
03864 fprintf( outFile, "%10d ", index[i]);
03865 }
03866 else j = i*num[1];
03867 for ( k=0; k<num[1]; ++k ) {
03868
03869 fprintf( outFile, "%10f ", vals[ j+k ]);
03870 }
03871 fprintf( outFile, "\n");
03872 }
03873 fclose(outFile);
03874 }
03875 SUMA_RETURNe;
03876 }
03877
03878
03879
03880
03881
03882
03883
03884
03885
03886
03887
03888
03889
03890
03891 float * SUMA_createColGradient( float *col, int numSeg, SUMA_Boolean allGvn ) {
03892
03893 int i, j, k, it;
03894 int numCol=0, numStdIncr, numColGvn;
03895 int *colRngInt=NULL, *colUsed=NULL, i_incr=0;
03896 int *bind_currCol=NULL, *distTOint=NULL, colIntArray[8];
03897 int dist_intTOrngBeg, dist_intTOrngEnd, *numColDiv=NULL, *tmpInt=NULL;
03898 float *colRng=NULL, color[8][3], *colSeg=NULL, *colIncr=NULL, *stdColIncr=NULL;
03899 float incR=0.0, incG=0.0, incB=0.0, temp, dist[2];
03900 SUMA_Boolean decr = NOPE, noGrad = NOPE;
03901 static char FuncName[]={"SUMA_createColGradient"};
03902
03903 SUMA_ENTRY;
03904 fprintf(SUMA_STDERR, "numSeg = %d\n", numSeg);
03905 if ( (col==NULL || numSeg<0) && allGvn) {
03906
03907
03908 allGvn = NOPE;
03909 }
03910
03911
03912
03913 if (col==NULL) {
03914
03915 colRngInt = SUMA_calloc(2, sizeof(int));
03916 colRng = SUMA_calloc(2, sizeof(int));
03917 colRng[0] = 0;
03918 colRng[1] = 7;
03919
03920
03921
03922
03923 }
03924 else {
03925 if ( col[0]<0 || col[1]>7 || col[0]<0 || col[1]>7) {
03926 fprintf(SUMA_STDERR, "\nError %s: Color Ranges not within 0->7. Exiting.\n", FuncName);
03927 exit(1);
03928 }
03929
03930
03931 if ( allGvn ) {
03932
03933 colRng = SUMA_calloc(numSeg, sizeof(float));
03934 for (i=0; i<numSeg; ++i) {
03935 colRng[i] = col[i]; }
03936 }
03937 else {
03938
03939 colRng = SUMA_calloc(2, sizeof(float));
03940 colRng[0] = col[0];
03941 colRng[1] = col[1];
03942
03943
03944 if ( col[1] < col[0] ) {
03945 decr = YUP;
03946 }
03947 }
03948 }
03949
03950
03951 color[0][0]=.75; color[0][1]=0; color[0][2]=0;
03952 color[1][0]=1; color[1][1]=.5; color[1][2]=0;
03953 color[2][0]=1; color[2][1]=1; color[2][2]=0;
03954 color[3][0]=0; color[3][1]=.75; color[3][2]=0;
03955 color[4][0]=0; color[4][1]=.75; color[4][2]=.75;
03956 color[5][0]=0; color[5][1]=0; color[5][2]=.75;
03957 color[6][0]=.5; color[6][1]=0; color[6][2]=.75;
03958 color[7][0]=.75; color[7][1]=0; color[7][2]=0;
03959
03960
03961 if (numSeg<0) numSeg = 700;
03962
03963
03964
03965 stdColIncr = SUMA_calloc( 71, sizeof(float));
03966 stdColIncr[0] = 0.0;
03967 colIntArray[0] = 0;
03968 for ( i=0; i<7; ++i) {
03969 colIntArray[i+1] = i+1;
03970 for ( j=1; j<11; ++j) {
03971 stdColIncr[ 10*i+j ] = stdColIncr[ 10*i ] + j*0.1;
03972 }
03973
03974 }
03975
03976
03977
03978
03979 if (allGvn) numColGvn = numSeg;
03980 else numColGvn = 2;
03981 if (!colRngInt)
03982 colRngInt = SUMA_calloc( 2*numColGvn, sizeof(int));
03983 bind_currCol = SUMA_calloc( 2, sizeof(int));
03984 distTOint = SUMA_calloc( numColGvn, sizeof(float));
03985
03986 for (i=0; i<numColGvn; ++i) {
03987
03988
03989 bind_currCol[0] = 0; bind_currCol[1] = 70;
03990 if ( !SUMA_binSearch( stdColIncr, colRng[i], bind_currCol ) ) {
03991 fprintf(SUMA_STDERR, "\nError %s: Failed in binary search !(%f < %f < %f). Exiting.\n\n", FuncName, stdColIncr[bind_currCol[0]], colRng[i], stdColIncr[bind_currCol[1]]);
03992 if (colRng) SUMA_free(colRng);
03993 if (stdColIncr) SUMA_free(stdColIncr);
03994 if (colRngInt) SUMA_free(colRngInt);
03995 if (bind_currCol) SUMA_free(bind_currCol);
03996 if (distTOint) SUMA_free(distTOint);
03997 exit(1);
03998 }
03999 if ( abs( stdColIncr[bind_currCol[0]]-colRng[i] ) < abs( stdColIncr[bind_currCol[1]]-colRng[i] ))
04000 colRng[i] = stdColIncr[bind_currCol[0]];
04001 else colRng[i] = stdColIncr[bind_currCol[1]];
04002
04003
04004
04005
04006 if ( abs( colRng[i] - (int)colRng[i] )<0.01 ) {
04007
04008 colRngInt[ 2*i ] = (int)colRng[i];
04009 colRngInt[ 2*i+1 ] = (int)colRng[i];
04010 }
04011 else {
04012 colRngInt[ 2*i ] = (int)colRng[i];
04013 colRngInt[ 2*i+1 ] = (int)colRng[i]+1;
04014 }
04015
04016
04017 distTOint[i] = bind_currCol[0]%10;
04018
04019 fprintf(SUMA_STDERR, "%d: %d < %f < %d dist = %d\n", i, colRngInt[2*i], colRng[i], colRngInt[2*i+1], distTOint[i]);
04020 }
04021
04022
04023
04024
04025 if ( !allGvn ) {
04026
04027
04028 if ( !decr ) numCol = colRngInt[3]-colRngInt[0] + 1;
04029 else numCol = colRngInt[1] - colRngInt[2] + 1;
04030
04031
04032 colUsed = SUMA_calloc(numCol, sizeof(int));
04033 colUsed[0] = colRngInt[0];
04034 for (i=1; i<numCol; ++i) {
04035 if ( !decr ) colUsed[i] = colUsed[i-1] + 1;
04036 else colUsed[i] = colUsed[i-1] - 1;
04037 }
04038
04039
04040 if ( !decr ) {
04041 dist_intTOrngBeg = distTOint[0];
04042 dist_intTOrngEnd = 10-distTOint[1];
04043 }
04044 else {
04045 dist_intTOrngBeg = 10-distTOint[0];
04046 dist_intTOrngEnd = distTOint[1];
04047 }
04048
04049
04050
04051
04052 if (numCol>3)
04053
04054 numStdIncr = 10*(numCol-1) - dist_intTOrngBeg - dist_intTOrngEnd;
04055 else if ( numCol==3 ) {
04056
04057 numStdIncr = 20 - dist_intTOrngBeg - dist_intTOrngEnd;
04058 }
04059 else {
04060
04061 numStdIncr = 10 - dist_intTOrngBeg - dist_intTOrngEnd;
04062 noGrad = YUP;
04063 }
04064
04065 if ( noGrad ) {
04066
04067
04068
04069 for (i=0; i<numSeg; ++i ) {
04070 k = 3*i;
04071 colSeg[k] = color[colRngInt[0]][0];
04072 colSeg[k+1] = color[colRngInt[0]][1];
04073 colSeg[k+2] = color[colRngInt[0]][2];
04074 }
04075 }
04076
04077 else {
04078
04079
04080
04081 numColDiv = SUMA_calloc( numCol-1, sizeof(int));
04082 numColDiv[0] = 10 - dist_intTOrngBeg;
04083 for (i=1; i<numCol-2; ++i) { numColDiv[i] = 10; }
04084 numColDiv[numCol-2] = 10 - dist_intTOrngEnd;
04085
04086
04087
04088
04089 colIncr = SUMA_calloc( 3*(numSeg-1)*(numStdIncr), sizeof(float));
04090
04091 k=0;
04092 for (i=0; i<numCol-1; ++i) {
04093
04094
04095 incR = (color[colUsed[i+1]][0] - color[colUsed[i]][0])/(10*(numSeg-1));
04096 incG = (color[colUsed[i+1]][1] - color[colUsed[i]][1])/(10*(numSeg-1));
04097 incB = (color[colUsed[i+1]][2] - color[colUsed[i]][2])/(10*(numSeg-1));
04098
04099
04100 for (j=0; j<numColDiv[i]*(numSeg-1); ++j) {
04101 colIncr[k++] = incR;
04102 colIncr[k++] = incG;
04103 colIncr[k++] = incB;
04104 }
04105 }
04106 }
04107
04108
04109
04110
04111 colSeg = SUMA_calloc( 3*numSeg, sizeof(float));
04112
04113 colSeg[0] = color[colUsed[0]][0];
04114 colSeg[1] = color[colUsed[0]][1];
04115 colSeg[2] = color[colUsed[0]][2];
04116 for (i=0; i<dist_intTOrngBeg; ++i) {
04117
04118
04119
04120 colSeg[0] = colSeg[0] + colIncr[0];
04121 colSeg[1] = colSeg[1] + colIncr[1];
04122 colSeg[2] = colSeg[2] + colIncr[2];
04123 }
04124
04125 i_incr = 0;
04126 fprintf(SUMA_STDERR, "numSeg = %d\n", numSeg);
04127 for (i=1; i<numSeg; ++i ) {
04128 k = 3*i;
04129
04130
04131 colSeg[k] = colSeg[k-3];
04132 colSeg[k+1] = colSeg[k-2];
04133 colSeg[k+2] = colSeg[k-1];
04134
04135
04136 for (j=0; j<numStdIncr; ++j) {
04137 colSeg[k] = colSeg[k] + colIncr[3*i_incr];
04138 colSeg[k+1] = colSeg[k+1] + colIncr[3*i_incr+1];
04139 colSeg[k+2] = colSeg[k+2] + colIncr[3*i_incr+2];
04140 ++i_incr;
04141 }
04142 }
04143 }
04144
04145 else {
04146
04147
04148 colSeg = SUMA_calloc( 3*numSeg, sizeof(float));
04149
04150 for (i=0; i<numSeg; ++i) {
04151
04152 k = 3*i;
04153 colSeg[k] = ((10-distTOint[i])/10)*color[colRngInt[2*i]][0] + (distTOint[i]/10)*color[colRngInt[2*i+1]][0];
04154 colSeg[k+1] = ((10-distTOint[i])/10)*color[colRngInt[2*i]][1] + (distTOint[i]/10)*color[colRngInt[2*i+1]][1];
04155 colSeg[k+2] = ((10-distTOint[i])/10)*color[colRngInt[2*i]][2] + (distTOint[i]/10)*color[colRngInt[2*i+1]][2];
04156 }
04157 }
04158
04159 SUMA_free(stdColIncr);
04160 SUMA_free(bind_currCol);
04161 SUMA_free(colRngInt);
04162 SUMA_free(distTOint);
04163 if (!allGvn) {
04164 SUMA_free(colUsed);
04165 SUMA_free(numColDiv);
04166 SUMA_free(colIncr);
04167 }
04168
04169 return colSeg;
04170 }
04171
04172
04173
04174
04175
04176
04177
04178
04179
04180
04181
04182
04183
04184
04185
04186
04187
04188 float * SUMA_assignColors( float *vals, float *cols, int numVal, int numCol, float *gradRng, float *valDiv ) {
04189
04190 int i, j, k;
04191 float *valCol=NULL;
04192 float min, max, segSize=0;
04193 static char FuncName[]={"SUMA_assignColors"};
04194
04195 SUMA_ENTRY;
04196 valCol = SUMA_calloc( 3*numVal, sizeof(float));
04197 valDiv = SUMA_calloc( numCol, sizeof(float));
04198
04199
04200 min = vals[0]; max = vals[0];
04201 for (i=0; i<numVal; ++i) {
04202 if (vals[i]<min) min = vals[i];
04203 else if (vals[i]>max) max = vals[i];
04204 }
04205
04206 if (gradRng==NULL) {
04207
04208 segSize = (max-min)/numCol;
04209
04210 for (i=0; i<numCol; ++i) {
04211 valDiv[i] = min+(i+1)*segSize;
04212 }
04213 }
04214 else {
04215
04216 segSize = (gradRng[1]-gradRng[0])/(numCol-2);
04217
04218 valDiv[0] = gradRng[0];
04219 valDiv[numCol-1] = max;
04220 for (i=1; i<numCol-1; ++i) {
04221 valDiv[i] = valDiv[0] + i*segSize;
04222 }
04223 }
04224
04225 for (i=0; i<numVal; ++i) {
04226
04227 j = 3*i;
04228 for (k=0; k<numCol; ++k) {
04229 if ( vals[i]<=valDiv[k] ) {
04230
04231 valCol[j] = cols[ 3*k ];
04232 valCol[j+1] = cols[ 3*k+1 ];
04233 valCol[j+2] = cols[ 3*k+2 ];
04234 break;
04235 }
04236 }
04237 }
04238 fprintf(SUMA_STDERR, "numCol = %d\n", numCol);
04239
04240 if (numCol<20) {
04241 fprintf(SUMA_STDERR, "COLOR RANGES:\n\t[%f, %f]\n", min, valDiv[0]);
04242 for (i=1; i<numCol; ++i) {
04243 fprintf(SUMA_STDERR, "\t(%f, %f]\n", valDiv[i-1], valDiv[i]);
04244 }
04245 fprintf(SUMA_STDERR, "\n");
04246 }
04247
04248 SUMA_free(valDiv);
04249
04250 return valCol;
04251 }
04252
04253
04254
04255
04256 SUMA_1dData *SUMA_Create_1dData (void)
04257 {
04258 static char FuncName[]={"SUMA_Create_1dData"};
04259 int i=0;
04260
04261 SUMA_1dData *data = NULL;
04262
04263 SUMA_ENTRY;
04264
04265 data = (SUMA_1dData *) SUMA_malloc (sizeof(SUMA_1dData));
04266 if (!data) {
04267 fprintf (SUMA_STDERR, "\nError %s: Failed to allocate for MI.\n", FuncName);
04268 SUMA_RETURN (NULL);
04269 }
04270
04271 data->nd_list = NULL;
04272 data->vxl_list = NULL;
04273 data->ijk_list = NULL;
04274 data->nvox_list = NULL;
04275 data->valArray = NULL;
04276
04277 SUMA_RETURN (data);
04278 }
04279
04280
04281
04282
04283 SUMA_Boolean SUMA_Free_1dData (SUMA_1dData *data)
04284 {
04285 static char FuncName[]={"SUMA_Free_1dData"};
04286
04287 SUMA_ENTRY;
04288
04289 if (!data) {
04290 SUMA_RETURN (YUP);
04291 }
04292 if (data->nd_list) SUMA_free (data->nd_list);
04293 if (data->vxl_list) SUMA_free (data->vxl_list);
04294 if (data->ijk_list) SUMA_free (data->ijk_list);
04295 if (data->nvox_list) SUMA_free (data->nvox_list);
04296 if (data->valArray) SUMA_free (data->valArray);
04297
04298 SUMA_free (data);
04299
04300 SUMA_RETURN (YUP);
04301 }
04302
04303
04304 #ifdef SUMA_AverageMaps_STAND_ALONE
04305
04306 void SUMA_AverageMaps_usage ()
04307
04308 {
04309 printf ("\nUsage: SUMA_AverageMaps <-maps mapFile index nodes valRange> [-cut cutoff] [-prefix fout]\n");
04310 printf ("\n\t-map: mapFile: 1D file name containing map\n\t index: index (begin with zero) of column to read\n\t nodes: index of node column in file (-1 if not included)\n\t valRange: ranges of interest for this column\n\n\t ex: -map file 2 -1 1.5 2 6 10.7\n\t reads third column from 'file', has no given indices \n\t and ranges of interest are 1.5->2 and 6->10.7\n\n\t Note: each file column requires its own '-map'\n\t\t max is 100\n");
04311 printf ("\n\t-cut: value ranges between which activations are not displayed.\n");
04312 printf ("\n\t-prefix: fout is prefix for output files\n\t\t(optional, default AvgMap)\n");
04313 printf ("\n\t Brenna D. Argall LBC/NIMH/NIH bargall@codon.nih.gov \n\t\t Mon February 24 146:23:42 EST 2003\n\n");
04314 exit (0);
04315 }
04316
04317
04318
04319
04320 int main (int argc, char *argv[])
04321 {
04322
04323 static char FuncName[]={"SUMA_AverageMaps-main"};
04324 int numMap, numVR, tmpNumVR, numCuts;
04325 int *clmn=NULL, *nodeClmn=NULL, *numRng=NULL;
04326 float *cutRng=NULL, *valsRng=NULL;
04327 float tmpValRng, tmp1, tmp2;
04328 SUMA_FileName *mapFiles=NULL;
04329 char fout[SUMA_MAX_DIR_LENGTH+SUMA_MAX_NAME_LENGTH];
04330 char avgFileNm[SUMA_MAX_DIR_LENGTH+SUMA_MAX_NAME_LENGTH];
04331 char *input;
04332 SUMA_Boolean brk, LocalHead=YUP, cut, tmpCut;
04333
04334 int kar, i, j=0, k;
04335 int mx_ndx, tmpNumVal, numNoCut, i_currRng;
04336 int *numVal=NULL, *ndx_list=NULL, *tempClmn=NULL, temp;
04337 float *tempValArray=NULL, *avgArray=NULL;
04338 int *avgIndex=NULL, *i_locInfo=NULL;
04339 SUMA_1dData **data=NULL;
04340
04341
04342 if (LocalHead) fprintf (SUMA_STDERR,"%s: Calling SUMA_Create_CommonFields ...\n", FuncName);
04343
04344 SUMAg_CF = SUMA_Create_CommonFields ();
04345 if (SUMAg_CF == NULL) {
04346 fprintf(SUMA_STDERR,"\nError %s: Failed in SUMA_Create_CommonFields\n", FuncName);
04347 exit(1);
04348 }
04349 if (LocalHead) fprintf (SUMA_STDERR,"%s: SUMA_Create_CommonFields Done.\n", FuncName);
04350
04351
04352
04353 if (argc < 2) {
04354 SUMA_AverageMaps_usage ();
04355 exit (1);
04356 }
04357
04358
04359 sprintf( fout, "%s", "AvgMap");
04360 mapFiles = (SUMA_FileName *)SUMA_calloc(100, sizeof(SUMA_FileName));
04361 clmn = SUMA_calloc(100, sizeof(int));
04362 nodeClmn = SUMA_calloc(100, sizeof(int));
04363 valsRng = SUMA_calloc( 1000, sizeof(float));
04364 cutRng = SUMA_calloc( 100, sizeof(float));
04365 numRng = SUMA_calloc( 100, sizeof(int));
04366 numMap = 0; numVR = 0; numCuts = 0;
04367 cut = NOPE;
04368 kar = 1;
04369 brk = NOPE;
04370 while (kar < argc) {
04371 if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
04372 SUMA_AverageMaps_usage ();
04373 exit (1);
04374 }
04375
04376 if (!brk && (strcmp(argv[kar], "-map") == 0 ))
04377 {
04378 if ( numMap >= 100 ) brk = YUP;
04379
04380 ++kar;
04381 if (kar >= argc) {
04382 fprintf (SUMA_STDERR, "need arguments after -map \n");
04383 exit (1);
04384 }
04385
04386
04387 mapFiles[numMap].FileName = argv[kar];
04388
04389
04390 ++kar;
04391 input = argv[kar];
04392 if ( strcmp(input, "-map")==0 || strcmp(input, "-cut")==0 || strcmp(input, "-prefix")==0 ) {
04393 fprintf(SUMA_STDERR, "\nError %s: Improper format for -map option. Exiting.\n", FuncName);
04394 exit(1);
04395 }
04396 else {
04397 clmn[ numMap ] = atoi(input);
04398 if (clmn[numMap]<0) {
04399 fprintf(SUMA_STDERR, "\nError %s: All column indices must be positive integers. Exiting.\n", FuncName);
04400 exit(1);
04401 }
04402 }
04403
04404
04405 ++kar;
04406 input = argv[kar];
04407 if ( strcmp(input, "-map")==0 || strcmp(input, "-cut")==0 || strcmp(input, "-prefix")==0 ) {
04408 fprintf(SUMA_STDERR, "\nError %s: Improper format for -map option. Exiting.\n", FuncName);
04409 exit(1);
04410 }
04411 else {nodeClmn[ numMap ] = atoi(input);
04412
04413
04414 tmpNumVR = 0;
04415 ++kar;
04416 input = argv[kar];
04417 if ( strcmp(input, "-map")==0 || strcmp(input, "-cut")==0 || strcmp(input, "-prefix")==0 ) {
04418
04419 numRng[numMap] = 0;
04420 }
04421 else {
04422 ++kar;
04423 while ( !(strcmp(input, "-map")==0) && !(strcmp(input, "-cut")==0) && !(strcmp(input, "-prefix")==0) && tmpNumVR<100 ) {
04424
04425 tmp1 = atof(argv[kar-1]);
04426 tmp2 = atof(argv[kar]);
04427
04428 if (tmp1<tmp2) {
04429 valsRng[ 2*numVR ] = tmp1;
04430 valsRng[ 2*numVR +1 ] = tmp2;
04431 }
04432 else {
04433 valsRng[ 2*numVR ] = tmp2;
04434 valsRng[ 2*numVR +1 ] = tmp1;
04435 }
04436 fprintf(SUMA_STDERR, "valRng = %f, %f\n", valsRng[2*numVR], valsRng[2*numVR+1]);
04437 ++kar; ++kar; ++numVR; ++tmpNumVR;
04438 if (kar<argc) input = argv[kar];
04439 else break;
04440 }
04441 numRng[numMap] = tmpNumVR;
04442 }
04443
04444 ++numMap;
04445 kar--;
04446 brk = YUP;
04447 }
04448 }
04449
04450 if (!brk && strcmp(argv[kar], "-cut") == 0)
04451 {
04452 ++kar; ++kar;
04453 input = argv[kar];
04454 if ( strcmp(input, "-map")==0 || strcmp(input, "-cut")==0 || strcmp(input, "-prefix")==0 ) {
04455 fprintf(SUMA_STDERR, "\nError %s: Improper format for -cut option (must come in pairs). Exiting.\n", FuncName);
04456 exit(1);
04457 }
04458 else {
04459 cutRng = SUMA_calloc( 100, sizeof(float));
04460 numCuts = 0;
04461 while ( !(strcmp(input, "-map")==0) && !(strcmp(input, "-cut")==0) && !(strcmp(input, "-prefix")==0) && numCuts<100 ) {
04462
04463 tmp1 = atof(argv[kar-1]);
04464 tmp2 = atof(argv[kar]);
04465
04466
04467 if (tmp1<tmp2) {
04468 cutRng[ 2*numCuts ] = tmp1;
04469 cutRng[ 2*numCuts +1 ] = tmp2;
04470 }
04471 else {
04472 cutRng[ 2*numCuts ] = tmp2;
04473 cutRng[ 2*numCuts +1 ] = tmp1;
04474 }
04475 fprintf(SUMA_STDERR, "cutRng = %f, %f\n", cutRng[2*numCuts], cutRng[2*numCuts+1]);
04476 ++kar; ++kar; ++numCuts;
04477 if (kar<argc) input = argv[kar];
04478 else break;
04479 }
04480 }
04481 kar--;
04482 brk = YUP;
04483 }
04484
04485
04486 if (!brk && strcmp(argv[kar], "-prefix") == 0)
04487 {
04488 kar ++;
04489 if (kar >= argc) {
04490 fprintf (SUMA_STDERR, "need argument after -prefix ");
04491 exit (1);
04492 }
04493 sprintf (fout, "%s", argv[kar]);
04494 brk = YUP;
04495 }
04496
04497
04498 if (!brk) {
04499 fprintf (SUMA_STDERR,"\nError %s: Option %s not understood. Try -help for usage\n", FuncName, argv[kar]);
04500 exit (1);
04501 }
04502 else {
04503 brk = NOPE;
04504 kar ++;
04505 }
04506 }
04507
04508
04509 sprintf( avgFileNm, "%s.1D", fout);
04510 if ( SUMA_filexists(avgFileNm)) {
04511 fprintf (SUMA_STDERR,"\nError %s: Output file %s exists.\nWill not overwrite.\n", FuncName, avgFileNm);
04512 exit(1);
04513 }
04514 fprintf(SUMA_STDERR, "\n");
04515
04516
04517
04518
04519 data = (SUMA_1dData **)SUMA_calloc( numMap, sizeof(SUMA_1dData));
04520 numVal = SUMA_calloc( numMap, sizeof(int));
04521 i_locInfo = SUMA_calloc( 6, sizeof(int));
04522 mx_ndx = 0;
04523
04524 for ( i=0; i<numMap; ++i ) {
04525
04526 numVal[i] = -1;
04527 tempClmn = SUMA_calloc( 100, sizeof(int));
04528 if( nodeClmn[i] >= 0) {
04529
04530 for ( j=0; j<6; ++j) {
04531 i_locInfo[j] = j;
04532 }
04533 tempClmn[0] = 2; tempClmn[1] = nodeClmn[i]; tempClmn[2] = clmn[i];
04534 }
04535 else {
04536
04537 for ( j=0; j<6; ++j) {
04538 i_locInfo[j] = -1;
04539 }
04540 tempClmn[0] = 1; tempClmn[1] = clmn[i];
04541 }
04542
04543 data[i] = (SUMA_1dData *)SUMA_Create_1dData();
04544 for (k=0; k<10; ++k) {
04545 fprintf(SUMA_STDERR, "Reading %c", mapFiles[i].FileName[k]);
04546 }
04547 fprintf(SUMA_STDERR, "\n");
04548
04549 SUMA_read1D( mapFiles[i].FileName, tempClmn, i_locInfo, data[i]);
04550
04551 if (data[i]->nd_list[data[i]->N_elem] > mx_ndx)
04552 mx_ndx = data[i]->nd_list[data[i]->N_elem];
04553
04554 SUMA_free(tempClmn);
04555 }
04556
04557
04558
04559
04560
04561 avgArray = SUMA_calloc( mx_ndx, sizeof(float));
04562 for (i=0; i<mx_ndx; ++i) {
04563 avgArray[i] = 0;
04564 i_currRng = 0;
04565 for ( j=0; j<numMap; ++j ) {
04566 if ( i <= data[j]->nd_list[data[j]->N_elem] ) {
04567
04568 if ( data[j]->nd_list[i] == i ) {
04569
04570
04571
04572 if ( numRng[j]==0 ) {
04573
04574 avgArray[i] = avgArray[i] + data[j]->valArray[i];
04575 }
04576 else {
04577
04578 for ( k=0; k<numRng[j]; ++k ) {
04579 if ( valsRng[ i_currRng+2*k ] < data[j]->valArray[i] &&
04580 valsRng[ i_currRng+2*k+1 ] > data[j]->valArray[i] ) {
04581
04582
04583
04584 avgArray[i] = avgArray[i] + data[j]->valArray[i];
04585 break;
04586 }
04587
04588 i_currRng = i_currRng + 2*numRng[j];
04589 }
04590 }
04591 }
04592 }
04593 }
04594 avgArray[i] = (float)avgArray[i] / (float)numMap;
04595 }
04596
04597
04598
04599 avgIndex = SUMA_calloc( mx_ndx, sizeof(int));
04600 if (cut) {
04601
04602 numNoCut=0;
04603 for (i=0; i<mx_ndx; ++i) {
04604 tmpCut = NOPE;
04605 for ( j=0; j<numCuts; ++j ) {
04606 if ( (avgArray[i] > cutRng[2*j]) ||
04607 (avgArray[i] < cutRng[2*j+1]) ) {
04608
04609 tmpCut = YUP;
04610 break;
04611 }
04612 }
04613 if ( tmpCut==NOPE ) {
04614
04615 avgIndex[numNoCut] = i;
04616 ++numNoCut;
04617 }
04618 }
04619 }
04620 else {
04621
04622 for ( i=0; i<mx_ndx; ++i ) { avgIndex[i] = i; }
04623 numNoCut = mx_ndx;
04624 }
04625
04626
04627
04628 fprintf (SUMA_STDERR, "%s: Now writing %s to disk ...\n", FuncName, avgFileNm);
04629 SUMA_write1D ( &numNoCut, avgArray, avgIndex, NULL, avgFileNm);
04630
04631
04632
04633 SUMA_free(mapFiles);
04634 SUMA_free(clmn);
04635 SUMA_free(avgArray);
04636 SUMA_free(avgIndex);
04637
04638 exit(0);
04639
04640 }
04641 #endif
04642
04643
04644
04645 #ifdef SUMA_GiveColor_STAND_ALONE
04646
04647 void SUMA_GiveColor_usage ()
04648
04649 {
04650 printf ("\nUsage: SUMA_GiveColor <-file fileNm index nodes> [-cr colRange] [-gr gradRange] [-seg numSeg] [-prefix fout]\n");
04651 printf ("\n\t-file: fileNm: 1D file name containing values\n\t index: index (begin with zero) of column to read\n\t nodes: index of node column in file (-1 if not included)\n");
04652 printf ("\n\n\t-cr: specifies roygbiv color range.\n\n\t Note: input 0->7 (roygbivr), low->high.\n\t 10 divisions per color accepted,\n\t so 3.4!=3 && 3.4!=4 but 3.4=3.42=3.446 etc\n\t Note: May specify each color explicitly (max 100).\n\n\t ex: -cr 1 3.5 has range of orange->green/blue (low->high).\n");
04653 printf ("\n\n\t-gr: value range over which the color gradient increments.\n\t input low->high\n\n\t ex: -gr 5 20 color gradient exists between values 5 and 20.\n");
04654 printf ("\n\n\t-seg: number of discrete segments of color range.\n\n\t Note: color range will not be continuous.\n");
04655 printf ("\n\t-prefix: prefix for output files (optional, default GivCol).\n");
04656 printf ("\n\t Brenna D. Argall LBC/NIMH/NIH bargall@codon.nih.gov \n\t\t Thursday March 20 14:23:42 EST 2003\n\n");
04657 exit (0);
04658 }
04659
04660
04661
04662
04663 int main (int argc, char *argv[])
04664 {
04665
04666 static char FuncName[]={"SUMA_GiveColor-main"};
04667 float *rngGrad=NULL, *rngCol=NULL, tmp1, tmp2;
04668 int clmn, nodeClmn, *numRng=NULL;
04669 char *input=NULL;
04670 SUMA_FileName file;
04671 char fout[SUMA_MAX_DIR_LENGTH+SUMA_MAX_NAME_LENGTH];
04672 char colFileNm[SUMA_MAX_DIR_LENGTH+SUMA_MAX_NAME_LENGTH];
04673 SUMA_Boolean brk, LocalHead=NOPE;
04674
04675 int i, j, kar, numSeg, numCol;
04676 int *columns=NULL, *i_locInfo=NULL;
04677 float high, low, *tmpArray=NULL;
04678 float *color=NULL, *colMap=NULL, *valArray=NULL;
04679 SUMA_1dData *data=NULL;
04680
04681
04682 if (LocalHead) fprintf (SUMA_STDERR,"%s: Calling SUMA_Create_CommonFields ...\n", FuncName);
04683
04684 SUMAg_CF = SUMA_Create_CommonFields ();
04685 if (SUMAg_CF == NULL) {
04686 fprintf(SUMA_STDERR,"\nError %s: Failed in SUMA_Create_CommonFields\n", FuncName);
04687 exit(1);
04688 }
04689 if (LocalHead) fprintf (SUMA_STDERR,"%s: SUMA_Create_CommonFields Done.\n", FuncName);
04690
04691
04692
04693 if (argc < 2) {
04694 SUMA_GiveColor_usage ();
04695 exit (1);
04696 }
04697
04698
04699 sprintf( fout, "%s", "GivCol");
04700 clmn = 0; nodeClmn = -1;
04701 numSeg = -1; numCol = 0;
04702 kar = 1;
04703 brk = NOPE;
04704
04705 while (kar < argc) {
04706 if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
04707 SUMA_GiveColor_usage ();
04708 exit (1);
04709 }
04710 if (!brk && (strcmp(argv[kar], "-file") == 0 ))
04711 {
04712 ++kar;
04713 if (kar >= argc) {
04714 fprintf (SUMA_STDERR, "need arguments after -file \n");
04715 exit (1);
04716 }
04717
04718
04719 file.FileName = argv[kar];
04720
04721
04722 ++kar;
04723 input = argv[kar];
04724 if ( strcmp(input, "-file")==0 || strcmp(input, "-cr")==0 || strcmp(input, "-gr")==0 || strcmp(input, "-seg")==0 || strcmp(input, "-prefix")==0 ) {
04725 fprintf(SUMA_STDERR, "\nError %s: Improper format for -file option. Exiting.\n", FuncName);
04726 exit(1);
04727 }
04728 else {
04729 clmn = atoi(input);
04730 if (clmn<0) {
04731 fprintf(SUMA_STDERR, "\nError %s: All column indices must be positive integers. Exiting.\n", FuncName);
04732 exit(1);
04733 }
04734 }
04735
04736
04737 ++kar;
04738 input = argv[kar];
04739 if ( strcmp(input, "-file")==0 || strcmp(input, "-cr")==0 || strcmp(input, "-gr")==0 || strcmp(input, "-seg")==0 || strcmp(input, "-prefix")==0 ) {
04740 fprintf(SUMA_STDERR, "\nError %s: Improper format for -map option. Exiting.\n", FuncName);
04741 exit(1);
04742 }
04743 else nodeClmn = atoi(input);
04744 brk = YUP;
04745 }
04746
04747
04748 if (!brk && strcmp(argv[kar], "-cr") == 0)
04749 {
04750 kar ++; kar++;
04751 if (kar >= argc) {
04752 fprintf (SUMA_STDERR, "need at least two arguments after -cr ");
04753 exit (1);
04754 }
04755 rngCol = SUMA_calloc( 100, sizeof(int));
04756 rngCol[0] = atof(argv[kar-1]);
04757 rngCol[1] = atof(argv[kar]);
04758 numCol = 2;
04759 kar ++;
04760 input = argv[kar];
04761
04762
04763 while ( !(strcmp(input, "-file")==0) && !(strcmp(input, "-cr")==0) && !(strcmp(input, "-gr")==0) && !(strcmp(input, "-seg")==0) && !(strcmp(input, "-prefix")==0) && numCol<100 ) {
04764
04765 rngCol[numCol++] = atof(input);
04766 kar ++;
04767 if (kar<argc) input = argv[kar];
04768 else break;
04769 }
04770 kar --;
04771 brk = YUP;
04772 }
04773
04774 if (!brk && strcmp(argv[kar], "-seg") == 0)
04775 {
04776 kar ++;
04777 if (kar >= argc) {
04778 fprintf (SUMA_STDERR, "need argument after -seg ");
04779 exit (1);
04780 }
04781 numSeg = atoi(argv[kar]);
04782 brk = YUP;
04783 }
04784
04785 if (!brk && strcmp(argv[kar], "-gr") == 0)
04786 {
04787 kar ++; kar++;
04788 if (kar >= argc) {
04789 fprintf (SUMA_STDERR, "need two arguments after -gr ");
04790 exit (1);
04791 }
04792 rngGrad = SUMA_calloc(2, sizeof(float));
04793 rngGrad[0] = atof(argv[kar-1]);
04794 rngGrad[1] = atof(argv[kar]);
04795 brk = YUP;
04796 }
04797
04798
04799 if (!brk && strcmp(argv[kar], "-prefix") == 0)
04800 {
04801 kar ++;
04802 if (kar >= argc) {
04803 fprintf (SUMA_STDERR, "need argument after -prefix ");
04804 exit (1);
04805 }
04806 sprintf (fout, "%s", argv[kar]);
04807 brk = YUP;
04808 }
04809
04810
04811 if (!brk) {
04812 fprintf (SUMA_STDERR,"\nError %s: Option %s not understood. Try -help for usage\n", FuncName, argv[kar]);
04813 exit (1);
04814 }
04815 else {
04816 brk = NOPE;
04817 kar ++;
04818 }
04819 }
04820
04821
04822 sprintf( colFileNm, "%s.col", fout);
04823 if ( SUMA_filexists(colFileNm) ) {
04824 fprintf (SUMA_STDERR,"\nError %s: Output file %s exists.\nWill not overwrite.\n\n", FuncName, colFileNm);
04825 exit(1);
04826 }
04827 if ( numCol!=2 && (numSeg!=numCol) && numCol!=0 ) {
04828 fprintf(SUMA_STDERR, "\nError %s: If colors given are not merely endpoints, exactly numSeg (%d, not %d) must be given. Exiting.\n", FuncName, numSeg, numCol);
04829 exit(1);
04830 }
04831
04832
04833
04834
04835 i_locInfo = SUMA_calloc( 6, sizeof(int));
04836
04837
04838 for ( j=1; j<6; ++j) {
04839 i_locInfo[j] = -1;
04840 }
04841
04842 if( nodeClmn >= 0) {
04843
04844 i_locInfo[0] = nodeClmn;
04845 columns = SUMA_calloc( 3, sizeof(int));
04846 columns[0] = 2; columns[1] = nodeClmn; columns[2] = clmn;
04847 }
04848 else {
04849
04850 i_locInfo[0] = -1;
04851 columns = SUMA_calloc( 2, sizeof(int));
04852 columns[0] = 1; columns[1] = clmn;
04853 }
04854
04855 data = (SUMA_1dData *)SUMA_Create_1dData();
04856 SUMA_read1D( file.FileName, columns, i_locInfo, data);
04857
04858
04859
04860 fprintf(SUMA_STDERR, "numseg = %d\n", numSeg);
04861 if ( numCol==numSeg )
04862 color = SUMA_createColGradient( rngCol, numSeg, YUP);
04863 else color = SUMA_createColGradient( rngCol, numSeg, NOPE);
04864
04865 if (numSeg<0) numSeg = 700;
04866 colMap = SUMA_assignColors( data->valArray, color, data->N_elem, numSeg, rngGrad, NULL);
04867
04868
04869
04870 fprintf (SUMA_STDERR, "\n%s: Now writing %s to disk ...\n\n", FuncName, colFileNm);
04871 SUMA_writeColorFile (colMap, data->N_elem, data->nd_list, colFileNm);
04872
04873
04874
04875 SUMA_free(columns);
04876 if (rngCol!=NULL) SUMA_free(rngCol);
04877 if (rngGrad!=NULL) SUMA_free(rngGrad);
04878 SUMA_Free_1dData(data);
04879 SUMA_free(i_locInfo);
04880 SUMA_free(color);
04881 SUMA_free(colMap);
04882
04883 exit(0);
04884
04885 }
04886 #endif
04887
04888
04889 #ifdef SUMA_Test_STAND_ALONE
04890
04891 int main (int argc, char *argv[])
04892 {
04893
04894 static char FuncName[]={"SUMA_Test-main"};
04895 float *valArray = NULL, *srtdArray = NULL;
04896 int num, i_curr, i=0;
04897
04898
04899 int *i_colm=NULL, *i_locInfo=NULL, N_Node;
04900 SUMA_1dData *data=NULL;
04901
04902 float *col=NULL, *rng=NULL, *points=NULL, *colors=NULL, *assgnCols=NULL;
04903
04904 char fileNm[1000];
04905 int *numW=NULL, j;
04906
04907 SUMA_Boolean LocalHead=NOPE;
04908
04909
04910 if (LocalHead) fprintf (SUMA_STDERR,"%s: Calling SUMA_Create_CommonFields ...\n", FuncName);
04911
04912 SUMAg_CF = SUMA_Create_CommonFields ();
04913 if (SUMAg_CF == NULL) {
04914 fprintf(SUMA_STDERR,"\nError %s: Failed in SUMA_Create_CommonFields\n", FuncName);
04915 exit(1);
04916 }
04917 if (LocalHead) fprintf (SUMA_STDERR,"%s: SUMA_Create_CommonFields Done.\n", FuncName);
04918
04919
04920
04921
04922
04923
04924
04925
04926
04927
04928
04929
04930
04931
04932
04933
04934
04935
04936 points = SUMA_calloc( 300, sizeof(float));
04937 for (i=0; i<300; ++i) {
04938 points[i] = i;
04939 }
04940 col = SUMA_calloc( 7, sizeof(float));
04941 for (i=0; i<7; ++i) {
04942 col[i] = i;
04943 }
04944
04945
04946
04947 rng = SUMA_calloc( 2, sizeof(float));
04948 rng[0] = 50; rng[1] = 250;
04949 fprintf(SUMA_STDERR, "before\n");
04950 colors = SUMA_createColGradient( col, 7, YUP );
04951 fprintf(SUMA_STDERR, "after\n");
04952 assgnCols = SUMA_assignColors( points, colors, 300, 7, rng, NULL );
04953 SUMA_writeColorFile (assgnCols, 300, NULL, "test_ccg.col");
04954
04955 SUMA_free(points);
04956 SUMA_free(col);
04957 SUMA_free(colors);
04958
04959
04960
04961
04962
04963
04964
04965
04966
04967
04968
04969
04970
04971
04972
04973
04974
04975
04976
04977
04978
04979
04980
04981
04982
04983
04984
04985
04986
04987
04988
04989
04990
04991
04992
04993
04994
04995
04996
04997
04998
04999
05000
05001
05002 exit(0);
05003
05004 }
05005 #endif
