1 #include "SUMA_suma.h"
2 
3 // drg - mac os 10.15 build needs explicit declaration
4 // this should be in glext.h. If warnings appear on 10.15 and higher Macs
5 // can remove this prototype
6 #ifdef DARWIN
7   extern void glWindowPos2s (GLshort x, GLshort y);
8 #endif
9 extern SUMA_CommonFields *SUMAg_CF;
10 extern SUMA_DO *SUMAg_DOv;
11 extern SUMA_SurfaceViewer *SUMAg_SVv;
12 extern int SUMAg_N_SVv;
13 extern int SUMAg_N_DOv;
14 
15 /* Volume rendering code, based on GLUT-3.7's advanced97/volume.c
16 See also SUMA_GLUT_volumedemo.c*/
17 
18 static GLfloat lightpos[4] = {150., 150., 150., 1.f};
19 
20 
21 // put the clip geometry (planes etc.) here for picking...
22 
SUMA_MoveCutplane(SUMA_VolumeObject * VO,int iplane,float d)23 int SUMA_MoveCutplane (SUMA_VolumeObject *VO, int iplane, float d)
24 {
25    static char FuncName[]={"SUMA_MoveCutplane"};
26 
27    SUMA_ENTRY;
28 
29    if (iplane < 0 || iplane > 5) {
30       SUMA_S_Err("Bad plane index");
31       SUMA_RETURN(0);
32    }
33 
34    VO->CutPlane[iplane][3] = VO->CutPlane[iplane][3]+d;
35    if (!SUMA_SetTextureClipPlaneSurface(VO, iplane)) {
36       SUMA_S_Err("Failed to set cutplane surface");
37       SUMA_RETURN(0);
38    }
39 
40    SUMA_RETURN(1);
41 }
42 
SUMA_adset_to_VE(SUMA_VolumeObject * VO,THD_3dim_dataset ** dsetp)43 SUMA_DSET *SUMA_adset_to_VE(SUMA_VolumeObject *VO, THD_3dim_dataset **dsetp)
44 {
45    static char FuncName[]={"SUMA_adset_to_VE"};
46    THD_3dim_dataset *dset=NULL;
47    THD_3dim_dataset *odset=NULL;
48    SUMA_DSET *sdset=NULL;
49    int n_VE=0, OverInd, OKdup=0, loc[2];
50    char orcode[6], *np=NULL, *dsetcmap=NULL, *forcode;
51    SUMA_ALL_DO *ado=(SUMA_ALL_DO *)VO;
52    SUMA_OVERLAYS *colplane = NULL, *curcolplane=NULL;
53    SUMA_Boolean SetupOverlay = YUP, MakeOverlayCurrent = YUP;
54    SUMA_Boolean LocalHead = NOPE;
55 
56    SUMA_ENTRY;
57 
58    dset = *dsetp;
59 
60    /* make sure data part of dset is loaded, not just the header          */
61    DSET_mallocize(dset); DSET_load(dset);
62    orcode[0] = ORIENT_typestr[dset->daxes->xxorient][0] ;
63    orcode[1] = ORIENT_typestr[dset->daxes->yyorient][0] ;
64    orcode[2] = ORIENT_typestr[dset->daxes->zzorient][0] ;
65    orcode[3] = '\0';
66    SUMA_LHv("dset orcode is %s\n", orcode);
67 
68    if ((forcode = getenv("SUMA_VO_Reorient")) &&
69        strcmp(forcode,"NO") && strcmp(forcode,"No") && strcmp(forcode,"no") ) {
70        if (!SUMA_ok_orstring(forcode)) {
71          SUMA_S_Err("Bad orientation string %s in env SUMA_VO_Reorient\n"
72                     "No reorienting done.", forcode);
73        } else if (strcmp(orcode,forcode) ) {
74          char sss[5];
75          SUMA_S_Note("Resampling %s from %s to %s, per user request.\n",
76                      DSET_HEADNAME(dset), orcode, forcode);
77          odset = r_new_resam_dset(dset, NULL, 0.0, 0.0, 0.0,
78                                   forcode, MRI_LINEAR, NULL, 1, 1);
79          sprintf(sss, ".%s",forcode);
80          np = SUMA_append_string(DSET_PREFIX(dset), sss);
81          EDIT_dset_items(  odset ,
82                          ADN_prefix      , np,
83                          ADN_none ) ;
84          tross_Copy_History( dset , odset ) ;
85          DSET_delete(dset); dset = odset; odset = NULL;
86          if (LocalHead && 0) {
87             SUMA_LH("Writing resampled dset");
88             DSET_write(dset);
89          }
90          SUMA_free(np); np = NULL;
91          *dsetp = dset;
92       }
93    }
94 
95    sdset = SUMA_afnidset2sumadset(&dset, 1, 1, 0);
96    if (dset) {
97       SUMA_S_Warn("Leakage! dset should be null by now...");
98    }
99 
100    /* Does this dset have a built in colormap?
101       If it does, then loadit into SCM */
102    if (!SUMA_Insert_Cmap_of_Dset(sdset)) {
103       SUMA_S_Err("Failed to insert Cmap");
104       SUMA_FreeDset(sdset); sdset = NULL;
105       SUMA_RETURN(NOPE);
106    }
107 
108    if (SetupOverlay) {
109       OverInd = -1;
110       SUMA_LH("Looking for pre-existing overlay on %s", SDSET_LABEL(sdset));
111       if ((colplane = SUMA_Fetch_OverlayPointerByDset (
112                               ado, sdset, &OverInd))) {
113          SUMA_LH("Col plane already present");
114          OKdup = 0;
115          if (colplane->OptScl->Clusterize) {
116             SUMA_S_Warn("Nothing implemented for volume clustering here,"
117                         "just going through the motions...");
118             colplane->OptScl->RecomputeClust = 1;
119          }
120       } else {
121          OKdup = 1;
122          SUMA_LH("Creating anew");
123          OverInd = SUMA_ADO_N_Overlays(ado);
124       }
125       if (!(colplane = SUMA_CreateOverlayPointer ( ADO_LABEL(ado),
126                                                 sdset, ADO_ID(ado), colplane))) {
127          SUMA_S_Err("Failed to create overlay");
128          SUMA_RETURN(NOPE);
129       }
130       colplane->isBackGrnd = NOPE;
131 
132       /* Add this plane to VO's Overlays */
133       SUMA_LH("VO has %d overlays", SUMA_ADO_N_Overlays(ado));
134       if (!SUMA_AddNewPlane (ado, colplane, SUMAg_DOv,
135                              SUMAg_N_DOv, OKdup)) {
136          SUMA_SL_Err("Failed in SUMA_AddNewPlane");
137          SUMA_FreeOverlayPointer(colplane);
138          if (!SUMA_DeleteDsetPointer(&sdset, SUMAg_CF->DsetList)) {
139             SUMA_S_Err("Failed to delete dset pointer");
140          }
141          SUMA_RETURN(NOPE);
142       }
143       SUMA_LH("VO now has %d overlays", SUMA_ADO_N_Overlays(ado));
144 
145       /* Match previous setting? */
146       curcolplane = SUMA_ADO_CurColPlane(ado);
147       if (SUMA_PreserveOverlaySettings(curcolplane, colplane)) {
148                                  /* attempt to preserve current situation */
149          SUMA_OVERLAYS *settingPlane = NULL;
150          settingPlane = curcolplane;
151          colplane->GlobalOpacity = settingPlane->GlobalOpacity;
152          colplane->ShowMode = settingPlane->ShowMode;
153          colplane->OptScl->BrightFact = settingPlane->OptScl->BrightFact;
154          colplane->OptScl->find = settingPlane->OptScl->find;
155          colplane->OptScl->tind = settingPlane->OptScl->tind;
156          colplane->OptScl->bind = settingPlane->OptScl->bind;
157          colplane->OptScl->UseThr = settingPlane->OptScl->UseThr;
158          colplane->OptScl->UseBrt = settingPlane->OptScl->UseBrt;
159          colplane->OptScl->ThrMode = settingPlane->OptScl->ThrMode;
160          colplane->OptScl->ThreshRange[0] =
161                                        settingPlane->OptScl->ThreshRange[0];
162          colplane->OptScl->ThreshRange[1] =
163                                        settingPlane->OptScl->ThreshRange[1];
164          colplane->OptScl->BrightRange[0] =
165                                        settingPlane->OptScl->BrightRange[0];
166          colplane->OptScl->BrightRange[1] =
167                                        settingPlane->OptScl->BrightRange[1];
168          colplane->OptScl->BrightMap[0] =
169                                        settingPlane->OptScl->BrightMap[0];
170          colplane->OptScl->BrightMap[1] =
171                                        settingPlane->OptScl->BrightMap[1];
172          colplane->SymIrange = settingPlane->SymIrange;
173          colplane->OptScl->IntRange[0] = settingPlane->OptScl->IntRange[0];
174          colplane->OptScl->IntRange[1] = settingPlane->OptScl->IntRange[1];
175          dsetcmap = NI_get_attribute(sdset->ngr,"SRT_use_this_cmap");
176          if (dsetcmap) {
177             SUMA_STRING_REPLACE(colplane->cmapname, dsetcmap);
178          } else {
179             SUMA_STRING_REPLACE(colplane->cmapname,
180                                 settingPlane->cmapname);
181          }
182          colplane->OptScl->Clusterize = settingPlane->OptScl->Clusterize;
183          colplane->OptScl->ClustOpt->AreaLim =
184             settingPlane->OptScl->ClustOpt->AreaLim;
185          colplane->OptScl->ClustOpt->DistLim =
186             settingPlane->OptScl->ClustOpt->DistLim;
187       } else {
188          /* set the opacity, index column and the range */
189          colplane->GlobalOpacity = YUP;
190          colplane->ShowMode = SW_SurfCont_DsetViewCol;
191          if (!OKdup) {/* only set this if first time creating plane*/
192             colplane->OptScl->BrightFact = 0.8;
193          }
194          colplane->OptScl->find = 0;
195          colplane->OptScl->tind = 0;
196          colplane->OptScl->bind = 0;
197          #if 0
198          SUMA_GetDsetColRange(sdset, 0, colplane->OptScl->IntRange, loc);
199          #else
200          colplane->OptScl->RangeUnits = SUMA_PERC_VALUE_UNITS;
201          colplane->OptScl->IntRange[0] = 2;
202          colplane->OptScl->IntRange[1] = 98;
203          colplane->OptScl->AutoIntRange = 0; /* turn of auto ranging
204                                  Otherwise SurfCont fields won't reflect
205                                  what will eventually get put into IntRange[]
206                                  in the ScaleToMap functions when the
207                                  SurfCont is first opened */
208          #endif
209          if (colplane->SymIrange) {
210             colplane->OptScl->IntRange[0] =
211                -fabs(SUMA_MAX_PAIR( colplane->OptScl->IntRange[0],
212                                     colplane->OptScl->IntRange[1]));
213             colplane->OptScl->IntRange[1] =
214                -colplane->OptScl->IntRange[0];
215          }
216 
217          /* stick a colormap onto that plane ? */
218          dsetcmap = NI_get_attribute(sdset->ngr,"SRT_use_this_cmap");
219          if (dsetcmap) {
220             SUMA_STRING_REPLACE(colplane->cmapname, dsetcmap);
221          } else {
222             /* don't worry, there's a default one */
223          }
224       }
225       if (colplane->OptScl->Clusterize)
226          colplane->OptScl->RecomputeClust = 1;
227       /* colorize the plane */
228       SUMA_LH("Colorizing Plane");
229       SUMA_ColorizePlane(colplane);
230 
231       /* SUMA_Show_ColorOverlayPlanes(&colplane, 1, 1); */
232 
233       /* set the new curColPlane to the newly loaded plane */
234       if (MakeOverlayCurrent) {
235          SUMA_X_SurfCont *SurfCont=SUMA_ADO_Cont(ado);
236          if (!SurfCont) {
237             SUMA_S_Err("OMG");
238          } else {
239             SurfCont->curColPlane = colplane;
240          }
241       }
242    }
243 
244    /* Add as new volume element */
245    n_VE = SUMA_VO_NumVE(VO);
246    VO->VE[n_VE] = (SUMA_VolumeElement*)SUMA_calloc(1,
247                                        sizeof(SUMA_VolumeElement));
248    VO->VE[n_VE]->dset_idcode_str = SUMA_copy_string(SDSET_ID(sdset));
249    if (!SUMA_InsertDsetPointer(&sdset, SUMAg_CF->DsetList, 0)) {
250       SUMA_S_Err("Failed to inset dset pointer. Replace is not enabled");
251       SUMA_RETURN(NOPE);
252    }
253 
254    if (!SUMA_VE_Set_Dims(VO->VE, n_VE)) {
255       SUMA_S_Err("Failed to set dims");
256       SUMA_RETURN(NOPE);
257    }
258 
259    /* Prep colorplane */
260    SUMA_LH("Prep texture");
261    if (!(VO->VE[n_VE]->texvec =
262          SUMA_VE_to_tex3d(VO->VE, n_VE, (byte)n_VE))) {
263       SUMA_S_Err("Failed in dset to text3d");
264       VO = SUMA_FreeVolumeObject(VO);
265       SUMA_RETURN(NOPE);
266    }
267    SUMA_LHv("Have slot %d,%s\n", n_VE, SUMA_VE_Headname(VO->VE, n_VE));
268 
269    /* Set the box limits */
270    SUMA_dset_extreme_corners(sdset,
271                            VO->VE[n_VE]->vo0, VO->VE[n_VE]->voN, 1);
272    SUMA_dset_extreme_corners(sdset,
273                            VO->VE[n_VE]->bo0, VO->VE[n_VE]->boN, 0);
274 
275 
276    SUMA_RETURN(sdset);
277 }
278 
SUMA_VE_to_tex3d(SUMA_VolumeElement ** VE,int iVE,byte col)279 GLubyte * SUMA_VE_to_tex3d(SUMA_VolumeElement **VE, int iVE, byte col)
280 {
281    static char FuncName[]={"SUMA_VE_to_tex3d"};
282    char *filename;
283    GLubyte *tex3ddata;
284    GLint max3dtexdims; /* maximum allowed 3d texture dimension */
285    GLint newval;
286    SUMA_DSET *sdset=NULL;
287    SUMA_Boolean LocalHead = NOPE;
288 
289    SUMA_ENTRY;
290 
291    if (!(sdset = SUMA_VE_dset(VE, iVE))) {
292       SUMA_S_Err("No volume found");
293       SUMA_RETURN(NULL);
294    }
295    if (!(tex3ddata =
296             (GLubyte *)SUMA_malloc(4*SUMA_VE_Nvox(VE, iVE)*sizeof(GLubyte)))) {
297       SUMA_S_Crit("Failed to allocate.");
298       SUMA_RETURN(NULL);
299    }
300    if (LocalHead) {
301       SUMA_LHv("Copying %d intensity in R, G, B, A \n",
302                SUMA_VE_Nvox(VE, iVE)*4);
303    }
304    #if 0
305    if (!SUMA_Colorize_dset_OBSOLETE(sdset, tex3ddata, col)) {
306       SUMA_S_Err("Failed to colorize VO");
307       SUMA_RETURN(NULL);
308    }
309    #endif
310 
311    SUMA_RETURN(tex3ddata);
312 }
313 
314 /* This function here is for illustrative purposes.
315    It may be too inefficient to have to allocate and
316    free SV for each colorizing operation.
317 
318 */
SUMA_Colorize_dset_OBSOLETE(SUMA_DSET * dset,byte * tex3ddata,byte colopt)319 SUMA_Boolean SUMA_Colorize_dset_OBSOLETE(SUMA_DSET *dset,
320                                 byte *tex3ddata, byte colopt)
321 {
322    static char FuncName[]={"SUMA_Colorize_dset_OBSOLETE"};
323    static SUMA_SCALE_TO_MAP_OPT *Opt=NULL;
324    static SUMA_COLOR_MAP *ColMap=NULL;
325    SUMA_COLOR_SCALED_VECT * SV= NULL;
326    float *floatvol=NULL;
327    byte *bytevol=NULL, am=0;
328    int i, j, i3;
329    int av=0;
330    SUMA_Boolean ans = YUP;
331    SUMA_Boolean LocalHead = NOPE;
332 
333    SUMA_ENTRY;
334 
335    /* setup some defaults for now. */
336    if (!Opt) { /* setup once, this struct should be part of VO, perhaps.
337                   Or part of sv structure, perhaps... */
338       Opt = SUMA_ScaleToMapOptInit();
339       Opt->alaAFNI=YUP;
340       if (LocalHead) {
341          SUMA_ShowScaleToMapOpt(Opt, NULL, 1);
342       }
343    }
344 
345    if (!ColMap) {
346       char *eee=getenv("SUMA_VO_ColorMap");
347       if (eee) {
348          if (!(ColMap = SUMA_FindNamedColMap(eee))) {
349             SUMA_S_Errv( "Colormap %s not found.\n"
350                          "Using bgyr64 instead.\n", eee);
351          }
352       } else {
353          eee = "bgyr64";
354       }
355       ColMap = SUMA_FindNamedColMap(eee);
356       if (!ColMap) {
357          SUMA_S_Errv("Could not get %s\n", eee);
358          SUMA_RETURN(NOPE);
359       }
360 
361    }
362    /* No need to colorize as was done in the days of olde.
363       Now colorization is handled in SUMA_Overlays_2_GLCOLAR4,
364       SUMA_ColorizePlane, and SUMA_ScaleToMap_Interactive */
365 
366    SUMA_RETURN(ans);
367 
368    /* Create temporary holding structure for colorized vectors */
369    if (!(SV = SUMA_Create_ColorScaledVect(SDSET_NVOX(dset), 0))) {
370       SUMA_S_Err("Failed to create SV");
371       ans = NOPE;      goto CLEANUP;
372    }
373 
374 
375    /* copy into tex3ddata */
376    if (!colopt) {
377       bytevol = (byte *)SUMA_calloc(SDSET_NVOX(dset), sizeof(byte));
378       if (!bytevol) {
379          SUMA_S_Err("Failed to allocate for bytevol");
380          ans = NOPE;      goto CLEANUP;
381       }
382       EDIT_coerce_scale_type(
383                   SDSET_VECLEN(dset) ,
384                   SDSET_BRICK_FACTOR(dset,0) ,
385                   SDSET_BRICK_TYPE(dset,0),
386                   SDSET_ARRAY(dset, 0) ,   /* input  */
387                   MRI_byte               , bytevol  ) ;
388       j=0;
389       for(i = 0; i < SDSET_NVOX(dset); i++) {
390          tex3ddata[j] = bytevol[i]; ++j;
391          tex3ddata[j] = bytevol[i]; ++j;
392          tex3ddata[j] = bytevol[i]; ++j;
393          tex3ddata[j] = bytevol[i]; ++j;
394       }
395       if (bytevol) SUMA_free(bytevol); bytevol=NULL;
396    } else {
397       /* put dset values in temporary floatvol float vector */
398       floatvol = (float *)SUMA_calloc(SDSET_NVOX(dset), sizeof(float));
399       if (!floatvol) {
400          SUMA_S_Err("Failed to allocate for floatvol");
401          ans = NOPE;      goto CLEANUP;
402       }
403       EDIT_coerce_scale_type(
404                   SDSET_NVOX(dset) ,
405                   SDSET_BRICK_FACTOR(dset,0) ,
406                   SDSET_BRICK_TYPE(dset,0),
407                   SDSET_ARRAY(dset, 0) ,   /* input  */
408                   MRI_float               , floatvol  ) ;
409       if (!SUMA_ScaleToMap_alaAFNI (floatvol, SDSET_NVOX(dset),
410                                  0.0, ColMap, Opt, SV)) {
411          SUMA_S_Err("Failed to colorize");
412          ans = NOPE;      goto CLEANUP;
413       }
414 
415       j=0;
416       for(i = 0; i < SDSET_NVOX(dset); i++) {
417          i3 = 3*i; am = 0;
418          tex3ddata[j] = (byte)(SV->cV[i3  ] * 255);
419                                    am = tex3ddata[j];  ++j;
420          tex3ddata[j] = (byte)(SV->cV[i3+1] * 255);
421             if (tex3ddata[j] > am) am = tex3ddata[j];  ++j;
422          tex3ddata[j] = (byte)(SV->cV[i3+2] * 255);
423             if (tex3ddata[j] > am) am = tex3ddata[j];  ++j;
424          if (SV->isMasked[i]) {
425             tex3ddata[j] = 0;
426          } else {
427             tex3ddata[j] = am;
428          }
429          ++j;
430       }
431 
432    }
433    CLEANUP:
434    if (SV) SUMA_Free_ColorScaledVect(SV); SV = NULL;
435    if (bytevol) SUMA_free(bytevol); bytevol = NULL;
436    if (floatvol) SUMA_free(floatvol); floatvol = NULL;
437 
438    SUMA_RETURN(ans);
439 }
440 
SUMA_CreateSphereList(void)441 void SUMA_CreateSphereList(void)
442 {
443    static char FuncName[]={"SUMA_CreateSphereList"};
444 
445    SUMA_ENTRY;
446 
447    SUMA_S_Note("Making sphere display list");
448    /* make a display list containing a sphere */
449    glNewList(1, GL_COMPILE);
450    {
451       static GLfloat lightpos[] = {150.f, 150.f, 150.f, 1.f};
452       static GLfloat material[] = {1.f, .5f, 1.f, 1.f};
453       GLUquadricObj *qobj = gluNewQuadric();
454       glPushAttrib(GL_LIGHTING_BIT);
455       glEnable(GL_LIGHTING);
456       glEnable(GL_LIGHT2);
457       glLightfv(GL_LIGHT2, GL_POSITION, lightpos);
458       glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, material);
459       gluSphere(qobj, 20.f, 20, 20);
460       gluDeleteQuadric(qobj);
461       glPopAttrib();
462    }
463    glEndList();
464    SUMA_RETURNe;
465 }
466 
SUMA_RecordEnablingState(SUMA_EnablingRecord * SER,char * Label)467 void SUMA_RecordEnablingState(SUMA_EnablingRecord *SER, char *Label)
468 {
469    static char FuncName[]={"SUMA_RecordEnablingState"};
470 
471    SUMA_ENTRY;
472 
473    if (!SER) {
474       SUMA_S_Err("NULL SER, how am I to record?");
475       SUMA_RETURNe;
476    }
477    snprintf(SER->Label,255, "%s", Label ? Label:"Unabeled");
478    SER->ALPHA_TEST = glIsEnabled(GL_ALPHA_TEST);
479    SER->DEPTH_TEST = glIsEnabled(GL_DEPTH_TEST);
480    SER->TEXTURE_3D_EXT = glIsEnabled(GL_TEXTURE_3D_EXT);
481    SER->TEXTURE_3D = glIsEnabled(GL_TEXTURE_3D);
482    SER->TEXTURE_2D = glIsEnabled(GL_TEXTURE_2D);
483    SER->TEXTURE_GEN_S = glIsEnabled(GL_TEXTURE_GEN_S);
484    SER->TEXTURE_GEN_T = glIsEnabled(GL_TEXTURE_GEN_T);
485    SER->TEXTURE_GEN_R = glIsEnabled(GL_TEXTURE_GEN_R);
486    SER->COLOR_MATERIAL = glIsEnabled(GL_COLOR_MATERIAL);
487    SER->CLIP_PLANE0 = glIsEnabled(GL_CLIP_PLANE0);
488    SER->CLIP_PLANE1 = glIsEnabled(GL_CLIP_PLANE1);
489    SER->CLIP_PLANE2 = glIsEnabled(GL_CLIP_PLANE2);
490    SER->CLIP_PLANE3 = glIsEnabled(GL_CLIP_PLANE3);
491    SER->CLIP_PLANE4 = glIsEnabled(GL_CLIP_PLANE4);
492    SER->CLIP_PLANE5 = glIsEnabled(GL_CLIP_PLANE5);
493    SER->LIGHTING = glIsEnabled(GL_LIGHTING);
494    SER->LIGHT0 = glIsEnabled(GL_LIGHT0);
495    SER->LIGHT1 = glIsEnabled(GL_LIGHT1);
496    SER->LIGHT2 = glIsEnabled(GL_LIGHT2);
497    SER->BLEND = glIsEnabled(GL_BLEND);
498    SER->LINE_SMOOTH = glIsEnabled(GL_LINE_SMOOTH);
499    /* SER-> = glIsEnabled(GL_); */
500    glGetFloatv(GL_CURRENT_COLOR, SER->CurCol);
501    glGetIntegerv(GL_COLOR_MATERIAL_PARAMETER, &(SER->ColMatParam));
502    glGetIntegerv(GL_COLOR_MATERIAL_FACE, &(SER->ColMatFace));
503 
504    SUMA_RETURNe;
505 }
506 
SUMA_RestoreEnablingState(SUMA_EnablingRecord * SER)507 void SUMA_RestoreEnablingState(SUMA_EnablingRecord *SER)
508 {
509    static char FuncName[]={"SUMA_RestoreEnablingState"};
510 
511    SUMA_ENTRY;
512    if (!SER) {
513       SUMA_S_Err("No pointer amigo");
514       SUMA_RETURNe;
515    }
516    if (SER->ALPHA_TEST) glEnable(GL_ALPHA_TEST);
517    else glDisable(GL_ALPHA_TEST);
518    if (SER->DEPTH_TEST) glEnable(GL_DEPTH_TEST);
519    else glDisable(GL_DEPTH_TEST);
520    if (SER->TEXTURE_3D_EXT) glEnable(GL_TEXTURE_3D_EXT);
521    else glDisable(GL_TEXTURE_3D_EXT);
522    if (SER->TEXTURE_3D) glEnable(GL_TEXTURE_3D);
523    else glDisable(GL_TEXTURE_3D);
524    if (SER->TEXTURE_2D) glEnable(GL_TEXTURE_2D);
525    else glDisable(GL_TEXTURE_2D);
526    if (SER->TEXTURE_GEN_S) glEnable(GL_TEXTURE_GEN_S);
527    else glDisable(GL_TEXTURE_GEN_S);
528    if (SER->TEXTURE_GEN_T) glEnable(GL_TEXTURE_GEN_T);
529    else glDisable(GL_TEXTURE_GEN_T);
530    if (SER->TEXTURE_GEN_R) glEnable(GL_TEXTURE_GEN_R);
531    else glDisable(GL_TEXTURE_GEN_R);
532    if (SER->CLIP_PLANE0) glEnable(GL_CLIP_PLANE0);
533    else glDisable(GL_CLIP_PLANE0);
534    if (SER->CLIP_PLANE1) glEnable(GL_CLIP_PLANE1);
535    else glDisable(GL_CLIP_PLANE1);
536    if (SER->CLIP_PLANE2) glEnable(GL_CLIP_PLANE2);
537    else glDisable(GL_CLIP_PLANE2);
538    if (SER->CLIP_PLANE3) glEnable(GL_CLIP_PLANE3);
539    else glDisable(GL_CLIP_PLANE3);
540    if (SER->CLIP_PLANE4) glEnable(GL_CLIP_PLANE4);
541    else glDisable(GL_CLIP_PLANE4);
542    if (SER->CLIP_PLANE5) glEnable(GL_CLIP_PLANE5);
543    else glDisable(GL_CLIP_PLANE5);
544    if (SER->LIGHTING) glEnable(GL_LIGHTING);
545    else glDisable(GL_LIGHTING);
546    if (SER->LIGHT0) glEnable(GL_LIGHT0);
547    else glDisable(GL_LIGHT0);
548    if (SER->LIGHT1) glEnable(GL_LIGHT1);
549    else glDisable(GL_LIGHT1);
550    if (SER->LIGHT2) glEnable(GL_LIGHT2);
551    else glDisable(GL_LIGHT2);
552    if (SER->BLEND) glEnable(GL_BLEND);
553    else glDisable(GL_BLEND);
554    if (SER->LINE_SMOOTH) glEnable(GL_LINE_SMOOTH);
555    else glDisable(GL_LINE_SMOOTH);
556    if (SER->COLOR_MATERIAL) glEnable(GL_COLOR_MATERIAL);
557    else glDisable(GL_COLOR_MATERIAL);
558 
559    /* For now, do not bother setting colors, etc. */
560 
561    /* if (SER->) glEnable();
562       else glDisable() */
563 
564    SUMA_RETURNe;
565 }
566 
SUMA_EnablingState_Info(SUMA_EnablingRecord * SERu)567 char *SUMA_EnablingState_Info(SUMA_EnablingRecord *SERu)
568 {
569    static char FuncName[]={"SUMA_EnablingState_Info"};
570    char *s=NULL;
571    SUMA_EnablingRecord SERl, *SER;
572    SUMA_STRING *SS=NULL;
573 
574    SUMA_ENTRY;
575 
576    SS = SUMA_StringAppend(NULL, NULL);
577    if (!SERu) {
578       SUMA_RecordEnablingState(&SERl, FuncName);
579       SER = &SERl;
580    }  else {
581       SER = SERu;
582    }
583 
584    SUMA_StringAppend_va(SS,"OpenGL State Record for %s\n", SER->Label);
585    SUMA_StringAppend_va(SS,"% 24s is %s\n",
586                        "GL_ALPHA_TEST", SER->ALPHA_TEST ? "+++":"---");
587    SUMA_StringAppend_va(SS,"% 24s is %s\n",
588                         "GL_DEPTH_TEST",SER->DEPTH_TEST ? "+++":"---");
589    SUMA_StringAppend_va(SS,"% 24s is %s\n",
590                "GL_TEXTURE_3D_EXT", SER->TEXTURE_3D_EXT ? "+++":"---");
591    SUMA_StringAppend_va(SS,"% 24s is %s\n",
592            "GL_TEXTURE_2D", SER->TEXTURE_2D ? "+++":"---");
593    SUMA_StringAppend_va(SS,"% 24s is %s\n",
594            "GL_TEXTURE_3D", SER->TEXTURE_3D ? "+++":"---");
595    SUMA_StringAppend_va(SS,"% 24s is %s\n",
596            "GL_TEXTURE_GEN_S", SER->TEXTURE_GEN_S ? "+++":"---");
597    SUMA_StringAppend_va(SS,"% 24s is %s\n",
598            "GL_TEXTURE_GEN_T", SER->TEXTURE_GEN_T ? "+++":"---");
599    SUMA_StringAppend_va(SS,"% 24s is %s\n",
600            "GL_TEXTURE_GEN_R", SER->TEXTURE_GEN_R ? "+++":"---");
601    SUMA_StringAppend_va(SS,"% 24s is %s\n",
602            "GL_CLIP_PLANE0", SER->CLIP_PLANE0 ? "+++":"---");
603    SUMA_StringAppend_va(SS,"% 24s is %s\n",
604            "GL_CLIP_PLANE1", SER->CLIP_PLANE1 ? "+++":"---");
605    SUMA_StringAppend_va(SS,"% 24s is %s\n",
606            "GL_CLIP_PLANE2", SER->CLIP_PLANE2 ? "+++":"---");
607    SUMA_StringAppend_va(SS,"% 24s is %s\n",
608            "GL_CLIP_PLANE3", SER->CLIP_PLANE3 ? "+++":"---");
609    SUMA_StringAppend_va(SS,"% 24s is %s\n",
610            "GL_CLIP_PLANE4", SER->CLIP_PLANE4 ? "+++":"---");
611    SUMA_StringAppend_va(SS,"% 24s is %s\n",
612            "GL_CLIP_PLANE5", SER->CLIP_PLANE5 ? "+++":"---");
613    SUMA_StringAppend_va(SS,"% 24s is %s\n",
614            "GL_LIGHTING", SER->LIGHTING ? "+++":"---");
615    SUMA_StringAppend_va(SS,"% 24s is %s\n",
616            "GL_COLOR_MATERIAL", SER->COLOR_MATERIAL ? "+++":"---");
617    SUMA_StringAppend_va(SS,"% 24s is %d\n",
618            "COLOR_MATERIAL_PARAMETER", SER->ColMatParam);
619    SUMA_StringAppend_va(SS,"% 24s is %d\n",
620            "COLOR_MATERIAL_FACE", SER->ColMatFace);
621    SUMA_StringAppend_va(SS,"% 24s is %.3f %.3f %.3f %.3f\n",
622            "CURRENT_COLOR",
623            SER->CurCol[0], SER->CurCol[1], SER->CurCol[2], SER->CurCol[3] );
624    SUMA_StringAppend_va(SS,"% 24s is %s\n",
625            "GL_LIGHT0", SER->LIGHT0 ? "+++":"---");
626    SUMA_StringAppend_va(SS,"% 24s is %s\n",
627            "GL_LIGHT1", SER->LIGHT1 ? "+++":"---");
628    SUMA_StringAppend_va(SS,"% 24s is %s\n",
629            "GL_LIGHT2", SER->LIGHT2 ? "+++":"---");
630    SUMA_StringAppend_va(SS,"% 24s is %s\n",
631            "GL_BLEND", SER->BLEND ? "+++":"---");
632    SUMA_StringAppend_va(SS,"% 24s is %s\n",
633            "GL_LINE_SMOOTH", SER->LINE_SMOOTH ? "+++":"---");
634 
635 /*
636    SUMA_StringAppend_va(SS,"% 24s is %s\n",
637            "GL_ ", SER-> ? "+++":"---");
638                         */
639    SUMA_SS2S(SS,s);
640 
641    SUMA_RETURN(s);
642 }
643 
SUMA_ShowEnablingState(SUMA_EnablingRecord * SER,FILE * out,char * preamble)644 void SUMA_ShowEnablingState(SUMA_EnablingRecord *SER, FILE *out,
645                             char *preamble) {
646    static char FuncName[]={"SUMA_ShowEnablingState"};
647    char *s=NULL;
648    SUMA_ENTRY;
649    if (!out) out = SUMA_STDOUT;
650 
651    s = SUMA_EnablingState_Info(SER);
652 
653    fprintf(out, "%s%s", preamble ? preamble:"", s);
654 
655    SUMA_free(s); s = NULL;
656 
657    SUMA_RETURNe;
658 }
659 
SUMA_DiffEnablingState_Info(SUMA_EnablingRecord * SERnew,SUMA_EnablingRecord * SERref)660 char *SUMA_DiffEnablingState_Info(SUMA_EnablingRecord *SERnew,
661                                   SUMA_EnablingRecord *SERref)
662 {
663    static char FuncName[]={"SUMA_DiffEnablingState_Info"};
664    char *s=NULL;
665    static SUMA_EnablingRecord SER, *SER_last=NULL;
666    SUMA_EnablingRecord now;
667    SUMA_STRING *SS=NULL;
668 
669    SUMA_ENTRY;
670 
671    if (!SERref) { /* Use last record */
672       if (!SER_last) { /* No last record */
673          SUMA_RecordEnablingState(&SER, "From Diff");
674          SER_last = &SER;
675       }
676       SERref = SER_last;
677    }
678    if (!SERnew) { /* Nothing given get current */
679       SUMA_RecordEnablingState(&now, "From Diff");
680       SERnew = &now;
681    }
682 
683    SS = SUMA_StringAppend(NULL, NULL);
684    SUMA_StringAppend_va(SS,"OpenGL State Diff: %s vs. %s\n",
685                         SERnew->Label, SERref->Label);
686    if (SERnew->ALPHA_TEST != SERref->ALPHA_TEST) {
687       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
688                               "GL_ALPHA_TEST",
689                               SERnew->ALPHA_TEST , SERref->ALPHA_TEST);
690    }
691    if (SERnew->DEPTH_TEST != SERref->DEPTH_TEST) {
692       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
693                               "GL_DEPTH_TEST",
694                               SERnew->DEPTH_TEST , SERref->DEPTH_TEST);
695    }
696    if (SERnew->TEXTURE_3D_EXT != SERref->TEXTURE_3D_EXT) {
697       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
698                               "GL_TEXTURE_3D_EXT",
699                               SERnew->TEXTURE_3D_EXT , SERref->TEXTURE_3D_EXT);
700    }
701    if (SERnew->TEXTURE_2D != SERref->TEXTURE_2D) {
702       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
703                               "GL_TEXTURE_2D",
704                               SERnew->TEXTURE_2D , SERref->TEXTURE_2D);
705    }
706    if (SERnew->TEXTURE_3D != SERref->TEXTURE_3D) {
707       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
708                               "GL_TEXTURE_3D",
709                               SERnew->TEXTURE_3D , SERref->TEXTURE_3D);
710    }
711    if (SERnew->TEXTURE_GEN_S != SERref->TEXTURE_GEN_S) {
712       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
713                               "GL_TEXTURE_GEN_S",
714                               SERnew->TEXTURE_GEN_S , SERref->TEXTURE_GEN_S);
715    }
716    if (SERnew->TEXTURE_GEN_T != SERref->TEXTURE_GEN_T) {
717       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
718                               "GL_TEXTURE_GEN_T",
719                               SERnew->TEXTURE_GEN_T , SERref->TEXTURE_GEN_T);
720    }
721    if (SERnew->TEXTURE_GEN_R != SERref->TEXTURE_GEN_R) {
722       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
723                               "GL_TEXTURE_GEN_R",
724                               SERnew->TEXTURE_GEN_R , SERref->TEXTURE_GEN_R);
725    }
726    if (SERnew->CLIP_PLANE0 != SERref->CLIP_PLANE0) {
727       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
728                               "GL_CLIP_PLANE0",
729                               SERnew->CLIP_PLANE0 , SERref->CLIP_PLANE0);
730    }
731    if (SERnew->CLIP_PLANE1 != SERref->CLIP_PLANE1) {
732       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
733                               "GL_CLIP_PLANE1",
734                               SERnew->CLIP_PLANE1 , SERref->CLIP_PLANE1);
735    }
736    if (SERnew->CLIP_PLANE2 != SERref->CLIP_PLANE2) {
737       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
738                               "GL_CLIP_PLANE2",
739                               SERnew->CLIP_PLANE2 , SERref->CLIP_PLANE2);
740    }
741    if (SERnew->ALPHA_TEST != SERref->ALPHA_TEST) {
742       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
743                               "GL_ALPHA_TEST",
744                               SERnew->ALPHA_TEST , SERref->ALPHA_TEST);
745    }
746    if (SERnew->CLIP_PLANE4 != SERref->CLIP_PLANE4) {
747       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
748                               "GL_CLIP_PLANE4",
749                               SERnew->CLIP_PLANE4 , SERref->CLIP_PLANE4);
750    }
751    if (SERnew->CLIP_PLANE5 != SERref->CLIP_PLANE5) {
752       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
753                               "GL_CLIP_PLANE5",
754                               SERnew->CLIP_PLANE5 , SERref->CLIP_PLANE5);
755    }
756    if (SERnew->LIGHTING != SERref->LIGHTING) {
757       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
758                               "GL_LIGHTING",
759                               SERnew->LIGHTING , SERref->LIGHTING);
760    }
761    if (SERnew->COLOR_MATERIAL != SERref->COLOR_MATERIAL) {
762       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
763                               "GL_COLOR_MATERIAL",
764                               SERnew->COLOR_MATERIAL , SERref->COLOR_MATERIAL);
765    }
766    if (SERnew->CurCol[0] != SERnew->CurCol[0] ||
767        SERnew->CurCol[1] != SERnew->CurCol[1] ||
768        SERnew->CurCol[2] != SERnew->CurCol[2] ||
769        SERnew->CurCol[3] != SERnew->CurCol[3]) {
770        SUMA_StringAppend_va(SS,
771                   "% 24s is %.3f %.3f %.3f %.3f vs %.3f %.3f %.3f %.3f\n",
772                   "CURRENT_COL",
773                   SERnew->CurCol[0], SERnew->CurCol[1],
774                   SERnew->CurCol[2], SERnew->CurCol[3],
775                   SERref->CurCol[0], SERref->CurCol[1],
776                   SERref->CurCol[2], SERref->CurCol[3] );
777    }
778    if (SERnew->ColMatParam != SERref->ColMatParam) {
779       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
780                               "COLOR_MATERIAL_PARAMETER",
781                               SERnew->ColMatParam , SERref->ColMatParam);
782    }
783    if (SERnew->ColMatFace != SERref->ColMatFace) {
784       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
785                               "COLOR_MATERIAL_FACE",
786                               SERnew->ColMatFace , SERref->ColMatFace);
787    }
788    if (SERnew->LIGHT0 != SERref->LIGHT0) {
789       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
790                               "LIGHT0",
791                               SERnew->LIGHT0 , SERref->LIGHT0);
792    }
793    if (SERnew->LIGHT1 != SERref->LIGHT1) {
794       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
795                               "LIGHT1",
796                               SERnew->LIGHT1 , SERref->LIGHT1);
797    }
798    if (SERnew->LIGHT2 != SERref->LIGHT2) {
799       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
800                               "LIGHT2",
801                               SERnew->LIGHT2 , SERref->LIGHT2);
802    }
803    if (SERnew->BLEND != SERref->BLEND) {
804       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
805                               "BLEND",
806                               SERnew->BLEND , SERref->BLEND);
807    }
808    if (SERnew->LINE_SMOOTH != SERref->LINE_SMOOTH) {
809       SUMA_StringAppend_va(SS,"% 24s is %d vs %d\n",
810                               "LINE_SMOOTH",
811                               SERnew->LINE_SMOOTH , SERref->LINE_SMOOTH);
812    }
813 
814 /*
815    SUMA_StringAppend_va(SS,"% 24s is %s\n",
816            "GL_ ", SER-> ? "+++":"---");
817                         */
818    SUMA_StringAppend_va(SS,"End of Diff.\n\n");
819    SUMA_SS2S(SS,s);
820 
821    /* Keep track of last visited */
822    SUMA_CopyEnablingState(SER_last, SERnew);
823    SUMA_RETURN(s);
824 }
825 
SUMA_DiffEnablingState(SUMA_EnablingRecord * SERnew,SUMA_EnablingRecord * SERref,FILE * out,char * preamble)826 void SUMA_DiffEnablingState(SUMA_EnablingRecord *SERnew,
827                             SUMA_EnablingRecord *SERref, FILE *out,
828                             char *preamble) {
829    static char FuncName[]={"SUMA_DiffEnablingState"};
830    char *s=NULL;
831    SUMA_ENTRY;
832    if (!out) out = SUMA_STDOUT;
833 
834    s = SUMA_DiffEnablingState_Info(SERnew, SERref);
835 
836    fprintf(out, "%s%s", preamble ? preamble:"", s);
837 
838    SUMA_free(s); s = NULL;
839 
840    SUMA_RETURNe;
841 }
842 
SUMA_CopyEnablingState(SUMA_EnablingRecord * SERnew,SUMA_EnablingRecord * SERref)843 int SUMA_CopyEnablingState(SUMA_EnablingRecord *SERnew,
844                            SUMA_EnablingRecord *SERref)
845 {
846    static char FuncName[]={"SUMA_CopyEnablingState"};
847 
848    SUMA_ENTRY;
849 
850    if (!SERnew || !SERref) SUMA_RETURN(NOPE);
851 
852    strcpy(SERnew->Label,SERref->Label);
853    SERnew->ALPHA_TEST = SERref->ALPHA_TEST ;
854    SERnew->DEPTH_TEST = SERref->DEPTH_TEST ;
855    SERnew->COLOR_MATERIAL = SERref->COLOR_MATERIAL ;
856    SERnew->TEXTURE_2D = SERref->TEXTURE_2D ;
857    SERnew->TEXTURE_3D = SERref->TEXTURE_3D ;
858    SERnew->TEXTURE_3D_EXT = SERref->TEXTURE_3D_EXT ;
859    SERnew->TEXTURE_GEN_S = SERref->TEXTURE_GEN_S ;
860    SERnew->TEXTURE_GEN_T = SERref->TEXTURE_GEN_T ;
861    SERnew->TEXTURE_GEN_R = SERref->TEXTURE_GEN_R ;
862    SERnew->CLIP_PLANE0 = SERref->CLIP_PLANE0 ;
863    SERnew->CLIP_PLANE1 = SERref->CLIP_PLANE1 ;
864    SERnew->CLIP_PLANE2 = SERref->CLIP_PLANE2 ;
865    SERnew->CLIP_PLANE3 = SERref->CLIP_PLANE3 ;
866    SERnew->CLIP_PLANE4 = SERref->CLIP_PLANE4 ;
867    SERnew->CLIP_PLANE5 = SERref->CLIP_PLANE5 ;
868    SERnew->LIGHTING = SERref->LIGHTING ;
869    SERnew->LIGHT0 = SERref->LIGHT0 ;
870    SERnew->LIGHT1 = SERref->LIGHT1 ;
871    SERnew->LIGHT2 = SERref->LIGHT2 ;
872    SERnew->BLEND = SERref->BLEND ;
873    SERnew->LINE_SMOOTH = SERref->LINE_SMOOTH ;
874    SUMA_COPY_VEC(SERref->CurCol,SERnew->CurCol,4,GLfloat, GLfloat);
875 
876    SUMA_RETURN(YUP);
877 }
878 
SUMA_dset_slice_corners(int slc,float * orig,float * del,int * nvox,float * corners)879 void SUMA_dset_slice_corners( int slc, float *orig, float *del,
880                               int *nvox, float *corners)
881 {
882    static char FuncName[]={"SUMA_dset_slice_corners"};
883    int kk=0;
884    SUMA_ENTRY;
885 
886    corners[kk]  = orig[0] + 0       * del[0]; ++kk;
887    corners[kk]  = orig[1] + 0       * del[1]; ++kk;
888    corners[kk]  = orig[2] + slc     * del[2]; ++kk;
889 
890    corners[kk]  = orig[0] + nvox[0] * del[0]; ++kk;
891    corners[kk]  = orig[1] + 0       * del[1]; ++kk;
892    corners[kk]  = orig[2] + slc     * del[2]; ++kk;
893 
894    corners[kk]  = orig[0] + nvox[0] * del[0]; ++kk;
895    corners[kk]  = orig[1] + nvox[1] * del[1]; ++kk;
896    corners[kk]  = orig[2] + slc     * del[2]; ++kk;
897 
898    corners[kk]  = orig[0] + 0       * del[0]; ++kk;
899    corners[kk]  = orig[1] + nvox[1] * del[1]; ++kk;
900    corners[kk]  = orig[2] + slc     * del[2]; ++kk;
901 
902    SUMA_RETURNe;
903 }
904 
SUMA_LoadVolDO(char * fname,SUMA_DO_CoordUnits coord_type,SUMA_VolumeObject ** VOp,byte PutVOinList)905 SUMA_Boolean SUMA_LoadVolDO (char *fname,
906                         SUMA_DO_CoordUnits coord_type, SUMA_VolumeObject **VOp,
907 			byte PutVOinList)
908 {
909    static char FuncName[]={"SUMA_LoadVolDO"};
910    SUMA_VolumeObject *VO=NULL;
911    THD_3dim_dataset *dset=NULL;
912 
913    SUMA_ENTRY;
914 
915    if (!fname) SUMA_RETURN(NOPE);
916    if (coord_type != SUMA_NORM_SCREEN_UNIT &&
917        coord_type != SUMA_WORLD) coord_type = SUMA_WORLD;
918 
919    if (!(dset = THD_open_dataset( fname ))) {
920       SUMA_S_Errv("Failed to open %s\n", fname);
921       SUMA_free(fname); fname = NULL;
922       SUMA_RETURN(NOPE);
923    }
924 
925    /* Create a DO from dset */
926    if (VOp) {
927       if (*VOp == NULL) {
928          if (!(VO = SUMA_CreateVolumeObject(fname))) {
929             SUMA_S_Err("Failed to create volume object");
930             if (dset) DSET_delete(dset);
931             SUMA_free(fname); fname = NULL;
932             SUMA_RETURN(NOPE);
933          }
934          *VOp = VO;
935       } else {
936          VO = *VOp;
937       }
938    } else {
939       if (!(VO = SUMA_CreateVolumeObject(fname))) {
940          SUMA_S_Err("Failed to create volume object");
941          if (dset) DSET_delete(dset);
942          SUMA_free(fname); fname = NULL;
943          SUMA_RETURN(NOPE);
944       }
945    }
946    /* put main dset into VO */
947    if (!(SUMA_AddDsetVolumeObject(VO, &dset))) {
948       SUMA_S_Err("Failed to add volume");
949       if (dset) DSET_delete(dset);
950       if (VO_N_VOLS(VO) == 0) {
951          VO = SUMA_FreeVolumeObject(VO);
952          if (VOp) *VOp = NULL;
953       }
954       SUMA_free(fname); fname = NULL;
955       SUMA_RETURN(NOPE);
956    }
957 
958    /* Set some display defaults */
959    {
960       SUMA_X_SurfCont *SurfCont = NULL;
961       SUMA_VOL_SAUX *VSaux=NULL;
962       VO->TexEnvMode = GL_REPLACE;
963       if ((SurfCont = SUMA_ADO_Cont((SUMA_ALL_DO *)VO)) &&
964           (VSaux = SUMA_ADO_VSaux((SUMA_ALL_DO *)VO))) {
965          /* Do the defaults, then modify per env */
966          VSaux->ShowAxSlc = 1;
967          SurfCont->Ax_slc->slice_num = (int)(SUMA_VO_N_Slices(VO, "Ax")/2.0);
968          SurfCont->Ax_slc->mont_inc = 1;
969 
970          VSaux->ShowSaSlc = 1;
971          SurfCont->Sa_slc->slice_num = (int)(SUMA_VO_N_Slices(VO, "Sa")/2.0);
972          SurfCont->Sa_slc->mont_num = 2;
973          SurfCont->Sa_slc->mont_inc =
974                         (int)SUMA_MAX_PAIR(SurfCont->Sa_slc->slice_num/2,1);
975          VSaux->ShowCoSlc = 0;
976          SurfCont->Co_slc->slice_num = (int)(SUMA_VO_N_Slices(VO, "Co")/2.0);
977 
978          VSaux->ShowVrSlc = 0;
979          VSaux->VrSelect = 0;
980          VSaux->SlicesAtCrosshair = 0;
981          /* Maybe params froms the env? */
982          SUMA_Set_VO_Slice_Params(SUMA_EnvVal("SUMA_VO_InitSlices"), VO);
983          /* And for selectable VR */
984          {
985             char *eee = getenv("SUMA_VrSelectable");
986             if (eee) {
987                if (strcmp(eee,"NO") == 0)  VSaux->VrSelect = NOPE;
988                else if (strcmp(eee,"YES") == 0) VSaux->VrSelect = YUP;
989                else {
990                   fprintf (SUMA_STDERR,
991                            "Warning %s:\n"
992                         "Bad value for environment variable SUMA_VrSelectable\n"
993                            "Assuming default of YES", FuncName);
994                   VSaux->VrSelect = YUP;
995                }
996             } else VSaux->VrSelect = YUP;
997          }
998       } else {
999          SUMA_S_Err("Failed to initialize volume display");
1000       }
1001    }
1002    if (PutVOinList) {
1003       /* Add VO into DO list */
1004       if (!SUMA_AddDO(SUMAg_DOv, &(SUMAg_N_DOv), (void *)VO,
1005                 	VO_type, coord_type)) {
1006 	 fprintf(SUMA_STDERR,"Error %s: Error Adding DO\n", FuncName);
1007 	 SUMA_RETURN(NOPE);
1008       }
1009    }
1010    SUMA_RETURN(YUP);
1011 }
1012 
SUMA_Set_VO_Slice_Params(char * params,SUMA_VolumeObject * VO)1013 int SUMA_Set_VO_Slice_Params(char *params, SUMA_VolumeObject *VO)
1014 {
1015    static char FuncName[]={"SUMA_Set_VO_Slice_Params"};
1016    NI_str_array *sar=NULL, *sub=NULL;
1017    int err=0, val, kk;
1018    float fval;
1019    SUMA_X_SurfCont *SurfCont = NULL;
1020    SUMA_VOL_SAUX *VSaux=NULL;
1021 
1022    SUMA_ENTRY;
1023    if (!params || params[0]=='\0') SUMA_RETURN(1);
1024    if (!VO) { SUMA_S_Err("NO VO"); SUMA_RETURN(0); }
1025    if (!(SurfCont = SUMA_ADO_Cont((SUMA_ALL_DO *)VO)) ||
1026        !(VSaux = SUMA_ADO_VSaux((SUMA_ALL_DO *)VO))) {
1027       SUMA_S_Err("Too early for this!");
1028       SUMA_RETURN(0);
1029    }
1030    if (!(sar = SUMA_NI_decode_string_list( params , ",; "))) {
1031       SUMA_S_Err("Huh?"); SUMA_RETURN(0);
1032    }
1033 
1034    /* Each string should be the form: variant:sn:num:inc */
1035    err = 0;
1036    VSaux->ShowAxSlc = 0;
1037    VSaux->ShowSaSlc = 0;
1038    VSaux->ShowCoSlc = 0;
1039    VSaux->ShowVrSlc = 0;
1040    for (kk=0; kk<sar->num && !err; ++kk) {
1041       if ((sub = SUMA_NI_decode_string_list( sar->str[kk] , ":")) &&
1042           (sub->num > 0)) {
1043                if (!strcasecmp(sub->str[0], "Ax") ||
1044                    !strcasecmp(sub->str[0], "hAx")) {
1045                   if (sub->str[0][0] == 'h') VSaux->ShowAxSlc = 0;
1046                   else VSaux->ShowAxSlc = 1;
1047                   if (sub->num > 1) {
1048                      fval = (float)strtod(sub->str[1], NULL);
1049                      if (fval > 0.0 && fval < 1.0) {
1050                         val = fval * (SUMA_VO_N_Slices(VO, "Ax")-1);
1051                      } else val = (int)fval;
1052                      if (val >= 0 &&
1053                          val < SUMA_VO_N_Slices(VO, "Ax")) {
1054                         SurfCont->Ax_slc->slice_num = val;
1055                      }
1056                   }
1057                   if (sub->num > 2) {
1058                      val = (int)strtod(sub->str[2], NULL);
1059                      if (val > 0 &&
1060                          val < SUMA_VO_N_Slices(VO, "Ax")) {
1061                          SurfCont->Ax_slc->mont_num = val;
1062                      }
1063                   }
1064                   if (sub->num > 3) {
1065                      fval = (float)strtod(sub->str[3], NULL);
1066                      if (fval > 0.0 && fval < 1.0) {
1067                         val = fval * (SUMA_VO_N_Slices(VO, "Ax")-1);
1068                      } else val = (int)fval;
1069                      if (val > 0 &&
1070                          val < SUMA_VO_N_Slices(VO, "Ax")) {
1071                          SurfCont->Ax_slc->mont_inc = val;
1072                      }
1073                   }
1074          } else if (!strcasecmp(sub->str[0], "Sa") ||
1075                     !strcasecmp(sub->str[0], "hSa")) {
1076                   if (sub->str[0][0] == 'h') VSaux->ShowSaSlc = 0;
1077                   else VSaux->ShowSaSlc = 1;
1078                   if (sub->num > 1) {
1079                      fval = (float)strtod(sub->str[1], NULL);
1080                      if (fval > 0.0 && fval < 1.0) {
1081                         val = fval * (SUMA_VO_N_Slices(VO, "Sa")-1);
1082                      } else val = (int)fval;
1083                      if (val >= 0 &&
1084                          val < SUMA_VO_N_Slices(VO, "Sa")) {
1085                         SurfCont->Sa_slc->slice_num = val;
1086                      }
1087                   }
1088                   if (sub->num > 2) {
1089                      val = (int)strtod(sub->str[2], NULL);
1090                      if (val > 0 &&
1091                          val < SUMA_VO_N_Slices(VO, "Sa")) {
1092                          SurfCont->Sa_slc->mont_num = val;
1093                      }
1094                   }
1095                   if (sub->num > 3) {
1096                      fval = (float)strtod(sub->str[3], NULL);
1097                      if (fval > 0.0 && fval < 1.0) {
1098                         val = fval * (SUMA_VO_N_Slices(VO, "Sa")-1);
1099                      } else val = (int)fval;
1100                      if (val > 0 &&
1101                          val < SUMA_VO_N_Slices(VO, "Sa")) {
1102                          SurfCont->Sa_slc->mont_inc = val;
1103                      }
1104                   }
1105          } else if (!strcasecmp(sub->str[0], "Co") ||
1106                     !strcasecmp(sub->str[0], "hCo")) {
1107                   if (sub->str[0][0] == 'h') VSaux->ShowCoSlc = 0;
1108                   else VSaux->ShowCoSlc = 1;
1109                   if (sub->num > 1) {
1110                      fval = (float)strtod(sub->str[1], NULL);
1111                      if (fval > 0.0 && fval < 1.0) {
1112                         val = fval * (SUMA_VO_N_Slices(VO, "Co")-1);
1113                      } else val = (int)fval;
1114                      if (val >= 0 &&
1115                          val < SUMA_VO_N_Slices(VO, "Co")) {
1116                         SurfCont->Co_slc->slice_num = val;
1117                      }
1118                   }
1119                   if (sub->num > 2) {
1120                      val = (int)strtod(sub->str[2], NULL);
1121                      if (val > 0 &&
1122                          val < SUMA_VO_N_Slices(VO, "Co")) {
1123                          SurfCont->Co_slc->mont_num = val;
1124                      }
1125                   }
1126                   if (sub->num > 3) {
1127                      fval = (float)strtod(sub->str[3], NULL);
1128                      if (fval > 0.0 && fval < 1.0) {
1129                         val = fval * (SUMA_VO_N_Slices(VO, "Co")-1);
1130                      } else val = (int)fval;
1131                      if (val > 0 &&
1132                          val < SUMA_VO_N_Slices(VO, "Co")) {
1133                          SurfCont->Co_slc->mont_inc = val;
1134                      }
1135                   }
1136          } else if (!strcasecmp(sub->str[0], "Vr") ||
1137                     !strcasecmp(sub->str[0], "hVr")) {
1138                   if (sub->str[0][0] == 'h') VSaux->ShowVrSlc = 0;
1139                   else VSaux->ShowVrSlc = 1;
1140          } else {
1141             SUMA_S_Err(
1142       "Slice variant %s not recognized for env SUMA_VO_InitSlices.\n"
1143       "Defaults will prevail.", sub->str[0]);
1144             err = 1;
1145             /* Just put the 'show' flags back where they were */
1146             VSaux->ShowAxSlc = 1;
1147             VSaux->ShowSaSlc = 1;
1148             VSaux->ShowCoSlc = 0;
1149             VSaux->ShowVrSlc = 0;
1150          }
1151 
1152       }
1153       sub = SUMA_free_NI_str_array(sub);
1154    }
1155    sar = SUMA_free_NI_str_array(sar);
1156 
1157    /* You must have something showing for now because otherwise you can't open
1158    a volume controller! */
1159    if (!VSaux->ShowAxSlc && !VSaux->ShowSaSlc &&
1160        !VSaux->ShowCoSlc && !VSaux->ShowVrSlc ) {
1161       SUMA_S_Note("For now, must force something to show");
1162       VSaux->ShowAxSlc = 1;
1163    }
1164 
1165    SUMA_RETURN(1);
1166 }
1167 
SUMA_Load3DTextureNIDOnel(NI_element * nel,SUMA_DO_CoordUnits defaultcoordtype)1168 SUMA_Boolean SUMA_Load3DTextureNIDOnel (NI_element *nel,
1169                                     SUMA_DO_CoordUnits defaultcoordtype)
1170 {
1171    static char FuncName[]={"SUMA_Load3DTextureNIDOnel"};
1172    char *fname=NULL;
1173    char *atr=NULL, stmp[128];
1174    int i = 0;
1175    SUMA_VolumeObject *VO=NULL;
1176    SUMA_DO_CoordUnits coord_type=SUMA_WORLD;
1177    SUMA_Boolean LocalHead = NOPE;
1178 
1179    SUMA_ENTRY;
1180 
1181    if (  !nel ||
1182          strcmp(nel->name,"3DTex")   )
1183       SUMA_RETURN(NOPE);
1184 
1185    if (NI_IS_STR_ATTR_EQUAL(nel, "read_status", "read")) SUMA_RETURN(YUP);
1186 
1187    NI_set_attribute(nel,"read_status","fail");
1188 
1189    if (! (fname = SUMA_copy_string(NI_get_attribute(nel,"filename"))) )
1190       SUMA_RETURN(NOPE);
1191    if (!SUMA_search_file(&fname, NULL)) { /* can't find it */
1192       SUMA_LH("Failed to find %s", fname);
1193       SUMA_free(fname); fname = NULL;
1194       SUMA_RETURN(NOPE);
1195    }
1196 
1197    /* Is there a particular coord type in nel? */
1198    if ((atr = NI_get_attribute(nel, "coord_type"))) {
1199       if ((coord_type = SUMA_CoordType(atr))
1200            == SUMA_COORD_TYPE_ERROR) {
1201          SUMA_S_Errv("Bad coord_type %s,"
1202                      "using default ",
1203                      atr);
1204          coord_type = defaultcoordtype;
1205       }
1206    }else {
1207       coord_type = defaultcoordtype;
1208    }
1209 
1210    if (!(SUMA_LoadVolDO(fname, coord_type, &VO,1))) {
1211       SUMA_S_Err("Failed to read %s", NI_get_attribute(nel,"filename"));
1212       SUMA_ifree(fname);
1213       SUMA_RETURN(NOPE);
1214    }
1215    SUMA_ifree(fname);
1216 
1217    /* any other dsets? */
1218    i = 0; sprintf(stmp,"overlay%d",i);
1219    while ((atr=NI_get_attribute(nel,stmp))) {
1220       SUMA_LHv("Loading %s\n", atr);
1221       if (! (fname = SUMA_copy_string(atr)) )
1222          SUMA_RETURN(NOPE);
1223       if (!SUMA_search_file(&fname, NULL)) { /* can't find it */
1224          SUMA_S_Errv("Failed to find %s\n", fname);
1225          SUMA_free(fname); fname = NULL;
1226          break;
1227       }
1228 
1229       if (!(SUMA_LoadVolDO( fname, coord_type, &VO,1))) {
1230          SUMA_S_Errv("Failed to open %s\n", fname);
1231          SUMA_free(fname); fname = NULL;
1232          break;
1233       }
1234 
1235       SUMA_LHv("Added %s\n", SUMA_VE_Headname(VO->VE, i+1));
1236       SUMA_free(fname); fname = NULL;
1237       ++i;
1238       sprintf(stmp,"overlay%d",i);
1239    }
1240 
1241 
1242    VO->TexEnvMode = SUMA_NIDO_TexEnvMode(nel, GL_REPLACE);
1243 
1244    /* store pointer copy of DO in nel */
1245    NI_SET_PTR(nel, "DO", VO);
1246    /* store idcode_str of DO in nel */
1247    NI_SET_STR(nel,"DO_idcode_str", VO->idcode_str);
1248 
1249    /* mark nel as read */
1250    NI_set_attribute(nel,"read_status","read");
1251 
1252    SUMA_RETURN(YUP);
1253 }
1254 
SUMA_VOof3DTextureNIDOnel(NI_element * nel)1255 SUMA_VolumeObject *SUMA_VOof3DTextureNIDOnel(NI_element *nel)
1256 {
1257    static char FuncName[]={"SUMA_VOof3DTextureNIDOnel"};
1258    SUMA_VolumeObject *VO = NULL;
1259    int ii;
1260    char *idcode_str=NULL;
1261    SUMA_Boolean LocalHead = NOPE;
1262 
1263    SUMA_ENTRY;
1264 
1265    if (!(idcode_str = NI_get_attribute(nel, "DO_idcode_str"))) {
1266       SUMA_S_Err("NULL nel DO_idcode_str");
1267       SUMA_RETURN(NULL);
1268    }
1269 
1270    for (ii=0; ii<SUMAg_N_DOv; ++ii) {
1271       if (SUMA_isVO(SUMAg_DOv[ii])) {
1272          VO = (SUMA_VolumeObject *)(SUMAg_DOv[ii].OP);
1273          if (!strcmp(VO->idcode_str, idcode_str)) SUMA_RETURN(VO);
1274       }
1275    }
1276 
1277    SUMA_LHv("DO for nel %s, %s, not found\n",
1278       nel->name, idcode_str);
1279 
1280    SUMA_RETURN(NULL);
1281 }
1282 
SUMA_3DTextureNIDOnelofVO(SUMA_VolumeObject * VO)1283 NI_element *SUMA_3DTextureNIDOnelofVO(SUMA_VolumeObject *VO)
1284 {
1285    static char FuncName[]={"SUMA_3DTextureNIDOnelofVO"};
1286    NI_element *nel = NULL;
1287 
1288    SUMA_ENTRY;
1289 
1290    SUMA_S_Err("Sorry, not implemented yet");
1291 
1292    SUMA_RETURN(nel);
1293 }
1294 
1295 /*!
1296    centers of first and last voxels in volume in RAI
1297 */
SUMA_dset_extreme_corners(SUMA_DSET * dset,float * mincorner,float * maxcorner,int voxcen)1298 void SUMA_dset_extreme_corners( SUMA_DSET *dset,
1299                                 float * mincorner, float *maxcorner,
1300                                 int voxcen)
1301 {
1302    static char FuncName[]={"SUMA_dset_extreme_corners"};
1303    float A[4][4], I[3], *v;
1304    int *dims;
1305    SUMA_Boolean LocalHead = NOPE;
1306 
1307    SUMA_ENTRY;
1308 
1309 
1310    if (mincorner) mincorner[0] = mincorner[1] = mincorner[2] = 0.0;
1311    if (maxcorner) maxcorner[0] = maxcorner[1] = maxcorner[2] = 0.0;
1312 
1313    if (!dset || !(v=SUMA_GetDatasetI2X(dset, A)) ||
1314        !(dims = SUMA_GetDatasetDimensions(dset)) ) {
1315       SUMA_S_Err("no valid ijk_to_dicom_real") ;
1316       SUMA_RETURNe;
1317    }
1318 
1319    if (mincorner) {
1320       if (voxcen) {
1321          mincorner[0] = A[0][3];
1322          mincorner[1] = A[1][3];
1323          mincorner[2] = A[2][3];
1324       } else {
1325          I[0] = I[1] = I[2] = -0.5;
1326          AFF44_MULT_I(mincorner, A, I);
1327       }
1328    }
1329    if (maxcorner) {
1330       if (voxcen) {
1331          I[0] = dims[0]-1;
1332          I[1] = dims[1]-1;
1333          I[2] = dims[2]-1;
1334          AFF44_MULT_I(maxcorner, A, I);
1335       } else {
1336          I[0] = dims[2]-0.5;
1337          I[1] = dims[1]-0.5;
1338          I[2] = dims[2]-0.5;
1339          AFF44_MULT_I(maxcorner, A, I);
1340       }
1341    }
1342 
1343    SUMA_RETURNe;
1344 }
1345 
1346 /*!
1347    Corners of box defining volume boundaries
1348 
1349    0,    0,    0
1350    Nx-1, 0,    0
1351    Nx-1, Ny-1, 0
1352    0,    Ny-1, 0
1353    0,    0,    Nz-1
1354    Nx-1, 0,    Nz-1
1355    Nx-1, Ny-1, Nz-1
1356    0,    Ny-1, Nz-1
1357 
1358 */
SUMA_dset_box_corners(SUMA_DSET * dset,float * corners,int voxcen)1359 SUMA_Boolean SUMA_dset_box_corners( SUMA_DSET *dset,
1360                                     float * corners, int voxcen)
1361 {
1362    static char FuncName[]={"SUMA_dset_box_corners"};
1363    float A[4][4], I[24], *v, *vi;
1364    int *dims, i;
1365    SUMA_Boolean LocalHead = NOPE;
1366 
1367    SUMA_ENTRY;
1368 
1369    if (!dset || !(v=SUMA_GetDatasetI2X(dset, A)) ||
1370        !(dims = SUMA_GetDatasetDimensions(dset)) ) {
1371       SUMA_S_Err("no valid ijk_to_dicom_real") ;
1372       SUMA_RETURN(NOPE);
1373    }
1374    if (!corners) {
1375       SUMA_S_Err("No return vehicle");
1376       SUMA_RETURN(NOPE);
1377    }
1378 
1379    /* Fill up the IJKs */
1380    i = 0;
1381    I[i++] = 0;         I[i++] = 0;         I[i++] = 0;
1382    I[i++] = dims[0]-1; I[i++] = 0;         I[i++] = 0;
1383    I[i++] = dims[0]-1; I[i++] = dims[1]-1; I[i++] = 0;
1384    I[i++] = 0;         I[i++] = dims[1]-1; I[i++] = 0;
1385    I[i++] = 0;         I[i++] = 0;         I[i++] = dims[2]-1;
1386    I[i++] = dims[0]-1; I[i++] = 0;         I[i++] = dims[2]-1;
1387    I[i++] = dims[0]-1; I[i++] = dims[1]-1; I[i++] = dims[2]-1;
1388    I[i++] = 0;         I[i++] = dims[1]-1; I[i++] = dims[2]-1;
1389 
1390    /* offset for non voxel center? */
1391    if (!voxcen) {
1392       for (i=0; i<24; ++i) {
1393          if (I[i] > 0.0f) I[i] += 0.5;
1394          else I[i] -= 0.5;
1395       }
1396    }
1397 
1398    for (i=0; i<24; i = i +3) {
1399       v = corners+i;
1400       vi = I+i;
1401       AFF44_MULT_I(v, A, vi);
1402    }
1403 
1404    SUMA_RETURN(YUP);
1405 }
1406 
1407 /* Older version, OK for cardinal axes only
1408    Use SUMA_dset_tex_slice_corners() */
SUMA_dset_tex_slice_corners_card(int slci,THD_3dim_dataset * dset,GLfloat * tcorners,GLfloat * corners,int dim,int voxcen)1409 void SUMA_dset_tex_slice_corners_card( int slci, THD_3dim_dataset *dset,
1410                               GLfloat *tcorners, GLfloat *corners,
1411                               int dim, int voxcen)
1412 {
1413    static char FuncName[]={"SUMA_dset_tex_slice_corners_card"};
1414    int kk=0;
1415    float orig[3] = { 0, 0, 0}, del[3] = { 0, 0, 0};
1416    int nvox[3] = { 0, 0, 0};
1417    int slcx, slcy, slcz=0;
1418    SUMA_Boolean LocalHead = NOPE;
1419 
1420    SUMA_ENTRY;
1421 
1422    if (!voxcen) {
1423       SUMA_LH("voxcen == 0 not supported in this old version");
1424    }
1425 
1426    orig[0] = dset->daxes->xxorg ;
1427    orig[1] = dset->daxes->yyorg ;
1428    orig[2] = dset->daxes->zzorg ;
1429    nvox[0] = DSET_NX(dset);
1430    nvox[1] = DSET_NY(dset);
1431    nvox[2] = DSET_NZ(dset);
1432    del[0] = dset->daxes->xxdel ;
1433    del[1] = dset->daxes->yydel ;
1434    del[2] = dset->daxes->zzdel ;
1435 
1436    switch (dim) {
1437       default:
1438          SUMA_S_Err("Bad dim value");
1439          SUMA_RETURNe;
1440       case 2:
1441          kk = 0;
1442     corners[kk] = orig[0] + 0       * del[0];
1443    tcorners[kk] = 0;                            ++kk;
1444     corners[kk] = orig[1] + 0       * del[1];
1445    tcorners[kk] = 0;                            ++kk;
1446     corners[kk] = orig[2] + slci    * del[2];
1447    tcorners[kk] = ((float)slci+0.5)/(float)nvox[2];++kk;
1448 
1449     corners[kk] = orig[0] + (nvox[0]-1) * del[0];
1450    tcorners[kk] = 1;                            ++kk;
1451     corners[kk] = orig[1] + 0       * del[1];
1452    tcorners[kk] = 0;                            ++kk;
1453     corners[kk] = orig[2] + slci    * del[2];
1454    tcorners[kk] = tcorners[2];                  ++kk;
1455 
1456 
1457     corners[kk] = orig[0] + (nvox[0]-1) * del[0];
1458    tcorners[kk] = 1;                            ++kk;
1459     corners[kk] = orig[1] + (nvox[1]-1) * del[1];
1460    tcorners[kk] = 1;                            ++kk;
1461     corners[kk] = orig[2] + slci    * del[2];
1462    tcorners[kk] = tcorners[2];                  ++kk;
1463 
1464     corners[kk] = orig[0] + 0       * del[0];
1465    tcorners[kk] = 0;                            ++kk;
1466     corners[kk] = orig[1] + (nvox[1]-1) * del[1];
1467    tcorners[kk] = 1;                            ++kk;
1468     corners[kk] = orig[2] + slci    * del[2];
1469    tcorners[kk] = tcorners[2];                  ++kk;
1470          break;
1471       case 1:
1472          kk = 0;
1473     corners[kk] = orig[0] + 0       * del[0];
1474    tcorners[kk] = 0;                            ++kk;
1475     corners[kk] = orig[1] + slci    * del[1];
1476    tcorners[kk] = ((float)slci+0.5)/(float)nvox[1];++kk;
1477     corners[kk] = orig[2] + 0       * del[2];
1478    tcorners[kk] = 0;                            ++kk;
1479 
1480     corners[kk] = orig[0] + (nvox[0]-1) * del[0];
1481    tcorners[kk] = 1;                            ++kk;
1482     corners[kk] = orig[1] + slci    * del[1];
1483    tcorners[kk] = tcorners[1];                  ++kk;
1484     corners[kk] = orig[2] + 0       * del[2];
1485    tcorners[kk] = 0;                            ++kk;
1486 
1487 
1488     corners[kk] = orig[0] + (nvox[0]-1) * del[0];
1489    tcorners[kk] = 1;                            ++kk;
1490     corners[kk] = orig[1] + slci    * del[1];
1491    tcorners[kk] = tcorners[1];                  ++kk;
1492     corners[kk] = orig[2] + (nvox[2]-1) * del[2];
1493    tcorners[kk] = 1;                            ++kk;
1494 
1495     corners[kk] = orig[0] + 0       * del[0];
1496    tcorners[kk] = 0;                            ++kk;
1497     corners[kk] = orig[1] + slci    * del[1];
1498    tcorners[kk] = tcorners[1];                  ++kk;
1499     corners[kk] = orig[2] + (nvox[2]-1) * del[2];
1500    tcorners[kk] = 1;                            ++kk;
1501          break;
1502       case 0:
1503          kk = 0;
1504     corners[kk] = orig[0] + slci    * del[0];
1505    tcorners[kk] = ((float)slci+0.5)/(float)nvox[0];++kk;
1506     corners[kk] = orig[1] + 0       * del[1];
1507    tcorners[kk] = 0;                            ++kk;
1508     corners[kk] = orig[2] + 0       * del[2];
1509    tcorners[kk] = 0;                            ++kk;
1510 
1511     corners[kk] = orig[0] + slci    * del[0];
1512    tcorners[kk] = tcorners[0];                  ++kk;
1513     corners[kk] = orig[1] + (nvox[1]-1) * del[1];
1514    tcorners[kk] = 1;                            ++kk;
1515     corners[kk] = orig[2] + 0       * del[2];
1516    tcorners[kk] = 0;                            ++kk;
1517 
1518 
1519     corners[kk] = orig[0] + slci    * del[0];
1520    tcorners[kk] = tcorners[0];                  ++kk;
1521     corners[kk] = orig[1] + (nvox[1]-1) * del[1];
1522    tcorners[kk] = 1;                            ++kk;
1523     corners[kk] = orig[2] + (nvox[2]-1) * del[2];
1524    tcorners[kk] = 1;                            ++kk;
1525 
1526     corners[kk] = orig[0] + slci    * del[0];
1527    tcorners[kk] = tcorners[0];                  ++kk;
1528     corners[kk] = orig[1] + 0       * del[1];
1529    tcorners[kk] = 0;                            ++kk;
1530     corners[kk] = orig[2] + (nvox[2]-1) * del[2];
1531    tcorners[kk] = 1;                            ++kk;
1532          break;
1533    }
1534 
1535    if (LocalHead) {
1536       SUMA_LHv("Slice %d, dim %d corners and textures\n", slci, dim);
1537       for (kk=0; kk<4; ++kk) {
1538          fprintf(SUMA_STDERR,
1539                   "c%d: %.3f   %.3f   %.3f\n"
1540                   "t%d: %.3f   %.3f   %.3f\n",
1541             kk, corners[3*kk],corners[3*kk+1],corners[3*kk+2],
1542             kk, tcorners[3*kk],tcorners[3*kk+1],tcorners[3*kk+2]);
1543 
1544       }
1545       fprintf(SUMA_STDERR,"\n");
1546    }
1547 
1548    SUMA_RETURNe;
1549 }
1550 
1551 /* Get slider value from texture corners,
1552    inverse of SUMA_dset_tex_slice_corners_gui*/
SUMA_dset_gui_slice_from_tex_slice(SUMA_VolumeElement ** VE,int ive,float * PlEq,int voxcen,char * variant,int * slider)1553 int SUMA_dset_gui_slice_from_tex_slice(SUMA_VolumeElement **VE, int ive,
1554                      float *PlEq, int voxcen,
1555                      char *variant,int *slider)
1556 {
1557    static char FuncName[]={"SUMA_dset_gui_slice_from_tex_slice"};
1558    char *orcode;
1559    int dim=0, nslc=0, *dims;
1560    float I[3]={0.0, 0.0, 0.0}, C[3]={0.0, 0.0, 0.0},
1561          Dir0[3] = {1, 0, 0}, Dir1[3] = {0, 1, 0}, Dir2[3] = {0, 0, 1},
1562          dd, d0, d1, d2;
1563    SUMA_DSET *dset=NULL;
1564    SUMA_Boolean LocalHead = NOPE;
1565 
1566    SUMA_ENTRY;
1567 
1568    if (ive < 0) ive = 0;
1569    if (!(dset = SUMA_VE_dset(VE, ive)) || !PlEq ||
1570        !(dims = SUMA_GetDatasetDimensions(dset))) {
1571       SUMA_S_Err("no dset or no variant") ;
1572       SUMA_RETURN(-1);
1573    }
1574 
1575    if (slider) *slider = -1;
1576 
1577    orcode = SUMA_Dset_orcode(dset);
1578    if (orcode[0] == 'X') { SUMA_S_Err("No orcode"); SUMA_RETURN(-1); }
1579 
1580    /* Take the normal and turn it to IJK land */
1581    AFF44_MULT_D(I, VE[ive]->X2I, PlEq);
1582    SUMA_UNITIZE_VEC(I,3);
1583 
1584    /* Find out which dim you're closest to */
1585    d0 = SUMA_MT_DOT(I, Dir0); dd = d0; dim = 0;
1586    d1 = SUMA_MT_DOT(I, Dir1);
1587    d2 = SUMA_MT_DOT(I, Dir2);
1588    if (SUMA_ABS(d1) > SUMA_ABS(dd)) {
1589       dim = 1; dd = d1;
1590    }
1591    if (SUMA_ABS(d2) > SUMA_ABS(dd)) {
1592       dim = 2; dd = d2;
1593    }
1594    SUMA_LH("PlEq: %f %f %f %f\n"
1595             "I  : %f %f %f\n"
1596             "Dots: %f %f %f, orcode %s, dim %d",
1597            PlEq[0], PlEq[1], PlEq[2], PlEq[3],
1598            I[0], I[1], I[2], d0, d1, d2, orcode, dim)
1599    if (variant) {
1600            if (orcode[dim] == 'I' || orcode[dim] == 'S') sprintf(variant,"Ax");
1601       else if (orcode[dim] == 'R' || orcode[dim] == 'L') sprintf(variant,"Sa");
1602       else if (orcode[dim] == 'A' || orcode[dim] == 'P') sprintf(variant,"Co");
1603    }
1604 
1605    /* Don't bother which slice number this is, just return the dim */
1606    SUMA_RETURN(dim);
1607 }
1608 
SUMA_dset_gui_slice_from_tex_slice_d(SUMA_VolumeElement ** VE,int ive,double * PlEq,int voxcen,char * variant,int * slider)1609 int SUMA_dset_gui_slice_from_tex_slice_d(SUMA_VolumeElement **VE, int ive,
1610                      double *PlEq, int voxcen,
1611                      char *variant,int *slider)
1612 {
1613    static char FuncName[]={"SUMA_dset_gui_slice_from_tex_slice_d"};
1614    float fv[4];
1615    if (!PlEq) return(-1);
1616    fv[0] = PlEq[0];
1617    fv[1] = PlEq[1];
1618    fv[2] = PlEq[2];
1619    fv[3] = PlEq[3];
1620    return(SUMA_dset_gui_slice_from_tex_slice(VE, ive, fv,
1621                                              voxcen, variant, slider));
1622 }
1623 
1624 /* Get texture corners from slider values
1625 \sa SUMA_dset_gui_slice_from_tex_slice*/
SUMA_dset_tex_slice_corners_gui(SUMA_VolumeElement ** VE,int ive,char * variant,int slider,GLfloat * tcorners,GLfloat * corners,GLfloat * slc_cen,float * PlEq,int voxcen)1626 int SUMA_dset_tex_slice_corners_gui(SUMA_VolumeElement **VE, int ive,
1627                                              char *variant,int slider,
1628                           GLfloat *tcorners, GLfloat *corners, GLfloat *slc_cen,
1629                           float *PlEq, int voxcen )
1630 {
1631    static char FuncName[]={"SUMA_dset_tex_slice_corners_gui"};
1632    char *orcode;
1633    int dim=0, nslc=0, *dims;
1634    float I[3]={0.0, 0.0, 0.0}, C[3]={0.0, 0.0, 0.0};
1635    SUMA_DSET *dset=NULL;
1636    SUMA_Boolean LocalHead = NOPE;
1637 
1638    SUMA_ENTRY;
1639 
1640    if (ive < 0) ive = 0;
1641    if (!(dset = SUMA_VE_dset(VE, ive)) || !variant||
1642        !(dims = SUMA_GetDatasetDimensions(dset))) {
1643       SUMA_S_Err("no dset or no variant") ;
1644       SUMA_RETURN(-1);
1645    }
1646 
1647    orcode = SUMA_Dset_orcode(dset);
1648    if (orcode[0] == 'X') { SUMA_S_Err("No orcode"); SUMA_RETURN(-1); }
1649 
1650    switch (variant[0]) {
1651       case 'A': /* axial slicing desired*/
1652          if (orcode[0] == 'I' || orcode[0] == 'S') {
1653             dim = 0; nslc = slider;
1654             if (orcode[0] == 'S') nslc = VE[ive]->Ni-1-slider;
1655             if (nslc < 0) nslc = 0;
1656             if (nslc >= VE[ive]->Ni) nslc = VE[ive]->Ni-1;
1657          } else if (orcode[1] == 'I' || orcode[1] == 'S') {
1658             dim = 1; nslc = slider;
1659             if (orcode[1] == 'S') nslc = VE[ive]->Nj-1-slider;
1660             if (nslc < 0) nslc = 0;
1661             if (nslc >= VE[ive]->Nj) nslc = VE[ive]->Nj-1;
1662          } else if (orcode[2] == 'I' || orcode[2] == 'S') {
1663             dim = 2; nslc = slider;
1664             if (orcode[2] == 'S') nslc = VE[ive]->Nk-1-slider;
1665             if (nslc < 0) nslc = 0;
1666             if (nslc >= VE[ive]->Nk) nslc = VE[ive]->Nk-1;
1667          }
1668          break;
1669       case 'S': /* sagittal slicing desired */
1670          if (orcode[0] == 'R' || orcode[0] == 'L') {
1671             dim = 0; nslc = slider;
1672             if (orcode[0] == 'L') nslc = VE[ive]->Ni-1-slider;
1673             if (nslc < 0) nslc = 0;
1674             if (nslc >= VE[ive]->Ni) nslc = VE[ive]->Ni-1;
1675          } else if (orcode[1] == 'R' || orcode[1] == 'L') {
1676             dim = 1; nslc = slider;
1677             if (orcode[1] == 'L') nslc = VE[ive]->Nj-1-slider;
1678             if (nslc < 0) nslc = 0;
1679             if (nslc >= VE[ive]->Nj) nslc = VE[ive]->Nj-1;
1680          } else if (orcode[2] == 'R' || orcode[2] == 'L') {
1681             dim = 2; nslc = slider;
1682             if (orcode[2] == 'L') nslc = VE[ive]->Nk-1-slider;
1683             if (nslc < 0) nslc = 0;
1684             if (nslc >= VE[ive]->Nk) nslc = VE[ive]->Nk-1;
1685          }
1686          break;
1687        case 'C': /* coronal slicing desired */
1688          if (orcode[0] == 'A' || orcode[0] == 'P') {
1689             dim = 0; nslc = slider;
1690             if (orcode[0] == 'P') nslc = VE[ive]->Ni-1-slider;
1691             if (nslc < 0) nslc = 0;
1692             if (nslc >= VE[ive]->Ni) nslc = VE[ive]->Ni-1;
1693          } else if (orcode[1] == 'A' || orcode[1] == 'P') {
1694             dim = 1; nslc = slider;
1695             if (orcode[1] == 'P') nslc = VE[ive]->Nj-1-slider;
1696             if (nslc < 0) nslc = 0;
1697             if (nslc >= VE[ive]->Nj) nslc = VE[ive]->Nj-1;
1698          } else if (orcode[2] == 'A' || orcode[2] == 'P') {
1699             dim = 2; nslc = slider;
1700             if (orcode[2] == 'P') nslc = VE[ive]->Nk-1-slider;
1701             if (nslc < 0) nslc = 0;
1702             if (nslc >= VE[ive]->Nk) nslc = VE[ive]->Nk-1;
1703          }
1704          break;
1705       case 'V': /* Volume rendering, do not do much */
1706          SUMA_S_Note("What to return here?");
1707          SUMA_RETURN(0);
1708          break;
1709       default:
1710          SUMA_S_Err("What gives? Variant '%s'", variant);
1711          SUMA_DUMP_TRACE("Bad variant trace");
1712          SUMA_RETURN(-1);
1713    }
1714    if (PlEq) { /* Equation of plane for slice in question */
1715       float mid_slc[3];
1716       /* Compute normal direction and slice center in index space */
1717       switch (dim) {
1718          case 0:
1719             I[0] = 1;            I[1] = 0;            I[2] = 0;
1720             C[0] = nslc;         C[1] = dims[1]/2.0;  C[2] = dims[2]/2.0;
1721             break;
1722          case 1:
1723             I[0] = 0;            I[1] = 1;            I[2] = 0;
1724             C[0] = dims[0]/2.0;  C[1] = nslc;         C[2] = dims[2]/2.0;
1725             break;
1726          case 2:
1727             I[0] = 0;            I[1] = 0;            I[2] = 1;
1728             C[0] = dims[0]/2.0;  C[1] = dims[1]/2.0;  C[2] = nslc;
1729             break;
1730       }
1731 
1732       AFF44_MULT_D(PlEq, VE[ive]->I2X, I); /* Compute normal in xyz */
1733 
1734       AFF44_MULT_I(mid_slc, VE[ive]->I2X, C); /* center of slice */
1735       SUMA_SHIFT_PLANE_TO_P(PlEq, mid_slc); /*go through mid of slice*/
1736    }
1737    if (tcorners || corners || slc_cen) {
1738       SUMA_dset_tex_slice_corners(nslc, dset, tcorners, corners, slc_cen,
1739                                   dim, voxcen);
1740    }
1741    SUMA_RETURN(dim);
1742 }
1743 
1744 /* Take an XYZ in RAI and return Ax, Sa, and Co GUI slider positions */
SUMA_XYZ_to_gui_slices(SUMA_VolumeElement ** VE,int ive,float * xyz,float * here)1745 float* SUMA_XYZ_to_gui_slices(SUMA_VolumeElement **VE, int ive,
1746                                     float *xyz, float *here)
1747 {
1748    static char FuncName[]={"SUMA_XYZ_to_gui_slices"};
1749    static float n[10][3];
1750    int icall=0;
1751    char *orcode;
1752    int dim=0, nslc=0, *dims;
1753    float I[3]={0.0, 0.0, 0.0}, C[3]={0.0, 0.0, 0.0};
1754    SUMA_DSET *dset=NULL;
1755    SUMA_Boolean LocalHead = NOPE;
1756 
1757    SUMA_ENTRY;
1758 
1759    if (icall > 9) icall=0;
1760    else ++icall;
1761 
1762    if (!here) here = (float *)(n[icall]);
1763    here[0] = here[1] = here[2] = -1;
1764 
1765    if (ive < 0) ive = 0;
1766    if (!xyz || !(dset = SUMA_VE_dset(VE, ive)) ||
1767        !(dims = SUMA_GetDatasetDimensions(dset)) ) {
1768       SUMA_S_Err("no dset or no variant") ;
1769       SUMA_RETURN(here);
1770    }
1771 
1772    orcode = SUMA_Dset_orcode(dset);
1773    if (orcode[0] == 'X') { SUMA_S_Err("No orcode"); SUMA_RETURN(here); }
1774 
1775    /* Change XYZ to I */
1776    AFF44_MULT_I(I, VE[ive]->X2I, xyz);
1777 
1778    /* Get the slider values for A, S, C, in this order */
1779    /* First the Ax slice number */
1780    dim = 0;
1781    if (orcode[0] == 'I' || orcode[0] == 'S') {
1782       if (orcode[0] == 'S') here[dim] = VE[ive]->Ni-1-SUMA_ROUND(I[0]);
1783       else here[dim] = SUMA_ROUND(I[0]);
1784       if (here[dim] < 0) here[dim] = 0;
1785       if (here[dim] >= VE[ive]->Ni) here[dim] = VE[ive]->Ni-1;
1786    } else if (orcode[1] == 'I' || orcode[1] == 'S') {
1787       if (orcode[1] == 'S') here[dim] = VE[ive]->Nj-1-SUMA_ROUND(I[1]);
1788       else here[dim] = SUMA_ROUND(I[1]);
1789       if (here[dim] < 0) here[dim] = 0;
1790       if (here[dim] >= VE[ive]->Nj) here[dim] = VE[ive]->Nj-1;
1791    } else if (orcode[2] == 'I' || orcode[2] == 'S') {
1792       if (orcode[2] == 'S') here[dim] = VE[ive]->Nk-1-SUMA_ROUND(I[2]);
1793       else here[dim] = SUMA_ROUND(I[2]);
1794       if (here[dim] < 0) here[dim] = 0;
1795       if (here[dim] >= VE[ive]->Nk) here[dim] = VE[ive]->Nk-1;
1796    }
1797    dim = 1;  /* Now look for Sag */
1798    if (orcode[0] == 'R' || orcode[0] == 'L') {
1799       if (orcode[0] == 'L') here[dim] = VE[ive]->Ni-1-SUMA_ROUND(I[0]);
1800       else here[dim] = SUMA_ROUND(I[0]);
1801       if (here[dim] < 0) here[dim] = 0;
1802       if (here[dim] >= VE[ive]->Ni) here[dim] = VE[ive]->Ni-1;
1803    } else if (orcode[1] == 'R' || orcode[1] == 'L') {
1804       if (orcode[1] == 'L') here[dim] = VE[ive]->Nj-1-SUMA_ROUND(I[1]);
1805       else here[dim] = SUMA_ROUND(I[1]);
1806       if (here[dim] < 0) here[dim] = 0;
1807       if (here[dim] >= VE[ive]->Nj) here[dim] = VE[ive]->Nj-1;
1808    } else if (orcode[2] == 'R' || orcode[2] == 'L') {
1809       if (orcode[2] == 'L') here[dim] = VE[ive]->Nk-1-SUMA_ROUND(I[2]);
1810       else here[dim] = SUMA_ROUND(I[2]);
1811       if (here[dim] < 0) here[dim] = 0;
1812       if (here[dim] >= VE[ive]->Nk) here[dim] = VE[ive]->Nk-1;
1813    }
1814    dim = 2;  /* Now look for Co */
1815    if (orcode[0] == 'A' || orcode[0] == 'P') {
1816       if (orcode[0] == 'P') here[dim] = VE[ive]->Ni-1-SUMA_ROUND(I[0]);
1817       else here[dim] = SUMA_ROUND(I[0]);
1818       if (here[dim] < 0) here[dim] = 0;
1819       if (here[dim] >= VE[ive]->Ni) here[dim] = VE[ive]->Ni-1;
1820    } else if (orcode[1] == 'A' || orcode[1] == 'P') {
1821       if (orcode[1] == 'P') here[dim] = VE[ive]->Nj-1-SUMA_ROUND(I[1]);
1822       else here[dim] = SUMA_ROUND(I[1]);
1823       if (here[dim] < 0) here[dim] = 0;
1824       if (here[dim] >= VE[ive]->Nj) here[dim] = VE[ive]->Nj-1;
1825    } else if (orcode[2] == 'A' || orcode[2] == 'P') {
1826       if (orcode[2] == 'P') here[dim] = VE[ive]->Nk-1-SUMA_ROUND(I[2]);
1827       else here[dim] = SUMA_ROUND(I[2]);
1828       if (here[dim] < 0) here[dim] = 0;
1829       if (here[dim] >= VE[ive]->Nk) here[dim] = VE[ive]->Nk-1;
1830    }
1831    SUMA_RETURN(here);
1832 }
1833 
1834 /* Set slices of VO at location xyz */
SUMA_VO_set_slices_XYZ(SUMA_VolumeObject * VOu,float * xyz)1835 SUMA_Boolean SUMA_VO_set_slices_XYZ(SUMA_VolumeObject *VOu, float *xyz)
1836 {
1837    static char FuncName[]={"SUMA_VO_set_slices_XYZ"};
1838    float *slices;
1839    int i;
1840    SUMA_VOL_SAUX *VSaux=NULL;
1841    SUMA_VolumeObject *VO=NULL;
1842    SUMA_Boolean LocalHead = NOPE;
1843 
1844    SUMA_ENTRY;
1845 
1846    if (!xyz) SUMA_RETURN(NOPE);
1847 
1848    for (i=0; i<SUMAg_N_DOv; ++i) {
1849       if (VOu) VO = VOu;
1850       else if (iDO_type(i) == VO_type) {
1851          VO = (SUMA_VolumeObject *)iDO_ADO(i);
1852       }else VO = NULL;
1853       VSaux = (SUMA_VOL_SAUX *)VDO_VSAUX(VO);
1854       if ((VO && VSaux->SlicesAtCrosshair) || VOu) { /* Do this if user
1855                                                supplies volume or
1856                                                if volumes request it */
1857          slices = SUMA_XYZ_to_gui_slices(VO->VE, 0, xyz, NULL);
1858          SUMA_LH("Jumping to slices Ax%f Sa%f Co%f",
1859                   slices[0], slices[1], slices[2]);
1860          SUMA_set_slice((SUMA_ALL_DO *)VO, "Ax",
1861                                   slices, "EngXYZ", 0);
1862          SUMA_set_slice((SUMA_ALL_DO *)VO, "Sa",
1863                                   slices+1, "EngXYZ", 0);
1864          SUMA_set_slice((SUMA_ALL_DO *)VO, "Co",
1865                                   slices+2, "EngXYZ", 0);
1866       }
1867       if (VOu) break;
1868    }
1869    SUMA_RETURN(YUP);
1870 }
1871 
SUMA_dset_tex_slice_corners(int slci,SUMA_DSET * dset,GLfloat * tcorners,GLfloat * corners,GLfloat * slc_cen,int dim,int voxcen)1872 void SUMA_dset_tex_slice_corners( int slci, SUMA_DSET *dset,
1873                            GLfloat *tcorners, GLfloat *corners, GLfloat *slc_cen,
1874                               int dim, int voxcen)
1875 {
1876    static char FuncName[]={"SUMA_dset_tex_slice_corners"};
1877    int kk=0;
1878    int slcx, slcy, slcz=0, *dims;
1879    float A[4][4], I[3], T[3];
1880    SUMA_Boolean LocalHead = NOPE;
1881 
1882    SUMA_ENTRY;
1883 
1884    if (tcorners) memset(tcorners, 0, 4*3*sizeof(float));
1885    if (corners) memset(corners, 0, 4*3*sizeof(float));
1886 
1887    if (!dset || !(SUMA_GetDatasetI2X(dset, A)) ||
1888        !(dims=SUMA_GetDatasetDimensions(dset)) ||
1889        !((tcorners && corners) || slc_cen) ) {
1890       SUMA_S_Err("no valid ijk_to_dicom_real or null input") ;
1891       SUMA_RETURNe;
1892    }
1893 
1894    if (tcorners && corners) {
1895       switch (dim) {
1896          default:
1897             SUMA_S_Err("Bad dim value");
1898             SUMA_RETURNe;
1899          case 2:
1900             kk = 0;
1901             I[0] = 0; I[1] = 0; I[2] = slci;
1902             T[0] = 0; T[1] = 0; T[2] = ((float)slci+0.5)/(float)dims[2];
1903             if (!voxcen) {I[0] -= 0.5; I[1] -= 0.5; }
1904             AFF44_MULT_I((corners+3*kk), A, I);
1905             SUMA_COPY_VEC(T, tcorners+3*kk, 3, float, GLfloat);
1906 
1907             kk = 1;
1908             I[0] = dims[0]-1; I[1] = 0; I[2] = slci;
1909             T[0] = 1;               T[1] = 0; T[2] = tcorners[2];
1910             if (!voxcen) {I[0] += 0.5; I[1] -= 0.5; }
1911             AFF44_MULT_I((corners+3*kk), A, I);
1912             SUMA_COPY_VEC(T, tcorners+3*kk, 3, float, GLfloat);
1913 
1914             kk = 2;
1915             I[0] = dims[0]-1; I[1] = dims[1]-1; I[2] = slci;
1916             T[0] = 1;               T[1] = 1;               T[2] = tcorners[2];
1917             if (!voxcen) {I[0] += 0.5; I[1] += 0.5; }
1918             AFF44_MULT_I((corners+3*kk), A, I);
1919             SUMA_COPY_VEC(T, tcorners+3*kk, 3, float, GLfloat);
1920 
1921             kk = 3;
1922             I[0] = 0; I[1] = dims[1]-1; I[2] = slci;
1923             T[0] = 0; T[1] = 1;               T[2] = tcorners[2];
1924             if (!voxcen) {I[0] -= 0.5; I[1] += 0.5; }
1925             AFF44_MULT_I((corners+3*kk), A, I);
1926             SUMA_COPY_VEC(T, tcorners+3*kk, 3, float, GLfloat);
1927             break;
1928 
1929          case 1:
1930             kk = 0;
1931             I[0] = 0; I[1] = slci;                                   I[2] = 0;
1932             T[0] = 0; T[1] = ((float)slci+0.5)/(float)dims[1]; T[2] = 0;
1933             if (!voxcen) {I[0] -= 0.5; I[2] -= 0.5; }
1934             AFF44_MULT_I((corners+3*kk), A, I);
1935             SUMA_COPY_VEC(T, tcorners+3*kk, 3, float, GLfloat);
1936 
1937             kk = 1;
1938             I[0] = dims[0]-1; I[1] = slci;        I[2] = 0;
1939             T[0] = 1;               T[1] = tcorners[1]; T[2] = 0;
1940             if (!voxcen) {I[0] += 0.5; I[2] -= 0.5; }
1941             AFF44_MULT_I((corners+3*kk), A, I);
1942             SUMA_COPY_VEC(T, tcorners+3*kk, 3, float, GLfloat);
1943 
1944             kk = 2;
1945             I[0] = dims[0]-1; I[1] = slci;        I[2] = dims[2]-1;
1946             T[0] = 1;               T[1] = tcorners[1]; T[2] = 1;
1947             if (!voxcen) {I[0] += 0.5; I[2] += 0.5; }
1948             AFF44_MULT_I((corners+3*kk), A, I);
1949             SUMA_COPY_VEC(T, tcorners+3*kk, 3, float, GLfloat);
1950 
1951             kk = 3;
1952             I[0] = 0; I[1] = slci;        I[2] = dims[2]-1;
1953             T[0] = 0; T[1] = tcorners[1]; T[2] = 1;
1954             if (!voxcen) {I[0] -= 0.5; I[2] += 0.5; }
1955             AFF44_MULT_I((corners+3*kk), A, I);
1956             SUMA_COPY_VEC(T, tcorners+3*kk, 3, float, GLfloat);
1957             break;
1958 
1959          case 0:
1960             kk = 0;
1961             I[0] = slci;                                   I[1] = 0; I[2] = 0;
1962             T[0] = ((float)slci+0.5)/(float)dims[0]; T[1] = 0; T[2] = 0;
1963             if (!voxcen) {I[1] -= 0.5; I[2] -= 0.5; }
1964             AFF44_MULT_I((corners+3*kk), A, I);
1965             SUMA_COPY_VEC(T, tcorners+3*kk, 3, float, GLfloat);
1966 
1967             kk = 1;
1968             I[0] = slci;        I[1] = dims[1]-1; I[2] = 0;
1969             T[0] = tcorners[0]; T[1] = 1;               T[2] = 0;
1970             if (!voxcen) {I[1] += 0.5; I[2] -= 0.5; }
1971             AFF44_MULT_I((corners+3*kk), A, I);
1972             SUMA_COPY_VEC(T, tcorners+3*kk, 3, float, GLfloat);
1973 
1974             kk = 2;
1975             I[0] = slci;        I[1] = dims[1]-1; I[2] = dims[2]-1;
1976             T[0] = tcorners[0]; T[1] = 1;               T[2] = 1;
1977             if (!voxcen) {I[1] += 0.5; I[2] += 0.5; }
1978             AFF44_MULT_I((corners+3*kk), A, I);
1979             SUMA_COPY_VEC(T, tcorners+3*kk, 3, float, GLfloat);
1980 
1981             kk = 3;
1982             I[0] = slci;        I[1] = 0; I[2] = dims[2]-1;
1983             T[0] = tcorners[0]; T[1] = 0; T[2] = 1;
1984             if (!voxcen) {I[1] -= 0.5; I[2] += 0.5; }
1985             AFF44_MULT_I((corners+3*kk), A, I);
1986             SUMA_COPY_VEC(T, tcorners+3*kk, 3, float, GLfloat);
1987             break;
1988       }
1989 
1990       if (LocalHead) {
1991          float cen[3] = {0.0, 0.0, 0.0};
1992          SUMA_LHv("Slice %d, dim %d %s corners and textures\n",
1993                   slci, dim,
1994                   voxcen ? "Voxel Center":"Slice Edge");
1995          for (kk=0; kk<4; ++kk) {
1996             fprintf(SUMA_STDERR,
1997                      "c%d: %.3f   %.3f   %.3f\n"
1998                      "t%d: %.3f   %.3f   %.3f\n",
1999                kk, corners[3*kk],corners[3*kk+1],corners[3*kk+2],
2000                kk, tcorners[3*kk],tcorners[3*kk+1],tcorners[3*kk+2]);
2001            cen[0] += corners[3*kk  ];
2002            cen[1] += corners[3*kk+1];
2003            cen[2] += corners[3*kk+2];
2004          }
2005          fprintf(SUMA_STDERR,"\n");
2006          cen[0] /= 4.0; cen[1] /= 4.0; cen[2] /= 4.0;
2007          fprintf(SUMA_STDERR, "Slice center RAI coord: %f %f %fmm\n",
2008                               cen[0], cen[1], cen[2]);
2009       }
2010    }
2011    if (slc_cen) {
2012       switch (dim) {
2013          default:
2014             SUMA_S_Err("Bad dim value");
2015             SUMA_RETURNe;
2016          case 2:
2017             I[0] = dims[0]/2.0; I[1] = dims[1]/2.0; I[2] = slci;
2018             if (!voxcen) {I[0] -= 0.5; I[1] -= 0.5; }
2019             AFF44_MULT_I((slc_cen), A, I);
2020             break;
2021          case 1:
2022             I[0] = dims[0]/2.0; I[1] = slci; I[2] = dims[2]/2.0;
2023             if (!voxcen) {I[0] -= 0.5; I[2] -= 0.5; }
2024             AFF44_MULT_I((slc_cen), A, I);
2025             break;
2026          case 0:
2027             I[0] = slci; I[1] = dims[1]/2.0; I[2] = dims[2]/2.0;
2028             if (!voxcen) {I[1] -= 0.5; I[2] -= 0.5; }
2029             AFF44_MULT_I((slc_cen), A, I);
2030             break;
2031       }
2032       SUMA_LH( "Slice center at RAI mm: %f %f %f",
2033                slc_cen[0], slc_cen[1], slc_cen[2]);
2034    }
2035 
2036    SUMA_RETURNe;
2037 }
2038 
SUMA_VO_SelectedSlice(SUMA_VolumeObject * vo,char * variant,float * scorners)2039 int SUMA_VO_SelectedSlice(SUMA_VolumeObject *vo, char *variant, float *scorners)
2040 {
2041    static char FuncName[]={"SUMA_VO_SelectedSlice"};
2042    SUMA_ALL_DO *ado=(SUMA_ALL_DO *)vo;
2043    SUMA_VOL_SAUX *VSaux = SUMA_ADO_VSaux(ado);
2044    int nslc = -1, dim, k;
2045    GLfloat slc_corners[12], slc_tcorners[12];
2046    SUMA_Boolean LocalHead = NOPE;
2047 
2048    SUMA_ENTRY;
2049 
2050    if (!( (VSaux = SUMA_ADO_VSaux(ado)) &&
2051           variant &&
2052           VSaux->PR &&
2053           VSaux->PR->iAltSel[SUMA_VOL_I] >= 0 &&
2054           VSaux->PR->iAltSel[SUMA_VOL_J] >= 0 &&
2055           VSaux->PR->iAltSel[SUMA_VOL_K] >= 0 ) ) {
2056       SUMA_RETURN(-1);
2057    }
2058 
2059    #if 0 /* Old way, does not work when in montage mode.
2060             Keep for the record */
2061    if ((dim = SUMA_dset_gui_slice_from_tex_slice_d(vo->VE, 0,
2062                                  VSaux->PR->dAltSel+SUMA_VOL_SLC_EQ0,
2063                                  0, variant,NULL))< 0) {
2064       SUMA_RETURN(-1);
2065    }
2066 
2067 
2068    nslc = VSaux->PR->iAltSel[SUMA_VOL_I+dim];
2069 
2070    SUMA_LH("Slice variant %s, dim %d. [%ld %ld %ld] -->%d",
2071             variant, dim,
2072             VSaux->PR->iAltSel[SUMA_VOL_I],
2073             VSaux->PR->iAltSel[SUMA_VOL_J],
2074             VSaux->PR->iAltSel[SUMA_VOL_K],
2075             nslc);
2076    if (nslc >= 0 && scorners) {
2077       SUMA_dset_tex_slice_corners_gui(vo->VE, 0, variant, nslc,
2078                           slc_tcorners, slc_corners,
2079                           NULL, NULL, 0 );
2080       for (k=0; k<12; ++k) scorners[k] = slc_corners[k];
2081    }
2082    #else
2083    SUMA_SlcCodeToVariant(VSaux->PR->iAltSel[SUMA_VOL_SLC_VARIANT],
2084                          variant);
2085    nslc = VSaux->PR->iAltSel[SUMA_VOL_SLC_NUM];
2086    SUMA_LH("Slice variant %s, [%ld %ld %ld] -->%d",
2087             variant,
2088             VSaux->PR->iAltSel[SUMA_VOL_I],
2089             VSaux->PR->iAltSel[SUMA_VOL_J],
2090             VSaux->PR->iAltSel[SUMA_VOL_K],
2091             nslc);
2092    if (nslc >= 0 && scorners &&
2093        VSaux->PR->iAltSel[SUMA_VOL_SLC_VARIANT] != SUMA_VR_VARIANT) {
2094       SUMA_dset_tex_slice_corners_gui(vo->VE, 0, variant, nslc,
2095                           slc_tcorners, slc_corners,
2096                           NULL, NULL, 0 );
2097       for (k=0; k<12; ++k) scorners[k] = slc_corners[k];
2098    }
2099    #endif
2100    SUMA_RETURN(nslc);
2101 }
2102 
2103 
SUMA_DrawVolumeDO_OLD(SUMA_VolumeObject * VO,SUMA_SurfaceViewer * sv)2104 SUMA_Boolean SUMA_DrawVolumeDO_OLD(SUMA_VolumeObject *VO, SUMA_SurfaceViewer *sv)
2105 {
2106    static char FuncName[]={"SUMA_DrawVolumeDO_OLD"};
2107    int i = 0, k = 0, j=0, ive=0;
2108    float iq[4]={0, 0, 0, 0}, vo0[3], voN[3];
2109    static int ipass=0, iplane = 0;
2110    GLfloat tex_corn[12] ;
2111    GLfloat slc_corn[12] ;
2112    GLfloat rotationMatrix[4][4], rt[4][4];
2113    GLboolean gl_dt, gl_bl;
2114    int ShowUnselected = 1;
2115    float tz = 0.0;
2116    static GLfloat init_rotationMatrix[4][4];
2117    static GLdouble dmatrix[16], init_mv_matrix[16];
2118    float* nlt; /*JB: temporary node list, because I do not want to type
2119                   "VO->SOcut[0]->NodeList" over and over...*/
2120    SUMA_Boolean LastTextureOnCutPlane=YUP;
2121    SUMA_Boolean LocalHead = NOPE;
2122 
2123    SUMA_ENTRY;
2124 
2125    if (!VO) SUMA_RETURN(NOPE);
2126    if (!sv) sv = &(SUMAg_SVv[0]);
2127 
2128    if (sv->DO_PickMode) {
2129       SUMA_LH("No need to draw volume in DO_PickMode");
2130       SUMA_RETURN(YUP);
2131    }
2132 
2133    if (!VO->Show) SUMA_RETURN(YUP);
2134 
2135    if (sv->PolyMode != SRM_Fill) {
2136       /* fill it up */
2137       glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
2138    }
2139 
2140    if (!VO->SOcut || !VO->SOcut[0]) SUMA_VO_InitCutPlanes(VO);
2141 
2142 
2143    if (0) {
2144       SUMA_LH( "Draw Clipping planes without modelview matrix enabled\n"
2145             "So planes are fixed in space.\n"
2146             "We're more used to creating cutting planes that \n"
2147             "are attached to the object. So should look into that\n"
2148             "instead. \n"
2149             "Note that the use of these clipping planes\n"
2150             "seem necessary when the volume is not padded with zeros.\n"
2151             "I don't quite know what that is, but I suspect it has to do\n"
2152             "with alphas being non zero at the tip and the placement of\n"
2153             "texture quads\n"  );
2154    }
2155    glClipPlane(GL_CLIP_PLANE0, VO->CutPlane[0]);
2156    glClipPlane(GL_CLIP_PLANE1, VO->CutPlane[1]);
2157    glClipPlane(GL_CLIP_PLANE2, VO->CutPlane[2]);
2158    glClipPlane(GL_CLIP_PLANE3, VO->CutPlane[3]);
2159    glClipPlane(GL_CLIP_PLANE4, VO->CutPlane[4]);
2160    glClipPlane(GL_CLIP_PLANE5, VO->CutPlane[5]); /* play with this
2161                                              one to change cuts*/
2162 
2163 
2164    /* Now we need to draw the sucker */
2165    SUMA_CHECK_GL_ERROR("OpenGL Error pre texture");
2166    glEnable(GL_TEXTURE_3D);
2167 
2168    /* enable clipping planes. */
2169    glEnable(GL_CLIP_PLANE0);
2170    glEnable(GL_CLIP_PLANE1);
2171    glEnable(GL_CLIP_PLANE2);
2172    glEnable(GL_CLIP_PLANE3);
2173    glEnable(GL_CLIP_PLANE4);
2174    glEnable(GL_CLIP_PLANE5);
2175 
2176    glTexEnvf(  GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE,
2177                VO->TexEnvMode   ); /* what happens if
2178                               there is color already on a vertex (I would likely
2179                               not need this for 3D textures...*/
2180 
2181    gl_dt = glIsEnabled(GL_DEPTH_TEST);
2182    gl_bl = glIsEnabled(GL_BLEND);
2183 
2184    ive = 0;
2185    while (VO->VE && VO->VE[ive]) {
2186       if (!VO->VE[ive]->texName) {
2187          if (!SUMA_CreateGL3DTexture(VO)) {
2188             SUMA_S_Err("Failed to create texture");
2189             SUMA_RETURN(NOPE);
2190          }
2191       }
2192 
2193       SUMA_LHv("About to bind texture %d for %s\n",
2194             VO->VE[ive]->texName[0], SUMA_VE_Headname(VO->VE, ive));
2195       glBindTexture(GL_TEXTURE_3D, VO->VE[ive]->texName[0]);
2196                                              /* make texName be current */
2197 
2198       /* Now generate the coordinates */
2199       SUMA_LHv( "About to generate polygons for dset %d (%p)\n",
2200                ive, SUMA_VE_dset(VO->VE, ive));
2201       SUMA_LHv( "About to generate polygons for dset %d (%s)\n",
2202                ive, SUMA_VE_Headname(VO->VE, ive));
2203       glShadeModel(GL_FLAT);     /* This should be reverted to sv's settings */
2204       if (!(gl_dt)) glEnable(GL_DEPTH_TEST);
2205       if (!(gl_bl)) glEnable(GL_BLEND);
2206 	   glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
2207       if (0 && LocalHead) {
2208          SUMA_GL_MAT_SHOW(GL_TEXTURE_MATRIX,"Tx PreSetup\n");
2209          SUMA_GL_MAT_SHOW(GL_MODELVIEW_MATRIX,"MV PreSetup\n");
2210       }
2211       SUMA_CHECK_GL_ERROR("OpenGL Error pre setup");
2212 
2213       /* The modelview matrix for drawing the quad vertices should remain the
2214       same, the texture matrix will reflect object rotations.
2215       IF you do want to see the 'slices', i.e. a stack of axial slices (for RAI
2216       dset), rather than volume rendering image, then allo the quads to be drawn
2217       with the scene's current modelview matrix and keep the texture matrix
2218       as the identity matrix.
2219       This slice rendering could be useful to show
2220       a triplet of slices (Ax, Sa, and Co), rather than the entire volume say.
2221       But that should
2222       be done as a separate nel, or perhaps a different way of drawing a dset
2223       within SUMA. One should also decide whether that display is better done
2224       as 3 2D textures, rather than storing an entire 3D texture and just drawing
2225       a few quads from it */
2226       if (!ipass) {
2227          /* store the Modelview_matrix to reuse later
2228          Later on, this should be created automatically, when
2229          dset is loaded. I am not sure if using whatever was there
2230          when dset was first loaded is ideal in terms of artifact reduction
2231          The important thing however, is that the quads are drawn with the
2232          same rotation matrix, translation should be applied normally.
2233          It is the texture matrix that is changed
2234          with mouse movements, thereby allowing the volume to appear to move
2235          with other objects in the scene */
2236 
2237          if (LocalHead) {
2238             SUMA_GL_MAT_SHOW(GL_MODELVIEW_MATRIX,"MV PreRecord\n");
2239          }
2240          glGetDoublev(GL_MODELVIEW_MATRIX, init_mv_matrix);
2241          if (LocalHead) {
2242             fprintf(stderr,"init_mv_matrix initialized to:\n");
2243             for (i=0;i<16;++i) fprintf(stderr,"%.3f\t", init_mv_matrix[i]);
2244             fprintf(stderr,"\n");
2245          }
2246          /* also store the rotation matrix at initialization, this is
2247          done to account for any prerotations already present on the scene */
2248          SUMA_build_rotmatrix(init_rotationMatrix,
2249                               sv->GVS[sv->StdView].currentQuat);
2250          SUMA_S_Note("The use of this ipass variable is silly and would\n"
2251                      "only work when we are loading one texture only once.\n"
2252                      "The parameters recorded here, should be part of nel,\n"
2253                      "at the initialization stage and not stored this way.\n"
2254                      "The same rant goes for silly clipping planes.\n"  );
2255          ++ipass;
2256       } else {
2257          if (0 && LocalHead) {
2258             SUMA_LH("init_mv_matrix recorded to be:");
2259             for (i=0;i<16;++i) fprintf(stderr,"%.3f\t", init_mv_matrix[i]);
2260             fprintf(stderr,"\n");
2261          }
2262       }
2263 
2264       /* first set the modelview to the one existing when the texture was
2265          first rendered*/
2266          glMatrixMode(GL_MODELVIEW); glPushMatrix(); /* get a new one*/
2267          glLoadIdentity(); /* set to identity*/
2268          /* Now apply constant Modelview for Quads*/
2269          glTranslatef ( sv->GVS[sv->StdView].translateVec[0],
2270                      sv->GVS[sv->StdView].translateVec[1], 0.0);
2271          glMultMatrixd(init_mv_matrix);
2272 
2273       /* now for the texture coordinates, these should be transformed by
2274          the inverse of what is done to the objects in the scene */
2275          glMatrixMode(GL_TEXTURE); glPushMatrix(); /* get a new one*/
2276          glLoadIdentity(); /* set to identity */
2277 
2278          /* Now apply the inverse rotation params to the texture coordinates */
2279          SUMA_build_rotmatrix(rotationMatrix, sv->GVS[sv->StdView].currentQuat);
2280 
2281          if (sv->ortho) {
2282             SUMA_TRANSP_4MATRIX(rotationMatrix); /* transpose = inverse
2283                                                 for rotation matrix*/
2284          } else {
2285             SUMA_INV_44ROTATIONMATRIX(rotationMatrix);
2286          }
2287          glTranslatef ( 0.5, 0.5, 0.5);
2288          glMultMatrixf(&rotationMatrix[0][0]);
2289          glMultMatrixf(&init_rotationMatrix[0][0]);  /* Need to multiply by
2290                                                         initial rotation applied
2291                                                         to the scene. (initial=at
2292                                                         the time quads were 1st
2293                                                         drawn*/
2294          glTranslatef (-0.5, -0.5, -0.5);
2295 
2296 
2297       if (0 && LocalHead) {
2298          SUMA_GL_MAT_SHOW(GL_TEXTURE_MATRIX,"Tx PreDraw\n");
2299          SUMA_GL_MAT_SHOW(GL_MODELVIEW_MATRIX,"MV PreDraw\n");
2300       }
2301 
2302       SUMA_CHECK_GL_ERROR("OpenGL Error pre draw");
2303       #if 1
2304       /* slice by slice drawing, doing the deed by hand*/
2305       for(i = 0; i < VO_NK(VO); i++) {
2306          glBegin(GL_QUADS);
2307             SUMA_dset_tex_slice_corners(i, SUMA_VE_dset(VO->VE,ive), tex_corn,
2308                                         slc_corn, NULL, 2, 0);
2309             for (k=0; k<4; ++k) {
2310                glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tex_corn[3*k+2]);
2311                      /* this one is affected by the Texture MatrixMode */
2312                glVertex3f(slc_corn[3*k], slc_corn[3*k+1], slc_corn[3*k+2]);
2313                      /* this one is affected by the Modelview matrixMode*/
2314             }
2315          glEnd();
2316       }
2317       #else /* automatic texture coordinate generation,
2318             same artifact, and more complicated splane and rplane stuff, leave it
2319            just in case I'll need it. */
2320        glEnable(GL_TEXTURE_GEN_S);
2321        glEnable(GL_TEXTURE_GEN_T);
2322        glEnable(GL_TEXTURE_GEN_R);
2323 
2324        glTexGeni(GL_S, GL_TEXTURE_GEN_MODE, GL_OBJECT_LINEAR);
2325        glTexGeni(GL_T, GL_TEXTURE_GEN_MODE, GL_OBJECT_LINEAR);
2326        glTexGeni(GL_R, GL_TEXTURE_GEN_MODE, GL_OBJECT_LINEAR);
2327 
2328        glTexGenfv(GL_S, GL_OBJECT_PLANE, splane);/*This ordering of the planes */
2329        glTexGenfv(GL_T, GL_OBJECT_PLANE, rplane);  /* resulted in proper axis */
2330        glTexGenfv(GL_R, GL_OBJECT_PLANE, tplane);  /*  alignment ...*/
2331 
2332 
2333 
2334        for(i = 0; i < DSET_NZ(VO->VE[ive]->dset); i++) {
2335          glBegin(GL_QUADS);
2336             SUMA_dset_tex_slice_corners(i,VO->VE[ive]->dset,
2337                                           tex_corn, slc_corn, NULL, 2, 0);
2338             for (k=0; k<4; ++k) {
2339                glVertex3f(slc_corn[3*k], slc_corn[3*k+1], slc_corn[3*k+2]);
2340             }
2341          glEnd();
2342       }
2343       glDisable(GL_TEXTURE_GEN_S);
2344       glDisable(GL_TEXTURE_GEN_T);
2345       glDisable(GL_TEXTURE_GEN_R);
2346       #endif
2347        SUMA_CHECK_GL_ERROR("OpenGL Error ddd");
2348 
2349       glMatrixMode(GL_TEXTURE); glPopMatrix(); /* pop matrix of GL_TEXTURE */
2350       glMatrixMode(GL_MODELVIEW);  glPopMatrix();/* and Modelview*/
2351 
2352       ++ive;
2353    }
2354 
2355    glFlush();
2356 
2357    /* disable the clipping planes */
2358    glDisable(GL_CLIP_PLANE0);
2359    glDisable(GL_CLIP_PLANE1);
2360    glDisable(GL_CLIP_PLANE2);
2361    glDisable(GL_CLIP_PLANE3);
2362    glDisable(GL_CLIP_PLANE4);
2363    glDisable(GL_CLIP_PLANE5);
2364 
2365 
2366    /* Here we create a texture on the cutplane from the last dset loaded
2367    At the moment, without this texture, nothing shows
2368    of the overlay volume */
2369    if (LastTextureOnCutPlane) {
2370       --ive; /* bring ive counter to last dset put into texture*/
2371       SUMA_dset_tex_slice_corners( 0, SUMA_VE_dset(VO->VE, ive),
2372                                    tex_corn, slc_corn, NULL, 2, 0);
2373 
2374       // Joachim says this is just wrong ...
2375       tz = 0.5+(-VO->CutPlane[0][3])/(float)VO_NK(VO);
2376 
2377       if (!gl_dt) glDisable(GL_DEPTH_TEST);
2378 
2379       /* If it were not for the slice textures shown here, then the overlay
2380          texture would not show up at all! */
2381       #if 0
2382       SUMA_LH("Texture on the slice, with triangles");
2383       glBegin(GL_TRIANGLES);
2384          k = 0;
2385          glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
2386             glVertex3f( slc_corn[3*k], slc_corn[3*k+1],
2387                                     -VO->CutPlane[0][3]); ++k;
2388          glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
2389             glVertex3f( slc_corn[3*k], slc_corn[3*k+1],
2390                                     -VO->CutPlane[0][3]); ++k;
2391          glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
2392             glVertex3f( slc_corn[3*k], slc_corn[3*k+1],
2393                                     -VO->CutPlane[0][3]);
2394 
2395          k = 0;
2396          glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
2397             glVertex3f( slc_corn[3*k], slc_corn[3*k+1],
2398                                     -VO->CutPlane[0][3]); k+=2;
2399          glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
2400             glVertex3f( slc_corn[3*k], slc_corn[3*k+1],
2401                                     -VO->CutPlane[0][3]); ++k;
2402          glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
2403             glVertex3f( slc_corn[3*k], slc_corn[3*k+1],
2404                                     -VO->CutPlane[0][3]);
2405 
2406       glEnd();
2407       #else
2408       SUMA_LH("Texture on the slice, QUADS?");
2409       glBegin(GL_QUADS);
2410          for (k=0; k<4; ++k) {
2411             glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
2412                   /* this one is affected by the Texture MatrixMode */
2413             glVertex3f(slc_corn[3*k], slc_corn[3*k+1],
2414                                           -VO->CutPlane[0][3]);
2415                   /* this one is affected by the Modelview matrixMode*/
2416          }
2417       glEnd();
2418       #endif
2419    }
2420       glDisable(GL_TEXTURE_3D);
2421    if (!gl_dt) glDisable(GL_DEPTH_TEST);
2422    if (!gl_bl) glDisable(GL_BLEND);
2423 
2424    glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
2425    glEnable(GL_COLOR_MATERIAL);
2426    for (iplane=0; iplane < 6; ++iplane) {
2427       if (VO->UseCutPlane[iplane]) {
2428          if (iplane == VO->SelectedCutPlane) glColor3f(1.0, 1.0, 1.0);
2429          else {
2430             if (ShowUnselected) {
2431                if (iplane==0 || iplane == 1) glColor3f(1.0, 0.0, 0.0);
2432                if (iplane==2 || iplane == 3) glColor3f(0.0, 1.0, 0.0);
2433                if (iplane==4 || iplane == 5) glColor3f(0.0, 0.0, 1.0);
2434             } else {
2435                continue;
2436             }
2437          }
2438          glBegin(GL_LINE_LOOP);
2439             nlt = VO->SOcut[iplane]->NodeList;
2440             glVertex3f( nlt[0],nlt[1],nlt[2] );
2441             glVertex3f( nlt[3],nlt[4],nlt[5] );
2442             glVertex3f( nlt[6],nlt[7],nlt[8] );
2443             glVertex3f( nlt[9],nlt[10],nlt[11] );
2444          glEnd();
2445       }
2446    }
2447    glDisable(GL_COLOR_MATERIAL);
2448 
2449    if (sv->PolyMode != SRM_Fill) {/* set fill mode back */
2450       SUMA_SET_GL_RENDER_MODE(sv->PolyMode);
2451    }
2452 
2453 
2454    SUMA_RETURN(YUP);
2455 
2456 }
2457 
SUMA_GET_VR_Slice_Pack(SUMA_VolumeObject * VO,SUMA_SurfaceViewer * sv)2458 SUMA_Boolean SUMA_GET_VR_Slice_Pack(SUMA_VolumeObject *VO,
2459                                     SUMA_SurfaceViewer *sv)
2460 {
2461    static char FuncName[]={"SUMA_GET_VR_Slice_Pack"};
2462    SUMA_VOL_SAUX *VSaux=NULL;
2463    SUMA_RENDERED_SLICE *rslc=NULL;
2464    SUMA_X_SurfCont *SurfCont = NULL;
2465    float *cen = NULL, Eq[4], *PlOff=NULL;
2466    int nn;
2467    int N_slc=150;
2468    SUMA_Boolean LocalHead = NOPE;
2469 
2470    SUMA_ENTRY;
2471 
2472    if (!VO || !(VSaux = SUMA_ADO_VSaux((SUMA_ALL_DO *)VO)) ||
2473        !(SurfCont = SUMA_ADO_Cont((SUMA_ALL_DO *)VO))) {
2474       SUMA_RETURN(NOPE);
2475    }
2476 
2477    if (  SurfCont->VR_fld->N_slice_num < 0 ||
2478          SurfCont->VR_fld->N_slice_num > 2000) {
2479       N_slc = 150;
2480    } else N_slc = (int) SurfCont->VR_fld->N_slice_num;
2481 
2482    cen = SUMA_VO_Grid_Center(VO, NULL);
2483    SUMA_ScreenPlane_WorldSpace(sv, cen, Eq);
2484    PlOff = (float *)SUMA_calloc(N_slc, sizeof(float));
2485    if (!PlOff || (SUMA_PlaneBoxSlice( sv->GVS[sv->StdView].ViewFrom, Eq,
2486                                       VO->VE[0]->bcorners, NULL, NULL,
2487                                       PlOff, N_slc) < 0)) {
2488       SUMA_S_Err("Failed to allocate or get %d slicing planes", N_slc);
2489       SUMA_ifree(PlOff);
2490       SUMA_RETURN(NOPE);
2491    }
2492    for (nn=0; nn<N_slc; ++nn) {
2493       rslc = (SUMA_RENDERED_SLICE *) SUMA_malloc(sizeof(SUMA_RENDERED_SLICE));
2494       rslc->Eq[0] = Eq[0]; rslc->Eq[1] = Eq[1]; rslc->Eq[2] = Eq[2];
2495       rslc->Eq[3] = PlOff[nn] ; sprintf(rslc->variant,"Vr");
2496       /* stick plane in list, last one to be rendered goes to top */
2497       SUMA_LH("Intersecting VR plane %f %f %f %f, on vol %s\n"
2498               "(origin %f %f %f)",
2499                rslc->Eq[0], rslc->Eq[1], rslc->Eq[2], rslc->Eq[3],
2500                SUMA_VE_Headname(VO->VE,0),
2501                VO->VE[0]->I2X[3][0],
2502                VO->VE[0]->I2X[3][1],VO->VE[0]->I2X[3][2]);
2503       dlist_ins_next(VSaux->vrslcl, dlist_head(VSaux->vrslcl), rslc);
2504    }
2505    SUMA_ifree(PlOff);
2506    SUMA_RETURN(YUP);
2507 }
2508 
SUMA_SlcVariantToCode(char * variant)2509 SUMA_VOL_REN_VARIANTS SUMA_SlcVariantToCode(char *variant)
2510 {
2511 
2512    static char FuncName[]={"SUMA_SlcVariantToCode"};
2513    if (!variant) {
2514       SUMA_S_Err("NULL variant");
2515       return(SUMA_ERR_VARIANT);
2516    }
2517    if (!strcmp(variant,"Ax")) return(SUMA_AX_VARIANT);
2518    if (!strcmp(variant,"Sa")) return(SUMA_SA_VARIANT);
2519    if (!strcmp(variant,"Co")) return(SUMA_CO_VARIANT);
2520    if (!strcmp(variant,"Vr")) return(SUMA_VR_VARIANT);
2521    SUMA_S_Err("Variant >%s< not recognized", variant);
2522    return(SUMA_ERR_VARIANT);
2523 }
2524 
SUMA_SlcCodeToVariant(SUMA_VOL_REN_VARIANTS v,char * variant)2525 void SUMA_SlcCodeToVariant(SUMA_VOL_REN_VARIANTS v, char *variant)
2526 {
2527    static char FuncName[]={"SUMA_SlcCodeToVariant"};
2528    SUMA_Boolean LocalHead = NOPE;
2529 
2530    variant[0] = '\0';
2531    switch (v) {
2532       default:
2533          SUMA_S_Err("Variant code %d unrecognized", v);
2534          SUMA_DUMP_TRACE("From here");
2535          variant[0] = '\0';
2536          break;
2537       case SUMA_ERR_VARIANT: /* Don't complain, happens if you remotely
2538                                 select a voxel when no slice is selected
2539                                 yet */
2540          SUMA_LH("Variant code %d is error flag", v);
2541          variant[0] = '\0';
2542          break;
2543       case SUMA_AX_VARIANT:
2544          variant[0] = 'A'; variant[1] = 'x'; variant[2] = '\0';
2545          break;
2546       case SUMA_SA_VARIANT:
2547          variant[0] = 'S'; variant[1] = 'a'; variant[2] = '\0';
2548          break;
2549       case SUMA_CO_VARIANT:
2550          variant[0] = 'C'; variant[1] = 'o'; variant[2] = '\0';
2551          break;
2552       case SUMA_VR_VARIANT:
2553          variant[0] = 'V'; variant[1] = 'r'; variant[2] = '\0';
2554          break;
2555    }
2556    return;
2557 }
2558 
2559 
SUMA_Get_Slice_Pack(SUMA_VolumeObject * VO,char * variant,SUMA_SurfaceViewer * sv)2560 SUMA_Boolean SUMA_Get_Slice_Pack(SUMA_VolumeObject *VO,
2561                                  char *variant, SUMA_SurfaceViewer *sv)
2562 {
2563    static char FuncName[]={"SUMA_Get_Slice_Pack"};
2564    SUMA_VOL_SAUX *VSaux=NULL;
2565    int i, ii1, ii0, iis, dim, *dims=NULL;
2566    SUMA_RENDERED_SLICE *rslc=NULL;
2567    SUMA_SLICE_FIELD *slc;
2568    SUMA_X_SurfCont *SurfCont = NULL;
2569    GLfloat slc_cen[6], scr_cen[6];
2570    SUMA_Boolean LocalHead = NOPE;
2571 
2572    SUMA_ENTRY;
2573 
2574    if (!VO || !(VSaux = SUMA_ADO_VSaux((SUMA_ALL_DO *)VO)) || !variant ||
2575        !(SurfCont = SUMA_ADO_Cont((SUMA_ALL_DO *)VO))) {
2576       SUMA_RETURN(NOPE);
2577    }
2578 
2579    switch (variant[0]) {
2580       case 'A':
2581          slc = SurfCont->Ax_slc;
2582          break;
2583       case 'S':
2584          slc = SurfCont->Sa_slc;
2585          break;
2586       case 'C':
2587          slc = SurfCont->Co_slc;
2588          break;
2589       case 'V': /* Volume Rendering */
2590          SUMA_RETURN(SUMA_GET_VR_Slice_Pack(VO, sv));
2591          break;
2592       default:
2593          slc = NULL;
2594          SUMA_S_Err("Bad variant");
2595          SUMA_RETURN(NOPE);
2596    }
2597 
2598    /* check for proper ordering direction
2599       (only needed if montage is in order ) */
2600    dim = SUMA_dset_tex_slice_corners_gui(VO->VE, 0, variant,
2601                               0,  /* most inferior slice */
2602                               NULL, NULL, slc_cen, NULL, 0);
2603    dim = SUMA_dset_tex_slice_corners_gui(VO->VE, 0, variant,
2604                               SUMA_VO_N_Slices(VO, variant)-1,/*superior++*/
2605                               NULL, NULL, slc_cen+3, NULL, 0);
2606    /* which of the two is closer to my nose? */
2607    SUMA_World2ScreenCoordsF(sv, 2, slc_cen, scr_cen, NULL, NOPE, NOPE);
2608 
2609 
2610    if (slc->mont_num > 1) {
2611       /* Form the slice planes. For now they are setup along the acquisition
2612          directions of the 0th VO->VE, but they need not be.
2613          Texture generation can then be carried out on any of the VEs
2614          regardless of their acquisition.                                */
2615       ii0 = slc->slice_num - ((slc->mont_num-1.0)/2.0)*slc->mont_inc;
2616       if (scr_cen[dim] > scr_cen[3+dim]) {
2617          iis = slc->mont_inc;
2618       } else {
2619          ii0 += (slc->mont_num-1)*slc->mont_inc;
2620          iis = -slc->mont_inc;
2621       }
2622       if (1) {
2623          i=0;
2624          while (i<slc->mont_num) {
2625             if (ii0 >= 0 && ii0 < SUMA_VO_N_Slices(VO, variant)) {
2626                SUMA_LH( "Top is closer bottom    %f %f %f --> scr z %f \n"
2627                         "              top       %f %f %f --> scr z %f \n",
2628                         slc_cen[0], slc_cen[1], slc_cen[2], scr_cen[2],
2629                         slc_cen[3], slc_cen[4], slc_cen[5], scr_cen[5]);
2630                rslc = (SUMA_RENDERED_SLICE *)
2631                               SUMA_malloc(sizeof(SUMA_RENDERED_SLICE));
2632                SUMA_dset_tex_slice_corners_gui(VO->VE, 0, variant, ii0,
2633                                                NULL, NULL, NULL, rslc->Eq, 0);
2634                rslc->slc_num = ii0; snprintf(rslc->variant, 15, "%s", variant);
2635                /* stick plane in list, last one rendered goes to top */
2636                SUMA_LH("Intersecting plane %f %f %f %f, on vol %s\n"
2637                        "(origin %f %f %f), variant %s",
2638                         rslc->Eq[0], rslc->Eq[1], rslc->Eq[2], rslc->Eq[3],
2639                         SUMA_VE_Headname(VO->VE,0),
2640                         VO->VE[0]->I2X[3][0],
2641                         VO->VE[0]->I2X[3][1],VO->VE[0]->I2X[3][2],
2642                         rslc->variant);
2643                dlist_ins_prev(VSaux->slcl, dlist_head(VSaux->slcl), rslc);
2644             }
2645             ++i;
2646             ii0 += iis;
2647          }
2648       } else {
2649          ii0 = slc->slice_num + slc->mont_num/2*slc->mont_inc;
2650          i=0;
2651          while (i<slc->mont_num) {
2652             if (ii0 >= 0 && ii0 < SUMA_VO_N_Slices(VO, variant)) {
2653                SUMA_LH( "Bot is closer bottom    %f %f %f --> scr z %f \n"
2654                         "              top       %f %f %f --> scr z %f \n",
2655                         slc_cen[0], slc_cen[1], slc_cen[2], scr_cen[2],
2656                         slc_cen[3], slc_cen[4], slc_cen[5], scr_cen[5]);
2657 
2658                rslc = (SUMA_RENDERED_SLICE *)
2659                               SUMA_malloc(sizeof(SUMA_RENDERED_SLICE));
2660                SUMA_dset_tex_slice_corners_gui(VO->VE, 0, variant, ii0,
2661                                                NULL, NULL, NULL, rslc->Eq, 0);
2662                rslc->slc_num = ii0; snprintf(rslc->variant, 15, "%s", variant);
2663                dlist_ins_prev(VSaux->slcl, dlist_head(VSaux->slcl), rslc);
2664             }
2665             ++i;
2666             ii0 -= slc->mont_inc;
2667          }
2668       }
2669    } else {
2670       SUMA_LH("Single slice %f %s", slc->slice_num, variant);
2671       rslc = (SUMA_RENDERED_SLICE *)
2672                            SUMA_malloc(sizeof(SUMA_RENDERED_SLICE));
2673       SUMA_dset_tex_slice_corners_gui(VO->VE, 0, variant, (int)slc->slice_num,
2674                                       NULL, NULL, NULL, rslc->Eq, 0);
2675       rslc->slc_num = (int)slc->slice_num;
2676       snprintf(rslc->variant, 15, "%s", variant);
2677       dlist_ins_prev(VSaux->slcl, dlist_head(VSaux->slcl), rslc);
2678 
2679    }
2680    SUMA_RETURN(YUP);
2681 }
2682 
SUMA_Count_All_VO_Textures(void)2683 int SUMA_Count_All_VO_Textures(void)
2684 {
2685    static char FuncName[]={"SUMA_Count_All_VO_Textures"};
2686    int i, j, c = 0;
2687    SUMA_ALL_DO *ado=NULL;
2688 
2689    for (i=0; i<SUMAg_N_DOv; ++i) {
2690       ado = iDO_ADO(i);
2691       if (ado->do_type == VO_type) {
2692          j=0;
2693          SUMA_VolumeObject *VO = (SUMA_VolumeObject *)ado;
2694          while (VO->VE && VO->VE[j]) {
2695             ++c;
2696             ++j;
2697          }
2698       }
2699    }
2700    return(c);
2701 }
2702 
2703 /*
2704    Check if texture has been loaded for a particular viewer
2705    NOTE: N_tex is only set is the texture was NOT loaded */
SUMA_SV_isTextureLoaded(SUMA_SurfaceViewer * sv,GLuint texName,int * N_tex)2706 SUMA_Boolean SUMA_SV_isTextureLoaded(SUMA_SurfaceViewer *sv,
2707                                      GLuint texName, int *N_tex)
2708 {
2709    static char FuncName[]={"SUMA_SV_isTextureLoaded"};
2710    int i=0;
2711 
2712    while (i<SUMA_MAX_DISPLAYABLE_OBJECTS && sv->LoadedTextures[i]!=-1) {
2713       if (sv->LoadedTextures[i] = (int)texName) return(YUP);
2714       ++i;
2715    }
2716    if (i == SUMA_MAX_DISPLAYABLE_OBJECTS && sv->LoadedTextures[i]!=-1) {
2717       SUMA_S_Warn("Looks like LoadedTextures is not plugged");
2718    }
2719    return(NOPE);
2720 }
2721 
SUMA_SV_Mark_Textures_Status(SUMA_SurfaceViewer * sv,char * MarkAs,SUMA_VolumeObject * VO,int j,int loadifneeded)2722 SUMA_Boolean SUMA_SV_Mark_Textures_Status(SUMA_SurfaceViewer *sv, char *MarkAs,
2723                                           SUMA_VolumeObject *VO, int j,
2724                                           int loadifneeded)
2725 {
2726    static char FuncName[]={"SUMA_SV_Mark_Textures_Status"};
2727    int N_tex = 0, i=0;
2728 
2729    SUMA_ENTRY;
2730 
2731    if (!sv || !MarkAs) {
2732       SUMA_RETURN(NOPE);
2733    }
2734    if (!strcmp(MarkAs, "unloaded_all")) {
2735       sv->LoadedTextures[0]=-1;
2736       SUMA_RETURN(YUP);
2737    } else if (!strcmp(MarkAs, "loaded_for_VO")) {
2738       if (!VO) SUMA_RETURN(NOPE);
2739       j = 0;
2740       while (VO->VE && VO->VE[j]) {
2741          if (!SUMA_SV_isTextureLoaded(sv, VO->VE[j]->texName[0], &N_tex)) {
2742             sv->LoadedTextures[N_tex] = VO->VE[j]->texName[0];
2743             sv->LoadedTextures[N_tex+1] = -1;
2744             if (loadifneeded) {
2745                SUMA_VE_LoadTexture(VO->VE, j);
2746             }
2747          }
2748          ++j;
2749       }
2750       SUMA_RETURN(YUP);
2751    } else if (!strcmp(MarkAs, "loaded_for_VO_one")) {
2752       if (!VO || j < 0 || !VO->VE || !VO->VE[j]) SUMA_RETURN(NOPE);
2753       if (!SUMA_SV_isTextureLoaded(sv, VO->VE[j]->texName[0], &N_tex)) {
2754          sv->LoadedTextures[N_tex] = VO->VE[j]->texName[0];
2755          sv->LoadedTextures[N_tex+1] = -1;
2756          if (loadifneeded) {
2757             SUMA_VE_LoadTexture(VO->VE, j);
2758          }
2759       }
2760       SUMA_RETURN(YUP);
2761    } else if (!strcmp(MarkAs, "loaded_all")) {
2762       SUMA_ALL_DO *ado=NULL;
2763       N_tex = 0;
2764       sv->LoadedTextures[N_tex]=-1;
2765       for (i=0; i<SUMAg_N_DOv; ++i) {
2766          ado = iDO_ADO(i);
2767          if (ado->do_type == VO_type) {
2768             VO = (SUMA_VolumeObject *)ado;
2769             j=0;
2770             while (VO->VE && VO->VE[j]) {
2771                sv->LoadedTextures[N_tex] = VO->VE[j]->texName[0];
2772                if (loadifneeded) {
2773                   SUMA_VE_LoadTexture(VO->VE, j);
2774                }
2775                ++N_tex;
2776                ++j;
2777             }
2778          }
2779       }
2780       sv->LoadedTextures[N_tex]=-1;
2781       SUMA_RETURN(YUP);
2782    } else {
2783       SUMA_S_Err("MarkAs %s not understood", MarkAs);
2784       SUMA_RETURN(NOPE);
2785    }
2786 }
2787 
SUMA_Draw_CIFTI_DO(SUMA_CIFTI_DO * CO,SUMA_SurfaceViewer * sv)2788 SUMA_Boolean SUMA_Draw_CIFTI_DO(SUMA_CIFTI_DO *CO,
2789                                SUMA_SurfaceViewer *sv)
2790 {
2791    static char FuncName[]={"SUMA_Draw_CIFTI_DO"};
2792    int i, pp=0;
2793    SUMA_ALL_DO *asdo=NULL;
2794    SUMA_Boolean LocalHead = NOPE;
2795 
2796    SUMA_ENTRY;
2797 
2798    SUMA_LH("We are not rendering CIFTI DOs directly anymore.\n"
2799       	   "Delete once we're sure we won't go down this route.\n"
2800 	   "CIFTI DOs are rendered by rendering their elementary\n"
2801 	   "DOs. And if you do go down that route, safer to refer\n"
2802 	   "to DOs by their IDs");
2803 
2804    SUMA_RETURN(YUP);
2805 
2806    for (i=0; i<CO->N_subdoms; ++i) {
2807       if (CO->subdoms_id[i]) {
2808       asdo = SUMA_CIFTI_subdom_ado(CO, i);
2809       switch (asdo->do_type) {
2810       	 case SO_type:
2811 	    if (!pp) SUMA_DrawMesh((SUMA_SurfaceObject *)asdo, sv);
2812 	    ++pp;
2813 	    break;
2814 	 case VO_type:
2815 	    SUMA_DrawVolumeDO((SUMA_VolumeObject*)asdo, sv);
2816 	    break;
2817 	 default:
2818 	    SUMA_S_Err("Not ready for subdom %d", asdo->do_type);
2819 	    break;
2820       }
2821       }
2822    }
2823 
2824 
2825    SUMA_RETURN(YUP);
2826 }
2827 
SUMA_DrawVolumeDO(SUMA_VolumeObject * VO,SUMA_SurfaceViewer * sv)2828 SUMA_Boolean SUMA_DrawVolumeDO(SUMA_VolumeObject *VO,
2829                                SUMA_SurfaceViewer *sv)
2830 {
2831    static char FuncName[]={"SUMA_DrawVolumeDO"};
2832    SUMA_ENTRY;
2833    if (!SUMA_DrawVolumeDO_slices(VO,sv)) {
2834       SUMA_S_Err("Failed to draw slices");
2835       SUMA_RETURN(NOPE);
2836    }
2837    if (!SUMA_DrawVolumeDO_3D(VO, sv)) {
2838       SUMA_S_Err("Failed to render volume");
2839       SUMA_RETURN(NOPE);
2840    }
2841    SUMA_RETURN(YUP);
2842 }
2843 
2844 
2845 /* Draw Volume Data, in slice mode for now */
SUMA_DrawVolumeDO_slices(SUMA_VolumeObject * VO,SUMA_SurfaceViewer * sv)2846 SUMA_Boolean SUMA_DrawVolumeDO_slices(SUMA_VolumeObject *VO,
2847                                     SUMA_SurfaceViewer *sv)
2848 {
2849    static char FuncName[]={"SUMA_DrawVolumeDO_slices"};
2850    int i = 0, k = 0, j=0, ive=0;
2851    float iq[4]={0, 0, 0, 0}, vo0[3], voN[3];
2852    static int ipass=0, iplane = 0;
2853    SUMA_RENDERED_SLICE *rslc=NULL;
2854    GLfloat tex_corn[18] ;
2855    GLfloat slc_corn[18] ;
2856    GLfloat rotationMatrix[4][4], rt[4][4];
2857    GLboolean gl_dt, gl_bl, gl_at;
2858    int ShowUnselected = 1, shmodel, nqd, ivelast;
2859    float tz = 0.0, I[3];
2860    static GLfloat init_rotationMatrix[4][4];
2861    static GLdouble dmatrix[16], init_mv_matrix[16];
2862    DListElmt *el=NULL;
2863    DList *st=NULL;
2864    SUMA_ALL_DO *ado = (SUMA_ALL_DO *)VO;
2865    SUMA_VOL_SAUX *VSaux = SUMA_ADO_VSaux(ado);
2866    SUMA_OVERLAYS *colp = NULL;
2867    SUMA_DSET *dset=NULL;
2868    float* nlt; /*JB: temporary node list, because I do not want to type
2869                   "VO->SOcut[0]->NodeList" over and over...*/
2870    SUMA_Boolean LastTextureOnCutPlane=YUP;
2871    SUMA_ATRANS_MODES trmode=SATM_ViewerDefault;
2872    SUMA_Boolean LocalHead = NOPE;
2873 
2874    SUMA_ENTRY;
2875 
2876    if (!VO) SUMA_RETURN(NOPE);
2877    if (!sv) sv = &(SUMAg_SVv[0]);
2878    if (!(colp = SUMA_ADO_CurColPlane(ado))) {
2879       SUMA_S_Err("Need colp here");
2880       SUMA_RETURN(NOPE);
2881    }
2882    if (sv->DO_PickMode) {
2883       SUMA_LH("No need to draw volume in DO_PickMode");
2884       SUMA_RETURN(YUP);
2885    }
2886 
2887    if (!VO->Show) SUMA_RETURN(YUP);
2888 
2889    if (!VO->SOcut || !VO->SOcut[0]) SUMA_VO_InitCutPlanes(VO);
2890 
2891    if (!VO->VE || !VO->VE[0]) {
2892       SUMA_S_Err("No elements?");
2893       SUMA_RETURN(NOPE);
2894    }
2895 
2896    /* setup slices per VO's standards, dictated by the first VE */
2897    dset = SUMA_VO_dset(VO);
2898    if (!(dset = SUMA_VO_dset(VO))) {
2899       SUMA_S_Err("No dset?");
2900       SUMA_RETURN(NOPE);
2901    }
2902 
2903    /* make sure all textures are ready */
2904    ive = 0;
2905    while (VO->VE && VO->VE[ive]) {
2906       if (!VO->VE[ive]->texName) {
2907          if (!SUMA_CreateGL3DTexture(VO)) {
2908             SUMA_S_Err("Failed to create texture");
2909             SUMA_RETURN(NOPE);
2910          }
2911          if (!SUMA_SV_Mark_Textures_Status(sv, "loaded_for_VO_one", VO, ive,0)) {
2912             SUMA_S_Err("Failed to mark texture as loaded");
2913             SUMA_RETURN(NOPE);
2914          }
2915       } else {
2916          /* Need to check if textures need to be reloaded.
2917          The actual loading needs to be redone only when a viewer has been
2918          closed and then open again. Not sure why that must be done but without
2919          a new load the textures go away after the viewer is reopened */
2920          if (!SUMA_SV_Mark_Textures_Status(sv, "loaded_for_VO_one", VO, ive, 1)){
2921             SUMA_S_Err("Failed to check or reload texture");
2922             SUMA_RETURN(NOPE);
2923          }
2924       }
2925       ++ive;
2926    }
2927 
2928 
2929    if (!SUMA_GLStateTrack( "new", &st, FuncName, NULL, NULL)) {
2930       SUMA_S_Err("Failed to create tracking list");
2931       SUMA_RETURN(NOPE);
2932    }
2933 
2934    if (sv->PolyMode != SRM_Fill) {
2935       /* fill it up */
2936       glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
2937    }
2938 
2939    /* Now we need to draw the sucker */
2940    SUMA_CHECK_GL_ERROR("OpenGL Error pre texture");
2941    glEnable(GL_TEXTURE_3D);
2942    if (!(gl_at = glIsEnabled(GL_ALPHA_TEST))) glEnable(GL_ALPHA_TEST);
2943    if (colp->AlphaThresh == 0.0f) glAlphaFunc(GL_ALWAYS, colp->AlphaThresh);
2944    else glAlphaFunc(GL_GREATER, colp->AlphaThresh);
2945                               /* Thresholded voxels alphas are set to 0 */
2946    glTexEnvf(  GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE,
2947                VO->TexEnvMode   ); /* what happens if
2948                               there is color already on a vertex (I would likely
2949                               not need this for 3D textures...*/
2950 
2951    gl_dt = glIsEnabled(GL_DEPTH_TEST);
2952    gl_bl = glIsEnabled(GL_BLEND);
2953 
2954    if (!(gl_dt)) glEnable(GL_DEPTH_TEST);
2955 
2956    glGetIntegerv(GL_SHADE_MODEL, &shmodel);
2957    if (shmodel != GL_FLAT)
2958       glShadeModel(GL_FLAT);
2959 
2960    trmode = VSaux->TransMode;
2961    if (trmode == SATM_ViewerDefault) {
2962       if ((trmode = SUMA_TransMode2ATransMode(sv->TransMode))
2963             <= SATM_ViewerDefault || trmode >= SATM_N_TransModes) {
2964          SUMA_S_Warn("Bad trans mode change from %d to %d",
2965                       sv->TransMode,trmode);
2966          trmode =  SATM_0;
2967       }
2968    }
2969    if (trmode <= SATM_ViewerDefault || trmode > SATM_N_TransModes) {
2970       SUMA_S_Warn("Bad trmode %d", trmode);
2971       trmode =  SATM_0;
2972    }
2973    if (trmode == SATM_ALPHA) {
2974       /* Setup blending options */
2975       /* See note below, but might be OK with few slices */
2976       if ((!gl_bl)) glEnable(GL_BLEND);
2977       glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
2978    } else {
2979       /* The safe way. Can't blend properly when showing
2980          slices, particularly stacks and
2981          multiple planes. Becomes a royal pain
2982          to render all in proper order*/
2983       if ((gl_bl)) glDisable(GL_BLEND);
2984       if (trmode > SATM_0) {
2985          SUMA_LHv("ATrans Mode %d\n", trmode );
2986          SUMA_SET_GL_TRANS_MODE(SUMA_ATransMode2TransMode(trmode), st, -1);
2987       }
2988    }
2989    /* empty list of rendered slices */
2990    while ((el = dlist_head(VSaux->slcl))) {
2991       dlist_remove(VSaux->slcl, el, (void **)&rslc);
2992       SUMA_Free_SliceListDatum((void *)rslc);
2993    }
2994 
2995    if (VSaux->ShowSaSlc && !SUMA_Get_Slice_Pack(VO, "Sa", sv)) {
2996       SUMA_S_Err("Failed to get Ax slice pack");
2997    }
2998 
2999    if (VSaux->ShowAxSlc && !SUMA_Get_Slice_Pack(VO, "Ax", sv)) {
3000       SUMA_S_Err("Failed to get Ax slice pack");
3001    }
3002 
3003 
3004    if (VSaux->ShowCoSlc && !SUMA_Get_Slice_Pack(VO, "Co", sv)) {
3005       SUMA_S_Err("Failed to get Ax slice pack");
3006    }
3007 
3008    SUMA_LH("Have %d slices to render", dlist_size(VSaux->slcl));
3009    SUMA_CHECK_GL_ERROR("OpenGL Error pre setup");
3010 
3011 
3012    /* Generate the textures for all VEs. Note, no blending of textures
3013       is done at the moment.
3014       If you have multiple VEs, you will need to render one slice at a
3015       time from each of the volumes. Blend them in the accumulate buffer
3016       and then return. Textures will need to be bound/unbound for each
3017       VE, each slice. Clunky but I won't bother with it until it proves
3018       too slow.*/
3019    if (SUMA_VO_NumVE(VO) > 1) {
3020       SUMA_S_Warn("Not ready to deal with multiple textures quite yet");
3021    }
3022 
3023    SUMA_LH("Have %d slices", dlist_size(VSaux->slcl));
3024    ivelast=-1;
3025    if ((el = dlist_tail(VSaux->slcl))) {
3026       do {
3027          rslc = (SUMA_RENDERED_SLICE *)el->data;
3028 
3029          ive = 0;
3030          while (VO->VE && VO->VE[ive]) {
3031             if (ive != ivelast) {
3032                glBindTexture(GL_TEXTURE_3D, VO->VE[ive]->texName[0]);
3033                ivelast = ive;
3034             }                                 /* make texName be current */
3035             /* compute plane intersection with VE[ive] */
3036             if ((nqd = SUMA_PlaneBoxIntersect( sv->GVS[sv->StdView].ViewFrom,
3037                                              rslc->Eq, VO->VE[ive]->bcorners,
3038                                              slc_corn)) > 2) {
3039                SUMA_LH("Have plane %f %f %f %f, %d pts on vol %s",
3040                            rslc->Eq[0], rslc->Eq[1], rslc->Eq[2], rslc->Eq[3],
3041                            nqd,  SUMA_VE_Headname(VO->VE,ive));
3042                glBegin(GL_POLYGON);
3043                   for (k=0; k<6; ++k) { /* draw all 6 points always, even
3044                                            when there are repetitions.
3045                                            Don't bother trimming to unique
3046                                            set unless this causes trouble */
3047                      /* change mm (edge coordinate to texture coords) */
3048                      AFF44_MULT_I(tex_corn, VO->VE[ive]->X2I, (slc_corn+3*k));
3049                      /* offset indices because slc_corn is on edge */
3050                      if (tex_corn[0] < 0) tex_corn[0] = 0;
3051                      else if (tex_corn[0] > VO->VE[ive]->Ni-1)
3052                         tex_corn[0] = VO->VE[ive]->Ni-1;
3053                      if (tex_corn[1] < 0) tex_corn[1] = 0;
3054                      else if (tex_corn[1] > VO->VE[ive]->Nj-1)
3055                         tex_corn[1] = VO->VE[ive]->Nj-1;
3056                      if (tex_corn[2] < 0) tex_corn[2] = 0;
3057                      else if (tex_corn[2] > VO->VE[ive]->Nk-1)
3058                         tex_corn[2] = VO->VE[ive]->Nk-1;
3059                      tex_corn[0] /= (float)(VO->VE[ive]->Ni-1);
3060                      tex_corn[1] /= (float)(VO->VE[ive]->Nj-1);
3061                      tex_corn[2] /= (float)(VO->VE[ive]->Nk-1);
3062                      glTexCoord3f(tex_corn[0],
3063                                   tex_corn[1], tex_corn[2]);
3064                            /* this one is affected by the Texture MatrixMode */
3065                      glVertex3f(slc_corn[3*k],
3066                                 slc_corn[3*k+1], slc_corn[3*k+2]);
3067                            /* this one is affected by the Modelview matrixMode*/
3068                   }
3069                glEnd();
3070             }
3071             if (ive > 0) {
3072                SUMA_S_Warn("Add blending here");
3073                /* Here is where you blend slice from this VE with the previous
3074                   This should work just fine as is actually, no need to blend
3075                   separately unless doing overlay on top always. In that case,
3076                   render each VE separately then blend results across VEs
3077                */
3078             }
3079             ++ive;
3080          }
3081 
3082          if (el != dlist_head(VSaux->slcl)) el = dlist_prev(el);
3083          else el = NULL;
3084       } while (el);
3085    }
3086    SUMA_CHECK_GL_ERROR("OpenGL Error ddd");
3087 
3088    glFlush();
3089    if (!gl_at) glDisable(GL_ALPHA_TEST);
3090 
3091    glDisable(GL_TEXTURE_3D);
3092 
3093    if (shmodel != GL_FLAT) glShadeModel(shmodel);
3094 
3095    if (sv->PolyMode != SRM_Fill) {/* set fill mode back */
3096       SUMA_SET_GL_RENDER_MODE(sv->PolyMode);
3097    }
3098 
3099    if (gl_dt) glEnable(GL_DEPTH_TEST);
3100    else glDisable(GL_DEPTH_TEST);
3101    if (gl_bl) glEnable(GL_BLEND);
3102    else glDisable(GL_BLEND);
3103 
3104    /* Now for the highlight */
3105    if (SUMA_SV_GetShowSelectedFaceSet(sv)) {
3106       int selslice = -1;
3107       float nlt[12];
3108       char variant[8]={"notset"};
3109       GLfloat No_Color[4] = { 0.0, 0.0, 0.0, 0.0 };
3110       selslice = SUMA_VO_SelectedSlice(VO, variant, nlt);
3111       if (selslice >= 0 && variant[0] != 'V') {
3112          SUMA_LH("Drawing %s Slice %d Selection Contour\n"
3113                  "%f %f %f --> %f %f %f ...\n",
3114                  variant, selslice,
3115                  nlt[0],nlt[1],nlt[2], nlt[3],nlt[4],nlt[5]);
3116          glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, No_Color);
3117          glColorMaterial(GL_FRONT, GL_EMISSION);
3118          glEnable(GL_COLOR_MATERIAL);
3119          glColor4f(0.25, 0.25, 0.25, 1.0);
3120          glBegin(GL_LINE_LOOP);
3121             glVertex3f( nlt[0],nlt[1],nlt[2] );
3122             glVertex3f( nlt[3],nlt[4],nlt[5] );
3123             glVertex3f( nlt[6],nlt[7],nlt[8] );
3124             glVertex3f( nlt[9],nlt[10],nlt[11] );
3125          glEnd();
3126          glColor4f(0, 0, 0, 0); /* Kill emission,        ZSS. March 2014
3127                         or risk hurting functions that expect it to be off. */
3128          glDisable(GL_COLOR_MATERIAL);
3129       } else {
3130          SUMA_LH("Either no selection or failed to find slice, "
3131                  "or VR intersection (variant=%s)", variant);
3132       }
3133    } else {
3134       SUMA_LH("Do not show selected faceset");
3135    }
3136 
3137    SUMA_LH("Undoing state changes, should fold ones above in here someday");
3138    SUMA_GLStateTrack("r", &st, FuncName, NULL, NULL);
3139 
3140    SUMA_RETURN(YUP);
3141 }
3142 
3143 
3144 /* Draw Volume Data, in 3D this time */
SUMA_DrawVolumeDO_3D(SUMA_VolumeObject * VO,SUMA_SurfaceViewer * sv)3145 SUMA_Boolean SUMA_DrawVolumeDO_3D(SUMA_VolumeObject *VO,
3146                                     SUMA_SurfaceViewer *sv)
3147 {
3148    static char FuncName[]={"SUMA_DrawVolumeDO_3D"};
3149    int i = 0, k = 0, j=0, ive=0;
3150    float iq[4]={0, 0, 0, 0}, vo0[3], voN[3];
3151    static int ipass=0, iplane = 0;
3152    SUMA_RENDERED_SLICE *rslc=NULL;
3153    GLfloat tex_corn[18] ;
3154    GLfloat slc_corn[18] ;
3155    GLfloat rotationMatrix[4][4], rt[4][4];
3156    GLboolean gl_dt, gl_bl, gl_at;
3157    int ShowUnselected = 1, shmodel, nqd, ivelast;
3158    DListElmt *el=NULL;
3159    DList *st=NULL;
3160    float tz = 0.0, I[3];
3161    static GLfloat init_rotationMatrix[4][4];
3162    static GLdouble dmatrix[16], init_mv_matrix[16];
3163    SUMA_ALL_DO *ado = (SUMA_ALL_DO *)VO;
3164    SUMA_VOL_SAUX *VSaux = SUMA_ADO_VSaux(ado);
3165    SUMA_OVERLAYS *colp = NULL;
3166    SUMA_DSET *dset=NULL;
3167    float* nlt; /*JB: temporary node list, because I do not want to type
3168                   "VO->SOcut[0]->NodeList" over and over...*/
3169    SUMA_Boolean LastTextureOnCutPlane=YUP;
3170    SUMA_ATRANS_MODES trmode=SATM_ViewerDefault;
3171    SUMA_Boolean LocalHead = NOPE;
3172 
3173    SUMA_ENTRY;
3174 
3175    if (!VO) SUMA_RETURN(NOPE);
3176    if (!sv) sv = &(SUMAg_SVv[0]);
3177    if (!(colp = SUMA_ADO_CurColPlane(ado))) {
3178       SUMA_S_Err("Need colp here");
3179       SUMA_RETURN(NOPE);
3180    }
3181    if (sv->DO_PickMode) {
3182       SUMA_LH("No need to draw volume in DO_PickMode");
3183       SUMA_RETURN(YUP);
3184    }
3185 
3186    if (!VO->Show) SUMA_RETURN(YUP);
3187 
3188    if (!VSaux->ShowVrSlc) SUMA_RETURN(YUP);
3189 
3190    if (!SUMA_GLStateTrack( "new", &st, FuncName, NULL, NULL)) {
3191       SUMA_S_Err("Failed to create tracking list");
3192       SUMA_RETURN(NOPE);
3193    }
3194 
3195 
3196    if (sv->PolyMode != SRM_Fill) {
3197       /* fill it up */
3198       glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
3199    }
3200 
3201    if (!VO->SOcut || !VO->SOcut[0]) SUMA_VO_InitCutPlanes(VO);
3202 
3203    /* Now we need to draw the sucker */
3204    SUMA_CHECK_GL_ERROR("OpenGL Error pre texture");
3205    glEnable(GL_TEXTURE_3D);
3206    /* You need a separate control for this ALPHA_TEST, perhaps.
3207       Enabling it, even at 0.1 threshold, causes ugly slice striping
3208       artifacts. It is possible we might need it for something later on...*/
3209    if (!(gl_at = glIsEnabled(GL_ALPHA_TEST))) glEnable(GL_ALPHA_TEST);
3210    #if 0
3211    if (colp->AlphaThresh == 0.0f) glAlphaFunc(GL_ALWAYS, colp->AlphaThresh);
3212    else glAlphaFunc(GL_GREATER, colp->AlphaThresh);
3213                               /* Thresholded voxels alphas are set to 0 */
3214    #else
3215    glAlphaFunc(GL_GREATER, 0.01); /* Might want to add a separate AlphaThresh
3216                                      for this. Usually you want it very small,
3217                                      or else you'll get ugly striping artifacts.
3218                                      But if you use the slice's alpha threshold,
3219                                      you'll need to reduce it low enough to make
3220                                      the volume look nice but that give the
3221                                      slices an ugly rim. */
3222    #endif
3223    glTexEnvf(  GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE,
3224                VO->TexEnvMode   ); /* what happens if
3225                               there is color already on a vertex (I would likely
3226                               not need this for 3D textures...*/
3227 
3228    gl_dt = glIsEnabled(GL_DEPTH_TEST);
3229    gl_bl = glIsEnabled(GL_BLEND);
3230 
3231    if (!(gl_dt)) glEnable(GL_DEPTH_TEST);
3232    if (!(gl_bl)) glEnable(GL_BLEND);
3233 
3234    { /* Transparency games */
3235    /* Beware of the trickery of transparency.
3236    For cheese cloth transp, rendering multiple objects in a transparent
3237    fashion is tricky.
3238    The bitmask autoshifting is not enough when many objects are in use.
3239    The stricking case is when you have two surfs and a volume, at 50% transp.
3240    If you shift for each object, then the first surface will get obscured by
3241    the volume (third object) because the twice shifted pattern will overlap
3242    with that of the first surf. One solution, now implemented, is to use the same
3243    mask shifting per class of objects. But then surfs won't be transparent
3244    between them (not a big deal). The other option is to set the transparency
3245    differently for different objects (i.e. not using the viewer-wide
3246    transparency). Other options might help under certain cases (see
3247    comment for function SUMA_StippleMask_shift().     ZSS March 13 2014    */
3248    trmode = VSaux->TransMode;
3249    if (trmode == SATM_ViewerDefault) {
3250       if ((trmode = SUMA_TransMode2ATransMode(sv->TransMode))
3251             <= SATM_ViewerDefault || trmode >= SATM_N_TransModes) {
3252          SUMA_S_Warn("Bad trans mode change from %d to %d",
3253                       sv->TransMode,trmode);
3254          trmode =  SATM_0;
3255       }
3256    }
3257    if (trmode <= SATM_ViewerDefault || trmode > SATM_N_TransModes) {
3258       SUMA_S_Warn("Bad trmode %d", trmode);
3259       trmode =  SATM_0;
3260    }
3261    if (trmode == SATM_ALPHA) {
3262       /* Setup blending options, override initial setting above */
3263       glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
3264    } else {
3265       /* The safe way. */
3266       if (trmode > SATM_0) {
3267          SUMA_LHv("ATrans Mode %d\n", trmode );
3268          SUMA_SET_GL_TRANS_MODE(SUMA_ATransMode2TransMode(trmode), st, VO_type);
3269       }
3270    }
3271    }
3272 
3273    if (!VO->VE || !VO->VE[0]) {
3274       SUMA_S_Err("No elements?");
3275       SUMA_RETURN(NOPE);
3276    }
3277 
3278    /* setup slices per VO's standards, dictated by the first VE */
3279    dset = SUMA_VO_dset(VO);
3280    if (!(dset = SUMA_VO_dset(VO))) {
3281       SUMA_S_Err("No dset?");
3282       SUMA_RETURN(NOPE);
3283    }
3284 
3285    /* make sure all textures are ready */
3286    ive = 0;
3287    while (VO->VE && VO->VE[ive]) {
3288       if (!VO->VE[ive]->texName) {
3289          if (!SUMA_CreateGL3DTexture(VO)) {
3290             SUMA_S_Err("Failed to create texture");
3291             SUMA_RETURN(NOPE);
3292          }
3293       }
3294       ++ive;
3295    }
3296    /* empty list of rendered slices */
3297    while ((el = dlist_head(VSaux->vrslcl))) {
3298       dlist_remove(VSaux->vrslcl, el, (void **)&rslc);
3299       SUMA_Free_SliceListDatum((void *)rslc);
3300    }
3301 
3302    /* Now create list of 3D slices */
3303 
3304    if (VSaux->ShowVrSlc && !SUMA_Get_Slice_Pack(VO, "Vr", sv)) {
3305       SUMA_S_Err("Failed to get VR slice pack");
3306    }
3307 
3308 
3309    SUMA_LH("Have %d slices to render", dlist_size(VSaux->vrslcl));
3310    /* Setup blending options */
3311       glGetIntegerv(GL_SHADE_MODEL, &shmodel);
3312       if (shmodel != GL_FLAT)
3313          glShadeModel(GL_FLAT);
3314 	   glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
3315       SUMA_CHECK_GL_ERROR("OpenGL Error pre setup");
3316 
3317 
3318    /* Generate the textures for all VEs. Note, no blending of textures
3319       is done at the moment.
3320       If you have multiple VEs, you will need to render one slice at a
3321       time from each of the volumes. Blend them in the accumulate buffer
3322       and then return. Textures will need to be bound/unbound for each
3323       VE, each slice. Clunky but I won't bother with it until it proves
3324       too slow.*/
3325    if (SUMA_VO_NumVE(VO) > 1) {
3326       SUMA_S_Warn("Not ready to deal with multiple textures quite yet");
3327    }
3328 
3329    SUMA_LH("Have %d slices", dlist_size(VSaux->vrslcl));
3330    ivelast=-1;
3331    if ((el = dlist_tail(VSaux->vrslcl))) {
3332       do {
3333          rslc = (SUMA_RENDERED_SLICE *)el->data;
3334 
3335          ive = 0;
3336          while (VO->VE && VO->VE[ive]) {
3337             if (ive != ivelast) {
3338                glBindTexture(GL_TEXTURE_3D, VO->VE[ive]->texName[0]);
3339                ivelast = ive;
3340             }                                 /* make texName be current */
3341             /* compute plane intersection with VE[ive] */
3342             if ((nqd = SUMA_PlaneBoxIntersect( sv->GVS[sv->StdView].ViewFrom,
3343                                              rslc->Eq, VO->VE[ive]->bcorners,
3344                                              slc_corn)) > 2) {
3345                SUMA_LH("Have plane %f %f %f %f, %d pts on vol %s",
3346                            rslc->Eq[0], rslc->Eq[1], rslc->Eq[2], rslc->Eq[3],
3347                            nqd,  SUMA_VE_Headname(VO->VE,ive));
3348                glBegin(GL_POLYGON);
3349                   for (k=0; k<6; ++k) { /* draw all 6 points always, even
3350                                            when there are repetitions.
3351                                            Don't bother trimming to unique
3352                                            set unless this causes trouble */
3353                      /* change mm (edge coordinate to texture coords) */
3354                      AFF44_MULT_I(tex_corn, VO->VE[ive]->X2I, (slc_corn+3*k));
3355                      /* offset indices because slc_corn is on edge */
3356                      if (tex_corn[0] < 0) tex_corn[0] = 0;
3357                      else if (tex_corn[0] > VO->VE[ive]->Ni-1)
3358                         tex_corn[0] = VO->VE[ive]->Ni-1;
3359                      if (tex_corn[1] < 0) tex_corn[1] = 0;
3360                      else if (tex_corn[1] > VO->VE[ive]->Nj-1)
3361                         tex_corn[1] = VO->VE[ive]->Nj-1;
3362                      if (tex_corn[2] < 0) tex_corn[2] = 0;
3363                      else if (tex_corn[2] > VO->VE[ive]->Nk-1)
3364                         tex_corn[2] = VO->VE[ive]->Nk-1;
3365                      tex_corn[0] /= (float)(VO->VE[ive]->Ni-1);
3366                      tex_corn[1] /= (float)(VO->VE[ive]->Nj-1);
3367                      tex_corn[2] /= (float)(VO->VE[ive]->Nk-1);
3368                      glTexCoord3f(tex_corn[0],
3369                                   tex_corn[1], tex_corn[2]);
3370                            /* this one is affected by the Texture MatrixMode */
3371                      glVertex3f(slc_corn[3*k],
3372                                 slc_corn[3*k+1], slc_corn[3*k+2]);
3373                            /* this one is affected by the Modelview matrixMode*/
3374                   }
3375                glEnd();
3376             }
3377             if (ive > 0) {
3378                SUMA_S_Warn("Add blending here");
3379                /* Here is where you blend slice from this VE with the previous
3380                   This should work just fine as is actually, no need to blend
3381                   separately unless doing overlay on top always. In that case,
3382                   render each VE separately then blend results across VEs
3383                */
3384             }
3385             ++ive;
3386          }
3387 
3388          if (el != dlist_head(VSaux->vrslcl)) el = dlist_prev(el);
3389          else el = NULL;
3390       } while (el);
3391    }
3392    SUMA_CHECK_GL_ERROR("OpenGL Error ddd");
3393 
3394    glFlush();
3395    if (!gl_at) glDisable(GL_ALPHA_TEST);
3396 
3397    glDisable(GL_TEXTURE_3D);
3398 
3399    if (shmodel != GL_FLAT) glShadeModel(shmodel);
3400 
3401    if (sv->PolyMode != SRM_Fill) {/* set fill mode back */
3402       SUMA_SET_GL_RENDER_MODE(sv->PolyMode);
3403    }
3404 
3405    if (gl_dt) glEnable(GL_DEPTH_TEST);
3406    else glDisable(GL_DEPTH_TEST);
3407    if (gl_bl) glEnable(GL_BLEND);
3408    else glDisable(GL_BLEND);
3409 
3410    /* This method should eventually replace the individual tests done above...*/
3411    SUMA_LH("Undoing state changes, should fold ones above in here someday");
3412    SUMA_GLStateTrack("r", &st, FuncName, NULL, NULL);
3413 
3414    SUMA_RETURN(YUP);
3415 }
3416 
3417 /* Draw Volume Data, exp version */
SUMA_DrawVolumeDO_exp(SUMA_VolumeObject * VO,SUMA_SurfaceViewer * sv)3418 SUMA_Boolean SUMA_DrawVolumeDO_exp(SUMA_VolumeObject *VO, SUMA_SurfaceViewer *sv)
3419 {
3420    static char FuncName[]={"SUMA_DrawVolumeDO_exp"};
3421    int i = 0, k = 0, j=0, ive=0;
3422    float iq[4]={0, 0, 0, 0}, vo0[3], voN[3];
3423    static int ipass=0, iplane = 0;
3424    SUMA_RENDERED_SLICE *rslc=NULL;
3425    GLfloat tex_corn[18] ;
3426    GLfloat slc_corn[18] ;
3427    GLfloat rotationMatrix[4][4], rt[4][4];
3428    GLboolean gl_dt, gl_bl;
3429    int ShowUnselected = 1, shmodel, nqd, ivelast;
3430    float tz = 0.0, I[3];
3431    static GLfloat init_rotationMatrix[4][4];
3432    static GLdouble dmatrix[16], init_mv_matrix[16];
3433    DListElmt *el=NULL;
3434    SUMA_ALL_DO *ado = (SUMA_ALL_DO *)VO;
3435    SUMA_VOL_SAUX *VSaux = SUMA_ADO_VSaux(ado);
3436    SUMA_DSET *dset=NULL;
3437    float* nlt; /*JB: temporary node list, because I do not want to type
3438                   "VO->SOcut[0]->NodeList" over and over...*/
3439    SUMA_Boolean LastTextureOnCutPlane=YUP;
3440    SUMA_Boolean LocalHead = NOPE;
3441 
3442    SUMA_ENTRY;
3443 
3444    #ifndef GL_VERSION_2_0
3445       /* GL must be old */
3446       SUMA_S_Err("Open GL < 2.0, glWindowPos2s() not yet supported (on all machines)");
3447       SUMA_RETURN(NOPE);
3448    #else
3449    if (!VO) SUMA_RETURN(NOPE);
3450    if (!sv) sv = &(SUMAg_SVv[0]);
3451 
3452    if (sv->DO_PickMode) {
3453       SUMA_LH("No need to draw volume in DO_PickMode");
3454       SUMA_RETURN(YUP);
3455    }
3456 
3457    if (!VO->Show) SUMA_RETURN(YUP);
3458 
3459    if (sv->PolyMode != SRM_Fill) {
3460       /* fill it up */
3461       glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
3462    }
3463 
3464    if (!VO->SOcut || !VO->SOcut[0]) SUMA_VO_InitCutPlanes(VO);
3465 
3466    /* Now we need to draw the sucker */
3467    SUMA_CHECK_GL_ERROR("OpenGL Error pre texture");
3468    glEnable(GL_TEXTURE_3D);
3469 
3470    glTexEnvf(  GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE,
3471                VO->TexEnvMode   ); /* what happens if
3472                               there is color already on a vertex (I would likely
3473                               not need this for 3D textures...*/
3474 
3475    gl_dt = glIsEnabled(GL_DEPTH_TEST);
3476    gl_bl = glIsEnabled(GL_BLEND);
3477 
3478    if (!(gl_dt)) glEnable(GL_DEPTH_TEST);
3479    if (!(gl_bl)) glEnable(GL_BLEND);
3480 
3481    if (!VO->VE || !VO->VE[0]) {
3482       SUMA_S_Err("No elements?");
3483       SUMA_RETURN(NOPE);
3484    }
3485 
3486    /* setup slices per VO's standards, dictated by the first VE */
3487    dset = SUMA_VO_dset(VO);
3488    if (!(dset = SUMA_VO_dset(VO))) {
3489       SUMA_S_Err("No dset?");
3490       SUMA_RETURN(NOPE);
3491    }
3492 
3493    /* make sure all textures are ready */
3494    ive = 0;
3495    while (VO->VE && VO->VE[ive]) {
3496       if (!VO->VE[ive]->texName) {
3497          if (!SUMA_CreateGL3DTexture(VO)) {
3498             SUMA_S_Err("Failed to create texture");
3499             SUMA_RETURN(NOPE);
3500          }
3501       }
3502       ++ive;
3503    }
3504    /* empty list of rendered slices */
3505    while ((el = dlist_head(VSaux->slcl))) {
3506       dlist_remove(VSaux->slcl, el, (void **)&rslc);
3507       SUMA_Free_SliceListDatum((void *)rslc);
3508    }
3509 
3510    if (VSaux->ShowSaSlc && !SUMA_Get_Slice_Pack(VO, "Sa", sv)) {
3511       SUMA_S_Err("Failed to get Ax slice pack");
3512    }
3513 
3514    if (VSaux->ShowAxSlc && !SUMA_Get_Slice_Pack(VO, "Ax", sv)) {
3515       SUMA_S_Err("Failed to get Ax slice pack");
3516    }
3517 
3518 
3519    if (VSaux->ShowCoSlc && !SUMA_Get_Slice_Pack(VO, "Co", sv)) {
3520       SUMA_S_Err("Failed to get Ax slice pack");
3521    }
3522 
3523    SUMA_LH("Have %d slices to render", dlist_size(VSaux->slcl));
3524    /* Setup blending options */
3525       glGetIntegerv(GL_SHADE_MODEL, &shmodel);
3526       if (shmodel != GL_FLAT)
3527          glShadeModel(GL_FLAT);
3528 	   glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
3529       SUMA_CHECK_GL_ERROR("OpenGL Error pre setup");
3530 
3531 
3532    /* Generate the textures for all VEs. Note, no blending of textures
3533       is done at the moment.
3534       If you have multiple VEs, you will need to render one slice at a
3535       time from each of the volumes. Blend them in the accumulate buffer
3536       and then return. Textures will need to be bound/unbound for each
3537       VE, each slice. Clunky but I won't bother with it until it proves
3538       too slow.*/
3539    if (SUMA_VO_NumVE(VO) > 1) {
3540       SUMA_S_Warn("Not ready to deal with multiple textures quite yet");
3541    }
3542 
3543    SUMA_LH("Have %d slices", dlist_size(VSaux->slcl));
3544    ivelast=-1;
3545    if ((el = dlist_tail(VSaux->slcl))) {
3546       do {
3547          rslc = (SUMA_RENDERED_SLICE *)el->data;
3548 
3549          ive = 0;
3550          while (VO->VE && VO->VE[ive]) {
3551             if (ive != ivelast) {
3552                glBindTexture(GL_TEXTURE_3D, VO->VE[ive]->texName[0]);
3553                ivelast = ive;
3554             }                                 /* make texName be current */
3555             /* compute plane intersection with VE[ive] */
3556             if ((nqd = SUMA_PlaneBoxIntersect( sv->GVS[sv->StdView].ViewFrom,
3557                                              rslc->Eq, VO->VE[ive]->bcorners,
3558                                              slc_corn)) > 2) {
3559                SUMA_LH("Have plane %f %f %f %f, %d pts on vol %s",
3560                            rslc->Eq[0], rslc->Eq[1], rslc->Eq[2], rslc->Eq[3],
3561                            nqd,  SUMA_VE_Headname(VO->VE,ive));
3562                glBegin(GL_POLYGON);
3563                   for (k=0; k<6; ++k) { /* draw all 6 points always, even
3564                                            when there are repetitions.
3565                                            Don't bother trimming to unique
3566                                            set unless this causes trouble */
3567                      /* change mm (edge coordinate to texture coords) */
3568                      AFF44_MULT_I(tex_corn, VO->VE[ive]->X2I, (slc_corn+3*k));
3569                      /* offset indices because slc_corn is on edge */
3570                      if (tex_corn[0] < 0) tex_corn[0] = 0;
3571                      else if (tex_corn[0] > VO->VE[ive]->Ni-1)
3572                         tex_corn[0] = VO->VE[ive]->Ni-1;
3573                      if (tex_corn[1] < 0) tex_corn[1] = 0;
3574                      else if (tex_corn[1] > VO->VE[ive]->Nj-1)
3575                         tex_corn[1] = VO->VE[ive]->Nj-1;
3576                      if (tex_corn[2] < 0) tex_corn[2] = 0;
3577                      else if (tex_corn[2] > VO->VE[ive]->Nk-1)
3578                         tex_corn[2] = VO->VE[ive]->Nk-1;
3579                      tex_corn[0] /= (float)(VO->VE[ive]->Ni-1);
3580                      tex_corn[1] /= (float)(VO->VE[ive]->Nj-1);
3581                      tex_corn[2] /= (float)(VO->VE[ive]->Nk-1);
3582                      glTexCoord3f(tex_corn[0],
3583                                   tex_corn[1], tex_corn[2]);
3584                            /* this one is affected by the Texture MatrixMode */
3585                      glVertex3f(slc_corn[3*k],
3586                                 slc_corn[3*k+1], slc_corn[3*k+2]);
3587                            /* this one is affected by the Modelview matrixMode*/
3588                   }
3589                glEnd();
3590             }
3591             #if 0
3592             {
3593                GLvoid *dbuf, *dlum;
3594                int ij; float ijmin, ijmax;
3595                GLfloat *uidb=NULL, *fff=NULL;
3596                /* grab depth buffer DOES NOT WORK*/
3597                dbuf = SUMA_grabPixels(5, sv->X->aWIDTH, sv->X->aHEIGHT);
3598                fff = (GLfloat *)dbuf;
3599                uidb = (GLfloat *)
3600                      SUMA_malloc(sv->X->aWIDTH*sv->X->aHEIGHT*sizeof(GLfloat));
3601                glReadPixels(0, 0, sv->X->aWIDTH, sv->X->aHEIGHT,
3602                             GL_DEPTH_COMPONENT, GL_FLOAT, uidb);
3603 
3604                ijmin = ijmax = (float)uidb[0];
3605                for (ij=0; ij<sv->X->aWIDTH*sv->X->aHEIGHT; ++ij) {
3606                   if (ijmin > uidb[ij]) {
3607                      ijmin = (float)uidb[ij];
3608                      fprintf(stderr,"min %f\n", (float)uidb[ij]);
3609                   }
3610                   if (ijmax < uidb[ij]) {
3611                      ijmax = (float)uidb[ij];
3612                      fprintf(stderr,"max %f\n", (float)uidb[ij]);
3613                   }
3614                }
3615                fprintf(stderr, "range %f %f\n", ijmin, ijmax);
3616                fprintf(stderr, "%f versus %f\n", fff[0], uidb[0]);
3617                /* grab luminance buffer */
3618                dlum = SUMA_grabPixels(1, sv->X->aWIDTH, sv->X->aHEIGHT);
3619                /* reset anything dark */
3620                /* rewrite depth buffer */
3621                SUMA_PixelsToDisk(sv, sv->X->aWIDTH, sv->X->aHEIGHT,
3622                                  uidb, 5, 1, "dbuf.jpg", 1, 1);
3623                SUMA_PixelsToDisk(sv, sv->X->aWIDTH, sv->X->aHEIGHT,
3624                                  dlum, 1, 1, "dlum.jpg", 1, 1);
3625                /* free buffers */
3626                SUMA_ifree(dbuf); SUMA_ifree(dlum);
3627             }
3628             #endif
3629             {
3630                GLfloat *dbuf;
3631                GLubyte *lbuf;
3632                int ij;
3633                float ijmin, ijmax;
3634                /* grab depth buffer DOES NOT WORK*/
3635                dbuf = (GLfloat *)
3636                      SUMA_malloc(sv->X->aWIDTH*sv->X->aHEIGHT*sizeof(GLfloat));
3637                glReadPixels(0, 0, sv->X->aWIDTH, sv->X->aHEIGHT,
3638                             GL_DEPTH_COMPONENT, GL_FLOAT, dbuf);
3639                ijmin = ijmax = (float)dbuf[0];
3640                for (ij=0; ij<sv->X->aWIDTH*sv->X->aHEIGHT; ++ij) {
3641                   if (ijmin > dbuf[ij]) {
3642                      ijmin = (float)dbuf[ij];
3643                      fprintf(stderr,"min %f\n", (float)dbuf[ij]);
3644                   }
3645                   if (ijmax < dbuf[ij]) {
3646                      ijmax = (float)dbuf[ij];
3647                      fprintf(stderr,"max %f\n", (float)dbuf[ij]);
3648                   }
3649                }
3650                fprintf(stderr, "depth range %f %f\n", ijmin, ijmax);
3651                /* grab luminance buffer */
3652                lbuf = (GLubyte *)SUMA_grabPixels(1,sv->X->aWIDTH,sv->X->aHEIGHT);
3653                ijmin = ijmax = (float)lbuf[0];
3654                for (ij=0; ij<sv->X->aWIDTH*sv->X->aHEIGHT; ++ij) {
3655                   if (ijmin > lbuf[ij]) {
3656                      ijmin = (float)lbuf[ij];
3657                      fprintf(stderr,"min %f\n", (float)lbuf[ij]);
3658                   }
3659                   if (ijmax < lbuf[ij]) {
3660                      ijmax = (float)lbuf[ij];
3661                      fprintf(stderr,"max %f\n", (float)lbuf[ij]);
3662                   }
3663                }
3664                fprintf(stderr, "luminance range %f %f\n", ijmin, ijmax);
3665                /* reset anything dark */
3666                for (ij=0; ij<sv->X->aWIDTH*sv->X->aHEIGHT; ++ij) {
3667                   if (lbuf[ij]<20) dbuf[ij]=-1.0; /* max it out */
3668                }
3669                /* rewrite depth buffer */
3670                glWindowPos2s(0,0); /* specify raster position */
3671                glDrawPixels(sv->X->aWIDTH, sv->X->aHEIGHT,
3672                             GL_DEPTH_COMPONENT, GL_FLOAT, dbuf);
3673                /* render last slice again */
3674                if ((nqd = SUMA_PlaneBoxIntersect( sv->GVS[sv->StdView].ViewFrom,
3675                                              rslc->Eq, VO->VE[ive]->bcorners,
3676                                              slc_corn)) > 2) {
3677                SUMA_LH("Have plane %f %f %f %f, %d pts on vol %s",
3678                            rslc->Eq[0], rslc->Eq[1], rslc->Eq[2], rslc->Eq[3],
3679                            nqd,  SUMA_VE_Headname(VO->VE,ive));
3680                glBegin(GL_POLYGON);
3681                   for (k=0; k<6; ++k) { /* draw all 6 points always, even
3682                                            when there are repetitions.
3683                                            Don't bother trimming to unique
3684                                            set unless this causes trouble */
3685                      /* change mm (edge coordinate to texture coords) */
3686                      AFF44_MULT_I(tex_corn, VO->VE[ive]->X2I, (slc_corn+3*k));
3687                      /* offset indices because slc_corn is on edge */
3688                      if (tex_corn[0] < 0) tex_corn[0] = 0;
3689                      else if (tex_corn[0] > VO->VE[ive]->Ni-1)
3690                         tex_corn[0] = VO->VE[ive]->Ni-1;
3691                      if (tex_corn[1] < 0) tex_corn[1] = 0;
3692                      else if (tex_corn[1] > VO->VE[ive]->Nj-1)
3693                         tex_corn[1] = VO->VE[ive]->Nj-1;
3694                      if (tex_corn[2] < 0) tex_corn[2] = 0;
3695                      else if (tex_corn[2] > VO->VE[ive]->Nk-1)
3696                         tex_corn[2] = VO->VE[ive]->Nk-1;
3697                      tex_corn[0] /= (float)(VO->VE[ive]->Ni-1);
3698                      tex_corn[1] /= (float)(VO->VE[ive]->Nj-1);
3699                      tex_corn[2] /= (float)(VO->VE[ive]->Nk-1);
3700                      glTexCoord3f(tex_corn[0],
3701                                   tex_corn[1], tex_corn[2]);
3702                            /* this one is affected by the Texture MatrixMode */
3703                      glVertex3f(slc_corn[3*k],
3704                                 slc_corn[3*k+1], slc_corn[3*k+2]);
3705                            /* this one is affected by the Modelview matrixMode*/
3706                   }
3707                glEnd();
3708                }
3709                /* junk it */
3710                SUMA_PixelsToDisk(sv, sv->X->aWIDTH, sv->X->aHEIGHT,
3711                                  dbuf, 5, 1, "dbuf.jpg", 1, 1);
3712                SUMA_PixelsToDisk(sv, sv->X->aWIDTH, sv->X->aHEIGHT,
3713                                  lbuf, 1, 1, "dlum.jpg", 1, 1);
3714                /* free buffers */
3715                SUMA_ifree(dbuf); SUMA_ifree(lbuf);
3716             }
3717             if (ive > 0) {
3718                SUMA_S_Warn("Add blending here");
3719                /* Here is where you blend slice from this VE with the previous
3720                   This should work just fine as is actually, no need to blend
3721                   separately unless doing overlay on top always. In that case,
3722                   render each VE separately then blend results across VEs
3723                */
3724             }
3725             ++ive;
3726          }
3727 
3728          if (el != dlist_head(VSaux->slcl)) el = dlist_prev(el);
3729          else el = NULL;
3730       } while (el);
3731    }
3732    SUMA_CHECK_GL_ERROR("OpenGL Error ddd");
3733 
3734    glFlush();
3735 
3736 
3737    /* Here we create a texture on the cutplane from the last dset loaded
3738    At the moment, without this texture, nothing shows
3739    of the overlay volume */
3740    if (0 && LastTextureOnCutPlane) {
3741       --ive; /* bring ive counter to last dset put into texture*/
3742       SUMA_dset_tex_slice_corners( 0, SUMA_VE_dset(VO->VE, ive),
3743                                    tex_corn, slc_corn, NULL, 2, 0);
3744 
3745       // Joachim says this is just wrong ...
3746       tz = 0.5+(-VO->CutPlane[0][3])/(float)VO_NK(VO);
3747 
3748       glEnable(GL_DEPTH_TEST);
3749 
3750       /* If it were not for the slice textures shown here, then the overlay
3751          texture would not show up at all! */
3752       #if 0
3753       SUMA_LH("Texture on the slice, with triangles");
3754       glBegin(GL_TRIANGLES);
3755          k = 0;
3756          glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
3757             glVertex3f( slc_corn[3*k], slc_corn[3*k+1],
3758                                     -VO->CutPlane[0][3]); ++k;
3759          glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
3760             glVertex3f( slc_corn[3*k], slc_corn[3*k+1],
3761                                     -VO->CutPlane[0][3]); ++k;
3762          glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
3763             glVertex3f( slc_corn[3*k], slc_corn[3*k+1],
3764                                     -VO->CutPlane[0][3]);
3765 
3766          k = 0;
3767          glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
3768             glVertex3f( slc_corn[3*k], slc_corn[3*k+1],
3769                                     -VO->CutPlane[0][3]); k+=2;
3770          glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
3771             glVertex3f( slc_corn[3*k], slc_corn[3*k+1],
3772                                     -VO->CutPlane[0][3]); ++k;
3773          glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
3774             glVertex3f( slc_corn[3*k], slc_corn[3*k+1],
3775                                     -VO->CutPlane[0][3]);
3776 
3777       glEnd();
3778       #else
3779       SUMA_LH("Texture on the slice, QUADS?");
3780       glBegin(GL_QUADS);
3781          for (k=0; k<4; ++k) {
3782             glTexCoord3f(tex_corn[3*k], tex_corn[3*k+1], tz);
3783                   /* this one is affected by the Texture MatrixMode */
3784             glVertex3f(slc_corn[3*k], slc_corn[3*k+1],
3785                                           -VO->CutPlane[0][3]);
3786                   /* this one is affected by the Modelview matrixMode*/
3787          }
3788       glEnd();
3789       #endif
3790       glDisable(GL_DEPTH_TEST);
3791    }
3792    glDisable(GL_TEXTURE_3D);
3793 
3794    if (!gl_bl) glDisable(GL_BLEND);
3795    if (shmodel != GL_FLAT) glShadeModel(shmodel);
3796    #if 0
3797    glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
3798    glEnable(GL_COLOR_MATERIAL);
3799    for (iplane=0; iplane < 6; ++iplane) {
3800       if (VO->UseCutPlane[iplane]) {
3801          if (iplane == VO->SelectedCutPlane) glColor3f(1.0, 1.0, 1.0);
3802          else {
3803             if (ShowUnselected) {
3804                if (iplane==0 || iplane == 1) glColor3f(1.0, 0.0, 0.0);
3805                if (iplane==2 || iplane == 3) glColor3f(0.0, 1.0, 0.0);
3806                if (iplane==4 || iplane == 5) glColor3f(0.0, 0.0, 1.0);
3807             } else {
3808                continue;
3809             }
3810          }
3811          glBegin(GL_LINE_LOOP);
3812             nlt = VO->SOcut[iplane]->NodeList;
3813             glVertex3f( nlt[0],nlt[1],nlt[2] );
3814             glVertex3f( nlt[3],nlt[4],nlt[5] );
3815             glVertex3f( nlt[6],nlt[7],nlt[8] );
3816             glVertex3f( nlt[9],nlt[10],nlt[11] );
3817          glEnd();
3818       }
3819    }
3820    glDisable(GL_COLOR_MATERIAL);
3821    #endif
3822 
3823    if (sv->PolyMode != SRM_Fill) {/* set fill mode back */
3824       SUMA_SET_GL_RENDER_MODE(sv->PolyMode);
3825    }
3826 
3827    if (gl_dt) glEnable(GL_DEPTH_TEST);
3828    else glDisable(GL_DEPTH_TEST);
3829    if (gl_bl) glEnable(GL_BLEND);
3830    else glDisable(GL_BLEND);
3831    SUMA_RETURN(YUP);
3832    #endif
3833 }
3834 
3835 
SUMA_VE_LoadTexture(SUMA_VolumeElement ** VE,int n)3836 SUMA_Boolean SUMA_VE_LoadTexture(SUMA_VolumeElement **VE, int n)
3837 {
3838    static char FuncName[]={"SUMA_VE_LoadTexture"};
3839    int i;
3840    SUMA_Boolean LocalHead = NOPE;
3841 
3842    SUMA_ENTRY;
3843 
3844    if (!VE || n < 0 || !VE[n]) {
3845       SUMA_S_Err("NULL input %p %d %p", VE, n, (VE && n >=0) ? VE[n]:NULL);
3846       SUMA_RETURN(NOPE);
3847    }
3848    glPixelStorei(GL_UNPACK_ALIGNMENT, 1); /* Have no padding at the
3849                                                 end of texel rows*/
3850    if (!VE[n]->texName) {
3851       VE[n]->texName = (GLuint *)SUMA_calloc(1,sizeof(GLuint));
3852       glGenTextures(1, VE[n]->texName); /* I just need 1 for now*/
3853    }
3854    if (!VE[n]->texvec) {
3855       SUMA_S_Err("NULL texvec!");
3856       SUMA_RETURN(NOPE);
3857    }
3858 
3859    glBindTexture(GL_TEXTURE_3D, VE[n]->texName[0]);
3860                         /* make texName be the current one
3861                            This will also create texture object
3862                            since this is the first time it is
3863                            called*/
3864 
3865    /* Should texture repeat or no ? - This does not need to be set with
3866       every draw call, it should be done at init time*/
3867    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP);
3868    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP);
3869    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP);
3870 
3871    /* How should magnification and reduction be handled? */
3872    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
3873    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
3874 
3875    SUMA_LHv("Storing texture for %s: %d %d %d\n",
3876             SUMA_VE_Headname(VE, n), SUMA_VE_Ni(VE, n),
3877             SUMA_VE_Nj(VE, n),  SUMA_VE_Nk(VE, n));
3878    /* And store the image poiner in question */
3879    glTexImage3D   (  GL_TEXTURE_3D,
3880                      0, /* texture level, highest resolution */
3881                      GL_RGBA, /* RGBA baby*/
3882                      SUMA_VE_Ni(VE,n),
3883                      SUMA_VE_Nj(VE,n),
3884                      SUMA_VE_Nk(VE,n),
3885 		               0, /* border is 0 wide. Might have to do borders and
3886                            split texture into two, if
3887                            volume is too big to fit into kangaroo's pouch */
3888                      GL_RGBA, GL_UNSIGNED_BYTE,
3889                      VE[n]->texvec);
3890    SUMA_RETURN(YUP);
3891 }
3892 
SUMA_VO_InitCutPlanes(SUMA_VolumeObject * VO)3893 SUMA_Boolean SUMA_VO_InitCutPlanes(SUMA_VolumeObject *VO)
3894 {
3895    static char FuncName[]={"SUMA_VO_InitCutPlanes"};
3896    int i,n;
3897    SUMA_Boolean LocalHead = NOPE;
3898 
3899    SUMA_ENTRY;
3900 
3901    if (!VO) SUMA_RETURN(NOPE);
3902 
3903       /* initialize clipping planes, assuming RAI.
3904    Even numbered planes cut above the coordinate
3905    Odd numbered planes cut below the coordinate.
3906    For a plane P, Objects at X,Y,Z where
3907       p[0]X+p[1]Y+p[2]Z+p[3]>=0 are displayed
3908    Those planes are applied in object space.
3909    */
3910 
3911 
3912    SUMA_LH("The cut off points should be interactive\n"
3913                "so this next assignment should be done \n"
3914                "with each drawing operation.\n"
3915                "These default values, based on First and Last Voxel\n"
3916                "should be saved in nel.\n"
3917                "This stuff is not used at the moment. Reconsider later");
3918 
3919    VO->CutPlane[0][3] = VO->VE[0]->voN[0]; /* Xmore */
3920    VO->CutPlane[1][3] = -VO->VE[0]->vo0[0]; /* Xless */
3921    VO->CutPlane[2][3] = VO->VE[0]->voN[1]; /* Ymore */
3922    VO->CutPlane[3][3] = -VO->VE[0]->vo0[1]; /* Yless */
3923    VO->CutPlane[4][3] = VO->VE[0]->voN[2]; /* Zmore */
3924    VO->CutPlane[5][3] = -VO->VE[0]->vo0[2]; /* Zless */
3925 
3926    VO->CutPlane[0][0] = 0.0;  /* Z (I/S) plane Clip below plane*/
3927    VO->CutPlane[0][1] = 0.0;
3928    VO->CutPlane[0][2] = 1.0;
3929    VO->CutPlane[0][3] = 0.8*VO->VE[0]->voN[2];
3930 
3931    VO->CutPlane[1][0] = 0.0;  /* Z (I/S) plane, clip above plane*/
3932    VO->CutPlane[1][1] = 0.0;
3933    VO->CutPlane[1][2] = -1.0;
3934    VO->CutPlane[1][3] = -0.8*VO->VE[0]->vo0[2];
3935 
3936    VO->CutPlane[2][0] = 0.0;  /* Y (A/P) plane, clip anterior plane*/
3937    VO->CutPlane[2][1] = 1.0;
3938    VO->CutPlane[2][2] = 0.0;
3939    VO->CutPlane[2][3] = 0.8*VO->VE[0]->voN[1];;
3940 
3941    VO->CutPlane[3][0] = 0.0;  /* Y (A/P) plane, clip posterior plane*/
3942    VO->CutPlane[3][1] = -1.0;
3943    VO->CutPlane[3][2] = 0.0;
3944    VO->CutPlane[3][3] = -0.8*VO->VE[0]->vo0[1];
3945 
3946    VO->CutPlane[4][0] = 1.0;  /* X (R/L) plane, clip right plane*/
3947    VO->CutPlane[4][1] = 0.0;
3948    VO->CutPlane[4][2] = 0.0;
3949    VO->CutPlane[4][3] = 0.8*VO->VE[0]->voN[0];
3950 
3951    VO->CutPlane[5][0] = -1.0;  /* X (R/L) plane, clip left plane*/
3952    VO->CutPlane[5][1] = 0.0;
3953    VO->CutPlane[5][2] = 0.0;
3954    VO->CutPlane[5][3] = -0.8*VO->VE[0]->vo0[0];
3955 
3956    for (i=0; i<6; ++i) {
3957       SUMA_SetTextureClipPlaneSurface(VO, i);
3958    }
3959 
3960    SUMA_RETURN(YUP);
3961 }
3962 
SUMA_CreateGL3DTexture(SUMA_VolumeObject * VO)3963 SUMA_Boolean SUMA_CreateGL3DTexture(SUMA_VolumeObject *VO)
3964 {
3965    static char FuncName[]={"SUMA_CreateGL3DTexture"};
3966    int i,n;
3967    SUMA_Boolean LocalHead = NOPE;
3968 
3969    SUMA_ENTRY;
3970 
3971    SUMA_LH("Initializing and creating texture\n"
3972            "So if I am handling the mixing, it is the \n"
3973            "color vector that ends up filling the texture map...");
3974    glPixelStorei(GL_UNPACK_ALIGNMENT, 1); /* Have no padding at the
3975                                                 end of texel rows*/
3976    if (0) {
3977    SUMA_S_Note(
3978          "Need to allow user to specify more than one texture.\n"
3979          "This would happen if each new nel got its own texName\n"
3980          "But this operation might fail if lots of textures are \n"
3981          "loaded, so I should guard against that.\n"
3982          "Regarding multiple textures, the simplest is to have\n"
3983          "anatomy and function in each texture and then render one\n"
3984          "after the other.\n"
3985          "Also, should check how and when textures are to be deleted\n"
3986          "and whether texture nel (and other NIDOs) are being properly\n"
3987          "disposed of when a DO is freed\n"
3988          "Lastly, is the issue of the dset header that I loath to \n"
3989          "duplicate into nel. I should just keep dset's header around\n"
3990          "somehow. For the moment, dset is one static variable in this\n"
3991          "file and that is clearly not appropriate. So have to deal with\n"
3992          "that problem too.\n");
3993    }
3994    n = 0;
3995    while (VO->VE && VO->VE[n]) {
3996       if (!VO->VE[n]->texName) {
3997          VO->VE[n]->texName = (GLuint *)SUMA_calloc(1,sizeof(GLuint));
3998          glGenTextures(1, VO->VE[n]->texName); /* I just need 1 for now*/
3999 
4000          if (!SUMA_VE_LoadTexture(VO->VE, n)) {
4001             SUMA_S_Err("Failed to load texture for %d", n);
4002             SUMA_RETURN(NOPE);
4003          }
4004       } else {
4005          SUMA_S_Note("Proably done already via SUMA_Overlays_2_GLCOLAR4's\n"
4006                      "call to SUMA_VE_LoadTexture. Does this function still \n"
4007                      "have a reason to exist?");
4008       }
4009       ++n;
4010    }
4011 
4012    if (!SUMA_VO_InitCutPlanes(VO)) {
4013       SUMA_S_Err("Failed to init cutplanes");
4014       SUMA_RETURN(NOPE);
4015    }
4016 
4017    SUMA_RETURN(YUP);
4018 }
4019 
SUMA_Draw3DTextureNIDOnel(NI_element * nel,SUMA_SurfaceObject * SO,SUMA_DO_CoordUnits default_coord_type,float * default_txcol,void * default_font,int default_node,SUMA_SurfaceViewer * sv)4020 SUMA_Boolean SUMA_Draw3DTextureNIDOnel (NI_element *nel,
4021                                     SUMA_SurfaceObject *SO,
4022                                     SUMA_DO_CoordUnits default_coord_type,
4023                                     float *default_txcol,
4024                                     void *default_font, int default_node,
4025                                     SUMA_SurfaceViewer *sv)
4026 {
4027    static char FuncName[]={"SUMA_Draw3DTextureNIDOnel"};
4028    SUMA_VolumeObject *VO=NULL;
4029    SUMA_Boolean LocalHead = NOPE;
4030 
4031    SUMA_ENTRY;
4032 
4033 
4034    if (!nel || strcmp(nel->name,"3DTex")) SUMA_RETURN(NOPE);
4035 
4036    if (NI_IS_STR_ATTR_EQUAL(nel,"read_status","fail")) {
4037       /* can't be read */
4038       SUMA_LH("read_status is fail, giving up");
4039       SUMA_RETURN(NOPE);
4040    }
4041 
4042 
4043    if (!NI_IS_STR_ATTR_EQUAL(nel,"read_status","read")) { /* read it */
4044       if (!SUMA_Load3DTextureNIDOnel(nel, default_coord_type)) {
4045          SUMA_LH("Failed to load 3d texture");
4046          SUMA_RETURN(NOPE);
4047       }
4048 
4049       if (!(VO = SUMA_VOof3DTextureNIDOnel(nel))) {
4050          SUMA_S_Err("Failed to find corresponding VO");
4051          SUMA_RETURN(NOPE);
4052       }
4053 
4054       if (!SUMA_CreateGL3DTexture(VO)) {
4055          SUMA_S_Err("Failed to create texture");
4056          SUMA_RETURN(NOPE);
4057       }
4058    } else {
4059       if (!(VO = SUMA_VOof3DTextureNIDOnel(nel))) {
4060          SUMA_S_Err("Failed to find corresponding VO");
4061          SUMA_RETURN(NOPE);
4062       }
4063    }
4064 
4065    if (!SUMA_DrawVolumeDO(VO, sv)) {
4066       SUMA_S_Err("Failed to draw volume");
4067       SUMA_RETURN(NOPE);
4068    }
4069 
4070 
4071 
4072 
4073    SUMA_RETURN(YUP);
4074 }
4075 
iPlane2Dim(int iplane)4076 int iPlane2Dim(int iplane)
4077 {
4078    if (iplane == 0 || iplane == 1) return(2);
4079    else if (iplane == 2 || iplane == 3) return(1);
4080    else if (iplane == 4 || iplane == 5) return(0);
4081    else return(-1);
4082 }
4083 
SUMA_SetTextureClipPlaneSurface(SUMA_VolumeObject * VO,int iplane)4084 SUMA_Boolean SUMA_SetTextureClipPlaneSurface(
4085                         SUMA_VolumeObject *VO, int iplane )
4086 {
4087    static char FuncName[]={"SUMA_SetTextureClipPlaneSurface"};
4088    SUMA_SurfaceObject *SO=NULL;
4089    static int nwarn = 0;
4090    int k;
4091    int *FaceSetList=NULL;
4092    float *NodeList=NULL;
4093    GLfloat tex_corn[12] ;
4094    GLfloat slc_corn[12] ;
4095    SUMA_DSET *dset=NULL;
4096    SUMA_NEW_SO_OPT *nsoopt=NULL;
4097    SUMA_Boolean LocalHead = NOPE;
4098 
4099    SUMA_ENTRY;
4100 
4101    if (!VO || !VO->VE || !(dset = SUMA_VE_dset(VO->VE,0))
4102        || iplane < 0 || iplane > 5) {
4103       SUMA_S_Err("Bad flow");
4104       SUMA_RETURN(NOPE);
4105    }
4106    if (!VO->SOcut[iplane]) {
4107       SUMA_LH("Creating surface for first time");
4108 
4109       NodeList = (float *)SUMA_calloc(3*4, sizeof(float));
4110             /* four nodes for a rectangle */
4111       FaceSetList = (int *)SUMA_calloc(2*3, sizeof(int));
4112             /* two triangles   */
4113 
4114       switch (iplane) {
4115          case 0:/* initialize for plane 0 */
4116             SUMA_dset_tex_slice_corners( 0, dset, tex_corn,
4117                                          slc_corn, NULL, iPlane2Dim(iplane), 0);
4118             k=0;
4119             NodeList[3*k] = slc_corn[3*k];
4120             NodeList[3*k+1] = slc_corn[3*k+1];
4121             NodeList[3*k+2] = -VO->CutPlane[iplane][3];
4122             k++;
4123             NodeList[3*k] = slc_corn[3*k];
4124             NodeList[3*k+1] = slc_corn[3*k+1];
4125             NodeList[3*k+2] = -VO->CutPlane[iplane][3];
4126             k++;
4127             NodeList[3*k] = slc_corn[3*k];
4128             NodeList[3*k+1] = slc_corn[3*k+1];
4129             NodeList[3*k+2] = -VO->CutPlane[iplane][3];
4130             k++;
4131             NodeList[3*k] = slc_corn[3*k];
4132             NodeList[3*k+1] = slc_corn[3*k+1];
4133             NodeList[3*k+2] = -VO->CutPlane[iplane][3];
4134             break;
4135          case 1:
4136             SUMA_dset_tex_slice_corners( 0, dset,  tex_corn,
4137                                          slc_corn, NULL, iPlane2Dim(iplane), 0);
4138             k=0;
4139             NodeList[3*k] = slc_corn[3*k];
4140             NodeList[3*k+1] = slc_corn[3*k+1];
4141             NodeList[3*k+2] = -VO->CutPlane[iplane][3];
4142             k++;
4143             NodeList[3*k] = slc_corn[3*k];
4144             NodeList[3*k+1] = slc_corn[3*k+1];
4145             NodeList[3*k+2] = -VO->CutPlane[iplane][3];
4146             k++;
4147             NodeList[3*k] = slc_corn[3*k];
4148             NodeList[3*k+1] = slc_corn[3*k+1];
4149             NodeList[3*k+2] = -VO->CutPlane[iplane][3];
4150             k++;
4151             NodeList[3*k] = slc_corn[3*k];
4152             NodeList[3*k+1] = slc_corn[3*k+1];
4153             NodeList[3*k+2] = -VO->CutPlane[iplane][3];
4154             break;
4155          case 2:
4156             SUMA_dset_tex_slice_corners( 0, dset,  tex_corn,
4157                                          slc_corn, NULL, iPlane2Dim(iplane), 0);
4158             k=0;
4159             NodeList[3*k] = slc_corn[3*k];
4160             NodeList[3*k+1] = -VO->CutPlane[iplane][3];
4161             NodeList[3*k+2] = slc_corn[3*k+2];
4162             k++;
4163             NodeList[3*k] = slc_corn[3*k];
4164             NodeList[3*k+1] = -VO->CutPlane[iplane][3];
4165             NodeList[3*k+2] = slc_corn[3*k+2];
4166             k++;
4167             NodeList[3*k] = slc_corn[3*k];
4168             NodeList[3*k+1] = -VO->CutPlane[iplane][3];
4169             NodeList[3*k+2] = slc_corn[3*k+2];
4170             k++;
4171             NodeList[3*k] = slc_corn[3*k];
4172             NodeList[3*k+1] = -VO->CutPlane[iplane][3];
4173             NodeList[3*k+2] = slc_corn[3*k+2];
4174             break;
4175          case 3:
4176             SUMA_dset_tex_slice_corners( 0, dset, tex_corn,
4177                                          slc_corn, NULL, iPlane2Dim(iplane), 0);
4178             k=0;
4179             NodeList[3*k] = slc_corn[3*k];
4180             NodeList[3*k+1] = -VO->CutPlane[iplane][3];
4181             NodeList[3*k+2] = slc_corn[3*k+2];
4182             k++;
4183             NodeList[3*k] = slc_corn[3*k];
4184             NodeList[3*k+1] = -VO->CutPlane[iplane][3];
4185             NodeList[3*k+2] = slc_corn[3*k+2];
4186             k++;
4187             NodeList[3*k] = slc_corn[3*k];
4188             NodeList[3*k+1] = -VO->CutPlane[iplane][3];
4189             NodeList[3*k+2] = slc_corn[3*k+2];
4190             k++;
4191             NodeList[3*k] = slc_corn[3*k];
4192             NodeList[3*k+1] = -VO->CutPlane[iplane][3];
4193             NodeList[3*k+2] = slc_corn[3*k+2];
4194             break;
4195          case 4:
4196             SUMA_dset_tex_slice_corners( 0, dset, tex_corn,
4197                                          slc_corn, NULL, iPlane2Dim(iplane), 0);
4198             k=0;
4199             NodeList[3*k] = -VO->CutPlane[iplane][3];
4200             NodeList[3*k+1] = slc_corn[3*k+1];
4201             NodeList[3*k+2] = slc_corn[3*k+2];
4202             k++;
4203             NodeList[3*k] = -VO->CutPlane[iplane][3];
4204             NodeList[3*k+1] = slc_corn[3*k+1];
4205             NodeList[3*k+2] = slc_corn[3*k+2];
4206             k++;
4207             NodeList[3*k] = -VO->CutPlane[iplane][3];
4208             NodeList[3*k+1] = slc_corn[3*k+1];
4209             NodeList[3*k+2] = slc_corn[3*k+2];
4210             k++;
4211             NodeList[3*k] = -VO->CutPlane[iplane][3];
4212             NodeList[3*k+1] = slc_corn[3*k+1];
4213             NodeList[3*k+2] = slc_corn[3*k+2];
4214             break;
4215          case 5:
4216             SUMA_dset_tex_slice_corners( 0, dset, tex_corn,
4217                                          slc_corn, NULL, iPlane2Dim(iplane), 0);
4218             k=0;
4219             NodeList[3*k] = -VO->CutPlane[iplane][3];
4220             NodeList[3*k+1] = slc_corn[3*k+1];
4221             NodeList[3*k+2] = slc_corn[3*k+2];
4222             k++;
4223             NodeList[3*k] = -VO->CutPlane[iplane][3];
4224             NodeList[3*k+1] = slc_corn[3*k+1];
4225             NodeList[3*k+2] = slc_corn[3*k+2];
4226             k++;
4227             NodeList[3*k] = -VO->CutPlane[iplane][3];
4228             NodeList[3*k+1] = slc_corn[3*k+1];
4229             NodeList[3*k+2] = slc_corn[3*k+2];
4230             k++;
4231             NodeList[3*k] = -VO->CutPlane[iplane][3];
4232             NodeList[3*k+1] = slc_corn[3*k+1];
4233             NodeList[3*k+2] = slc_corn[3*k+2];
4234             break;
4235          default:
4236             SUMA_S_Errv("Bad iplane %d\n", iplane);
4237             SUMA_RETURN(NOPE);
4238       }
4239       FaceSetList[0] = 0;
4240       FaceSetList[1] = 1;
4241       FaceSetList[2] = 2;
4242 
4243       FaceSetList[3] = 0;
4244       FaceSetList[4] = 2;
4245       FaceSetList[5] = 3;
4246 
4247       nsoopt = SUMA_NewNewSOOpt();
4248       VO->SOcut[iplane] = SUMA_NewSO( &NodeList, 4, &FaceSetList, 2, nsoopt );
4249       /* free nsoopt here */
4250       nsoopt = SUMA_FreeNewSOOpt(nsoopt);
4251    }
4252 
4253    NodeList = VO->SOcut[iplane]->NodeList;
4254    FaceSetList = VO->SOcut[iplane]->FaceSetList;
4255 
4256 
4257    /* adjust depending on which plane is called for */
4258    switch (iplane) {
4259       case 0:
4260          SUMA_dset_tex_slice_corners( 0, dset, tex_corn,
4261                                       slc_corn, NULL, iPlane2Dim(iplane), 0);
4262          k=0;
4263          NodeList[3*k] = slc_corn[3*k];
4264          NodeList[3*k+1] = slc_corn[3*k+1];
4265          NodeList[3*k+2] = -VO->CutPlane[iplane][3];
4266          k++;
4267          NodeList[3*k] = slc_corn[3*k];
4268          NodeList[3*k+1] = slc_corn[3*k+1];
4269          NodeList[3*k+2] = -VO->CutPlane[iplane][3];
4270          k++;
4271          NodeList[3*k] = slc_corn[3*k];
4272          NodeList[3*k+1] = slc_corn[3*k+1];
4273          NodeList[3*k+2] = -VO->CutPlane[iplane][3];
4274          k++;
4275          NodeList[3*k] = slc_corn[3*k];
4276          NodeList[3*k+1] = slc_corn[3*k+1];
4277          NodeList[3*k+2] = -VO->CutPlane[iplane][3];
4278 
4279          break;
4280       case 1:
4281          SUMA_dset_tex_slice_corners( 0, dset, tex_corn,
4282                                       slc_corn, NULL, iPlane2Dim(iplane), 0);
4283          k=0;
4284          NodeList[3*k] = slc_corn[3*k];
4285          NodeList[3*k+1] = slc_corn[3*k+1];
4286          NodeList[3*k+2] = VO->CutPlane[iplane][3];
4287          k++;
4288          NodeList[3*k] = slc_corn[3*k];
4289          NodeList[3*k+1] = slc_corn[3*k+1];
4290          NodeList[3*k+2] = VO->CutPlane[iplane][3];
4291          k++;
4292          NodeList[3*k] = slc_corn[3*k];
4293          NodeList[3*k+1] = slc_corn[3*k+1];
4294          NodeList[3*k+2] = VO->CutPlane[iplane][3];
4295          k++;
4296          NodeList[3*k] = slc_corn[3*k];
4297          NodeList[3*k+1] = slc_corn[3*k+1];
4298          NodeList[3*k+2] = VO->CutPlane[iplane][3];
4299          break;
4300       case 2:
4301          SUMA_dset_tex_slice_corners( 0, dset, tex_corn,
4302                                       slc_corn, NULL, iPlane2Dim(iplane), 0);
4303          k=0;
4304          NodeList[3*k] = slc_corn[3*k];
4305          NodeList[3*k+1] = -VO->CutPlane[iplane][3];
4306          NodeList[3*k+2] = slc_corn[3*k+2];
4307          k++;
4308          NodeList[3*k] = slc_corn[3*k];
4309          NodeList[3*k+1] = -VO->CutPlane[iplane][3];
4310          NodeList[3*k+2] = slc_corn[3*k+2];
4311          k++;
4312          NodeList[3*k] = slc_corn[3*k];
4313          NodeList[3*k+1] = -VO->CutPlane[iplane][3];
4314          NodeList[3*k+2] = slc_corn[3*k+2];
4315          k++;
4316          NodeList[3*k] = slc_corn[3*k];
4317          NodeList[3*k+1] = -VO->CutPlane[iplane][3];
4318          NodeList[3*k+2] = slc_corn[3*k+2];
4319          break;
4320       case 3:
4321          SUMA_dset_tex_slice_corners( 0, dset, tex_corn,
4322                                       slc_corn, NULL, iPlane2Dim(iplane), 0);
4323          k=0;
4324          NodeList[3*k] = slc_corn[3*k];
4325          NodeList[3*k+1] = VO->CutPlane[iplane][3];
4326          NodeList[3*k+2] = slc_corn[3*k+2];
4327          k++;
4328          NodeList[3*k] = slc_corn[3*k];
4329          NodeList[3*k+1] = VO->CutPlane[iplane][3];
4330          NodeList[3*k+2] = slc_corn[3*k+2];
4331          k++;
4332          NodeList[3*k] = slc_corn[3*k];
4333          NodeList[3*k+1] = VO->CutPlane[iplane][3];
4334          NodeList[3*k+2] = slc_corn[3*k+2];
4335          k++;
4336          NodeList[3*k] = slc_corn[3*k];
4337          NodeList[3*k+1] = VO->CutPlane[iplane][3];
4338          NodeList[3*k+2] = slc_corn[3*k+2];
4339          break;
4340       case 4:
4341          SUMA_dset_tex_slice_corners( 0, dset, tex_corn,
4342                                       slc_corn, NULL, iPlane2Dim(iplane), 0);
4343          k=0;
4344          NodeList[3*k] = -VO->CutPlane[iplane][3];
4345          NodeList[3*k+1] = slc_corn[3*k+1];
4346          NodeList[3*k+2] = slc_corn[3*k+2];
4347          k++;
4348          NodeList[3*k] = -VO->CutPlane[iplane][3];
4349          NodeList[3*k+1] = slc_corn[3*k+1];
4350          NodeList[3*k+2] = slc_corn[3*k+2];
4351          k++;
4352          NodeList[3*k] = -VO->CutPlane[iplane][3];
4353          NodeList[3*k+1] = slc_corn[3*k+1];
4354          NodeList[3*k+2] = slc_corn[3*k+2];
4355          k++;
4356          NodeList[3*k] = -VO->CutPlane[iplane][3];
4357          NodeList[3*k+1] = slc_corn[3*k+1];
4358          NodeList[3*k+2] = slc_corn[3*k+2];
4359          break;
4360       case 5:
4361          SUMA_dset_tex_slice_corners( 0, dset, tex_corn,
4362                                       slc_corn, NULL, iPlane2Dim(iplane), 0);
4363          k=0;
4364          NodeList[3*k] =  VO->CutPlane[iplane][3];
4365          NodeList[3*k+1] = slc_corn[3*k+1];
4366          NodeList[3*k+2] = slc_corn[3*k+2];
4367          k++;
4368          NodeList[3*k] =  VO->CutPlane[iplane][3];
4369          NodeList[3*k+1] = slc_corn[3*k+1];
4370          NodeList[3*k+2] = slc_corn[3*k+2];
4371          k++;
4372          NodeList[3*k] =  VO->CutPlane[iplane][3];
4373          NodeList[3*k+1] = slc_corn[3*k+1];
4374          NodeList[3*k+2] = slc_corn[3*k+2];
4375          k++;
4376          NodeList[3*k] =  VO->CutPlane[iplane][3];
4377          NodeList[3*k+1] = slc_corn[3*k+1];
4378          NodeList[3*k+2] = slc_corn[3*k+2];
4379          break;
4380       default:
4381          SUMA_S_Err("Should not be here");
4382          break;
4383    }
4384 
4385    SUMA_RETURN(YUP);
4386 }
4387 
SUMA_TextureClipPlaneSurfaces(int * N_SOv)4388 SUMA_SurfaceObject ** SUMA_TextureClipPlaneSurfaces(int *N_SOv)
4389 {
4390    static char FuncName[]={"SUMA_TextureClipPlaneSurfaces"};
4391    SUMA_SurfaceObject **SOv=NULL;
4392    SUMA_VolumeObject *VO=NULL;
4393    int ii=0, jj, kk;
4394    SUMA_Boolean LocalHead = NOPE;
4395 
4396 
4397    SUMA_ENTRY;
4398 
4399    *N_SOv=0;
4400    for (ii=0; ii<SUMAg_N_DOv; ++ii) {
4401       if (SUMA_isVO(SUMAg_DOv[ii])) *N_SOv = *N_SOv+6;
4402    }
4403 
4404    SOv = (SUMA_SurfaceObject **)SUMA_calloc(*N_SOv,
4405                                  sizeof(SUMA_SurfaceObject *));
4406    kk = 0;
4407    for (ii=0; ii<SUMAg_N_DOv; ++ii) {
4408       if (SUMA_isVO(SUMAg_DOv[ii])) {
4409          VO = (SUMA_VolumeObject *)(SUMAg_DOv[ii].OP);
4410          for (jj=0; jj<6; ++jj) {
4411             if (VO->UseCutPlane[jj]) {
4412                SOv[kk] = VO->SOcut[jj]; ++kk;
4413             }
4414          }
4415       }
4416    }
4417 
4418    *N_SOv = kk;
4419 
4420    SUMA_RETURN(SOv);
4421 }
4422 
SUMA_VolumeObjectOfClipPlaneSurface(SUMA_SurfaceObject * SO)4423 SUMA_VolumeObject *SUMA_VolumeObjectOfClipPlaneSurface(SUMA_SurfaceObject *SO)
4424 {
4425    static char FuncName[]={"SUMA_VolumeObjectOfClipPlaneSurface"};
4426    SUMA_VolumeObject *VO=NULL, *VOr = NULL;
4427    int ii=0, jj, kk;
4428    SUMA_Boolean LocalHead = NOPE;
4429 
4430    SUMA_ENTRY;
4431 
4432    VO = NULL; VOr=NULL;
4433    for (ii=0; ii<SUMAg_N_DOv; ++ii) {
4434       if (SUMA_isVO(SUMAg_DOv[ii])) {
4435          VO = (SUMA_VolumeObject *)(SUMAg_DOv[ii].OP);
4436          for (jj=0; jj<6; ++jj) {
4437             if (VO->SOcut[jj] == SO) {
4438                if (!VOr) VOr = VO;
4439                else {
4440                   SUMA_S_Err("Found more than one VO for SO");
4441                   SUMA_RETURN(NULL);
4442                }
4443             }
4444          }
4445       }
4446    }
4447 
4448    SUMA_RETURN(VOr);
4449 }
4450 
SUMA_VO_SlicesAtCrosshair(SUMA_VolumeObject * VO)4451 int SUMA_VO_SlicesAtCrosshair(SUMA_VolumeObject *VO)
4452 {
4453    SUMA_VOL_SAUX *VSaux = VDO_VSAUX(VO);
4454    if (VSaux) return(VSaux->SlicesAtCrosshair);
4455    return(0);
4456 }
4457 
4458