00001 #ifndef SUMA_BRAINWRAP_INCLUDED
00002 #define SUMA_BRAINWRAP_INCLUDED
00003 
00004 typedef enum { SUMA_3dSS_NO_PAUSE = 0, SUMA_3dSS_DEMO_PAUSE, SUMA_3dSS_INTERACTIVE } SUMA_3DSS_MODES;
00005 
00006 float SUMA_LoadPrepInVol (SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, SUMA_SurfaceObject **SOhull);
00007 int SUMA_Find_IminImax (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, int ni, 
00008                         float *MinMax, float *MinMax_dist , float *MinMax_over, float *MinMax_over_dist,
00009                         float *Means, float *undershish, float *overshish, int *dvecind_under, int *dvecind_over, int ShishMax);
00010 int SUMA_SkullMask (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, SUMA_COMM_STRUCT *cs);
00011 int SUMA_StretchToFitLeCerveau (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, SUMA_COMM_STRUCT *cs);
00012 byte *SUMA_FindVoxelsInSurface_SLOW (SUMA_SurfaceObject *SO, SUMA_VOLPAR *VolPar, int *N_inp, int fillhole) ;
00013 short *SUMA_SurfGridIntersect (SUMA_SurfaceObject *SO, float *NodeIJKlist, SUMA_VOLPAR *VolPar, int *N_inp, int fillhole, THD_3dim_dataset *fillholeset);
00014 short *SUMA_FindVoxelsInSurface (SUMA_SurfaceObject *SO, SUMA_VOLPAR *VolPar, int *N_inpnt, int  fillhole, THD_3dim_dataset *fillholeset) ;
00015 int SUMA_Reposition_Touchup(SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, float limtouch, SUMA_COMM_STRUCT *cs) ;
00016 float *SUMA_Suggest_Touchup(SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, float limtouch, SUMA_COMM_STRUCT *cs, int *N_touch);
00017 float *SUMA_Suggest_Touchup_Grad(SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, float limtouch, SUMA_COMM_STRUCT *cs, int *N_touch);
00018 int SUMA_DidUserQuit(void);
00019 EDIT_options *SUMA_BlankAfniEditOptions(void);
00020 /*!
00021    SUMA_WRAP_BRAIN_SMOOTH(niter, bufp1, bufp2);
00022    \brief a chunk used in two places in SUMA_BrainWrap.
00023    Requires two float pointers than can be null on the first call
00024    but must be freed at the end 
00025 */
00026 #define SUMA_WRAP_BRAIN_SMOOTH_NN(niter, dsmooth, refNodeList){ \
00027    SUMA_SurfaceObject m_SOref;   \
00028    int m_in;   \
00029    float *m_a, m_P2[2][3], m_U[3], m_Un, m_Rref, m_R, m_Dr, m_Dn;  \
00030    if (!refNodeList) {  \
00031       refNodeList = (float *)SUMA_malloc(sizeof(float)*SO->N_Node*3);   \
00032       if (!refNodeList) { SUMA_SL_Crit("Failed to allocate for refNodeList!"); SUMA_RETURN(NOPE); }   \
00033    }  \
00034    SUMA_SO_RADIUS(SO, m_Rref);   /* get the radius before shrinking */\
00035    memcpy((void*)refNodeList, (void *)SO->NodeList, SO->N_Node * 3 * sizeof(float)); /* copy original surface coords */ \
00036    dsmooth = SUMA_NN_GeomSmooth( SO, niter, SO->NodeList, 3, SUMA_ROW_MAJOR, dsmooth, cs);    \
00037    memcpy((void*)SO->NodeList, (void *)dsmooth, SO->N_Node * 3 * sizeof(float)); /* copy smoothed surface coords */  \
00038    SUMA_RECOMPUTE_NORMALS(SO);   \
00039    /* matching area fails for some reason. */\
00040    if (0) { \
00041       /* just grow surface by 3 mm outward, fails for meshes of different size and different degrees of smoothing*/  \
00042       for (m_in=0; m_in < SO->N_Node; ++m_in) { \
00043          m_a = &(SO->NodeList[3*m_in]); \
00044          SUMA_UNIT_VEC(SO->Center, m_a, m_U, m_Un);\
00045          SUMA_POINT_AT_DISTANCE_NORM(m_U, SO->Center, (m_Un+2), m_P2);   \
00046          SO->NodeList[3*m_in] = m_P2[0][0]; SO->NodeList[3*m_in+1] = m_P2[0][1]; SO->NodeList[3*m_in+2] = m_P2[0][2];   \
00047       }  \
00048    }  else {   /* an approximate resizing based on radius */\
00049       SUMA_SO_RADIUS(SO, m_R);   /* get the radius after shrinking */\
00050       m_Dr = ( m_Rref - m_R ) / m_Rref;   \
00051       for (m_in=0; m_in < SO->N_Node; ++m_in) { \
00052          m_a = &(SO->NodeList[3*m_in]); \
00053          SUMA_UNIT_VEC(SO->Center, m_a, m_U, m_Un);\
00054          m_Dn = m_Dr*m_Un + m_Un;   \
00055          if (m_Un) { \
00056             SUMA_POINT_AT_DISTANCE_NORM(m_U, SO->Center, m_Dn, m_P2);  \
00057             SO->NodeList[3*m_in] = m_P2[0][0]; SO->NodeList[3*m_in+1] = m_P2[0][1]; SO->NodeList[3*m_in+2] = m_P2[0][2];   \
00058          }  \
00059       }  \
00060    }  \
00061 }
00062 
00063 /* Area matching is failing for some reason (negative area patches ?) */
00064 #define SUMA_WRAP_BRAIN_SMOOTH_MATCHAREA(niter, dsmooth, refNodeList){ \
00065    SUMA_SurfaceObject m_SOref;   \
00066    if (!refNodeList) {  \
00067       refNodeList = (float *)SUMA_malloc(sizeof(float)*SO->N_Node*3);   \
00068       if (!refNodeList) { SUMA_SL_Crit("Failed to allocate for refNodeList!"); SUMA_RETURN(NOPE); }   \
00069    }  \
00070    memcpy((void*)refNodeList, (void *)SO->NodeList, SO->N_Node * 3 * sizeof(float)); /* copy original surface coords */ \
00071    dsmooth = SUMA_NN_GeomSmooth( SO, niter, SO->NodeList, 3, SUMA_ROW_MAJOR, dsmooth, cs);    \
00072    memcpy((void*)SO->NodeList, (void *)dsmooth, SO->N_Node * 3 * sizeof(float)); /* copy smoothed surface coords */  \
00073    m_SOref.PolyArea = NULL; m_SOref.NodeDim = 3; m_SOref.FaceSetDim = 3;   \
00074    m_SOref.NodeList = refNodeList; m_SOref.N_Node = SO->N_Node;   \
00075    m_SOref.FaceSetList = SO->FaceSetList; m_SOref.N_FaceSet = SO->N_FaceSet;  \
00076    SUMA_EquateSurfaceAreas(SO, &m_SOref, 0.1, cs);  \
00077 }
00078 
00079 /*!
00080    this one is too wimpy to get rid of folding 
00081 */
00082 #define SUMA_WRAP_BRAIN_SMOOTH_TAUB(niter, dsmooth, refNodeList){ \
00083       dsmooth = SUMA_Taubin_Smooth( SO, NULL,   \
00084                                     0.6307, -.6732, SO->NodeList, \
00085                                     niter, 3, SUMA_ROW_MAJOR, dsmooth, cs);   \
00086       memcpy((void*)SO->NodeList, (void *)dsmooth, SO->N_Node * 3 * sizeof(float)); \
00087 }
00088 
00089 
00090 
00091 /*!
00092    Stop using this monster, use SUMA_Reposition_Touchup
00093 */
00094 #define SUMA_REPOSITION_TOUCHUP(limtouch){  \
00095       int m_in, m_N_troub = 0, m_cond1=0, m_cond2=0, m_cond3 = 0;   \
00096       float m_MinMax_over[2], m_MinMax[2], m_MinMax_dist[2], m_MinMax_over_dist[2], m_Means[3], m_tb; \
00097       float m_U[3], m_Un, *m_a, m_P2[2][3], *m_norm, m_shft; \
00098       float *m_touchup=NULL, **m_wgt=NULL, *m_dsmooth=NULL; \
00099       m_touchup = (float *)SUMA_calloc(SO->N_Node, sizeof(float));   \
00100       if (!m_touchup) { SUMA_SL_Crit("Failed to allocate"); exit(1); }  \
00101       for (m_in=0; m_in<SO->N_Node; ++m_in) {   \
00102          SUMA_Find_IminImax(SO, Opt, m_in,  m_MinMax, m_MinMax_dist, m_MinMax_over, m_MinMax_over_dist, m_Means, NULL, NULL, NULL, NULL, 0); \
00103          /* Shift the node outwards:   \
00104             0- the minimum location of the minimum over the node is not 0
00105             AND
00106             1- if the minimum over the node is closer than the minimum below the node  or the internal minimum is more than 1mm away \
00107                and the internal min is masked to 0 or the minimum below the node just sucks, > tb and the minimum over the node is < tb\
00108                The second part of the condition is meant to better model a finger of cortex that has enough csf on the inner surface \
00109                to have it be set to 0 by SpatNorm this dark csf would keep the brain surface from reaching the outer surface \
00110             AND   \
00111             2- if the maximum over the node is much higher than the maximum below the node (fat on skull), \
00112                it must be beyond the location of the minimum over the node   \
00113             AND   \
00114             3- The new position of the node is not much larger than the current value
00115          */ \
00116          m_tb = (m_MinMax[1] - Opt->t2) * 0.5 + Opt->t2; \
00117          m_cond1 = (    (m_MinMax_over_dist[0] < m_MinMax_dist[0]) || (m_MinMax_dist[0] > 1 && m_MinMax[0] >= m_MinMax_over[0] ) \
00118                      || (m_MinMax[0] > m_tb && m_MinMax_over[0] < m_MinMax[0]) ); \
00119          if (m_MinMax_over[1] > m_MinMax[1] && m_MinMax_over[1] > 0.9 * Opt->t98 && (m_MinMax_over_dist[1] < m_MinMax_over_dist[0]) ) m_cond2 = 0;  \
00120          else m_cond2 = 1; \
00121          if (m_MinMax_over[0] > 1.2 * m_Means[0]) m_cond3 = 0;\
00122          else m_cond3 = 1; \
00123          if (Opt->NodeDbg == m_in) {   \
00124                m_a = &(SO->NodeList[3*m_in]);   \
00125                fprintf(SUMA_STDERR, "%s: Debug during touchup for node %d:\n"   \
00126                                     " MinMax     =[%.1f,%.1f] MinMax_dist     = [%.2f,%.2f] \n"   \
00127                                     " MinMax_over=[%.1f,%.1f] MinMax_over_dist= [%.2f,%.2f] \n"   \
00128                                     " Val at node %.1f, mean below %.1f, mean above %.1f\n"  \
00129                                     " zalt= %f, r/2 = %f\n"    \
00130                                     " Conditions: %f %d %d %d\n",    \
00131                        FuncName, m_in, \
00132                        m_MinMax[0], m_MinMax[1], m_MinMax_dist[0], m_MinMax_dist[1], \
00133                        m_MinMax_over[0], m_MinMax_over[1], m_MinMax_over_dist[0], m_MinMax_over_dist[1], \
00134                        m_Means[0], m_Means[1], m_Means[2],   \
00135                        m_a[2] - SO->Center[2], Opt->r/2, \
00136                        m_MinMax_over_dist[0] , m_cond1 , m_cond2, m_cond3); \
00137          }  \
00138          if (m_MinMax_over_dist[0] && m_cond1 && m_cond2 && m_cond3) {   \
00139                m_a = &(SO->NodeList[3*m_in]);   \
00140                /* reposition nodes */  \
00141                if (0 && LocalHead) { fprintf(SUMA_STDERR, "%s: Suggest repair for node %d (better min above):\n"   \
00142                                     " MinMax     =[%.1f,%.1f] MinMax_dist     = [%.2f,%.2f] \n"   \
00143                                     " MinMax_over=[%.1f,%.1f] MinMax_over_dist= [%.2f,%.2f] \n"   \
00144                                     " Val at node %.1f, mean below %.1f, mean above %.1f\n"  \
00145                                     " zalt= %f, r/2 = %f\n",    \
00146                        FuncName, m_in, \
00147                        m_MinMax[0], m_MinMax[1], m_MinMax_dist[0], m_MinMax_dist[1], \
00148                        m_MinMax_over[0], m_MinMax_over[1], m_MinMax_over_dist[0], m_MinMax_over_dist[1], \
00149                        m_Means[0], m_Means[1], m_Means[2],   \
00150                        m_a[2] - SO->Center[2], Opt->r/2); } \
00151                {  /* Used to insist on value and top half if (m_MinMax_over_dist[0] && (m_a[2] - SO->Center[2] > 0))  */\
00152                                              /* to avoid going into eye balls */\
00153                                              /* Then added !SUMA_IS_EYE_ZONE(m_a,SO->Center) to avoid eyes */  \
00154                                              /* but that is no longer necessary with fancier expansion conditions. */\
00155                   m_touchup[m_in] = m_MinMax_over_dist[0]; \
00156                }  \
00157             ++m_N_troub;   \
00158          }  \
00159          /* nodes in the front of the brain that are sitting in the midst of a very bright area \
00160          should be moved further down until they hit the brain vs non brain area*/ \
00161       }  \
00162       if (LocalHead) fprintf (SUMA_STDERR,"%s: ********************* %d troubled nodes found\n", FuncName, m_N_troub); \
00163       if (1){ \
00164          /* fix the big bumps: Bad for brain surface that has lots of bump on top, like 201, 
00165          that kind of smoothing is better performed on an inner model of the skull .... Someday ....*/\
00166          if (LocalHead) fprintf (SUMA_STDERR,"%s: ********************* Now filtering...\n", FuncName); \
00167          m_wgt = SUMA_Chung_Smooth_Weights(SO);   \
00168          if (!m_wgt) { \
00169             SUMA_SL_Err("Failed to compute weights.\n"); \
00170             exit(1); \
00171          }  \
00172          m_dsmooth = SUMA_Chung_Smooth (SO, m_wgt, 200, 20, m_touchup, 1, SUMA_COLUMN_MAJOR, NULL, ps->cs);   \
00173          if (LocalHead) fprintf (SUMA_STDERR,"%s: ********************* filtering done.\n", FuncName);   \
00174       } else { \
00175          m_dsmooth = m_touchup; m_touchup = NULL;  \
00176       }  \
00177       /* add the changes */   \
00178       for (m_in=0; m_in<SO->N_Node; ++m_in) {   \
00179          m_a = &(SO->NodeList[3*m_in]);   \
00180          if (1 || !SUMA_IS_EYE_ZONE(m_a,SO->Center)) { \
00181             if (m_a[2] - SO->Center[2] > 10 )  m_shft = m_touchup[m_in]; /* no smooth baby for top where bumpy sulci can occur */ \
00182             else m_shft = m_dsmooth[m_in];   \
00183             if (m_shft) { \
00184                m_a = &(SO->NodeList[3*m_in]);   \
00185                m_norm = &(SO->NodeNormList[3*m_in]);  \
00186                SUMA_POINT_AT_DISTANCE(m_norm, m_a, SUMA_MIN_PAIR(m_shft, limtouch), m_P2);   \
00187                SO->NodeList[3*m_in] = m_P2[0][0]; SO->NodeList[3*m_in+1] = m_P2[0][1]; SO->NodeList[3*m_in+2] = m_P2[0][2];   \
00188                if (0 && LocalHead) fprintf(SUMA_STDERR,"%s: Acting on node %d because it is in top half, boosting by %f\n", \
00189                         FuncName, m_in, SUMA_MIN_PAIR(m_shft, limtouch));   \
00190             }\
00191          }  \
00192       }  \
00193       if (m_touchup) SUMA_free(m_touchup); m_touchup = NULL;   \
00194       if (m_dsmooth) SUMA_free(m_dsmooth); m_dsmooth = NULL;   \
00195       if (m_wgt) SUMA_free2D ((char **)m_wgt, SO->N_Node); m_wgt = NULL;   \
00196 }
00197 
00198 #define SUMA_IS_EYE_ZONE(n,c) ( ( ( (n)[1] - (c)[1] ) < -10 && ( (n)[2] - (c)[2] ) < 0 ) ? 1 : 0 ) 
00199 
00200 
00201 /*!
00202    finds the mean coordinate of the neighboring nodes (not including node i)
00203 */
00204 #define SUMA_MEAN_NEIGHB_COORD(SO, i, mnc) { \
00205    int m_j, m_in; \
00206    mnc[0] = 0; mnc[1] = 0; mnc[2] = 0; \
00207    for (m_j = 0; m_j < SO->FN->N_Neighb[i]; ++m_j) {  \
00208       m_in = 3*SO->FN->FirstNeighb[i][m_j]; \
00209       mnc[0] += SO->NodeList[m_in]; \
00210       mnc[1] += SO->NodeList[m_in+1]; \
00211       mnc[2] += SO->NodeList[m_in+2]; \
00212    }  \
00213    mnc[0] = mnc[0] / (SO->FN->N_Neighb[i]); mnc[1] = mnc[1] / (SO->FN->N_Neighb[i]);  mnc[2] = mnc[2] / (SO->FN->N_Neighb[i]); \
00214 }
00215 /*!
00216    find the mean segment length 
00217 */
00218 #define SUMA_MEAN_SEGMENT_LENGTH(SO, l) { \
00219    int m_cnt=0, m_j;   \
00220    float *m_p1, *m_p2, m_d=0.0, m_v[3];  \
00221    l = 0.0; \
00222    for (m_j=0; m_j < SO->EL->N_EL; ++m_j) {  \
00223       if (SO->EL->ELps[m_j][2] >= 0) { \
00224          m_p1 = &(SO->NodeList[3*SO->EL->EL[m_j][0]]); m_p2 = &(SO->NodeList[3*SO->EL->EL[m_j][1]]); \
00225          SUMA_UNIT_VEC(m_p1, m_p2, m_v, m_d);   \
00226          l += m_d;   \
00227          ++m_cnt; \
00228       }  \
00229    }  \
00230    if (m_cnt) l /= m_cnt; \
00231 }
00232 
00233 /*!
00234    find the earliest max step, quit looking beyond a threshold
00235 */
00236 #define SUMA_MAX_SHISH_JUMP(vec, diffmax, i_diffmax, thresh_clip){  \
00237    int m_kk = 1, m_thresh = 0;   \
00238    float m_diff;  \
00239    m_diff = 0; diffmax = 0;  i_diffmax = 0;\
00240    while (vec[m_kk] >=0 && m_kk < ShishMax && !m_thresh) {   \
00241       m_diff = vec[m_kk] - vec[m_kk - 1]; \
00242       if (m_diff > diffmax) {diffmax = m_diff; i_diffmax = m_kk-1; }  \
00243       if (diffmax > thresh_clip) m_thresh = 1;   \
00244       ++ m_kk; \
00245    }  \
00246 }
00247 
00248 /*!
00249    find the earliest max negative step
00250 */
00251 #define SUMA_MAX_NEG_SHISH_JUMP(vec, diffmax, i_diffmax){  \
00252    int m_kk = 1, m_thresh = 0;   \
00253    float m_diff;  \
00254    m_diff = 0; diffmax = 0;  i_diffmax = 0;\
00255    while (vec[m_kk] >=0 && m_kk < ShishMax && !m_thresh) {   \
00256       m_diff = vec[m_kk] - vec[m_kk - 1]; \
00257       if (m_diff < diffmax) {\
00258          diffmax = m_diff; i_diffmax = m_kk-1; \
00259       }  \
00260       ++ m_kk; \
00261    }  \
00262 }
00263 /*!
00264    find the number of negative jumps above a certain threshold 
00265 */
00266 #define SUMA_MIN_SHISH_JUMP(vec, diffmin, i_diffmin, thresh_clip, N_neg){  \
00267    int m_kk = 1;   \
00268    float m_diff;  \
00269    m_diff = 0; diffmin = 0;  i_diffmin = 0; N_neg=0;\
00270    while (vec[m_kk] >=0 && m_kk < ShishMax ) {   \
00271       m_diff = vec[m_kk] - vec[m_kk - 1]; \
00272       if (m_diff < diffmin) { diffmin = m_diff; i_diffmax = m_kk-1; }  \
00273       if (SUMA_ABS(diffmin) > thresh_clip) ++N_neg;   \
00274       ++ m_kk; \
00275    }  \
00276 }
00277 
00278 /*! \brief fill a shish d units lond starting at xyz and moving along direction U
00279    
00280       \param XYZ coordinates of starting point
00281       \param U the famed direction vector
00282       \param in_vol AFNI volume structure
00283       \param dvec data vector 
00284       \param stp size of the step to take 
00285       \param n_stp number of steps to take
00286       \param nxx, nxy number of voxels in the x and x*y directions
00287       \param shish a vector to store dvec's values as it crosses them 
00288       \param N_shishmax (int) maximum number of values allowed in shish
00289 */
00290 
00291 #define SUMA_ONE_SHISH_PLEASE(nl, U, in_vol, dvec, stp, n_stp, nxx, nxy, shish, N_shishmax){\
00292    int m_istep, m_nind;   \
00293    static THD_fvec3 m_ncoord, m_ndicom; \
00294    static THD_ivec3 m_nind3;  \
00295    static double m_jmp, m_td[3];  \
00296    m_istep = 0; \
00297    m_jmp = 0.0;  \
00298    while (m_istep <= n_stp) {   \
00299       m_td[0] = m_jmp * U[0]; m_td[1] = m_jmp * U[1]; m_td[2] = m_jmp * U[2]; \
00300       /* get 1d coord of point */   \
00301       m_ndicom.xyz[0] = nl[0] + m_td[0] ; m_ndicom.xyz[1] = nl[1]+ m_td[1]; m_ndicom.xyz[2] = nl[2]+ m_td[2];  \
00302       m_ncoord = THD_dicomm_to_3dmm(in_vol, m_ndicom);   \
00303       m_nind3 = THD_3dmm_to_3dind(in_vol, m_ncoord);  \
00304       m_nind = m_nind3.ijk[0] + m_nind3.ijk[1] * nxx + m_nind3.ijk[2] * nxy;  \
00305       if (m_istep < N_shishmax) shish[m_istep] = dvec[m_nind];  \
00306       else break; \
00307       m_jmp += stp;  \
00308       ++m_istep;  \
00309    }  \
00310    if (m_istep < N_shishmax) shish[m_istep] = -1;  \
00311 }
00312 
00313 
00314 #endif