{   \
   SUMA_OFFSET_LL_DATUM * m_dat = NULL; \
   DListElmt *m_Elm = NULL;   \
   do {  \
     if (m_Elm) m_Elm = m_Elm->next;   \
     else m_Elm =  dlist_head(list); \
     m_dat = (SUMA_OFFSET_LL_DATUM *)m_Elm->data; \
   } while(m_dat->layer != LayInd && m_Elm != dlist_tail(list));   \
   if (m_dat->layer != LayInd) Elm = NULL;  \
   else Elm = m_Elm; \
}

{  \
   SUMA_OFFSET_LL_DATUM * m_dat = NULL; \
   DListElmt *m_Elm = NULL;   \
   do {  \
     if (m_Elm) m_Elm = m_Elm->next;   \
     else m_Elm =  dlist_head(list); \
     m_dat = (SUMA_OFFSET_LL_DATUM *)m_Elm->data; \
   } while(m_dat->ni != n_jne && m_Elm != dlist_tail(list));   \
   if (m_dat->ni != n_jne) Elm = NULL;  \
   else Elm = m_Elm; \
}

00819 {
00820    static char FuncName[]={"SUMA_AddNodeToLayer"};
00821    static SUMA_Boolean LocalHead = NOPE;
00822    
00823    /* is this a new layer */
00824    if (LayInd > OffS->N_layers) { /* error */
00825       SUMA_SL_Err("LayInd > OffS->N_layers. This should not be!");
00826       SUMA_RETURN(NOPE);
00827    } else if (LayInd == OffS->N_layers) { /* need a new one */
00828       SUMA_LH("Adding layer");
00829       OffS->N_layers += 1;
00830       OffS->layers = (SUMA_NODE_NEIGHB_LAYER *) SUMA_realloc(OffS->layers, OffS->N_layers*sizeof(SUMA_NODE_NEIGHB_LAYER));
00831       OffS->layers[LayInd].N_AllocNodesInLayer = 200;
00832       OffS->layers[LayInd].NodesInLayer = (int *) SUMA_malloc(OffS->layers[LayInd].N_AllocNodesInLayer * sizeof(int));
00833       OffS->layers[LayInd].N_NodesInLayer = 0;
00834    }
00835    
00836    OffS->layers[LayInd].N_NodesInLayer += 1;
00837    /* do we need to reallocate for NodesInLayer ? */
00838    if (OffS->layers[LayInd].N_NodesInLayer ==  OffS->layers[LayInd].N_AllocNodesInLayer) { /* need more space */
00839       SUMA_LH("reallocating neighbors");
00840       OffS->layers[LayInd].N_AllocNodesInLayer += 200;
00841       OffS->layers[LayInd].NodesInLayer = (int *) SUMA_realloc (OffS->layers[LayInd].NodesInLayer, OffS->layers[LayInd].N_AllocNodesInLayer * sizeof(int));
00842    }
00843    
00844    OffS->layers[LayInd].NodesInLayer[OffS->layers[LayInd].N_NodesInLayer - 1] = n;
00845    
00846    SUMA_RETURN(YUP); 
00847 }

06585 {
06586    static char FuncName[]={"SUMA_Allocate_SPI"};
06587    int i;
06588    SUMA_SURF_PLANE_INTERSECT *SPI = NULL;
06589    
06590    SUMA_ENTRY;
06591    
06592    SPI = (SUMA_SURF_PLANE_INTERSECT *) SUMA_malloc(sizeof(SUMA_SURF_PLANE_INTERSECT));
06593    if (!SPI) {
06594       fprintf (SUMA_STDERR, "Error %s: Could not allocate for SPI\n", FuncName);
06595       SUMA_RETURN (SPI);
06596    }
06597    
06598    SPI->IntersEdges = (int *) SUMA_calloc (SO->EL->N_EL, sizeof(int)); /* allocate for the max imaginable*/
06599    SPI->IntersNodes = (float *) SUMA_calloc (3 * SO->EL->N_EL, sizeof(float));
06600    SPI->isEdgeInters = (SUMA_Boolean *) SUMA_calloc (SO->EL->N_EL, sizeof(SUMA_Boolean));
06601    SPI->IntersTri = (int *) SUMA_calloc (SO->N_FaceSet, sizeof(int));
06602    SPI->isNodeInMesh = (SUMA_Boolean *) SUMA_calloc (SO->N_Node, sizeof(SUMA_Boolean));
06603    SPI->isTriHit = (SUMA_Boolean *) SUMA_calloc (SO->N_FaceSet, sizeof(SUMA_Boolean));
06604 
06605    if (!SPI->IntersEdges || !SPI->IntersTri || !SPI->IntersNodes || !SPI->isTriHit || !SPI->isEdgeInters)
06606       {
06607          fprintf (SUMA_STDERR, "Error %s: Could not allocate \n", FuncName);
06608          SUMA_RETURN (SPI);
06609       }
06610    
06611    SPI->N_IntersEdges = 0;
06612    SPI->N_IntersTri = 0;
06613    SPI->N_NodesInMesh = 0;  
06614    SUMA_RETURN (SPI);
06615 }

00624 {
00625    static char FuncName[] = {"SUMA_ApplyAffine"};
00626    float **XYZo, **Mr, **XYZn, D[3];
00627    int i, i3, idbg = 0;
00628    SUMA_Boolean LocalHead = NOPE;
00629    
00630    SUMA_ENTRY;
00631    
00632    if (!NodeList || N_Node <=0) { 
00633       SUMA_SL_Err("Bad Entries.\n");
00634       SUMA_RETURN(NOPE);
00635    }
00636    
00637    Mr = (float **)SUMA_allocate2D(3, 3, sizeof(float));
00638    XYZn = (float **)SUMA_allocate2D(3, 1, sizeof(float));
00639    XYZo = (float **)SUMA_allocate2D(3, 1, sizeof(float));
00640    
00641    SUMA_LH("Forming Mr");
00642    Mr[0][0] = M[0][0]; Mr[0][1] = M[0][1]; Mr[0][2] = M[0][2]; 
00643    Mr[1][0] = M[1][0]; Mr[1][1] = M[1][1]; Mr[1][2] = M[1][2]; 
00644    Mr[2][0] = M[2][0]; Mr[2][1] = M[2][1]; Mr[2][2] = M[2][2];
00645    D[0] = M[0][3]; D[1] = M[1][3]; D[2] = M[2][3];
00646    
00647    SUMA_LH("Transforming");
00648    if (LocalHead ) {
00649       i3 = 3*idbg;
00650       fprintf (SUMA_STDERR,"In: %f %f %f\n", NodeList[i3], NodeList[i3+1], NodeList[i3+2]);
00651    }
00652    for (i=0; i< N_Node; ++i) {
00653       i3 = 3 * i;
00654       if (!center) {
00655          XYZo[0][0] = NodeList[i3]; XYZo[1][0] = NodeList[i3+1]; XYZo[2][0] = NodeList[i3+2];
00656       } else {
00657          XYZo[0][0] = NodeList[i3] - center[0]; XYZo[1][0] = NodeList[i3+1] - center[1]; XYZo[2][0] = NodeList[i3+2] - center[2];
00658       }   
00659 
00660       SUMA_MULT_MAT(Mr, XYZo, XYZn, 3, 3, 1, float,float,float);
00661       
00662       if (!center) { 
00663          NodeList[i3] = XYZn[0][0]+D[0]; NodeList[i3+1] = XYZn[1][0]+D[1]; NodeList[i3+2] = XYZn[2][0]+D[2]; 
00664       } else {
00665          NodeList[i3] = XYZn[0][0]+D[0] + center[0]; NodeList[i3+1] = XYZn[1][0]+D[1]+ center[1]; NodeList[i3+2] = XYZn[2][0]+D[2]+ center[2]; 
00666       }
00667       
00668    }
00669    if (LocalHead ) {
00670       i3 = 3*idbg;
00671       fprintf (SUMA_STDERR,"Out: %f %f %f\n", NodeList[i3], NodeList[i3+1], NodeList[i3+2]);
00672    }
00673    SUMA_LH("Done");
00674    
00675    SUMA_free2D((char**)Mr, 3);
00676    SUMA_free2D((char**)XYZn, 3);
00677    SUMA_free2D((char**)XYZo, 3);
00678    
00679    SUMA_RETURN(YUP);
00680 }

01434 {
01435    static char FuncName[]={"SUMA_AreaDiff"};
01436    double da, *fp, Dr, A;
01437    static int ncall=0;
01438    int i;
01439    static double Rref = 0.0, Aref = 0.0;
01440    SUMA_SurfaceObject *SO, *SOref;
01441    SUMA_COMM_STRUCT *cs=NULL;
01442    SUMA_AreaDiffDataStruct *fdata = (SUMA_AreaDiffDataStruct*)fvdata ;
01443    SUMA_Boolean LocalHead = NOPE;
01444 
01445    SUMA_ENTRY;
01446    
01447    if (!fdata) {
01448       SUMA_LH("Reset");
01449       Rref = 0.0; Aref = 0.0;
01450       ncall = 0;
01451       SUMA_RETURN(0.0);
01452    }
01453    
01454    SO = fdata->SO;
01455    SOref = fdata->SOref;
01456    cs = fdata->cs;
01457    
01458    if (!ncall) {
01459       SUMA_LH("Initializing, calculating Aref and Rref");
01460       Aref = fdata->Aref;
01461       Rref = fdata->Rref;
01462       if (LocalHead) { fprintf(SUMA_STDERR,"%s: Reference volume = %f, radius = %f \n", FuncName, Aref, Rref); }
01463       if (cs->Send) { /* send the first monster (it's SOref "in SUMA" that's being modified on the fly) */
01464          if (!SUMA_SendToSuma (SOref, cs, (void *)SO->NodeList, SUMA_NODE_XYZ, 1)) {
01465          SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted.");
01466          }
01467       }
01468    }
01469    
01470    A = SUMA_NewAreaAtRadius(SO, r, Rref, fdata->tmpList);
01471    da = Aref - A; /* the area difference */
01472    if (LocalHead) {
01473       fprintf(SUMA_STDERR,"%s: Call %d, A = %f, Aref = %f, da = %f\n", FuncName,ncall, A, Aref, da);
01474    }
01475       
01476    /* need an update ? */
01477    if (cs->Send) { /* send the update (it's SOref "in SUMA" that's being modified on the fly) */
01478       if (!SUMA_SendToSuma (SOref, cs, (void *)fdata->tmpList, SUMA_NODE_XYZ, 1)) {
01479       SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted.");
01480       }
01481    }
01482 
01483    ++ncall;
01484    
01485    SUMA_RETURN(da);
01486 }

06336 {
06337    static char FuncName[]={"SUMA_AssignTriBranch"};
06338    int *IntersEdgesCopy = NULL, N_IntersEdgesCopy, i_Branch, E1, kedge, i, 
06339          N_iBranch[SUMA_MAX_BRANCHES], NBlist[SUMA_MAX_BRANCHES], iBranch = 0, 
06340          N_Branch, Bcnt, ilist, j, ivisit, *VisitationOrder, TriCheck;
06341    SUMA_Boolean Local_IntersEdgesCopy = NOPE;
06342    int *TriBranch = NULL; /*!< Vector of SO->EL->N_EL / 3 elements but with only N_IntersTri meaningful values. If TriBranch[j] = b
06343                      then triangle j in SO->FaceSet is a member of Branch b. Branch numbering starts at 1 and may not be consecutive. 
06344                      A branch is a collection of connected triangles and may form a closed loop */
06345    SUMA_TRI_BRANCH *Bv = NULL;
06346    SUMA_Boolean LocalHead = NOPE;
06347    
06348    SUMA_ENTRY;
06349 
06350 
06351    i_Branch = 0;
06352    
06353    /* make a copy of SPI->N_IntersEdges since this vector will be modified and might be needed later, 
06354       might make that optional later and just copy pointers if IntersEdgesCopy is not needed elsewhere.
06355       IntersEdgesCopy flag is used to decide on freeing IntersEdgesCopy or not.*/
06356 
06357    VisitationOrder = (int *)SUMA_calloc (SO->N_FaceSet, sizeof (int)); /* keeps track of the order in which triangles are visited. This is used for creating branches with the proper sequence */
06358    TriBranch = (int *)SUMA_calloc (SO->EL->N_EL / 3,  sizeof(int));
06359    
06360    if (!VisitationOrder || !TriBranch) {
06361       fprintf (SUMA_STDERR, "Error %s: Failed to allocate.\n", FuncName);
06362       if (TriBranch) SUMA_free(TriBranch);
06363       if (VisitationOrder) SUMA_free(VisitationOrder);
06364       SUMA_RETURN (NULL);   
06365    }
06366    
06367    N_IntersEdgesCopy = SPI->N_IntersEdges;
06368    if (DoCopy) {
06369       IntersEdgesCopy = (int *) SUMA_calloc (N_IntersEdgesCopy, sizeof (int));
06370       Local_IntersEdgesCopy = YUP;
06371       for (i=0; i < N_IntersEdgesCopy; ++i) {
06372          IntersEdgesCopy[i] = SPI->IntersEdges[i];
06373       }
06374    }else {
06375       Local_IntersEdgesCopy = NOPE;
06376       IntersEdgesCopy = SPI->IntersEdges;
06377    }
06378 
06379    if (!IntersEdgesCopy) {
06380      fprintf (SUMA_STDERR, "Error %s: Failed to allocate for or receive IntersEdgesCopy.\n", FuncName);
06381      if (TriBranch) SUMA_free(TriBranch);
06382      if (VisitationOrder) SUMA_free(VisitationOrder);
06383      SUMA_RETURN (NULL);
06384    }
06385 
06386    ivisit = 0;
06387    while (N_IntersEdgesCopy) {
06388    
06389       if (!i_Branch && Nx >= 0) {
06390          /* start from edge containing Nx, this is only done at the starting point (i_Branch = 0) */
06391          E1 = -1;
06392          i=0;
06393          while (i < N_IntersEdgesCopy && E1 < 0) {
06394             if ( (SO->EL->EL[IntersEdgesCopy[i]][0] == Nx) || (SO->EL->EL[IntersEdgesCopy[i]][1] == Nx) ) {
06395                E1 = IntersEdgesCopy[i];
06396                kedge = i;
06397             }  
06398             ++i;
06399          }
06400       }else {
06401          /* no starting orders, start from any decent edge */
06402          /* find an edge with one hosting triangle */
06403          E1 = SUMA_Find_Edge_Nhost (SO->EL, IntersEdgesCopy, N_IntersEdgesCopy, &kedge, 1);
06404       }      
06405 
06406       if (E1 < 0) { /* no such edge found, take first edge in InInter */
06407             kedge = 0;
06408             E1 = IntersEdgesCopy[kedge];
06409             if (LocalHead) fprintf (SUMA_STDERR, "%s: No 1 host edge edge found.\n", FuncName);
06410       }else {
06411             if (LocalHead) fprintf (SUMA_STDERR, "%s: Found edge.\n", FuncName);
06412       }
06413       
06414       /* remove this edge from the list */
06415       --(N_IntersEdgesCopy);
06416       if (LocalHead) fprintf (SUMA_STDERR, "%s: kedge = %d, N_IntersEdgesCopy = %d.\n", FuncName, kedge, N_IntersEdgesCopy);
06417       IntersEdgesCopy[kedge] = IntersEdgesCopy[N_IntersEdgesCopy];
06418 
06419       /* start a new i_Branch - All i_Branch indices must be > 0*/
06420       ++i_Branch;   
06421       if (i_Branch > SUMA_MAX_BRANCHES-1) {
06422          fprintf (SUMA_STDERR, "Error %s: No more than %d branches allowed.\n", FuncName, SUMA_MAX_BRANCHES);
06423          SUMA_RETURN (NULL); 
06424       } 
06425       
06426       /* mark the triangle containing E1 */
06427       if (LocalHead) fprintf (SUMA_STDERR, "%s: Marking triangle %d with branch %d.\n", FuncName, SO->EL->ELps[E1][1], i_Branch);
06428       TriBranch[SO->EL->ELps[E1][1]] = i_Branch;
06429       VisitationOrder[ivisit] = SO->EL->ELps[E1][1]; ++ivisit;
06430       
06431       if (LocalHead) fprintf (SUMA_STDERR, "%s: Called recursive SUMA_Mark_Tri.\n", FuncName);
06432       if (!SUMA_Mark_Tri (SO->EL, E1, i_Branch, TriBranch, IntersEdgesCopy, &(N_IntersEdgesCopy), VisitationOrder, &ivisit)) {
06433          fprintf (SUMA_STDERR, "Error %s: Failed in SUMA_Mark_Tri.\n", FuncName);
06434       }
06435       if (LocalHead) fprintf (SUMA_STDERR, "%s: Returned from recursive SUMA_Mark_Tri.\n", FuncName);
06436 
06437       /* repeat till all edges are used up */
06438    }
06439 
06440    if (Local_IntersEdgesCopy) {
06441       if (LocalHead) fprintf (SUMA_STDERR, "%s: freeing IntersEdgesCopy.\n", FuncName);
06442       SUMA_free(IntersEdgesCopy); 
06443    }else {
06444       /* also change N_IntersEdges */
06445       SPI->N_IntersEdges = N_IntersEdgesCopy;
06446    }
06447 
06448    /* SUMA_disp_dvect (TriBranch, SO->N_FaceSet);  */
06449 
06450    N_Branch = i_Branch;
06451 
06452    /* determine the number of branch elements to allocate for - IDIOT PROOF, doing it in the recursive function SUMA_Mark_Tri was annoying*/
06453    for (i=0; i <= N_Branch; ++i) N_iBranch[i] = 0; /* remember, Branch numbering starts at 1 */
06454    if (LocalHead) fprintf (SUMA_STDERR, "%s: Searching all %d intersected triangles.\n", FuncName, SPI->N_IntersTri);
06455    Bcnt = 0;
06456    for (i=0; i < SO->N_FaceSet; ++i) {
06457       if (TriBranch[i]) {
06458          /* fprintf (SUMA_STDERR, "%d:%d\t", TriBranch[i], N_iBranch[TriBranch[i]]); */
06459          ++Bcnt;
06460          N_iBranch[TriBranch[i]] = N_iBranch[TriBranch[i]] + 1;
06461       }
06462    }
06463    
06464    #if 0
06465       fprintf (SUMA_STDERR, "Values in N_iBranch, idiot proof:\n");
06466       SUMA_disp_dvect (N_iBranch, N_Branch+1);
06467       fprintf (SUMA_STDERR, "\n");
06468    #endif
06469    
06470    if (LocalHead) fprintf (SUMA_STDERR, "%s: Found %d triangles belonging to a branch out of %d intersected triangles.\n", FuncName, Bcnt, SPI->N_IntersTri);
06471      
06472    /* Now you want to create a vector of N_Branches to represent the intersection */
06473    Bv = (SUMA_TRI_BRANCH *) SUMA_malloc (sizeof(SUMA_TRI_BRANCH)*(N_Branch+1)); /* you should only need N_Branch, but that 1 won't hurt ...*/
06474    if (!Bv) {
06475       fprintf (SUMA_STDERR, "Error %s: Could not allocate for Bv.\n", FuncName);
06476       SUMA_RETURN (NULL);
06477    }
06478    
06479    /* initialize allocated Bv elements */
06480    for (i=0; i<= N_Branch; ++i) { /* You have allocated for N_Branch+1*/
06481       Bv[i].list = NULL;
06482       Bv[i].N_list = 0;
06483    } 
06484    
06485    Bcnt = 0;
06486    for (i=0; i<= N_Branch; ++i) { /* Branch numbering starts at 1 */
06487       if (N_iBranch[i]) {
06488          /* something in that branch, allocate and initialize*/
06489          if (LocalHead) fprintf (SUMA_STDERR, "%s: Allocating for %d elements, Old Branch %d, New Branch %d.\n", FuncName, N_iBranch[i], i, Bcnt); 
06490          Bv[Bcnt].list = (int *) SUMA_calloc (N_iBranch[i]+1, sizeof(int));
06491          Bv[Bcnt].N_list = N_iBranch[i];
06492          Bv[Bcnt].iBranch = Bcnt;
06493          NBlist[i] = Bcnt; /* store new indexing for Branches */
06494          ++Bcnt;
06495       }
06496       
06497    }
06498    
06499    /* store the total number of used branches */
06500    *BranchCount = Bcnt;
06501    
06502    /* now fill up the branches*/
06503    if (LocalHead) fprintf (SUMA_STDERR, "%s: Filling up branches...\n", FuncName);
06504    for (i=0; i <= N_Branch; ++i) N_iBranch[i] = 0; /* now use this vector as a counter for the filling at each new branch index */
06505    for (i=0; i < SPI->N_IntersTri; ++i) { /* only go over visited triangles */
06506       TriCheck = TriBranch[VisitationOrder[i]]; /* returns the branch number of triangle VisitationOrder[i] */
06507       if (TriCheck) {
06508          Bcnt = NBlist[TriCheck]; /* get the new branch number from the original (old) branch number */
06509          #if 0
06510          fprintf (SUMA_STDERR,"%s: Tricheck = %d\n", FuncName, TriCheck); */
06511          if (Bcnt >= *BranchCount) {
06512             fprintf (SUMA_STDERR, "\aError %s: BranchCount = %d <= Bcnt = %d.\n", FuncName, *BranchCount, Bcnt);
06513          }
06514          if (N_iBranch[Bcnt] >= Bv[Bcnt].N_list) {
06515             fprintf (SUMA_STDERR, "\aError %s: Bcnt = %d. N_iBranch[Bcnt] = %d >= Bv[Bcnt].N_list = %d\n", FuncName, Bcnt, N_iBranch[Bcnt], Bv[Bcnt].N_list);
06516          }
06517          #endif
06518          Bv[Bcnt].list[N_iBranch[Bcnt]] = VisitationOrder[i]; /* store the index of the visited triangle in that branch */
06519          N_iBranch[Bcnt] += 1; /* store the number of elements in that branch */
06520       }
06521    }
06522    
06523    if (LocalHead) fprintf (SUMA_STDERR, "%s: freeing ...\n", FuncName);
06524    if (VisitationOrder) SUMA_free(VisitationOrder);
06525    if (TriBranch) SUMA_free(TriBranch);
06526    SUMA_RETURN (Bv);    
06527 }

01548                                                                                                                           {
01549    static char FuncName[]={"SUMA_BinaryZeroSearch"};
01550    int Niter;
01551    double x, fx;
01552    SUMA_Boolean done;
01553    SUMA_Boolean LocalHead = NOPE;
01554 
01555    SUMA_ENTRY;
01556    
01557    if (Nitermax < 0) Nitermax = 1000;
01558 
01559    x = 0.0;
01560    Niter = 0;
01561    done = NOPE;
01562    while(!done && Niter < Nitermax) {
01563       x = (a+b)/2.0;
01564       fx = (*f)(x, fdata);
01565       if (LocalHead) fprintf(SUMA_STDERR,"%s: %d\ta=%.4f\tb=%.4f\tx=%.4f\tfx=%.4f\n", 
01566                                           FuncName, Niter, a, b, x, fx);
01567       if (fx < 0) a = x;
01568       else b = x;
01569       if (fabs(fx) < tol) done = YUP;
01570       ++Niter;
01571    }
01572    
01573    if (!done) {
01574       SUMA_SL_Warn(  "Reached iteration limit\n"
01575                      "without converging.\n");
01576    }
01577    
01578    SUMA_RETURN(x);
01579 }

00946 {
00947    static char FuncName[]={"SUMA_CalcNeighbDist"};
00948    float **DistFirstNeighb=NULL, *a, *b;
00949    int i, j;
00950    static SUMA_Boolean LocalHead = NOPE;
00951    
00952    SUMA_ENTRY;
00953    
00954    if (!SO) { SUMA_RETURN(NULL); }
00955    if (!SO->FN) { SUMA_RETURN(NULL); }
00956    
00957    DistFirstNeighb = (float **)SUMA_allocate2D(SO->FN->N_Node, SO->FN->N_Neighb_max, sizeof(float));
00958    if (!DistFirstNeighb) {
00959       SUMA_SL_Crit("Failed to allocate for DistFirstNeighb");
00960       SUMA_RETURN(NULL);
00961    }
00962    for (i=0; i < SO->FN->N_Node; ++i) {
00963       a = &(SO->NodeList[3*SO->FN->NodeId[i]]);
00964       for (j=0; j < SO->FN->N_Neighb[i]; ++j) {
00965          b = &(SO->NodeList[3*SO->FN->FirstNeighb[i][j]]);
00966          SUMA_SEG_LENGTH(a, b, DistFirstNeighb[i][j]);
00967          if (SO->FN->NodeId[i] == 5 && SO->FN->FirstNeighb[i][j] == 133092) {
00968             fprintf (SUMA_STDERR, "%f %f %f\n%f %f %f\n%f\n", 
00969                SO->NodeList[3*SO->FN->NodeId[i]], SO->NodeList[3*SO->FN->NodeId[i]+1], SO->NodeList[3*SO->FN->NodeId[i]+2],
00970                SO->NodeList[3*SO->FN->FirstNeighb[i][j]], SO->NodeList[3*SO->FN->FirstNeighb[i][j]+1], 
00971                SO->NodeList[3*SO->FN->FirstNeighb[i][j]+2], DistFirstNeighb[i][j]);
00972          }
00973       }
00974    }
00975    
00976    SUMA_RETURN (DistFirstNeighb);
00977 }

03072 {
03073    static char FuncName[]={"SUMA_Chung_Smooth"};
03074    float *fout_final = NULL, *fbuf=NULL, *fin=NULL, *fout=NULL, *fin_next = NULL;
03075    float delta_time, fp, dfp, fpj, minfn=0.0, maxfn=0.0;
03076    int n , k, j, niter, vnk, os;
03077    SUMA_Boolean LocalHead = NOPE;
03078    
03079    SUMA_ENTRY;
03080 
03081    if (N_iter % 2) {
03082       SUMA_SL_Err("N_iter must be an even number\n");
03083       SUMA_RETURN(NULL);
03084    }
03085    
03086    if (!SO || !wgt || !fin_orig) {
03087       SUMA_SL_Err("NULL SO or wgt or fin_orig\n");
03088       SUMA_RETURN(NULL);
03089    }
03090    if (!SO->FN) {
03091       SUMA_SL_Err("NULL SO->FN\n");
03092       SUMA_RETURN(NULL);
03093    }
03094    
03095    if (vpn < 1) {
03096       SUMA_SL_Err("vpn < 1\n");
03097       SUMA_RETURN(NULL);
03098    }  
03099    
03100    if (fout_final_user == fin_orig) {
03101       SUMA_SL_Err("fout_final_user == fin_orig");
03102       SUMA_RETURN(NULL);
03103    }
03104    
03105    if (!fout_final_user) { /* allocate for output */
03106       fout_final = (float *)SUMA_calloc(SO->N_Node * vpn, sizeof(float));
03107       if (!fout_final) {
03108          SUMA_SL_Crit("Failed to allocate for fout_final\n");
03109          SUMA_RETURN(NULL);
03110       }
03111    }else {
03112       fout_final = fout_final_user; /* pre-allocated */
03113    }
03114    
03115    /* allocate for buffer */
03116    fbuf = (float *)SUMA_calloc(SO->N_Node * vpn, sizeof(float));
03117    if (!fbuf) {
03118       SUMA_SL_Crit("Failed to allocate for fbuf\n");
03119       SUMA_RETURN(NULL);
03120    }
03121    
03122    
03123    if (cs->Send) { /* send the first monster */
03124       if (!SUMA_SendToSuma (SO, cs, (void *)fin_orig, SUMA_NODE_RGBAb, 1)) {
03125          SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted.");
03126       }
03127    }
03128    
03129    fin_next = fin_orig;
03130    delta_time= (FWHM * FWHM)/(16*N_iter*log(2));
03131    switch (d_order) {
03132       case SUMA_COLUMN_MAJOR:
03133          for (niter=0; niter < N_iter; ++niter) {
03134             if ( niter % 2 ) { /* odd */
03135                fin = fin_next; /* input from previous output buffer */
03136                fout = fout_final; /* results go into final vector */
03137                fin_next = fout_final; /* in the next iteration, the input is from fout_final */
03138             } else { /* even */
03139                /* input data is in fin_new */
03140                fin = fin_next;
03141                fout = fbuf; /* results go into buffer */
03142                fin_next = fbuf; /* in the next iteration, the input is from the buffer */
03143             }
03144             for (k=0; k < vpn; ++k) {
03145                os = SO->N_Node*k;   /* node value indexing offset */
03146                for (n=0; n < SO->N_Node; ++n) {
03147                   vnk = n+os; /* index of kth value at node n */
03148                   fp = fin[vnk]; /* kth value at node n */
03149                   dfp = 0.0;
03150                   if (SO->FN->N_Neighb[n]) minfn = maxfn = fin[SO->FN->FirstNeighb[n][0]+os];
03151                   for (j=0; j < SO->FN->N_Neighb[n]; ++j) {
03152                      fpj = fin[SO->FN->FirstNeighb[n][j]+os]; /* value at jth neighbor of n */
03153                      if (fpj < minfn) minfn = fpj;
03154                      if (fpj > maxfn) maxfn = fpj;
03155                      dfp += wgt[n][j] * (fpj - fp); 
03156                   }/* for j*/
03157                   fout[vnk] = fin[vnk] + delta_time * dfp;
03158                   if (fout[vnk] < minfn) fout[vnk] = minfn;
03159                   if (fout[vnk] > maxfn) fout[vnk] = maxfn;
03160                }/* for n */ 
03161                  
03162                if (cs->Send) {
03163                   if (!SUMA_SendToSuma (SO, cs, (void *)fout, SUMA_NODE_RGBAb, 1)) {
03164                      SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted.");
03165                   }
03166                }
03167             } /* for k */
03168          }/* for niter */
03169          break;
03170       case SUMA_ROW_MAJOR:
03171          SUMA_SL_Err("Row Major not implemented");
03172          SUMA_RETURN(NULL);
03173          break;
03174       default:
03175          SUMA_SL_Err("Bad Major, very bad.\n");
03176          SUMA_RETURN(NULL);
03177          break;
03178    }
03179    
03180    if (fbuf) SUMA_free(fbuf); fbuf = NULL;
03181                
03182    SUMA_RETURN(fout);
03183 }

02096 {
02097    static char FuncName[]={"SUMA_Chung_Smooth_Weights"};
02098    float **wgt=NULL, *coord_nbr=NULL, *cotan=NULL, *tfp=NULL;
02099    float dv[3], p[3], q[3];
02100    float area, area_p, area_q, dot_p, dot_q;
02101    int i, j, k, n, j3p1, j3m1, n3, j3=0, nj, nj3, i_dbg;
02102    SUMA_Boolean LocalHead = NOPE;
02103    
02104    SUMA_ENTRY;
02105 
02106    if (!SO) {
02107       SUMA_SL_Err("Null SO");
02108       SUMA_RETURN(NULL);
02109    }
02110    if (!SO->FN) {
02111       SUMA_SL_Err("Null SO->FN");
02112       SUMA_RETURN(NULL);
02113    }
02114    
02115    i_dbg = -1; /* index of debugging node, set to -1 for no debugging */
02116    /* in the demo mesh that Moo Chung gave me, 
02117    Node 17 has the following neighbors in ccw order:
02118    231 230 250 261 239 236 - 231 230 250 261 239 236 
02119    Set i_dbg = 17 and turn on LocalHead to get a confirmation of this
02120    by this function*/
02121    /* implement the non-parametric weight estimation method */
02122    wgt = (float **)SUMA_allocate2D(SO->N_Node, SO->FN->N_Neighb_max, sizeof(float));  /* vector of node weights */
02123    coord_nbr = (float *)SUMA_malloc((SO->FN->N_Neighb_max + 2) * sizeof(float) * 3); /* vector of neighboring node coordinates */ 
02124    cotan = (float *)SUMA_malloc(SO->FN->N_Neighb_max * sizeof(float)); 
02125    if (!wgt || !coord_nbr || !cotan) {
02126       SUMA_SL_Crit("Failed to allocate for wgt &/|coord_nbr &/|cotan");
02127       SUMA_RETURN(NULL);
02128    }
02129    
02130    for (n=0; n < SO->N_Node; ++n) {
02131       n3 = 3 * n;
02132       /* translate the coordinates of the neighboring nodes to make n be the origin */
02133       for (j=0; j<SO->FN->N_Neighb[n]; ++j) {
02134          j3 = 3 * (j+1);
02135          nj = SO->FN->FirstNeighb[n][j]; nj3 = 3 * nj;
02136          coord_nbr[j3] = SO->NodeList[nj3] - SO->NodeList[n3];
02137          coord_nbr[j3+1] = SO->NodeList[nj3+1] - SO->NodeList[n3+1];
02138          coord_nbr[j3+2] = SO->NodeList[nj3+2] - SO->NodeList[n3+2];
02139       }  /* for j */
02140       /* padd with last neighbor at the very beginning and 1st neighbor at the end 
02141          in matlab: coord_nbr = [coord_nbr(:,last_neighb) coord_nbr coord_nbr(:,first_neighb)];
02142       */   
02143       for (k=0; k < 3; ++k) coord_nbr[k] = coord_nbr[j3+k];  
02144       j3 = 3 * ( SO->FN->N_Neighb[n] + 1);
02145       for (k=0; k < 3; ++k) coord_nbr[j3+k] = coord_nbr[3+k];
02146       if (LocalHead && n == i_dbg) { SUMA_disp_vect (coord_nbr, 3 * (SO->FN->N_Neighb[n] + 2)) ;  }
02147       
02148       
02149       area = 0.0;
02150       for (j=1; j<=SO->FN->N_Neighb[n]; ++j) { 
02151          j3 = 3 * j; j3p1 = 3 * (j+1); j3m1 = 3 * (j-1);
02152          for (k=0; k < 3; ++k) dv[k] = coord_nbr[j3p1+k] - coord_nbr[j3+k]; 
02153          tfp = &(coord_nbr[j3p1]);
02154          dot_p = SUMA_MT_DOT (tfp, dv);
02155          SUMA_MT_CROSS(p, tfp, dv);
02156          for (k=0; k < 3; ++k) dv[k] = coord_nbr[j3m1+k] - coord_nbr[j3+k]; 
02157          tfp = &(coord_nbr[j3m1]);
02158          dot_q = SUMA_MT_DOT (tfp, dv);
02159          SUMA_MT_CROSS(q, tfp, dv);
02160          
02161          SUMA_NORM(area_p, p); 
02162          SUMA_NORM(area_q, q); 
02163          
02164          cotan[j-1] = dot_p/area_p + dot_q/area_q;
02165          area += area_p/2.0;
02166          if (LocalHead && n == i_dbg) {
02167             fprintf (SUMA_STDERR,"[%d->%d] area_p, area_q = %f, %f\n",
02168                                   n, SO->FN->FirstNeighb[n][j-1],
02169                                   area_p / 2.0,  area_q / 2.0);
02170          }
02171       }
02172       
02173       for (j=0; j<SO->FN->N_Neighb[n]; ++j) {
02174          wgt[n][j] = cotan[j]/area;
02175       }
02176       if (LocalHead && n == i_dbg) {
02177          fprintf (SUMA_STDERR,"%s: Weight Results for neighbors of %d (matlab node %d):\n",
02178                               FuncName, n, n+1);
02179          SUMA_disp_vect (wgt[n], SO->FN->N_Neighb[n]);
02180       }
02181    }  /* for n */
02182 
02183    /* free local variables */
02184    if (coord_nbr) SUMA_free(coord_nbr); coord_nbr = NULL;
02185    if (cotan) SUMA_free(cotan); cotan = NULL;
02186 
02187    SUMA_RETURN(wgt);
02188 }

00133 {
00134    static char FuncName[]={"SUMA_CreateVTI"};
00135    SUMA_VTI *vti = NULL;
00136    
00137    SUMA_ENTRY;
00138    if (!N_TriIndex) {
00139       SUMA_SL_Err("Nothing to do !");
00140       SUMA_RETURN(vti);
00141    }
00142    
00143    vti = (SUMA_VTI *)SUMA_malloc(sizeof(SUMA_VTI));
00144    vti->N_TriIndex = N_TriIndex;
00145    if (TriIndex ) {
00146       vti->TriIndex = TriIndex;
00147    }else {
00148       /* create empty copy */
00149       vti->TriIndex = (int *)SUMA_calloc(N_TriIndex, sizeof(int));
00150       if (!vti->TriIndex) {
00151          SUMA_SL_Crit("Failed to allocate for vti->TriIndex");
00152          SUMA_RETURN(NULL);
00153       }
00154    }
00155    vti->N_IntersectedVoxels = (int *)SUMA_calloc(N_TriIndex, sizeof(int));
00156    vti->IntersectedVoxels = (int **)SUMA_calloc(N_TriIndex, sizeof(int*));
00157    if (!vti->N_IntersectedVoxels || !vti->IntersectedVoxels) {
00158          SUMA_SL_Crit("Failed to allocate for vti's innerds");
00159          SUMA_RETURN(NULL);
00160    }
00161 
00162    SUMA_RETURN(vti);
00163 }

06838 {
06839    static char FuncName[] = {"SUMA_Dijkstra"};
06840    SUMA_Boolean LocalHead = NOPE;
06841    float *L = NULL, Lmin = -1.0, le = 0.0, DT_DIJKSTRA;
06842    int i, iw, iv, v, w, N_Neighb, *Path = NULL;
06843    struct  timeval  start_time;
06844    SUMA_DIJKSTRA_PATH_CHAIN *DC = NULL, *DCi, *DCp;
06845    SUMA_Boolean Found = NOPE;
06846    /* variables for method 2 */
06847    int N_Lmins, *vLmins, *vLocInLmins, iLmins, ReplacingNode, ReplacedNodeLocation;
06848    float *Lmins; 
06849    
06850    
06851    SUMA_ENTRY;
06852    
06853    *Lfinal = -1.0;
06854    *N_Path = 0;
06855    
06856    /* make sure Both Nx and Ny exist in isNodeInMesh */
06857    if (!isNodeInMesh[Nx]) {
06858       fprintf (SUMA_STDERR,"\aError %s: Node %d (Nx) is not in mesh.\n", FuncName, Nx);
06859       SUMA_RETURN (NULL);
06860    }  
06861    if (!isNodeInMesh[Ny]) {
06862       fprintf (SUMA_STDERR,"\aError %s: Node %d (Ny) is not in mesh.\n", FuncName, Ny);
06863       SUMA_RETURN (NULL);
06864    }
06865 
06866    if (!SO->FN) {
06867       fprintf (SUMA_STDERR, "Error %s: SO does not have FN structure.\n", FuncName);
06868       SUMA_RETURN (NULL);
06869    }
06870 
06871    if (LocalHead) {
06872       /* Start timer for next function */
06873       SUMA_etime(&start_time,0);      
06874    }
06875    
06876    /* allocate for chain */
06877    DC = (SUMA_DIJKSTRA_PATH_CHAIN *) SUMA_malloc (sizeof(SUMA_DIJKSTRA_PATH_CHAIN) * SO->N_Node);
06878    if (!DC) {
06879       fprintf (SUMA_STDERR, "Error %s: Could not allocate. \n", FuncName);
06880       SUMA_RETURN (NULL);
06881    }
06882    
06883    switch (Method_Number) {
06884    
06885       case 0:  /* Method 0, Brute force */
06886          /* allocate for vertices labels */
06887          L = (float *) SUMA_calloc (SO->N_Node, sizeof (float));
06888          if (!L) {
06889             fprintf (SUMA_STDERR, "Error %s: Failed to allocate.\n", FuncName);
06890             SUMA_free(DC);
06891             SUMA_RETURN (NULL);
06892          }
06893 
06894          /* label all vertices with very large numbers, initialize path previous pointers to null */
06895          for (i=0; i < SO->N_Node; ++i) {
06896             L[i] = LARGE_NUM;   
06897             DC[i].Previous = NULL;
06898          }
06899          /* label starting vertex with 0 */
06900          L[Nx] = 0.0;
06901          Lmin = 0.0;
06902          v = Nx;
06903          *Lfinal = -1.0;
06904          /* initialize path at Nx */
06905          DC[Nx].Previous = NULL;
06906          DC[Nx].node = Nx;
06907          DC[Nx].le = 0.0;
06908          DC[Nx].order = 0;
06909          *N_Path = 0;
06910          /* Brute force method */
06911          do {
06912             /* find v in Mesh / L(v) is minimal */
06913             /* this sucks up a lot of time because it is searching the entire set of SO->N_Node instead of the one that was intersected only.
06914             This can be sped up, considerably */
06915             SUMA_MIN_LOC_VEC(L, SO->N_Node, Lmin, v);   /* locates and finds the minimum of L, nodes not in mesh will keep their large values and will not be picked*/
06916             if (!isNodeInMesh[v]) {
06917                fprintf (SUMA_STDERR, "\aERROR %s: Dijkstra derailed. v = %d, Lmin = %f\n. Try another point.", FuncName, v, Lmin);
06918                SUMA_free (L);
06919                SUMA_free(DC);
06920                SUMA_RETURN (NULL); 
06921             }
06922             if (v == Ny) {
06923                if (LocalHead) fprintf (SUMA_STDERR, "%s: Done.\n", FuncName);
06924                *Lfinal = L[v];
06925                Found = YUP;
06926             } else {
06927                N_Neighb = SO->FN->N_Neighb[v];
06928                for (i=0; i < N_Neighb; ++i) {
06929                   w = SO->FN->FirstNeighb[v][i];
06930                   if (isNodeInMesh[w]) {
06931                      iw = 3*w;
06932                      iv = 3*v;
06933                      le = sqrt ( (SO->NodeList[iw] - SO->NodeList[iv]) * (SO->NodeList[iw] - SO->NodeList[iv]) +
06934                                  (SO->NodeList[iw+1] - SO->NodeList[iv+1]) * (SO->NodeList[iw+1] - SO->NodeList[iv+1]) +
06935                                  (SO->NodeList[iw+2] - SO->NodeList[iv+2]) * (SO->NodeList[iw+2] - SO->NodeList[iv+2]) );
06936                      if (L[w] > L[v] + le ) {
06937                         L[w] = L[v] + le;  
06938                         /* update the path */
06939                         DCp = &(DC[v]); /* previous path */
06940                         DC[w].Previous = (void *) DCp;
06941                         DC[w].le = le;
06942                         DC[w].node = w;
06943                         DC[w].order = DCp->order + 1;
06944                      } 
06945                   }
06946                }
06947 
06948                /* remove node v from isNodeInMesh and reset their distance value to a very large one, 
06949                   this way you do not have to reinitialize this variable. */
06950                isNodeInMesh[v] = NOPE;
06951                *N_isNodeInMesh -= 1;
06952                L[v] = LARGE_NUM; 
06953                Found = NOPE;
06954             }
06955          } while (*N_isNodeInMesh > 0 && !Found);
06956 
06957          if (!Found) {
06958             fprintf (SUMA_STDERR, "Error %s: No more nodes in mesh, failed to reach target.\n", FuncName);
06959             SUMA_free (L);
06960             SUMA_free(DC);
06961             SUMA_RETURN (NULL);
06962          }else {
06963             if (LocalHead) fprintf (SUMA_STDERR, "%s: Path between Nodes %d and %d is %f.\n", FuncName, Nx, Ny, *Lfinal);
06964          }
06965 
06966 
06967          if (LocalHead) {
06968             /* stop timer */
06969             DT_DIJKSTRA = SUMA_etime(&start_time,1);
06970             fprintf (SUMA_STDERR, "%s: Method 1- Elapsed time in function %f seconds.\n", FuncName, DT_DIJKSTRA);
06971          }
06972 
06973          SUMA_free(L);
06974          break;
06975 
06976       case 1:  /********* Method 1- faster minimum searching *******************/
06977          if (LocalHead) {
06978             /* Start timer for next function */
06979             SUMA_etime(&start_time,0);      
06980          }
06981 
06982          /* allocate for vertices labels and minimums vectors*/
06983          L = (float *) SUMA_calloc (SO->N_Node, sizeof (float));        /* L[i] = distance to a node i*/
06984          Lmins = (float *) SUMA_calloc (SO->N_Node, sizeof (float));    /* Lmins = vector containing minimum calculated distances to node */
06985          vLmins = (int *) SUMA_calloc (SO->N_Node, sizeof (int));       /* vLmins[i] = index (into L) of the node having a distance Lmins[i] */
06986          vLocInLmins = (int *) SUMA_calloc (SO->N_Node, sizeof (int));  /* vLocInLmin[j] = index (into Lmins) of a node having index j (into L) */
06987 
06988          if (!L || !Lmins || !vLmins || !vLocInLmins) {
06989             fprintf (SUMA_STDERR, "Error %s: Failed to allocate.\n", FuncName);
06990             SUMA_RETURN (NULL);
06991          }
06992 
06993          /* label all vertices with very large numbers and initialize vLocInLmins to -1*/
06994          for (i=0; i < SO->N_Node; ++i) {
06995             L[i] = LARGE_NUM;
06996             Lmins[i] = LARGE_NUM;   
06997             vLocInLmins[i] = -1;            
06998             DC[i].Previous = NULL;
06999          }
07000 
07001          /* label starting vertex with 0 */
07002          L[Nx] = 0.0;
07003          *Lfinal = -1.0;
07004 
07005          /* initialize values of vectors used to keep track of minimum values of L and their corresponding nodes */
07006          Lmins[0] = 0.0;
07007          vLmins[0] = Nx;
07008          vLocInLmins[Nx] = 0;
07009          N_Lmins = 1;
07010 
07011          /* initialize path at Nx */
07012          DC[Nx].Previous = NULL;
07013          DC[Nx].node = Nx;
07014          DC[Nx].le = 0.0;
07015          DC[Nx].order = 0;
07016          *N_Path = 0;
07017          
07018          /* method with efficient tracking of minimum */
07019          if (LocalHead) fprintf (SUMA_STDERR, "%s: about to MIN_LOC ....N_isNodeInMesh = %d\n", FuncName, *N_isNodeInMesh);
07020          do {
07021             /* find v in Mesh / L(v) is minimal */
07022             SUMA_MIN_LOC_VEC(Lmins, N_Lmins, Lmin, iLmins);   /* locates the minimum value in Lmins vector */
07023             v = vLmins[iLmins];   /* get the node for this Lmin value */
07024             if (!isNodeInMesh[v]) {
07025                fprintf (SUMA_STDERR, "\aERROR %s: Dijkstra derailed. v = %d, Lmin = %f\n. Try another point.", FuncName, v, Lmin);
07026                SUMA_free (L);
07027                SUMA_free (Lmins);
07028                SUMA_free(vLmins);
07029                SUMA_free(vLocInLmins);
07030                SUMA_free(DC);
07031                SUMA_RETURN (NULL);
07032             }
07033             #ifdef LOCALDEBUG
07034                fprintf (SUMA_STDERR, "%s: Node v = %d.\n", FuncName, v);
07035             #endif
07036             if (v == Ny) {
07037                if (LocalHead) fprintf (SUMA_STDERR, "%s: Done.\n", FuncName);
07038                *Lfinal = L[v];
07039                Found = YUP;
07040             } else {
07041                N_Neighb = SO->FN->N_Neighb[v];
07042                for (i=0; i < N_Neighb; ++i) {
07043                   w = SO->FN->FirstNeighb[v][i];
07044                   if (isNodeInMesh[w]) {
07045                      iw = 3*w;
07046                      iv = 3*v;
07047                      le = sqrt ( (SO->NodeList[iw] - SO->NodeList[iv]) * (SO->NodeList[iw] - SO->NodeList[iv]) +
07048                                  (SO->NodeList[iw+1] - SO->NodeList[iv+1]) * (SO->NodeList[iw+1] - SO->NodeList[iv+1]) +
07049                                  (SO->NodeList[iw+2] - SO->NodeList[iv+2]) * (SO->NodeList[iw+2] - SO->NodeList[iv+2]) );
07050                      if (L[w] > L[v] + le ) {
07051                         #ifdef LOCALDEBUG
07052                            fprintf (SUMA_STDERR, "%s: L[%d]=%f > L[%d] = %f + le = %f.\n", FuncName, w, L[w], v, L[v], le);
07053                         #endif
07054                         L[w] = L[v] + le; 
07055                         /* update the path */
07056                         DCp = &(DC[v]); /* previous path */
07057                         DC[w].Previous = (void *) DCp;
07058                         DC[w].le = le;
07059                         DC[w].node = w;
07060                         DC[w].order = DCp->order + 1;
07061                         
07062                         if (vLocInLmins[w] < 0) { 
07063                            #ifdef LOCALDEBUG
07064                               fprintf (SUMA_STDERR, "%s: adding entry for w = %d - First Hit. \n", FuncName, w);
07065                            #endif
07066                            Lmins[N_Lmins] = L[w]; /* add this value to Lmins vector */
07067                            vLmins[N_Lmins] = w; /* store the node for this Lmins value */
07068                            vLocInLmins[w] = N_Lmins; /* store where that node is represented in Lmins */
07069                            ++N_Lmins;  /* increment N_Lmins */  
07070                         } else {
07071                            #ifdef LOCALDEBUG
07072                               fprintf (SUMA_STDERR, "%s: modifying entry for w = %d  Second Hit.\n", FuncName, w); */
07073                            #endif
07074                            Lmins[vLocInLmins[w]] = L[w]; /* update value for Lmins */
07075                         }                        
07076                      }else {
07077                         #ifdef LOCALDEBUG
07078                            fprintf (SUMA_STDERR, "%s: L[%d]=%f < L[%d] = %f + le = %f.\n", FuncName, w, L[w], v, L[v], le); */
07079                         #endif
07080                      } 
07081                   }
07082                }
07083 
07084                /* remove node v from isNodeInMesh and reset their distance value to a very large one, 
07085                   this way you do not have to reinitialize this variable. */
07086                isNodeInMesh[v] = NOPE;
07087                *N_isNodeInMesh -= 1;
07088                L[v] = LARGE_NUM; 
07089                Found = NOPE;
07090 
07091                /* also remove the values (by swapping it with last element) for this node from Lmins */
07092                #ifdef LOCALDEBUG
07093                   {
07094                      int kkk;
07095                      fprintf (SUMA_STDERR,"Lmins\tvLmins\tvLocInLmins\n");
07096                      for (kkk=0; kkk < N_Lmins; ++kkk) fprintf (SUMA_STDERR,"%f\t%d\t%d\n", Lmins[kkk], vLmins[kkk], vLocInLmins[vLmins[kkk]] );
07097                
07098                   }
07099                #endif
07100                
07101                if (vLocInLmins[v] >= 0) { /* remove its entry if there is one */
07102                   #ifdef LOCALDEBUG
07103                      fprintf (SUMA_STDERR, "%s: removing node v = %d. N_Lmins = %d\n", FuncName,  v, N_Lmins);
07104                   #endif
07105                   --N_Lmins;
07106                   ReplacingNode = vLmins[N_Lmins];
07107                   ReplacedNodeLocation = vLocInLmins[v];
07108                   Lmins[vLocInLmins[v]] = Lmins[N_Lmins];
07109                   vLmins[vLocInLmins[v]] = vLmins[N_Lmins];
07110                   vLocInLmins[ReplacingNode] = ReplacedNodeLocation;
07111                   vLocInLmins[v] = -1;
07112                   Lmins[N_Lmins] = LARGE_NUM; 
07113                }
07114             }
07115          } while (*N_isNodeInMesh > 0 && !Found);
07116 
07117          if (!Found) {
07118             fprintf (SUMA_STDERR, "Error %s: No more nodes in mesh, failed to reach target %d. NLmins = %d\n", FuncName, Ny, N_Lmins);
07119             SUMA_free (L);
07120             SUMA_free (Lmins);
07121             SUMA_free(vLmins);
07122             SUMA_free(vLocInLmins);
07123             SUMA_free(DC);
07124             SUMA_RETURN (NULL);
07125          }else {
07126             if (LocalHead) fprintf (SUMA_STDERR, "%s: Path between Nodes %d and %d is %f.\n", FuncName, Nx, Ny, *Lfinal);
07127          }
07128 
07129 
07130          if (LocalHead) {
07131             /* stop timer */
07132             DT_DIJKSTRA = SUMA_etime(&start_time,1);
07133             fprintf (SUMA_STDERR, "%s: Method 2- Elapsed time in function %f seconds.\n", FuncName, DT_DIJKSTRA);
07134          }
07135 
07136          SUMA_free(L);
07137          SUMA_free(Lmins);
07138          SUMA_free(vLmins);
07139          SUMA_free(vLocInLmins);
07140          break;   /********** Method 1- faster minimum searching **************/
07141       default: 
07142          fprintf (SUMA_STDERR, "Error %s: No such method (%d).\n", FuncName, Method_Number);
07143          if (DC) SUMA_free(DC);
07144          SUMA_RETURN (NULL);
07145          break;
07146    }
07147    
07148    /* now reconstruct the path */
07149    *N_Path = DC[Ny].order+1;
07150    Path = (int *) SUMA_calloc (*N_Path, sizeof(int));
07151    if (!Path) {
07152       fprintf (SUMA_STDERR, "Error %s: Failed to allocate.\n", FuncName);
07153       if (DC) SUMA_free(DC);
07154       SUMA_RETURN (NULL);
07155    }
07156    
07157    DCi = &(DC[Ny]);
07158    iv = *N_Path - 1;
07159    Path[iv] = Ny;
07160    if (iv > 0) {
07161       do {
07162          --iv;
07163          DCp = (SUMA_DIJKSTRA_PATH_CHAIN *) DCi->Previous;
07164          Path[iv] = DCp->node;
07165          DCi = DCp;
07166       } while (DCi->Previous);
07167    }
07168    
07169    if (iv != 0) {
07170       fprintf (SUMA_STDERR, "Error %s: iv = %d. This should not be.\n", FuncName, iv);
07171    }  
07172    
07173    SUMA_free(DC);
07174    SUMA_RETURN (Path);
07175 }

01701 {
01702    static char FuncName[]={"SUMA_EquateSurfaceAreas"};
01703    int iter, i, iter_max, ndiv;
01704    double a, b, d;
01705    SUMA_AreaDiffDataStruct fdata;
01706    SUMA_Boolean LocalHead = NOPE;
01707    
01708    SUMA_ENTRY;
01709 
01710    if (!SO || !SOref) { SUMA_SL_Err("NULL surfaces"); SUMA_RETURN(NOPE); }
01711    if (SO->N_Node != SOref->N_Node || SO->N_FaceSet != SOref->N_FaceSet) { SUMA_SL_Err("Surfaces not isotopic"); SUMA_RETURN(NOPE); }
01712    
01713    if (LocalHead) {
01714       fprintf(SUMA_STDERR, "%s:\n"
01715                            " SO    Center: %f, %f, %f\n"
01716                            " SOref Center: %f, %f, %f\n"
01717                            , FuncName, 
01718                            SO->Center[0], SO->Center[1], SO->Center[2],
01719                            SOref->Center[0], SOref->Center[1], SOref->Center[2]);  
01720    }
01721      
01722    /* fill up fdata */
01723    fdata.SO = SO; fdata.SOref = SOref; fdata.cs = cs;
01724    fdata.tmpList = (float *)SUMA_malloc(SOref->NodeDim * SOref->N_Node * sizeof(float));
01725    if (!fdata.tmpList) {
01726       SUMA_SL_Err("Failed to allocate");
01727       SUMA_RETURN(0);
01728    }
01729 
01730    if (!SUMA_GetAreaDiffRange(&fdata, &a, &b)) {
01731       SUMA_SL_Err("Failed to get range");
01732       SUMA_RETURN(NOPE);
01733    }
01734    
01735    if (LocalHead) {
01736       fprintf(SUMA_STDERR,"%s:\na = %f\tb=%f\n", FuncName, a, b);
01737    }      
01738    SUMA_BinaryZeroSearch(a, b, SUMA_AreaDiff, &fdata, 500, tol);  
01739    
01740    /* now make the new node list be SO's thingy*/
01741    SUMA_free(SO->NodeList); SO->NodeList = fdata.tmpList; fdata.tmpList = NULL;
01742        
01743    SUMA_RETURN(YUP);
01744 }

01943 {
01944    static char FuncName[]={"SUMA_EquateSurfaceSize"};
01945    int i=0, j=0, cnt = 0, istrt, istp;
01946    struct timeval start_time, start_time_all;
01947    float etime_GetOffset, etime_GetOffset_all, ave_dist= 0.0, dj = 0.0, ave_dist_ref= 0.0, *a=NULL;
01948    float P2[2][3], U[3], Un;
01949    SUMA_GET_OFFSET_STRUCT *OffS = NULL;
01950    SUMA_Boolean LocalHead = NOPE;
01951    
01952    SUMA_ENTRY;
01953 
01954    if (!SO || !SOref) { SUMA_SL_Err("NULL surfaces"); SUMA_RETURN(NOPE); }
01955    if (SO->N_Node != SOref->N_Node || SO->N_FaceSet != SOref->N_FaceSet) { SUMA_SL_Err("Surfaces not isotopic"); SUMA_RETURN(NOPE); }
01956    
01957    if (LocalHead) {
01958       fprintf(SUMA_STDERR, "%s:\n"
01959                            " SO    Center: %f, %f, %f\n"
01960                            " SOref Center: %f, %f, %f\n"
01961                            " max_off = %f\n", FuncName, 
01962                            SO->Center[0], SO->Center[1], SO->Center[2],
01963                            SOref->Center[0], SOref->Center[1], SOref->Center[2],
01964                            max_off);  
01965    }
01966       
01967    OffS = SUMA_Initialize_getoffsets (SOref->N_Node);
01968    #ifdef FROM_THIS_NODE
01969    istrt = FROM_THIS_NODE;
01970    istp = TO_THIS_NODE+1;
01971    #else
01972    istrt = 0;
01973    istp = SOref->N_Node;
01974    #endif
01975    for (i =istrt ; i<istp; ++i) {
01976       if (i == 0) {
01977          SUMA_etime(&start_time,0);
01978       }
01979       SUMA_getoffsets2 (i, SOref, max_off, OffS, NULL, 0);
01980       /* find average distance between nodes within offset and center of surface */
01981       a = &(SOref->NodeList[3*i]); SUMA_SEG_LENGTH(a, SOref->Center, ave_dist_ref); 
01982       cnt = 1;
01983       #ifdef FROM_THIS_NODE
01984             fprintf(SUMA_STDERR, "%s: Considering the following %d neighbors to:\n"
01985                                  "i=%d; [%f, %f, %f]; d[%d] = %f\n", FuncName, OffS->N_Nodes,
01986                                     i, SOref->NodeList[3*i], SOref->NodeList[3*i+1], SOref->NodeList[3*i+2], 
01987                                     cnt, ave_dist_ref);
01988       #endif
01989       for (j=0; j<OffS->N_Nodes; ++j)
01990       {
01991          
01992          if (i!=j && OffS->LayerVect[j] >= 0 && OffS->OffVect[j] <= max_off)
01993          {
01994             
01995             a = &(SOref->NodeList[3*j]); SUMA_SEG_LENGTH(a, SOref->Center, dj);
01996             ave_dist_ref += dj;
01997             ++cnt;
01998             #ifdef FROM_THIS_NODE
01999                fprintf(SUMA_STDERR, ""
02000                                  "j=%d; [%f, %f, %f]; d[%d] = %f\n",
02001                                     j, SOref->NodeList[3*j], SOref->NodeList[3*j+1], SOref->NodeList[3*j+2], 
02002                                     cnt, dj);
02003             #endif
02004          }
02005       }
02006       ave_dist_ref /=  (float)cnt;
02007       /* move node i to the reference average location
02008       Do not travel along normals, you should travel along
02009       radial direction Center-->node*/
02010       a = &(SO->NodeList[3*i]); SUMA_UNIT_VEC(SO->Center, a, U, Un);
02011       if (Un) {
02012          SUMA_POINT_AT_DISTANCE_NORM(U, SO->Center, ave_dist_ref, P2);
02013          SO->NodeList[3*i] = P2[0][0]; SO->NodeList[3*i+1] = P2[0][1]; SO->NodeList[3*i+2] = P2[0][2];
02014       } else {
02015             SUMA_SL_Err("Identical points!\n"
02016                         "No coordinates modified");
02017       }
02018       
02019       if (LocalHead) {
02020          if (! (i%999)) {
02021             a = &(SO->NodeList[3*i]);
02022             SUMA_SEG_LENGTH(a, SOref->Center, dj);
02023             fprintf(SUMA_STDERR, "%s:\n"
02024                            "node i=%d, avg_dist_ref = %f\ncnt = %d\n"
02025                            "Check on P2: New dist =%f ?=? %f\n", 
02026                            FuncName, i, ave_dist_ref, cnt, dj, ave_dist_ref);
02027             etime_GetOffset = SUMA_etime(&start_time,1);
02028             fprintf(SUMA_STDERR, "%s: Search to %f mm took %f seconds for %d nodes.\n"
02029                   "Projected completion time: %f minutes\n",
02030                   FuncName, max_off, etime_GetOffset, i+1,
02031                   etime_GetOffset * SO->N_Node / 60.0 / (i+1));
02032          }
02033       }
02034       if (! (i%99) && cs) {
02035          if (cs->Send) { /* send the first monster (it's SOref  "in SUMA" that's being modified on the fly*/
02036             if (!SUMA_SendToSuma (SOref, cs, (void *)SO->NodeList, SUMA_NODE_XYZ, 1)) {
02037             SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted.");
02038             }
02039          }
02040       }
02041       
02042       #ifdef FROM_THIS_NODE
02043       {
02044          FILE *fid=NULL;
02045          char *outname=NULL, tmp[20];
02046          int ii;
02047          if (cs->Send) { /* send the first monster (it's SOref that's being modified on the fly*/
02048             if (!SUMA_SendToSuma (SOref, cs, (void *)SO->NodeList, SUMA_NODE_XYZ, 1)) {
02049             SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted.");
02050             }
02051          }
02052          sprintf(tmp,"offset_n%d", FROM_THIS_NODE);
02053          outname = SUMA_Extension("", ".1D", YUP);
02054          outname = SUMA_append_replace_string(outname, "offset.1D", "", 1);
02055          fid = fopen(outname, "w"); free(outname); outname = NULL;
02056          if (!fid) {
02057             SUMA_SL_Err("Could not open file for writing.\nCheck file permissions, disk space.\n");
02058          } else {
02059             fprintf (fid,"#Column 1 = Node index\n"
02060                          "#column 2 = Neighborhood layer\n"
02061                          "#Column 3 = Distance from node %d\n", 99);
02062             for (ii=0; ii<SO->N_Node; ++ii) {
02063                if (OffS->LayerVect[ii] >= 0) {
02064                   fprintf(fid,"%d\t%d\t%f\n", ii, OffS->LayerVect[ii], OffS->OffVect[ii]);
02065                }
02066             }
02067             fclose(fid);
02068          }
02069          { int jnk; fprintf(SUMA_STDOUT,"Pausing, next node is %d...", i+1); jnk = getchar(); fprintf(SUMA_STDOUT,"\n"); }
02070       }
02071       #endif
02072              
02073       /* recycle offsets structure */
02074       SUMA_Recycle_getoffsets (OffS);
02075       
02076    }   
02077    
02078    /* offsets are to be freed */
02079    SUMA_Free_getoffsets(OffS); OffS = NULL;
02080    
02081    SUMA_RETURN(YUP);
02082 }

01757 {
01758    static char FuncName[]={"SUMA_EquateSurfaceVolumes"};
01759    int iter, i, iter_max, ndiv;
01760    double a, b, d;
01761    SUMA_VolDiffDataStruct fdata;
01762    SUMA_Boolean LocalHead = NOPE;
01763    
01764    SUMA_ENTRY;
01765 
01766    if (!SO || !SOref) { SUMA_SL_Err("NULL surfaces"); SUMA_RETURN(NOPE); }
01767    if (SO->N_Node != SOref->N_Node || SO->N_FaceSet != SOref->N_FaceSet) { SUMA_SL_Err("Surfaces not isotopic"); SUMA_RETURN(NOPE); }
01768    
01769    if (LocalHead) {
01770       fprintf(SUMA_STDERR, "%s:\n"
01771                            " SO    Center: %f, %f, %f\n"
01772                            " SOref Center: %f, %f, %f\n"
01773                            , FuncName, 
01774                            SO->Center[0], SO->Center[1], SO->Center[2],
01775                            SOref->Center[0], SOref->Center[1], SOref->Center[2]);  
01776    }
01777      
01778    /* fill up fdata */
01779    fdata.SO = SO; fdata.SOref = SOref; fdata.cs = cs;
01780    fdata.tmpList = (float *)SUMA_malloc(SOref->NodeDim * SOref->N_Node * sizeof(float));
01781    if (!fdata.tmpList) {
01782       SUMA_SL_Err("Failed to allocate");
01783       SUMA_RETURN(0);
01784    }
01785 
01786    if (!SUMA_GetVolDiffRange(&fdata, &a, &b)) {
01787       SUMA_SL_Err("Failed to get range");
01788       SUMA_RETURN(NOPE);
01789    }
01790    
01791    if (LocalHead) {
01792       fprintf(SUMA_STDERR,"%s:\na = %f\tb=%f\n", FuncName, a, b);
01793    }      
01794    SUMA_BinaryZeroSearch(a, b, SUMA_VolDiff, &fdata, 500, tol);  
01795    
01796    /* now make the new node list be SO's thingy*/
01797    SUMA_free(SO->NodeList); SO->NodeList = fdata.tmpList; fdata.tmpList = NULL;
01798        
01799    SUMA_RETURN(YUP);
01800 }

00482 {
00483    static char FuncName[]={"SUMA_FillToVoxelMask"};
00484    byte *isin = NULL, *visited=NULL;
00485    DList*candlist=NULL;
00486    DListElmt *dothiselm=NULL;
00487    int dothisvoxel, itmp;
00488    int nl[50], N_n, in ,neighb, nijk, i, j, k, nij;
00489    SUMA_Boolean LocalHead = NOPE;
00490    
00491    SUMA_ENTRY;
00492    
00493    *N_in = 0;
00494    
00495    if (!ijkmask) {
00496       SUMA_SL_Err("Nothing to do");
00497       SUMA_RETURN(NULL);
00498    }
00499    if (ijkmask[ijkseed]) {
00500       SUMA_SL_Err("Seed is on mask. Bad business.");
00501       SUMA_RETURN(NULL);
00502    }
00503    
00504    nij = ni * nj;
00505    nijk = ni * nj * nk;
00506    
00507    if (LocalHead) {
00508       SUMA_1D_2_3D_index (ijkseed, i, j, k, ni, nij);
00509       fprintf(SUMA_STDERR,"%s:\nSeed is %d %d %d\n", FuncName, i, j, k); 
00510    }
00511    candlist = (DList*)SUMA_malloc(sizeof(DList));
00512    visited = (byte *)SUMA_calloc(nijk, sizeof(byte));
00513    if (!visited || !candlist) {
00514       SUMA_SL_Crit("Failed to allocate for visited or candlist");  
00515       SUMA_RETURN(NULL);  
00516    }
00517    
00518    if (usethisisin) isin = usethisisin;
00519    else {
00520       isin = (byte *)SUMA_calloc(nijk, sizeof(byte));
00521       if (!isin) {
00522          SUMA_SL_Crit("Failed to allocate");
00523          SUMA_RETURN(NULL);
00524       }
00525    }
00526    
00527    dothisvoxel = ijkseed;
00528    dlist_init(candlist, NULL);
00529    
00530    
00531    isin[dothisvoxel] = 1; ++(*N_in); /* Add voxel to cluster */
00532    visited[dothisvoxel] = 1;  
00533    dlist_ins_next(candlist, dlist_tail(candlist), (void *)dothisvoxel); /* Add voxel as next candidate*/
00534    
00535    while (dlist_size(candlist)) {
00536       /* find neighbors in its vicinity */
00537       dothiselm = dlist_head(candlist); dothisvoxel = (int) dothiselm->data;
00538       N_n = SUMA_VoxelNeighbors (dothisvoxel, ni, nj, nk, SUMA_VOX_NEIGHB_FACE, nl);
00539       /* remove node from candidate list */
00540       dlist_remove(candlist, dothiselm, (void*)&itmp);
00541       /* search to see if any are to be assigned */
00542       for (in=0; in<N_n; ++in) { 
00543          neighb = nl[in];
00544          if (!ijkmask[neighb]) {
00545             isin[neighb] = 1; ++(*N_in); /* Add voxel to cluster */
00546             /* mark it as a candidate if it has not been visited as a candidate before */
00547             if (!visited[neighb]) {
00548                dlist_ins_next(candlist, dlist_tail(candlist), (void *)neighb);
00549                visited[neighb] = 1;   
00550             }
00551          }
00552       }
00553    }
00554    
00555    if (visited) SUMA_free(visited); visited = NULL;
00556    if (candlist) { dlist_destroy(candlist); SUMA_free(candlist); candlist  = NULL; }
00557 
00558    
00559    SUMA_RETURN(isin);
00560 }

06788 {
06789    static char FuncName[]={"SUMA_Find_Edge_Nhost"};
06790    
06791    SUMA_ENTRY;
06792 
06793    for (*i=0; *i < N_IsInter; ++(*i)) {
06794       if (EL->ELps[IsInter[*i]][2] == Nhost) SUMA_RETURN (IsInter[*i]);
06795    }
06796    
06797    SUMA_RETURN (-1);
06798 
06799 }

02700 {
02701    static char FuncName[]={"SUMA_FormNeighbOffset"};
02702    int i, ii, il, jl, noffs, ShowNode = 4;
02703    SUMA_GET_OFFSET_STRUCT *OffS = NULL;
02704    struct  timeval start_time;
02705    float etime_GetOffset, mean_N_Neighb, dist, dist_norm;   
02706    SUMA_OFFSET_STRUCT *OffS_out=NULL;
02707    SUMA_Boolean LocalHead = NOPE;
02708    
02709    SUMA_ENTRY;
02710    if (!SO) { SUMA_SL_Err("NULL SO"); SUMA_RETURN(NULL); }
02711    if (!SO->FN) {
02712       SUMA_SL_Err("NULL SO->FN");
02713       SUMA_RETURN(NULL);
02714    }
02715    OffS_out = (SUMA_OFFSET_STRUCT *)SUMA_malloc(SO->N_Node * sizeof(SUMA_OFFSET_STRUCT));
02716    
02717    SUMA_etime(&start_time,0);
02718    
02719    OffS = SUMA_Initialize_getoffsets (SO->N_Node);
02720    mean_N_Neighb = 0;
02721    dist_norm = 1.1 * OffsetLim;
02722    for (i=0; i < SO->N_Node; ++i) {
02723       SUMA_getoffsets2 (i, SO, OffsetLim, OffS, NULL, 0);
02724       /* Now store all the relevant info in OffS in OffS_out[i] */
02725       OffS_out[i].N_Neighb = 0; 
02726       for (il=1; il<OffS->N_layers; ++il) {
02727          OffS_out[i].N_Neighb += OffS->layers[il].N_NodesInLayer;
02728       }
02729       OffS_out[i].Neighb_ind = (int *)SUMA_malloc(OffS_out[i].N_Neighb * sizeof(int));
02730       OffS_out[i].Neighb_dist = (float *)SUMA_malloc(OffS_out[i].N_Neighb * sizeof(float));
02731       mean_N_Neighb += OffS_out[i].N_Neighb;
02732       noffs = 0;
02733       for (il=1; il<OffS->N_layers; ++il) {
02734          for (jl=0; jl<OffS->layers[il].N_NodesInLayer; ++jl) {
02735             OffS_out[i].Neighb_ind[noffs] = OffS->layers[il].NodesInLayer[jl]; 
02736             #if 1
02737             /* don't play fancy with the weights here */
02738             OffS_out[i].Neighb_dist[noffs] = OffS->OffVect[OffS_out[i].Neighb_ind[noffs]];
02739             #else
02740             dist = OffS->OffVect[OffS_out[i].Neighb_ind[noffs]];
02741             if (dist > OffsetLim) OffS_out[i].Neighb_dist[noffs] = 0;
02742             else OffS_out[i].Neighb_dist[noffs] = (dist ); 
02743             #endif
02744             ++noffs;
02745          }
02746       }
02747       
02748       /* Show me the offsets for one node*/
02749       if (0) {
02750          if (i == ShowNode) {
02751             FILE *fid=NULL;
02752             char *outname=NULL;
02753             outname = SUMA_Extension("SomethingOffset", ".1D", YUP);
02754             outname = SUMA_append_replace_string(outname, "offset.1D", "", 1);
02755             fid = fopen(outname, "w"); free(outname); outname = NULL;
02756             if (!fid) {
02757                SUMA_SL_Err("Could not open file for writing.\nCheck file permissions, disk space.\n");
02758             } else {
02759                fprintf (fid,"#Column 1 = Node index\n"
02760                             "#column 2 = Neighborhood layer\n"
02761                             "#Column 3 = Distance from node %d\n", ShowNode);
02762                for (ii=0; ii<SO->N_Node; ++ii) {
02763                   if (OffS->LayerVect[ii] >= 0) {
02764                      fprintf(fid,"%d\t%d\t%f\n", ii, OffS->LayerVect[ii], OffS->OffVect[ii]);
02765                   }
02766                }
02767                fclose(fid);
02768             }
02769          }
02770       }
02771                
02772       if (i == 99) {
02773          etime_GetOffset = SUMA_etime(&start_time,1);
02774          fprintf(SUMA_STDERR, "%s: Search to %f mm took %f seconds for %d nodes.\n"
02775                               "Projected completion time: %f minutes\n"
02776                               "Projected memory need for structure %f MB", 
02777                               FuncName, OffsetLim, etime_GetOffset, i+1, 
02778                               etime_GetOffset * SO->N_Node / 60.0 / (i+1),
02779                               (mean_N_Neighb / (i+1) * 8 + 12)* SO->N_Node/1000000.0);
02780       }
02781       SUMA_Recycle_getoffsets (OffS);
02782    }
02783    SUMA_Free_getoffsets(OffS); OffS = NULL;
02784    
02785    etime_GetOffset = SUMA_etime(&start_time,1);
02786    fprintf(SUMA_STDERR, "%s: Search to %f mm took %f seconds for %d nodes.\n"
02787                         "Mean number of neighbors per node: %f\n",
02788                         FuncName, OffsetLim, etime_GetOffset, SO->N_Node, mean_N_Neighb / SO->N_Node);
02789    
02790    SUMA_RETURN(OffS_out);           
02791 } 

00859 {
00860    static char FuncName[]={"SUMA_Free_getoffsets"};
00861    int i = 0;
00862    static SUMA_Boolean LocalHead = NOPE;
00863    
00864    SUMA_ENTRY;
00865    
00866    if (!OffS) SUMA_RETURN(NULL);
00867    
00868    if (OffS->layers) {
00869       for (i=0; i< OffS->N_layers; ++i) if (OffS->layers[i].NodesInLayer) SUMA_free(OffS->layers[i].NodesInLayer);
00870       SUMA_free(OffS->layers);
00871    }
00872    
00873    if (OffS->OffVect) SUMA_free(OffS->OffVect);
00874    if (OffS->LayerVect) SUMA_free(OffS->LayerVect);
00875    SUMA_free(OffS); OffS = NULL;
00876    
00877    SUMA_RETURN(NULL);
00878 }

02799 {
02800    static char FuncName[]={"SUMA_free_NeighbOffset"};
02801    int i;
02802    SUMA_ENTRY;
02803 
02804    if (!SO) {
02805       SUMA_S_Err("NULL SO!");
02806       SUMA_RETURN(NULL);
02807    }
02808    if (!OffS_out) SUMA_RETURN(NULL);
02809    for (i=0; i < SO->N_Node; ++i) {
02810       OffS_out[i].N_Neighb = 0;
02811       if (OffS_out[i].Neighb_dist) SUMA_free(OffS_out[i].Neighb_dist); OffS_out[i].Neighb_dist = NULL;
02812       if (OffS_out[i].Neighb_ind) SUMA_free(OffS_out[i].Neighb_ind); OffS_out[i].Neighb_ind = NULL;
02813    }
02814    SUMA_free(OffS_out); 
02815    SUMA_RETURN(NULL);
02816 }

01142 {
01143    static char FuncName[]={"SUMA_Free_Offset_ll_Datum"};
01144    SUMA_OFFSET_LL_DATUM *dt;
01145    
01146    SUMA_ENTRY;
01147    
01148    if (data) {
01149       dt = (SUMA_OFFSET_LL_DATUM *)data; 
01150       SUMA_free(dt);
01151    }
01152    
01153    SUMA_RETURNe;
01154 }   

06621 {
06622    static char FuncName[]={"SUMA_free_SPI"};
06623    
06624    SUMA_ENTRY;
06625    
06626    if (!SPI) SUMA_RETURNe;
06627    if (SPI->IntersTri) SUMA_free(SPI->IntersTri);
06628    if (SPI->IntersNodes) SUMA_free(SPI->IntersNodes);
06629    if (SPI->IntersEdges) SUMA_free(SPI->IntersEdges);
06630    if (SPI->isNodeInMesh) SUMA_free(SPI->isNodeInMesh); 
06631    if (SPI->isTriHit) SUMA_free (SPI->isTriHit);
06632    if (SPI->isEdgeInters) SUMA_free (SPI->isEdgeInters);
06633      
06634    if (SPI) SUMA_free(SPI);
06635    
06636    SUMA_RETURNe;
06637 }

06559 {
06560 
06561    static char FuncName[]={"SUMA_free_STB"};
06562    int i;
06563    
06564    SUMA_ENTRY;
06565    
06566    for (i=0; i < N_Bv; ++i) {
06567       if (Bv[i].list) SUMA_free(Bv[i].list);
06568    }
06569    if (Bv) SUMA_free(Bv);
06570    
06571    SUMA_RETURNe;
06572     
06573 }

05299 {
05300    static char FuncName[]={"SUMA_freePatch"};
05301    
05302    SUMA_ENTRY;
05303    
05304    
05305    if (Patch->FaceSetIndex) SUMA_free(Patch->FaceSetIndex);
05306    if (Patch->FaceSetList) SUMA_free(Patch->FaceSetList);
05307    if (Patch->nHits) SUMA_free(Patch->nHits);
05308    if (Patch) SUMA_free(Patch);
05309    SUMA_RETURN(YUP);
05310    
05311 }

00166 {
00167    static char FuncName[]={"SUMA_FreeVTI"};
00168    int i;
00169    
00170    SUMA_ENTRY;
00171    
00172    if (!vti) SUMA_RETURN(NULL);
00173    if (vti->TriIndex) SUMA_free(vti->TriIndex);
00174    if (vti->IntersectedVoxels) {
00175       for (i=0; i<vti->N_TriIndex; ++i) {
00176          if (vti->IntersectedVoxels[i]) free(vti->IntersectedVoxels[i]);
00177       }
00178       SUMA_free(vti->IntersectedVoxels);
00179    }
00180    if (vti->N_IntersectedVoxels) SUMA_free(vti->N_IntersectedVoxels);    
00181    SUMA_free(vti);
00182      
00183    SUMA_RETURN(NULL);
00184 }     

07485 {  static char FuncName[]={"SUMA_FromIntEdgeToIntEdge"};
07486    int Tri2 = 0, cnt, Incident[5], N_Incident;
07487    float dx, dy, dz;
07488    SUMA_Boolean Found, LocalHead = NOPE;
07489    
07490    SUMA_ENTRY;
07491 
07492    if (Tri < 0 || E1 < 0 || E2 < 0) {
07493       fprintf (SUMA_STDERR, "Error %s: Tri (%d) or E1 (%d) or E2 (%d) is negative!\n", FuncName, Tri, E1, E2);
07494       SUMA_RETURN (NOPE);
07495    }
07496    
07497    
07498    dx = (SPI->IntersNodes[3*E2] - SPI->IntersNodes[3*E1]);
07499    dy = (SPI->IntersNodes[3*E2+1] - SPI->IntersNodes[3*E1+1]);
07500    dz = (SPI->IntersNodes[3*E2+2] - SPI->IntersNodes[3*E1+2]);
07501    if (LocalHead) {
07502       fprintf (SUMA_STDERR, "%s: Entered - Tri %d, E1 %d [%d %d], E2 %d [%d %d]\n\tdx = %f dy = %f dz = %f\n", 
07503          FuncName, Tri, E1, EL->EL[E1][0], EL->EL[E1][1], E2, EL->EL[E2][0], EL->EL[E2][1], dx, dy, dz);
07504    }
07505    *d += sqrt( dx * dx + dy * dy + dz * dz);
07506    
07507    if (*d > dmax) {
07508       /* path already longer than Dijkstra path, no need to search further in this direction, get out with this d value */
07509       fprintf (SUMA_STDERR, "%s: Path longer than dmax. Returning.\n", FuncName);
07510       SUMA_RETURN (YUP);
07511    }
07512    
07513    if (EL->EL[E2][0] == Ny || EL->EL[E2][1] == Ny) {
07514       fprintf (SUMA_STDERR, "%s: Found Ny, d = %f\n", FuncName, *d);
07515       if (!Visited[Tri]) {
07516          /* add triangle to path */
07517          tPath[*N_tPath] = Tri;
07518          ++*N_tPath;
07519       }
07520       SUMA_RETURN (YUP);
07521    } else if (Visited[Tri]) {
07522       fprintf (SUMA_STDERR, "Error %s: Triangle %d already visited.\n",FuncName, Tri); 
07523       SUMA_RETURN (NOPE);
07524    }
07525      
07526    /* mark triangle as visited */
07527    if (LocalHead) fprintf (SUMA_STDERR, "%s: Marking triangle %d and adding %dth element to tPath.\n", FuncName, Tri, *N_tPath);
07528    Visited[Tri] = YUP;
07529    
07530    /* add triangle to path */
07531    tPath[*N_tPath] = Tri;
07532    ++*N_tPath;
07533    
07534    /* now get the second intersected triangle, incident to E2 */
07535    if (LocalHead) fprintf (SUMA_STDERR, "%s: Searching for triangles incident to E2 %d.\n", FuncName, E2);
07536    if (!SUMA_Get_Incident(EL->EL[E2][0], EL->EL[E2][1], EL, Incident, &N_Incident, 1)) {
07537       fprintf (SUMA_STDERR,"Error %s: Failed to get Incident triangles.\n", FuncName);
07538       SUMA_RETURN (NOPE);
07539    }
07540    
07541    /* find Tri2 such that Tri2 != Tri and Tri2 is an intersected triangle */
07542    cnt = 0;
07543    Found = NOPE;
07544    while (cnt < N_Incident && !Found) {
07545       if (SPI->isTriHit[Incident[cnt]] && Incident[cnt] != Tri && !Visited[Incident[cnt]]) {
07546          Found = YUP;
07547          Tri2 = Incident[cnt];
07548       }
07549       ++cnt;
07550    }
07551    
07552    if (!Found) {
07553       fprintf (SUMA_STDERR,"Error %s: Could not find next triangle.\n", FuncName);
07554       SUMA_RETURN (NOPE);
07555    }
07556    
07557    Tri = Tri2;
07558    E1 = E2;
07559    
07560    /* now find the new E2 */
07561    if (LocalHead) fprintf (SUMA_STDERR, "%s: Finding new E2.\n", FuncName);
07562    
07563    if (!SUMA_isSameEdge (EL, EL->Tri_limb[Tri][0], E1) && SPI->isEdgeInters[EL->Tri_limb[Tri][0]]) {
07564       E2 = EL->Tri_limb[Tri][0];
07565    }else if (!SUMA_isSameEdge (EL, EL->Tri_limb[Tri][1], E1) && SPI->isEdgeInters[EL->Tri_limb[Tri][1]]) {
07566       E2 = EL->Tri_limb[Tri][1];
07567    }else if (!SUMA_isSameEdge (EL, EL->Tri_limb[Tri][2], E1) && SPI->isEdgeInters[EL->Tri_limb[Tri][2]]) {
07568       E2 = EL->Tri_limb[Tri][2];
07569    }else {
07570       fprintf (SUMA_STDERR,"Error %s: No E2 found.\n", FuncName);
07571       SUMA_RETURN (NOPE);
07572    }
07573    
07574    /* call the same function again */
07575    if (!SUMA_FromIntEdgeToIntEdge (Tri, E1, E2, EL, SPI, Ny, Visited, d, dmax, tPath, N_tPath)) {
07576       fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_FromIntEdgeToIntEdge.\n", FuncName);
07577       SUMA_RETURN (NOPE);
07578    }
07579    
07580    SUMA_RETURN (YUP);
07581 }

01593 {
01594    static char FuncName[]={"SUMA_GetAreaDiffRange"};
01595    double a = 0.0, b = 0.0, nat=0, nbt=0, An, Bn;
01596    SUMA_Boolean LocalHead = NOPE;
01597 
01598    SUMA_ENTRY;
01599 
01600    /* decide on segment range */
01601    fdata->Aref = fabs((double)SUMA_Mesh_Area(fdata->SOref, NULL, -1));
01602    SUMA_SO_RADIUS(fdata->SOref, fdata->Rref);
01603    fdata->A = fabs((double)SUMA_Mesh_Area(fdata->SO, NULL, -1));
01604    SUMA_SO_RADIUS(fdata->SO, fdata->R);
01605 
01606    /* a very simple range setting. might very well fail at times */
01607    if (fdata->Aref - fdata->A < 0) { 
01608       a = fdata->R; /* choose 'a' such that Aref - A is < 0, Choose b later*/
01609       An = fdata->A;
01610       b = fdata->Rref;
01611       do {
01612          SUMA_LH("Looking for b");
01613          b *= 1.3;  /* choose A that Aref - A is > 0 */
01614          Bn = SUMA_NewAreaAtRadius(fdata->SO, b, fdata->Rref, fdata->tmpList);
01615          ++nbt;
01616       } while ( fdata->Aref - Bn < 0 && nbt < 200);
01617    }else{
01618       b = fdata->R; /* choose 'b' such that Aref - A is > 0, Choose a later*/
01619       Bn = fdata->A;
01620       a = fdata->Rref;
01621       do {
01622          SUMA_LH("Looking for a");
01623          a *= 0.7;
01624          An = SUMA_NewAreaAtRadius(fdata->SO, a, fdata->Rref, fdata->tmpList);
01625          ++nat;
01626       } while ( fdata->Aref -  An> 0 && nat < 200);
01627    }
01628 
01629    *ap = a; *bp = b;
01630 
01631    if (nat >= 200 || nbt >= 200) {
01632       SUMA_SL_Err("Failed to find segment.");
01633       SUMA_RETURN(NOPE);
01634    }
01635 
01636    if (LocalHead) {
01637       fprintf (SUMA_STDERR,"%s:\nChosen range is [%f %f] with Areas [%f %f], reference Area is %f\n", 
01638             FuncName, a, b, An, Bn, fdata->Aref);
01639    }
01640    
01641    SUMA_RETURN(YUP); 
01642 }

05341 {
05342    static char FuncName[]={"SUMA_GetContour"};
05343    SUMA_EDGE_LIST * SEL=NULL;
05344    SUMA_PATCH *Patch = NULL;
05345    int i, Tri, Tri1, Tri2, sHits;
05346    SUMA_CONTOUR_EDGES *CE = NULL;
05347    SUMA_Boolean *isNode=NULL;
05348    SUMA_Boolean LocalHead = NOPE;
05349    
05350    SUMA_ENTRY;
05351 
05352    *N_ContEdges = -1;
05353    
05354    /* get the Node member structure if needed*/
05355    if (!SO->MF) {
05356       SUMA_SLP_Err("Member FaceSet not created.\n");
05357       SUMA_RETURN(CE);
05358    }  
05359 
05360    /* create a flag vector of which node are in Nodes */
05361    isNode = (SUMA_Boolean *) SUMA_calloc(SO->N_Node, sizeof(SUMA_Boolean));
05362    if (!isNode) {
05363       SUMA_SLP_Crit("Failed to allocate for isNode");
05364       SUMA_RETURN(CE);
05365    }
05366    
05367    for (i=0; i < N_Node; ++i) isNode[Nodes[i]] = YUP;
05368   
05369    if (UseThisPatch) {
05370       SUMA_LH("Using passed patch");
05371       Patch = UseThisPatch;
05372    } else { 
05373       SUMA_LH("Creating patch");
05374       switch (ContourMode) {
05375          case 0:
05376             Patch = SUMA_getPatch (Nodes, N_Node, SO->FaceSetList, SO->N_FaceSet, SO->MF, 2);
05377             break;
05378          case 1:
05379             Patch = SUMA_getPatch (Nodes, N_Node, SO->FaceSetList, SO->N_FaceSet, SO->MF, 1);
05380             break;
05381          default:
05382             SUMA_SL_Err("Bad contour mode"); SUMA_RETURN(NULL);
05383             break;
05384       }
05385    }
05386    if (LocalHead) SUMA_ShowPatch (Patch,NULL);
05387    
05388    if (Patch->N_FaceSet) {
05389       SEL = SUMA_Make_Edge_List_eng (Patch->FaceSetList, Patch->N_FaceSet, SO->N_Node, SO->NodeList, 0, NULL);
05390    
05391       if (0 && LocalHead) SUMA_Show_Edge_List (SEL, NULL);
05392       /* allocate for maximum */
05393       CE = (SUMA_CONTOUR_EDGES *) SUMA_malloc(SEL->N_EL * sizeof(SUMA_CONTOUR_EDGES));
05394       if (!CE) {
05395          SUMA_SLP_Crit("Failed to allocate for CE");
05396          SUMA_RETURN(CE);
05397       }
05398    
05399       switch (ContourMode) {
05400          case 0: /* a pretty contour, edges used here may not be the outermost of the patch */
05401             /* edges that are part of unfilled triangles are good */
05402             i = 0;
05403             *N_ContEdges = 0;
05404             while (i < SEL->N_EL) {
05405                if (SEL->ELps[i][2] == 2) {
05406                   Tri1 = SEL->ELps[i][1];
05407                   Tri2 = SEL->ELps[i+1][1];
05408                   sHits = Patch->nHits[Tri1] + Patch->nHits[Tri2];
05409                   if (sHits == 5 || sHits == 4) { /* one tri with 3 hits and one with 2 hits or 2 Tris with 2 hits each */
05410                      /* Pick edges that are part of only one triangle with three hits */
05411                                                  /* or two triangles with two hits */
05412                      /* There's one more condition, both nodes have to be a part of the original list */
05413                      if (isNode[SEL->EL[i][0]] && isNode[SEL->EL[i][1]]) {
05414                         CE[*N_ContEdges].n1 = SEL->EL[i][0];
05415                         CE[*N_ContEdges].n2 = SEL->EL[i][1];
05416                         ++ *N_ContEdges;
05417 
05418                         if (LocalHead) {
05419                            fprintf (SUMA_STDERR,"%s: Found edge made up of nodes [%d %d]\n",
05420                               FuncName, SEL->EL[i][0], SEL->EL[i][1]);
05421                         }
05422                      }
05423                   }
05424                }
05425 
05426                if (SEL->ELps[i][2] > 0) {
05427                   i += SEL->ELps[i][2];
05428                } else {
05429                   i ++;
05430                }
05431             }
05432             break;
05433          case 1: /* outermost contour, not pretty, but good for getting the outermost edge */
05434             i = 0;
05435             *N_ContEdges = 0;
05436             while (i < SEL->N_EL) {
05437                if (SEL->ELps[i][2] == 1) {
05438                   CE[*N_ContEdges].n1 = SEL->EL[i][0];
05439                   CE[*N_ContEdges].n2 = SEL->EL[i][1];
05440                   ++ *N_ContEdges;
05441                   if (LocalHead) {
05442                            fprintf (SUMA_STDERR,"%s: Found edge made up of nodes [%d %d]\n",
05443                               FuncName, SEL->EL[i][0], SEL->EL[i][1]);
05444                   }
05445                }
05446                if (SEL->ELps[i][2] > 0) {
05447                   i += SEL->ELps[i][2];
05448                } else {
05449                   i ++;
05450                }
05451             }
05452             break;
05453          default:
05454             SUMA_SL_Err("Bad ContourMode");
05455             SUMA_RETURN(NULL);
05456             break;            
05457       }
05458       
05459       /* Now reallocate */
05460       if (! *N_ContEdges) {
05461          SUMA_free(CE); CE = NULL;
05462          SUMA_RETURN(CE);
05463       }else {
05464          CE = (SUMA_CONTOUR_EDGES *) SUMA_realloc (CE, *N_ContEdges * sizeof(SUMA_CONTOUR_EDGES));
05465          if (!CE) {
05466             SUMA_SLP_Crit("Failed to reallocate for CE");
05467             SUMA_RETURN(CE);
05468          }
05469       }
05470          
05471       SUMA_free_Edge_List (SEL); 
05472    }
05473    
05474    if (!UseThisPatch) {
05475       SUMA_freePatch (Patch); 
05476    }
05477    Patch = NULL;
05478    
05479    SUMA_free(isNode);
05480    
05481    SUMA_RETURN(CE);
05482 }

02564 {
02565    static char FuncName[]={"SUMA_GetOffset2Offset"};
02566    int il, jl, noffs;
02567    SUMA_Boolean LocalHead = NOPE;
02568    
02569    SUMA_ENTRY;
02570    
02571    if (!GOS || !OS) {
02572       SUMA_SL_Err("NULL input"); SUMA_RETURN(NOPE);
02573    }
02574    
02575    OS->N_Neighb = 0; 
02576    for (il=1; il<GOS->N_layers; ++il) {
02577       OS->N_Neighb += GOS->layers[il].N_NodesInLayer;
02578    }
02579    OS->Neighb_ind = (int *)SUMA_malloc(OS->N_Neighb * sizeof(int));
02580    OS->Neighb_dist = (float *)SUMA_malloc(OS->N_Neighb * sizeof(float));
02581    if (!OS->Neighb_ind || !OS->Neighb_dist) {
02582       SUMA_SL_Crit("Failed to allocate.");
02583       SUMA_RETURN(NOPE);
02584    }
02585    
02586    noffs = 0;
02587    for (il=1; il<GOS->N_layers; ++il) {
02588       for (jl=0; jl<GOS->layers[il].N_NodesInLayer; ++jl) {
02589          OS->Neighb_ind[noffs] = GOS->layers[il].NodesInLayer[jl]; 
02590          OS->Neighb_dist[noffs] = GOS->OffVect[OS->Neighb_ind[noffs]];
02591          ++noffs;
02592       }
02593    }
02594    
02595    SUMA_RETURN(YUP);
02596 }

00683 {
00684    static char FuncName[]={"SUMA_getoffsets"};
00685    int i, ni, iseg;
00686    float Off_tmp;
00687    SUMA_Boolean Visit = NOPE;
00688    static SUMA_Boolean LocalHead = NOPE;
00689    
00690    SUMA_ENTRY;
00691    
00692    #if DoCheck
00693    if (!SO->FN || !SO->EL) {
00694       SUMA_SL_Err("SO->FN &/| SO->EL are NULL.\n");
00695       SUMA_RETURN(NOPE);
00696    }
00697    #endif
00698    
00699    #if DBG
00700    if (LocalHead) fprintf(SUMA_STDERR,"%s: Working node %d, %d neighbs. lim = %f\n", 
00701                                     FuncName, n, SO->FN->N_Neighb[n], lim);
00702    #endif
00703    
00704    for (i=0; i < SO->FN->N_Neighb[n]; ++i) {
00705       ni = SO->FN->FirstNeighb[n][i]; /* for notational sanity */
00706       iseg = SUMA_FindEdge (SO->EL, n, SO->FN->FirstNeighb[n][i]);
00707       #if DoCheck
00708       if (iseg < 0) {
00709          SUMA_SL_Err("Failed to find segment");
00710          SUMA_RETURN(NOPE);
00711       }
00712       #endif
00713       
00714       Off_tmp = Off[n] + SO->EL->Le[iseg];   /* that is the distance from n (original n) to ni along
00715                                                 that particular path */
00716                                              
00717       Visit = NOPE;
00718       if (Off[ni] < 0 || Off_tmp < Off[ni]) { /* Distance improvement, visit/revist that node */
00719          if (Off_tmp < lim) { /* only record if less than lim */
00720             Visit = YUP;
00721             Off[ni] = Off_tmp;
00722          } 
00723       } 
00724       
00725       #if DBG
00726       if (LocalHead) fprintf(SUMA_STDERR,"%s: %d --> %d. Visit %d, Current %f, Old %f\n", 
00727          FuncName, n, ni, Visit, Off_tmp, Off[ni]);
00728       #endif
00729       
00730       #if 0
00731          { int jnk; fprintf(SUMA_STDOUT,"Pausing ..."); jnk = getchar(); fprintf(SUMA_STDOUT,"\n"); }
00732       #endif
00733 
00734       if (Visit) { /* a new node has been reached with an offset less than limit, go down that road */
00735          if (!SUMA_getoffsets (ni, SO, Off, lim))  {
00736             SUMA_SL_Err("Failed in SUMA_getoffsets");
00737             SUMA_RETURN (NOPE);
00738          }
00739       }
00740    }
00741 
00742    SUMA_RETURN(YUP);
00743 }

01051 {
01052    static char FuncName[]={"SUMA_getoffsets2"};
01053    int LayInd, il, n_il, n_jne, k, n_prec = -1, n_k, jne, iseg=0;
01054    float Off_tmp, Seg, *a, *b, minSeg, SegPres; /*! *** SegPres added Jul 08 04, ZSS bug before ... */
01055    SUMA_Boolean Visit = NOPE;
01056    SUMA_Boolean AllDone = NOPE;
01057    static SUMA_Boolean LocalHead = NOPE;
01058    
01059    SUMA_ENTRY;
01060    
01061    if (!OffS) {
01062       SUMA_SL_Err("NULL OffS");
01063       SUMA_RETURN(NOPE);
01064    }
01065    
01066    /* setup 0th layer */
01067    OffS->OffVect[n] = 0.0;   /* n is at a distance 0.0 from itself */
01068    OffS->LayerVect[n] = 0;   /* n is on the zeroth layer */
01069    OffS->layers[0].N_NodesInLayer = 1;
01070    OffS->layers[0].NodesInLayer[0] = n;
01071    if (CoverThisNode) { 
01072       if (CoverThisNode[n]) {
01073          CoverThisNode[n] = 0; --N_CoverThisNode;
01074       }
01075    }
01076    LayInd = 1;  /* index of next layer to build */
01077    AllDone = NOPE;
01078    while (!AllDone) {
01079       
01080       AllDone = YUP; /* assume that this would be the last layer */
01081       for (il=0; il < OffS->layers[LayInd - 1].N_NodesInLayer; ++il) { /* go over all nodes in previous layer */
01082          n_il =  OffS->layers[LayInd - 1].NodesInLayer[il]; /* node from previous layer */
01083          for (jne=0; jne < SO->FN->N_Neighb[n_il]; ++jne) { /* go over all the neighbours of node n_il */
01084             n_jne = SO->FN->FirstNeighb[n_il][jne];        /* node that is an immediate neighbor to n_il */
01085             if (OffS->LayerVect[n_jne] < 0) { /* node is not assigned to a layer yet */
01086                OffS->LayerVect[n_jne] =  LayInd;    /* assign new layer index to node */
01087                OffS->OffVect[n_jne] = 0.0;          /* reset its distance from node n */
01088                SUMA_AddNodeToLayer (n_jne, LayInd, OffS);   /* add the node to the nodes in the layer */
01089                minSeg = 100000.0;
01090                n_prec = -1; 
01091                Seg = 0.0;
01092                SegPres = 0.0;
01093                for (k=0; k < SO->FN->N_Neighb[n_jne]; ++k) { /* calculate shortest distance of node to any precursor */  
01094                   n_k = SO->FN->FirstNeighb[n_jne][k];
01095                   if (OffS->LayerVect[n_k] == LayInd - 1) { /* this neighbor is a part of the previous layer, good */
01096                      if (n_prec < 0) n_prec = SO->FN->FirstNeighb[n_jne][0];
01097                      a = &(SO->NodeList[3*n_k]); b = &(SO->NodeList[3*n_jne]);
01098                      /* this is the slow part, too many redundant computations. 
01099                         Computation time is cut by a factor > 2 if Seg was set to a constant
01100                         However, attempts at accessing pre-calculated segment lengths
01101                         proved to be slower. See Comments in function help*/
01102                      SUMA_SEG_LENGTH_SQ (a, b, Seg);                    
01103                      if (OffS->OffVect[n_prec] + Seg < minSeg) {
01104                         minSeg = Seg + OffS->OffVect[n_prec];
01105                         SegPres = Seg;
01106                         n_prec = n_k;
01107                      }
01108                   }
01109                }/* for k */
01110                
01111                if (n_prec < 0) { /* bad news */
01112                   SUMA_SL_Crit("No precursor found for node.");
01113                   OffS = SUMA_Free_getoffsets (OffS);
01114                   SUMA_RETURN(NOPE);
01115                } else {
01116                   OffS->OffVect[n_jne] = OffS->OffVect[n_prec] + sqrt(SegPres); SegPres = 0.0;
01117                   if (!CoverThisNode) {
01118                      if (OffS->OffVect[n_jne] < lim) { /* must go at least one more layer */
01119                         AllDone = NOPE;
01120                      }
01121                   } else {
01122                      if (CoverThisNode[n_jne]) {
01123                         CoverThisNode[n_jne] = 0; --N_CoverThisNode;
01124                      }
01125                      if (N_CoverThisNode > 0) {
01126                         AllDone = NOPE;
01127                      }
01128                   }
01129                }
01130             } /* node not already in layer */
01131             
01132          } /* for jne */
01133       
01134       } /* for il */    
01135       ++LayInd;
01136    } /* while AllDone */
01137    
01138    SUMA_RETURN(YUP);
01139 }

01193 {
01194    static char FuncName[]={"SUMA_getoffsets_ll"};
01195    int LayInd, il, n_il, n_jne, k, n_prec = -1, n_k, jne, iseg=0;
01196    float Off_tmp, Seg, *a, *b, minSeg, SegPres; /*! *** SegPres added Jul 08 04, ZSS bug before ... */
01197    SUMA_Boolean Visit = NOPE;
01198    SUMA_Boolean AllDone = NOPE;
01199    SUMA_OFFSET_LL_DATUM * n_dat = NULL, *dat = NULL, *dat_nk = NULL, *dat_prec = NULL, *dat_ne=NULL;
01200    DList *list = NULL;
01201    DListElmt *elm = NULL, *elm_prec = NULL, *elm_ne=NULL, *elm_nk=NULL;
01202    static SUMA_Boolean LocalHead = NOPE;
01203    
01204    SUMA_ENTRY;
01205    
01206    
01207    /* create the list */
01208    SUMA_LH("Initializing list ...");
01209    list = (DList *)SUMA_malloc(sizeof(DList));
01210    dlist_init(list, SUMA_Free_Offset_ll_Datum);
01211    
01212    /* setup 0th layer */
01213    SUMA_LH("New OffsetDatum");
01214    n_dat = SUMA_New_Offset_ll_Datum(n, 0);
01215    n_dat->off = 0.0;   /* n is at a distance 0.0 from itself */
01216    dlist_ins_next(list, dlist_tail(list), (void*)n_dat);
01217    
01218    if (CoverThisNode) { 
01219       if (CoverThisNode[n]) {
01220          CoverThisNode[n] = 0; --N_CoverThisNode;
01221       }
01222    }
01223    LayInd = 1;  /* index of next layer to build */
01224    AllDone = NOPE;
01225    while (!AllDone) {
01226       AllDone = YUP; /* assume that this would be the last layer */
01227       elm = NULL;
01228          do {
01229             if (!elm) { SUMA_BEGINNING_OF_LAYER(list, (LayInd-1), elm); }
01230             else elm = elm->next;
01231             if (!elm) {
01232                SUMA_SL_Err("Could not find beginning of layer!");
01233                SUMA_RETURN(NULL);
01234             }
01235             dat = (SUMA_OFFSET_LL_DATUM *)elm->data;
01236             if (dat->layer == LayInd -1) {
01237                n_il = dat->ni;
01238                for (jne=0; jne < SO->FN->N_Neighb[n_il]; ++jne) { /* go over all the neighbours of node n_il */
01239                   n_jne = SO->FN->FirstNeighb[n_il][jne];        /* node that is an immediate neighbor to n_il */
01240                   SUMA_FIND_ELMENT_FOR_NODE(list, n_jne, elm_ne);
01241                   if (!elm_ne) { /* node not in any layer */
01242                      dat_ne = SUMA_New_Offset_ll_Datum(n_jne, LayInd); /* create an element for it */
01243                      dat_ne->off = 0.0;
01244                      dlist_ins_next(list, dlist_tail(list), (void*)dat_ne);
01245                      minSeg = 100000.0;
01246                      n_prec = -1; 
01247                      Seg = 0.0;
01248                      SegPres = 0.0;
01249                      for (k=0; k < SO->FN->N_Neighb[n_jne]; ++k) { /* calculate shortest distance of node to any precursor */  
01250                         n_k = SO->FN->FirstNeighb[n_jne][k];
01251                         SUMA_FIND_ELMENT_FOR_NODE(list, n_k, elm_nk); 
01252                         if (n_prec < 0 && elm_nk) { 
01253                            n_prec = n_k; elm_prec = elm_nk; 
01254                            dat_prec = (SUMA_OFFSET_LL_DATUM *)elm_prec->data;
01255                         }
01256                         if (elm_nk) {
01257                            dat_nk = (SUMA_OFFSET_LL_DATUM *)elm_nk->data;
01258                            if (dat_nk->layer == LayInd - 1) { /* this neighbor is a part of the previous layer, good */
01259                               a = &(SO->NodeList[3*n_k]); b = &(SO->NodeList[3*n_jne]);
01260                               /* this is the slow part, too many redundant computations. 
01261                                  Computation time is cut by a factor > 2 if Seg was set to a constant
01262                                  However, attempts at accessing pre-calculated segment lengths
01263                                  proved to be slower. See Comments in function help*/
01264                               SUMA_SEG_LENGTH_SQ (a, b, Seg);                    
01265                               if (dat_prec->off + Seg < minSeg) {
01266                                  minSeg = Seg + dat_prec->off;
01267                                  SegPres = Seg;
01268                                  n_prec = n_k;
01269                                  elm_prec = elm_nk;
01270                                  dat_prec = dat_nk;
01271                               }
01272                            }
01273                         } /* if elm_nk */
01274                      }/* for k */
01275                      if (n_prec < 0) { /* bad news */
01276                         SUMA_SL_Crit("No precursor found for node.");
01277                         SUMA_RETURN(NULL);
01278                      } else {
01279                         dat_ne->off = dat_prec->off + sqrt(SegPres); SegPres = 0.0;
01280                         if (!CoverThisNode) {
01281                            if (dat_ne->off < lim) { /* must go at least one more layer */
01282                               AllDone = NOPE;
01283                            }
01284                         } else {
01285                            if (CoverThisNode[n_jne]) {
01286                               CoverThisNode[n_jne] = 0; --N_CoverThisNode;
01287                            }
01288                            if (N_CoverThisNode > 0) {
01289                               AllDone = NOPE;
01290                            }
01291                         }
01292                      }
01293                   } /* if elm_ne */
01294                } /* for jne */
01295             } /* dat->layer == LayInd */
01296          }  while (dat->layer == (LayInd-1) && elm != dlist_tail(list));
01297       
01298       ++LayInd;
01299    } /* while AllDone */
01300    
01301    SUMA_RETURN(list);
01302 }

05226 {
05227    int * BeenSelected;
05228    int i, j, node, ip, ip2, NP;
05229    SUMA_PATCH *Patch;
05230    static char FuncName[]={"SUMA_getPatch"};
05231    
05232    SUMA_ENTRY;
05233    
05234    NP = 3;
05235    BeenSelected = (int *)SUMA_calloc (N_Full_FaceSetList, sizeof(int));
05236    Patch = (SUMA_PATCH *)SUMA_malloc(sizeof(SUMA_PATCH));
05237    
05238    if (!BeenSelected || !Patch) {
05239       fprintf (SUMA_STDERR,"Error %s: Could not allocate for BeenSelected or patch.\n", FuncName);
05240       SUMA_RETURN(NULL);
05241    }
05242    /* find out the total number of facesets these nodes are members of */
05243    Patch->N_FaceSet = 0; /* total number of facesets containing these nodes */
05244    for (i=0; i < N_Nodes; ++i) {
05245       node = NodesSelected[i];
05246       for (j=0; j < Memb->N_Memb[node]; ++j) {
05247          if (!BeenSelected[Memb->NodeMemberOfFaceSet[node][j]]) {
05248             /* this faceset has not been selected, select it */
05249             ++ Patch->N_FaceSet;
05250          }
05251          ++ BeenSelected[Memb->NodeMemberOfFaceSet[node][j]];
05252       }   
05253    }
05254    
05255    /* now load these facesets into a new matrix */
05256    
05257    Patch->FaceSetList = (int *) SUMA_calloc (Patch->N_FaceSet * 3, sizeof(int));
05258    Patch->FaceSetIndex = (int *) SUMA_calloc (Patch->N_FaceSet, sizeof(int));
05259    Patch->nHits = (int *) SUMA_calloc (Patch->N_FaceSet, sizeof(int));
05260    
05261    if (!Patch->FaceSetList || !Patch->FaceSetIndex || !Patch->nHits) {
05262       fprintf (SUMA_STDERR,"Error %s: Could not allocate for Patch->FaceSetList || Patch_FaceSetIndex.\n", FuncName);
05263       SUMA_RETURN(NULL);
05264    }
05265    
05266    j=0;
05267    for (i=0; i < N_Full_FaceSetList; ++i) {
05268       if (BeenSelected[i] >= MinHits) {
05269          Patch->nHits[j] = BeenSelected[i];
05270          Patch->FaceSetIndex[j] = i;
05271          ip = NP * j;
05272          ip2 = NP * i;
05273          Patch->FaceSetList[ip] = Full_FaceSetList[ip2];
05274          Patch->FaceSetList[ip+1] = Full_FaceSetList[ip2+1];
05275          Patch->FaceSetList[ip+2] = Full_FaceSetList[ip2+2];
05276          ++j;
05277       }
05278    }
05279    
05280    /* reset the numer of facesets because it might have changed given the MinHits condition,
05281    It won't change if MinHits = 1.
05282    It's OK not to change the allocated space as long as you are using 1D arrays*/
05283    Patch->N_FaceSet = j;
05284    
05285    if (BeenSelected) SUMA_free(BeenSelected);
05286    
05287    SUMA_RETURN(Patch);   
05288 }

01650 {
01651    static char FuncName[]={"SUMA_GetVolDiffRange"};
01652    double a = 0.0, b = 0.0, nat=0, nbt=0;
01653    SUMA_Boolean LocalHead = NOPE;
01654 
01655    SUMA_ENTRY;
01656 
01657    /* decide on segment range */
01658    fdata->Vref = fabs((double)SUMA_Mesh_Volume(fdata->SOref, NULL, -1));
01659    SUMA_SO_RADIUS(fdata->SOref, fdata->Rref);
01660    fdata->V = fabs((double)SUMA_Mesh_Volume(fdata->SO, NULL, -1));
01661    SUMA_SO_RADIUS(fdata->SO, fdata->R);
01662 
01663    /* a very simple range setting. might very well fail at times */
01664    if (fdata->Vref - fdata->V < 0) { 
01665       a = fdata->R; /* choose 'a' such that Vref - V is < 0, Choose b later*/
01666       b = fdata->Rref;
01667       do {
01668          SUMA_LH("Looking for b");
01669          b *= 1.3; ++nbt; /* choose b that Vref - V is > 0 */
01670       } while ( fdata->Vref - SUMA_NewVolumeAtRadius(fdata->SO, b, fdata->Rref, fdata->tmpList) < 0 && nbt < 200);
01671    }else{
01672       b = fdata->R; /* choose 'b' such that Vref - V is > 0, Choose a later*/
01673       a = fdata->Rref;
01674       do {
01675          SUMA_LH("Looking for a");
01676          a *= 0.7; ++nat;
01677       } while ( fdata->Vref - SUMA_NewVolumeAtRadius(fdata->SO, a, fdata->Rref, fdata->tmpList) > 0 && nat < 200);
01678    }
01679 
01680    *ap = a; *bp = b;
01681 
01682    if (nat >= 200 || nbt >= 200) {
01683       SUMA_SL_Err("Failed to find segment.");
01684       SUMA_RETURN(NOPE);
01685    }
01686 
01687    SUMA_RETURN(YUP); 
01688 }

00201 {
00202    static char FuncName[]={"SUMA_GetVoxelsIntersectingTriangle"};
00203    int ti, nx, ny, nz, nxy, nxyz, N_inbox, n1, n2, n3, nt, nt3, nijk, nf;
00204    int N_alloc, N_realloc, en, *voxelsijk=NULL, N_voxels1d = 0, *voxels1d = NULL;
00205    int *TriIndex=NULL, N_TriIndex;
00206    float dxyz[3];
00207    float tol_dist = 0; /* A way to fatten up the shell a little bit. Set to 0 if no fat is needed */
00208    float *p1, *p2, *p3, min_v[3], max_v[3], p[3], dist;
00209    FILE *fp=NULL;
00210    SUMA_Boolean LocalHead = YUP;
00211    
00212    SUMA_ENTRY;
00213       
00214    if (SO->FaceSetDim != 3 || SO->NodeDim != 3) {
00215       SUMA_SL_Err("SO->FaceSetDim != 3 || SO->NodeDim != 3"); 
00216       SUMA_RETURN(NULL);
00217    }
00218    if (!vti) {
00219       SUMA_SL_Err("vti must be non NULL");
00220       SUMA_RETURN(NULL);
00221    }
00222    if (vti->N_TriIndex <= 0) {
00223       SUMA_SL_Err("vti must be initialized");
00224       SUMA_RETURN(NULL);
00225    }
00226    
00227    TriIndex = vti->TriIndex;
00228    N_TriIndex = vti->N_TriIndex;
00229    
00230    nx = VolPar->nx; ny = VolPar->ny; nz = VolPar->nz; nxy = nx * ny; nxyz = nx * ny * nz;
00231    
00232    if (LocalHead) {
00233       fp = fopen("SUMA_GetVoxelsIntersectingTriangle.1D","w");
00234       if (fp) fprintf(fp, "# Voxels from %s that intersect the triangles \n", VolPar->filecode);
00235    }   
00236    /* cycle through all triangles and find voxels that intersect them */
00237    N_alloc = 2000; /* expected maximum number of voxels in triangle's bounding box */
00238    N_realloc = 0;
00239    voxelsijk = (int *)SUMA_malloc(sizeof(int)*N_alloc*3);
00240    if (!voxelsijk) { SUMA_SL_Crit("Failed to Allocate!"); SUMA_RETURN(NULL);  }   
00241    dxyz[0] = VolPar->dx; dxyz[1] = VolPar->dy; dxyz[2] = VolPar->dz;
00242    for (ti=0; ti<N_TriIndex; ++ti) {
00243       if (LocalHead) fprintf(SUMA_STDERR,"%s: Now processing %dth triangle\n", FuncName, ti);
00244       nf = TriIndex[ti];
00245       if (LocalHead) fprintf(SUMA_STDERR,"\t\tindexed %d, \n", nf);
00246       n1 = SO->FaceSetList[SO->FaceSetDim*nf]; n2 = SO->FaceSetList[SO->FaceSetDim*nf+1]; n3 = SO->FaceSetList[SO->FaceSetDim*nf+2];
00247       if (LocalHead) fprintf(SUMA_STDERR,"\t\tmade up of  nodes %d, %d. %d . \n", n1, n2, n3);
00248       /* find the bounding box of the triangle */
00249       p1 = &(NodeIJKlist[3*n1]); p2 = &(NodeIJKlist[3*n2]); p3 = &(NodeIJKlist[3*n3]); 
00250       SUMA_TRIANGLE_BOUNDING_BOX(p1, p2, p3, min_v, max_v);
00251       
00252       /* quick check of preallocate size of voxelsijk */
00253       en =((int)(max_v[0] - min_v[0] + 5) * (int)(max_v[1] - min_v[1] + 5) * (int)(max_v[2] - min_v[2] + 5)); 
00254       if ( en > N_alloc) {
00255          ++N_realloc; if (N_realloc > 5) { SUMA_SL_Warn("Reallocating, increase limit to improve speed.\nEither triangles too large or grid too small"); }
00256          N_alloc = 2*en;
00257          voxelsijk = (int *)SUMA_realloc(voxelsijk, 3*N_alloc*sizeof(int));
00258          if (!voxelsijk) { SUMA_SL_Crit("Failed to Allocate!"); SUMA_RETURN(NULL); }
00259       } 
00260       /* find the list of voxels inhabiting this box */
00261       N_inbox = 0;
00262       if (!SUMA_VoxelsInBox(voxelsijk, &N_inbox, min_v, max_v)) {
00263          SUMA_SL_Err("Unexpected error!"); SUMA_RETURN(NULL); 
00264       }
00265       if (!N_inbox) { SUMA_SL_Err("Unexpected error, no voxels in box!"); SUMA_RETURN(NULL);  }
00266       if (N_inbox >= N_alloc) { SUMA_SL_Err("Allocation trouble!"); SUMA_RETURN(NULL);  }
00267       if (LocalHead) fprintf(SUMA_STDERR,"\t\t%d nodes in box\n", N_inbox);
00268       
00269       /* allocate for 1D indices of voxels intersecting the triangle */
00270       if (voxels1d) {
00271          SUMA_SL_Err("NULL pointer expected here");
00272          SUMA_RETURN(NULL);
00273       }
00274       if (LocalHead) fprintf(SUMA_STDERR,"\t\tShit man, %d nodes in box\n", N_inbox);
00275       voxels1d = (int *)malloc(N_inbox * sizeof(int)); /* Too many SUMA_mallocs keep it simple here, this function is called many many times */
00276       if (LocalHead) fprintf(SUMA_STDERR,"\t\tWTF man, %d nodes in box\n", N_inbox);
00277       N_voxels1d=0;
00278       if (!voxels1d) {
00279          SUMA_SL_Crit("Failed to allocate voxels1d");
00280          SUMA_RETURN(NULL);
00281       }
00282       /* mark these voxels as inside the business */
00283       if (LocalHead) fprintf(SUMA_STDERR,"%s:\t\tabout to process %d voxels\n", FuncName, N_inbox);
00284       for (nt=0; nt < N_inbox; ++nt) {
00285          nt3 = 3*nt;
00286          if (voxelsijk[nt3] < nx &&  voxelsijk[nt3+1] < ny &&  voxelsijk[nt3+2] < nz) {
00287             nijk = SUMA_3D_2_1D_index(voxelsijk[nt3], voxelsijk[nt3+1], voxelsijk[nt3+2], nx , nxy);  
00288             { 
00289                /* what side of the plane is this voxel on ? */
00290                p[0] = (float)voxelsijk[nt3]; p[1] = (float)voxelsijk[nt3+1]; p[2] = (float)voxelsijk[nt3+2]; 
00291                SUMA_DIST_FROM_PLANE(p1, p2, p3, p, dist);
00292                /* Does voxel contain any of the nodes ? */
00293                if (tol_dist && SUMA_ABS(dist) < tol_dist) dist = tol_dist; /* Fatten the representation a little bit 
00294                                                              There are holes in the mask created using
00295                                                              the condition below alone. I am not sure
00296                                                              why that is the case but whatever gap there
00297                                                              is in one plane, results from a thick line in
00298                                                              the other. Don't know if that is a bug yet
00299                                                              or an effect of discretization. At any rate
00300                                                              it should not affect what I plan to do with
00301                                                              this. Could the bug be in 
00302                                                              SUMA_isVoxelIntersect_Triangle?*/
00303                if (!(SUMA_IS_STRICT_NEG(VolPar->Hand * dist))) { /* voxel is outside (along normal) */
00304                   /* does this triangle actually intersect this voxel ?*/
00305                   if (SUMA_isVoxelIntersect_Triangle (p, dxyz, p1, p2, p3)) {
00306                      /* looks good, store it */
00307                      if (LocalHead) fprintf(SUMA_STDERR,"nt %d, N_voxels1d %d\n", nt, N_voxels1d);
00308                      voxels1d[N_voxels1d] = nijk; ++N_voxels1d;
00309                      if (fp) fprintf(fp, "%d %d %d\n", voxelsijk[nt3], voxelsijk[nt3+1], voxelsijk[nt3+2]);
00310                   } 
00311                }
00312                
00313             }
00314          }
00315       }
00316       /* store the results */
00317       vti->IntersectedVoxels[ti] = voxels1d; 
00318       vti->N_IntersectedVoxels[ti] = N_voxels1d;
00319       voxels1d = NULL; N_voxels1d = 0; 
00320    }
00321    
00322    if (LocalHead) {
00323       if (fp) fclose(fp); fp = NULL;
00324    } 
00325    SUMA_RETURN(vti);
00326 }

00759 {
00760    static char FuncName[]={"SUMA_Initialize_getoffsets"};
00761    int i;
00762    SUMA_GET_OFFSET_STRUCT *OffS = NULL;
00763    
00764    SUMA_ENTRY;
00765    
00766    if (N_Node <= 0) {
00767       SUMA_SL_Err("Bad values for N_Node");
00768       SUMA_RETURN (OffS);
00769    }
00770    
00771    OffS = (SUMA_GET_OFFSET_STRUCT *)SUMA_malloc(sizeof(SUMA_GET_OFFSET_STRUCT));
00772    if (!OffS) {
00773       SUMA_SL_Err("Failed to allocate for OffS");
00774       SUMA_RETURN (OffS);
00775    }
00776    
00777    OffS->OffVect = (float *) SUMA_malloc(N_Node * sizeof(float));
00778    OffS->LayerVect = (int *) SUMA_malloc(N_Node * sizeof(int));
00779    OffS->N_Nodes = N_Node;
00780    
00781    if (!OffS->LayerVect || !OffS->OffVect) {
00782       SUMA_SL_Err("Failed to allocate for OffS->LayerVect &/| OffS->OffVect");
00783       SUMA_free(OffS);
00784       SUMA_RETURN (OffS);
00785    }
00786    
00787    /* initialize vectors */
00788    for (i=0; i< N_Node; ++i) {
00789       OffS->OffVect[i] = 0.0;
00790       OffS->LayerVect[i] = -1;
00791    }
00792    
00793    /* add a zeroth layer for node n */
00794    OffS->N_layers = 1;
00795    OffS->layers = (SUMA_NODE_NEIGHB_LAYER *) SUMA_malloc(OffS->N_layers * sizeof(SUMA_NODE_NEIGHB_LAYER));
00796    OffS->layers[0].N_AllocNodesInLayer = 1;
00797    OffS->layers[0].NodesInLayer = (int *) SUMA_malloc(OffS->layers[0].N_AllocNodesInLayer * sizeof(int));
00798    OffS->layers[0].N_NodesInLayer = 0;   
00799    
00800    SUMA_RETURN (OffS);
00801    
00802 }

07286 {
07287    static char FuncName[]={"SUMA_IntersectionStrip"};
07288    int *tPath1 = NULL, *tPath2 = NULL, Incident[50], N_Incident, Nx = -1, 
07289       Ny = -1, Tri = -1, Tri1 = -1, istart, n2 = -1, n3 = -1, E1, E2, cnt, N_tPath1, N_tPath2;
07290    float d1, d2;
07291    SUMA_Boolean *Visited = NULL, Found, LocalHead = NOPE;
07292    
07293    SUMA_ENTRY;
07294    
07295    /* find the edge containing the 1st 2 nodes of the Dijkstra path */
07296    /* find the triangle that contains the edge formed by the 1st 2 nodes of the Dijkstra path and is intersected by the plane*/
07297    Tri1 = -1;
07298    if (LocalHead) {
07299       fprintf (SUMA_STDERR, "%s: Looking for a triangle containing nodes [%d %d].\n", FuncName, nPath[0], nPath[1]);
07300    }
07301    
07302    Found = SUMA_Get_Incident(nPath[0], nPath[1], SO->EL, Incident, &N_Incident, 1);
07303    if (!Found) {
07304       /* no such triangle, get a triangle that contains nPath[0] and is intersected */
07305       fprintf (SUMA_STDERR, "%s: No triangle contains nodes [%d %d].\n", FuncName, nPath[0], nPath[1]);
07306       if (nPath[0] == SO->N_Node - 1) {
07307          fprintf (SUMA_STDERR, "Warning %s: 1st node is last node of surface, traversing path backwards.\n", FuncName);
07308          Nx = nPath[N_nPath - 1];
07309          Ny = nPath[0];
07310       }else {
07311          Nx = nPath[0];
07312          Ny = nPath[N_nPath - 1];
07313       }
07314       istart = SO->EL->ELloc[Nx];
07315       /* find an edge containing the first node and belonging to an intersected triangle */
07316       Found = NOPE;
07317       while (SO->EL->EL[istart][0] == Nx && !Found) {
07318          Tri = SO->EL->ELps[istart][1];
07319          if (SPI->isTriHit[Tri]) {
07320             Found = YUP;
07321             Tri1 = Tri;
07322          }
07323          ++istart;
07324       } 
07325    }else {
07326       Nx = nPath[0];
07327       Ny = nPath[N_nPath - 1];
07328       
07329       /* find which of these triangles was intersected */
07330       if (LocalHead) {
07331          fprintf (SUMA_STDERR, "%s: Found %d triangles containing nodes [%d %d].\n", FuncName, N_Incident, nPath[0], nPath[1]);
07332          for (cnt = 0; cnt < N_Incident; ++cnt) fprintf (SUMA_STDERR, "%d isHit %d\n", Incident[cnt], SPI->isTriHit[Incident[cnt]]);
07333          fprintf (SUMA_STDERR, "\n"); 
07334       }
07335       Found = NOPE;
07336       cnt = 0;
07337       while (cnt < N_Incident && !Found) {
07338          if (SPI->isTriHit[Incident[cnt]]) {
07339             Found = YUP;
07340             Tri1 = Incident[cnt];
07341          }
07342          ++cnt;
07343       }
07344    }
07345    
07346    if (!Found) {
07347       fprintf (SUMA_STDERR, "Error %s: Starting Edge could not be found.\n", FuncName);
07348       SUMA_RETURN (NULL);
07349    }else if (LocalHead) {
07350       fprintf (SUMA_STDERR, "%s: Starting with triangle %d.\n", FuncName, Tri1);
07351    }
07352 
07353    /* found starting triangle edge, begin with side 1 */
07354    if (SO->FaceSetList[3*Tri1] == Nx) {
07355       n2 = SO->FaceSetList[3*Tri1+1];
07356       n3 = SO->FaceSetList[3*Tri1+2];
07357    } else if (SO->FaceSetList[3*Tri1+1] == Nx) {
07358       n2 = SO->FaceSetList[3*Tri1];
07359       n3 = SO->FaceSetList[3*Tri1+2];
07360    } else if (SO->FaceSetList[3*Tri1+2] == Nx) {
07361       n2 = SO->FaceSetList[3*Tri1];
07362       n3 = SO->FaceSetList[3*Tri1+1];
07363    } else {
07364       fprintf (SUMA_STDERR, "Error %s: Triangle %d does not contain Nx %d.\n", FuncName, Tri1, Nx);
07365       SUMA_RETURN (NULL);
07366    }  
07367    
07368    
07369    
07370    E1 = SUMA_FindEdgeInTri (SO->EL, Nx, n2, Tri1);
07371    if (!SPI->isEdgeInters[E1]) {
07372       E1 = SUMA_FindEdgeInTri (SO->EL, Nx, n3, Tri1);
07373    }
07374    /* now choose E2 such that E2 is also intersected */
07375    if (!SUMA_isSameEdge (SO->EL, SO->EL->Tri_limb[Tri1][0], E1) && SPI->isEdgeInters[SO->EL->Tri_limb[Tri1][0]]) {
07376       E2 = SO->EL->Tri_limb[Tri1][0];
07377    }else if (!SUMA_isSameEdge (SO->EL, SO->EL->Tri_limb[Tri1][1], E1) && SPI->isEdgeInters[SO->EL->Tri_limb[Tri1][1]]) {
07378       E2 = SO->EL->Tri_limb[Tri1][1];
07379    }else if (!SUMA_isSameEdge (SO->EL, SO->EL->Tri_limb[Tri1][2], E1) && SPI->isEdgeInters[SO->EL->Tri_limb[Tri1][2]]) {
07380       E2 = SO->EL->Tri_limb[Tri1][2];
07381    }else {
07382       fprintf (SUMA_STDERR,"Error %s: No E2 found.\n", FuncName);
07383       SUMA_RETURN (NULL);
07384    }
07385 
07386    Visited = (SUMA_Boolean *) SUMA_calloc (SO->N_FaceSet, sizeof(SUMA_Boolean));
07387    tPath1 = (int *) SUMA_calloc (SO->N_FaceSet, sizeof(int));
07388    if (!Visited || !tPath1) {
07389       fprintf (SUMA_STDERR, "Error %s: Failed to allocate.\n", FuncName);
07390       if (Visited) SUMA_free(Visited);
07391       if (tPath2) SUMA_free(tPath1); 
07392       SUMA_RETURN (NULL);
07393    }
07394    
07395    N_tPath1 = 0;
07396    if (!SUMA_FromIntEdgeToIntEdge (Tri1, E1, E2, SO->EL, SPI, Ny, Visited, &d1, dmax, tPath1, &N_tPath1)) {
07397       fprintf (SUMA_STDERR, "Error %s: Failed in SUMA_FromIntEdgeToIntEdge.\n", FuncName);
07398       if (Visited) SUMA_free(Visited);
07399       if (tPath2) SUMA_free(tPath1);      
07400       SUMA_RETURN (NULL);
07401    }
07402    
07403    if (LocalHead) {
07404       fprintf (SUMA_STDERR, "%s: Found a distance of %f.\n\n\n", FuncName, d1);
07405    }
07406    
07407    /* Now try going in the other direction, E2->E1 */
07408    cnt = E2;
07409    E2 = E1;
07410    E1 = cnt;
07411    
07412    /* reset the values of Visited */
07413    for (cnt=0; cnt < SO->N_FaceSet; ++cnt) if (Visited[cnt]) Visited[cnt] = NOPE;
07414    
07415    tPath2 = (int *) SUMA_calloc (SO->N_FaceSet, sizeof(int));
07416    if (!Visited || !tPath2) {
07417       fprintf (SUMA_STDERR, "Error %s: Failed to allocate.\n", FuncName);
07418       if (Visited) SUMA_free(Visited);
07419       if (tPath1) SUMA_free(tPath1);
07420       if (tPath2) SUMA_free(tPath2);
07421       SUMA_RETURN (NULL);
07422    }
07423 
07424    N_tPath2 = 0;
07425    if (!SUMA_FromIntEdgeToIntEdge (Tri1, E1, E2, SO->EL, SPI, Ny, Visited, &d2, dmax, tPath2, &N_tPath2)) {
07426       fprintf (SUMA_STDERR, "Error %s: Failed in SUMA_FromIntEdgeToIntEdge.\n", FuncName);
07427       if (Visited) SUMA_free(Visited);
07428       if (tPath1) SUMA_free(tPath1);
07429       if (tPath2) SUMA_free(tPath2);
07430       SUMA_RETURN (NULL);
07431    }
07432    
07433    if (Visited) SUMA_free(Visited);
07434    
07435    if (LocalHead) {
07436       fprintf (SUMA_STDERR, "%s: Found a distance of %f.\n", FuncName, d2);
07437    }
07438    
07439    if (d2 < d1) {
07440       *N_tPath = N_tPath2;
07441       *dinters = d2;
07442       if (tPath1) SUMA_free(tPath1);
07443       SUMA_RETURN (tPath2);
07444    } else {
07445       *dinters = d1;
07446       *N_tPath = N_tPath1;
07447       if (tPath2) SUMA_free(tPath2);
07448       SUMA_RETURN (tPath1);
07449    }
07450    
07451 }

07243 {
07244    static char FuncName[]={"SUMA_isSameEdge"};
07245    
07246    SUMA_ENTRY;
07247 
07248    if (EL->EL[E1][0] == EL->EL[E2][0] && EL->EL[E1][1] == EL->EL[E2][1]) {
07249       SUMA_RETURN (YUP);
07250    } else {
07251       SUMA_RETURN (NOPE);
07252    }
07253    
07254 }