00001 #include "SUMA_suma.h"
00002 
00003 extern SUMA_CommonFields *SUMAg_CF;
00004 extern SUMA_DO *SUMAg_DOv;
00005 extern SUMA_SurfaceViewer *SUMAg_SVv;
00006 extern int SUMAg_N_SVv; 
00007 extern int SUMAg_N_DOv;  
00008 
00009 
00010 SUMA_M2M_STRUCT *SUMA_NewM2M(char *SO1_id, int N_SO1_nodes, char *SO2_id)
00011 {
00012    static char FuncName[]={"SUMA_NewM2M"};
00013    SUMA_M2M_STRUCT *M2M=NULL;
00014    SUMA_ENTRY;
00015    
00016    if (!SO1_id || !SO2_id) SUMA_RETURN(M2M);
00017    
00018    M2M = (SUMA_M2M_STRUCT*)SUMA_malloc(sizeof(SUMA_M2M_STRUCT));
00019    
00020    M2M->M1Nn = N_SO1_nodes;
00021    M2M->M1n = (int*)SUMA_calloc(M2M->M1Nn, sizeof(int));
00022    M2M->M2t_M1n = (int*)SUMA_calloc(M2M->M1Nn, sizeof(int));
00023    M2M->M2Nne_M1n = (int*)SUMA_calloc(M2M->M1Nn, sizeof(int));
00024    M2M->M2ne_M1n = (int**)SUMA_calloc(M2M->M1Nn, sizeof(int*));
00025    M2M->M2pb_M1n = (float *)SUMA_calloc(2*M2M->M1Nn, sizeof(float));
00026    M2M->M2p_M1n = (float *)SUMA_calloc(3*M2M->M1Nn, sizeof(float));
00027    M2M->PD = (double *)SUMA_calloc(M2M->M1Nn, sizeof(double));
00028    M2M->M2we_M1n = (double**)SUMA_calloc(M2M->M1Nn, sizeof(double*));
00029    if (!M2M->M1n || !M2M->M2t_M1n || !M2M->M2Nne_M1n || !M2M->M2ne_M1n || !M2M->M2we_M1n) {
00030       SUMA_SL_Crit("Failed to allocate");
00031       SUMA_RETURN(NULL);
00032    }
00033    
00034    M2M->M1_IDcode = SUMA_copy_string(SO1_id);
00035    M2M->M2_IDcode = SUMA_copy_string(SO2_id);
00036    
00037    
00038    SUMA_RETURN(M2M);
00039 }  
00040 
00041 SUMA_M2M_STRUCT *SUMA_FreeM2M(SUMA_M2M_STRUCT *M2M)
00042 {
00043    static char FuncName[]={"SUMA_FreeM2M"};
00044    int i;
00045    
00046    SUMA_ENTRY;
00047    
00048    if (!M2M) SUMA_RETURN(NULL);
00049    if (M2M->M2we_M1n) {
00050       for (i=0; i<M2M->M1Nn; ++i) { 
00051          if (M2M->M2we_M1n[i]) {
00052             SUMA_free(M2M->M2we_M1n[i]);
00053             M2M->M2we_M1n[i] = NULL;
00054          }
00055       }
00056       SUMA_free(M2M->M2we_M1n);
00057       M2M->M2we_M1n = NULL;
00058    }
00059    if (M2M->M2ne_M1n) {
00060       for (i=0; i<M2M->M1Nn; ++i) { 
00061          if (M2M->M2ne_M1n[i]) {
00062             SUMA_free(M2M->M2ne_M1n[i]);
00063             M2M->M2ne_M1n[i] = NULL;
00064          }
00065       }
00066       SUMA_free(M2M->M2ne_M1n);
00067       M2M->M2ne_M1n = NULL;
00068    }
00069    if (M2M->M1n) SUMA_free(M2M->M1n); M2M->M1n = NULL;
00070    if (M2M->M2t_M1n) SUMA_free(M2M->M2t_M1n); M2M->M2t_M1n= NULL;
00071    if (M2M->M2Nne_M1n) SUMA_free(M2M->M2Nne_M1n); M2M->M2Nne_M1n = NULL;
00072    if (M2M->M2pb_M1n) SUMA_free(M2M->M2pb_M1n); M2M->M2pb_M1n = NULL;
00073    if (M2M->M2p_M1n) SUMA_free(M2M->M2p_M1n); M2M->M2p_M1n = NULL;
00074    if (M2M->PD) SUMA_free(M2M->PD); M2M->PD = NULL;
00075    if (M2M->M1_IDcode) SUMA_free(M2M->M1_IDcode); M2M->M1_IDcode = NULL;
00076    if (M2M->M2_IDcode) SUMA_free(M2M->M2_IDcode); M2M->M2_IDcode = NULL;
00077 
00078    SUMA_free(M2M);
00079    SUMA_RETURN(NULL);     
00080 } 
00081   
00082 char *SUMA_M2M_node_Info (SUMA_M2M_STRUCT *M2M, int node)
00083 {
00084    static char FuncName[]={"SUMA_M2M_node_Info"};
00085    char *s = NULL;
00086    SUMA_STRING *SS = NULL;
00087    int i, found, j;
00088    
00089 
00090    SS = SUMA_StringAppend(NULL, NULL);
00091       
00092    if (!M2M) { SS = SUMA_StringAppend(SS,"NULL M2M"); goto CLEAN_RETURN; }
00093    
00094    if (M2M->M1_IDcode) { SS = SUMA_StringAppend_va(SS, "M1_IDcode %s\n", M2M->M1_IDcode); }
00095    else { SS = SUMA_StringAppend_va(SS, "M1_IDcode is NULL\n"); }
00096    if (M2M->M2_IDcode) { SS = SUMA_StringAppend_va(SS, "M2_IDcode %s\n", M2M->M2_IDcode); }
00097    else { SS = SUMA_StringAppend_va(SS, "M2_IDcode is NULL\n"); }
00098   
00099    i = 0; found = 0;
00100    while (i < M2M->M1Nn && !found) {
00101       if (M2M->M1n[i] == node) {
00102          found = 1;
00103       } else ++i;
00104    }
00105    
00106    if (!found) { SS = SUMA_StringAppend_va (SS, "Node %d not found in M2M->M1n", node); goto CLEAN_RETURN; }
00107    
00108    SS = SUMA_StringAppend_va (SS, "Mapping results for node %d (n1) of mesh 1 (M1):\n", M2M->M1n[i]);
00109    SS = SUMA_StringAppend_va (SS, "Index of triangle (t2) in mesh 2 (M2) hosting n1: %d\n", M2M->M2t_M1n[i]);
00110    SS = SUMA_StringAppend_va (SS, "Projection coordinates in t2 (%f,%f,%f)\n", M2M->M2p_M1n[3*i], M2M->M2p_M1n[3*i+1], M2M->M2p_M1n[3*i+2]);
00111    SS = SUMA_StringAppend_va (SS, "Projection barycentric coordinates in t2 (%g,%g)\n", M2M->M2pb_M1n[2*i], M2M->M2pb_M1n[2*i+1]);
00112    SS = SUMA_StringAppend_va (SS, "Projection distance of n1 onto t2 is: %g\n", M2M->PD[i]);
00113    SS = SUMA_StringAppend_va (SS, "Number of nodes (n2) in M2 considered neighbors to n1: %d\n", M2M->M2Nne_M1n[i]);
00114    SS = SUMA_StringAppend_va (SS, "n2   \tw2weight\n");
00115    for (j=0; j< M2M->M2Nne_M1n[i]; ++j) {
00116       SS = SUMA_StringAppend_va (SS, "%s\t%g\n", MV_format_fval2(M2M->M2ne_M1n[i][j], 5), M2M->M2we_M1n[i][j]);
00117    }
00118    
00119    
00120    CLEAN_RETURN:
00121    SUMA_SS2S(SS,s);
00122    
00123    SUMA_RETURN(s);
00124 } 
00125 /*!
00126    \brief a function to find the mapping on nodes from SO1's mesh M1 to SO2's mesh M2
00127    The mapping method is Nearest Neighbor, each node nj considered from M1, will get for neighbors 
00128    the nodes forming the triangle in M2 that nj projects to.
00129    
00130    \param SO1 (SUMA_SurfaceObject *) surface 1
00131    \param SO2 (SUMA_SurfaceObject *) surface 2
00132    \param NL_1 (int *) list of node indices to consider from surface 1 
00133                IF NULL, then all nodes of surface 1 are considered
00134    \param N_NL_1 (int) number of values in NL_1 (= SO1->N_Node if NL_1 is NULL)
00135    \param PD_1 (float *) vector of SO1->N_Node triplets specifying the direction
00136                          of the projection of nodes in NL_1. 
00137                          PD_1[3*NL_1[j]], PD_1[3*NL_1[j]+1], PD_1[3*NL_1[j]+2] are the projection directions
00138                          of node NL_1[j].
00139                          If NULL then the projection direction for node NL_1[j] is:
00140                          SO1->NodeNormList[3*NL_1[j]], SO1->NodeNormList[3*NL_1[j]+1], SO1->NodeNormList[3*NL_1[j]+2] is used
00141    \param dlim (float) maximum distance to search in each direction, along the projection vector. 
00142                         If 0, then a default of 100.0 is used.
00143    \return M2M (SUMA_M2M_STRUCT *) Mesh to Mesh mapping structure, see SUMA_SurfaceToSurface.h for details
00144 */ 
00145 SUMA_M2M_STRUCT *SUMA_GetM2M_NN( SUMA_SurfaceObject *SO1, SUMA_SurfaceObject *SO2,
00146                                  int *oNL_1, int N_NL_1, float *PD_1, float dlim, 
00147                                  int NodeDbg)
00148 {
00149    static char FuncName[]={"SUMA_GetM2M_NN"};
00150    SUMA_M2M_STRUCT *M2M = NULL;
00151    int *NL_1 = NULL;
00152    int j, id, id2, nnt, k, nj, t3, j3;
00153    float *triNode0, *triNode1, *triNode2, *hit;
00154    float delta_t;
00155    double *wv, wgt[3], weight_tot; 
00156    float P0[3], P1[3], P2[3], N0[3];
00157    float Points[2][3];
00158    SUMA_MT_INTERSECT_TRIANGLE *MTI = NULL;
00159    struct timeval tt; 
00160    SUMA_Boolean LocalHead = NOPE;
00161    
00162    SUMA_ENTRY;
00163    
00164    SUMA_etime (&tt, 0);
00165    
00166    if (!SO1 || !SO2) { SUMA_SL_Err("NULL input"); goto CLEAN_EXIT; }
00167    if (!oNL_1) N_NL_1 = SO1->N_Node;
00168    if (N_NL_1 < 1) { SUMA_SL_Err("No nodes to consider"); goto CLEAN_EXIT; }
00169    if (dlim <= 0) dlim = 100.0; 
00170    /* start filling M2M */
00171    M2M = SUMA_NewM2M(SO1->idcode_str, N_NL_1, SO2->idcode_str);
00172    if (!M2M) { SUMA_SL_Crit("Failed to create M2M"); goto CLEAN_EXIT; }
00173    
00174    /* fill up M2M->M1n */
00175    if (!oNL_1) { for (j=0; j<N_NL_1; ++j) M2M->M1n[j]= j; }
00176    else { for (j=0; j<N_NL_1; ++j) M2M->M1n[j]= oNL_1[j]; }
00177    
00178    if (!PD_1) PD_1 = SO1->NodeNormList;
00179    
00180    for (j = 0; j < M2M->M1Nn; ++j) {
00181       j3    = 3 * j;
00182       nj    = M2M->M1n[j];
00183       id    = SO1->NodeDim * nj;
00184       P0[0] = SO1->NodeList[id];
00185       P0[1] = SO1->NodeList[id+1];
00186       P0[2] = SO1->NodeList[id+2];
00187 
00188       N0[0] = PD_1[id  ];
00189       N0[1] = PD_1[id+1];
00190       N0[2] = PD_1[id+2];
00191       
00192       SUMA_POINT_AT_DISTANCE(N0, P0, dlim, Points);
00193       
00194       P1[0] = Points[0][0];
00195       P1[1] = Points[0][1];
00196       P1[2] = Points[0][2];
00197       P2[0] = Points[1][0];
00198       P2[1] = Points[1][1];
00199       P2[2] = Points[1][2];
00200  
00201       /* now determine the distance along normal */
00202       MTI = SUMA_MT_intersect_triangle(P0, P1, SO2->NodeList, SO2->N_Node, SO2->FaceSetList, SO2->N_FaceSet, MTI);
00203       if (LocalHead) fprintf(SUMA_STDERR,"%s: number of hits for node %d : %d\n", FuncName, nj, MTI->N_hits);  
00204       if (MTI->N_hits ==0) {
00205          if (LocalHead) fprintf(SUMA_STDERR, "%s: Could not find hit for node %d in either direction.\n", FuncName, nj);
00206          M2M->M2Nne_M1n[j] = 0;
00207          M2M->M2t_M1n[j] = -1;
00208          M2M->PD[j] = 0.0;
00209          M2M->M2pb_M1n[2*j  ] = -1.0;
00210          M2M->M2pb_M1n[2*j+1] = -1.0;
00211          j3 = 3*j; hit = &(M2M->M2p_M1n[j3]);
00212          hit[0] = -1.0;
00213          hit[1] = -1.0;
00214          hit[2] = -1.0;      
00215       } else {
00216          if (LocalHead) {
00217              for (k = 0; k < MTI->N_el; k++) {
00218                if (MTI->isHit[k] == YUP) fprintf(SUMA_STDERR, "%s: hit %d: %f (%f, %f)\n",FuncName, k, MTI->t[k], MTI->u[k], MTI->v[k]);
00219             }
00220          }
00221          M2M->M2t_M1n[j] = MTI->ifacemin;
00222          M2M->PD[j] = MTI->t[MTI->ifacemin];
00223          
00224          /* create the neighborhood list */
00225          M2M->M2Nne_M1n[j] = 3; /* always three neighbors for NN method, first node always the closest! */
00226          M2M->M2ne_M1n[j] = (int *)SUMA_malloc(M2M->M2Nne_M1n[j]*sizeof(int));
00227                                          *(M2M->M2ne_M1n[j]  ) = MTI->inodemin;   /*! That is the closest node */
00228          nnt = (MTI->inodeminlocal+1)%3; *(M2M->M2ne_M1n[j]+1) = SO2->FaceSetList[3*MTI->ifacemin+nnt]; /*index of second node in t2 */
00229          nnt = (MTI->inodeminlocal+2)%3; *(M2M->M2ne_M1n[j]+2) = SO2->FaceSetList[3*MTI->ifacemin+nnt]; /*index of third node in t2 */
00230 
00231          /* Now for the weights of each neighboring node */
00232          M2M->M2we_M1n[j] = (double *)SUMA_malloc(M2M->M2Nne_M1n[j]*sizeof(double));
00233 
00234          /* store the hit location */
00235          j3 = 3*j; hit = &(M2M->M2p_M1n[j3]);
00236          hit[0] = MTI->P[0];
00237          hit[1] = MTI->P[1];
00238          hit[2] = MTI->P[2];
00239          if (M2M->M1n[j] == NodeDbg) {
00240             fprintf(SUMA_STDERR, "%s: Hit coords for node %d of M1: \n"
00241                                  "%f %f %f\n"
00242                                  "%f %f %f\n", FuncName, M2M->M1n[j], hit[0], hit[1], hit[2], 
00243                                  M2M->M2p_M1n[j3], M2M->M2p_M1n[j3+1], M2M->M2p_M1n[j3+2]);
00244          }
00245          /* store the barycentric (u,v) location of intersection */
00246          M2M->M2pb_M1n[2*j  ] = MTI->u[MTI->ifacemin];
00247          M2M->M2pb_M1n[2*j+1] = MTI->v[MTI->ifacemin];
00248 
00249          /**determine weights which are the barycetric corrdinates of the intersection node
00250             The intersected triangle is formed by the 1st three nodes stored in M2ne_M1n
00251             RESPECT THE ORDER in M2ne_M1n */
00252             /* get pointers to x,y,z of each node of intersected triangle*/
00253             t3 = 3*MTI->ifacemin;
00254             triNode0 = &(SO2->NodeList[ 3*M2M->M2ne_M1n[j][0] ]);
00255             triNode1 = &(SO2->NodeList[ 3*M2M->M2ne_M1n[j][1] ]);
00256             triNode2 = &(SO2->NodeList[ 3*M2M->M2ne_M1n[j][2] ]);
00257 
00258          SUMA_TRI_AREA( ((MTI->P)), triNode1, triNode2, wgt[0] ); 
00259          SUMA_TRI_AREA( ((MTI->P)), triNode0, triNode2, wgt[1] ); 
00260          SUMA_TRI_AREA( ((MTI->P)), triNode0, triNode1, wgt[2] ); 
00261 
00262          weight_tot =  wgt[0] + wgt[1] + wgt[2];
00263 
00264          wv = M2M->M2we_M1n[j];
00265          if (weight_tot) {
00266             wv[0] = wgt[0] / weight_tot;
00267             wv[1] = wgt[1] / weight_tot;
00268             wv[2] = wgt[2] / weight_tot;
00269          }else { /* some triangles have zero area in FreeSurfer surfaces */
00270             wv[0] = wv[1] =  wv[2] = 1.0/3.0;
00271          }
00272       }
00273 
00274       
00275       if (!(j%500) && j) {
00276          delta_t = SUMA_etime(&tt, 1);
00277          fprintf (SUMA_STDERR, " [%d]/[%d] %.2f/100%% completed. Dt = %.2f min done of %.2f min total\r" ,  j, N_NL_1, (float)j / N_NL_1 * 100, delta_t/60, delta_t/j * N_NL_1/60);
00278          if (LocalHead) {
00279             char *s = NULL;
00280             s = SUMA_M2M_node_Info(M2M, M2M->M1n[j]);
00281             fprintf(SUMA_STDERR,"\n***\n%s\n***\n", s); 
00282             SUMA_free(s); s = NULL;
00283          }
00284       }
00285 
00286    }
00287 
00288    if (LocalHead) {
00289       delta_t = SUMA_etime(&tt, 1);
00290       fprintf (SUMA_STDERR, " [%d]/[%d] %.2f/100%% completed. Dt = %.2f min done of %.2f min total\r" ,  j, N_NL_1, (float)j / N_NL_1 * 100, delta_t/60, delta_t/j * N_NL_1/60);
00291       fprintf (SUMA_STDERR, "\n");
00292    }
00293 
00294    CLEAN_EXIT:
00295    if (MTI) MTI = SUMA_Free_MT_intersect_triangle(MTI); 
00296 
00297    
00298    SUMA_RETURN(M2M);
00299 }
00300 
00301 /*!
00302    \brief A function to interpolate data from one mesh onto another
00303    \param M2M (SUMA_M2M_STRUCT *) 
00304    \param far_data (float *) Data from mesh 2. The vector in far_data can 
00305                              represent an nvec * ncol matrix of values stored
00306                              in row major or column major order. Think of each
00307                              column as a separate sub-brick.
00308    \param ncol (int) number of columns in far_data
00309    \param nrow (int) number of rows in far_data (this number must be equal to the number of nodes on mesh 2!)
00310    \param d_order (SUMA_INDEXING_ORDER) SUMA_ROW_MAJOR, i.e. xyz xyz xyz xyz
00311                                        SUMA_COLUMN_MAJOR, i.e. xxxx yyyy zzzz
00312    \param useCloset (int) 1 means use only data form the closest node
00313                     0 means use data from all neighbors in M2M
00314    \return dt (float *) interpolation of far_data from mesh 2 (M2) onto nodes of M1
00315                         dt is ncol*M2M->M1Nn in the same order as d_order
00316 */
00317 float *SUMA_M2M_interpolate(SUMA_M2M_STRUCT *M2M, float *far_data, int ncol, int nrow, SUMA_INDEXING_ORDER d_order, int useClosest )
00318 {
00319    static char FuncName[]={"SUMA_M2M_interpolate"};
00320    int j, k, i, nk, nkid, njid, N_k, nj;
00321    float *dt=NULL;
00322    SUMA_Boolean LocalHead = NOPE;
00323 
00324    SUMA_ENTRY;
00325    
00326    if (!M2M || !far_data) {
00327       SUMA_SL_Err("NULL input");
00328       SUMA_RETURN(dt);
00329    }
00330    /* allocation */
00331    dt = (float *)SUMA_calloc(M2M->M1Nn*ncol, sizeof(float));
00332    if (!dt) { SUMA_SL_Crit("Failed to allocate"); SUMA_RETURN(dt); }
00333 
00334    /* here we go */
00335    if (d_order == SUMA_ROW_MAJOR) {
00336       if (!useClosest) {
00337          SUMA_LH("Using all neighbors, ROW MAJOR interpolation");
00338          for (j=0; j<M2M->M1Nn; ++j) {
00339             nj = M2M->M1n[j]; /* node on M1 */
00340             njid = j*ncol; /* ROW MAJOR BABY */
00341             N_k = M2M->M2Nne_M1n[j]; 
00342             for (i=0; i<ncol; ++i) { /* for each column */
00343                dt[njid+i] = 0.0;
00344                for (k=0; k<N_k; ++k) { /* for each neighbor */
00345                   nk = M2M->M2ne_M1n[j][k]; 
00346                   nkid = nk * ncol; /* ROW MAJOR BABY */
00347                   dt[njid+i] += far_data[nkid+i] * M2M->M2we_M1n[j][k]; 
00348                }
00349             }
00350          }   
00351       } else {
00352          SUMA_LH("Using immediate neighbor, ROW MAJOR interpolation");
00353          k = 0; /* just the closest neighbor  */
00354          for (j=0; j<M2M->M1Nn; ++j) {
00355             nj = M2M->M1n[j]; /* node on M1 */
00356             njid = j*ncol; /* ROW MAJOR BABY */
00357             for (i=0; i<ncol; ++i) { /* for each column */
00358                dt[njid+i] = 0.0;      
00359                if (M2M->M2Nne_M1n[j]) { /* Some nodes have no neighbors! */
00360                   /* k = 0, set above, just the closest neighbor  */
00361                   nk = M2M->M2ne_M1n[j][k]; 
00362                   nkid = nk * ncol; /* ROW MAJOR BABY */
00363                   dt[njid+i] += far_data[nkid+i];
00364                }
00365             }
00366          }
00367       }
00368    } else if (d_order == SUMA_COLUMN_MAJOR) {
00369       if (!useClosest) {
00370          SUMA_LH("Using all neighbors, COLUMN MAJOR interpolation");
00371          for (i=0; i<ncol; ++i) { /* for each column */
00372             for (j=0; j<M2M->M1Nn; ++j) { /* for each node on M1 */
00373                nj = M2M->M1n[j]; 
00374                njid = j+i*M2M->M1Nn; /* index of nj's ith column entry into dt, COLUMN MAJOR BABY */
00375                dt[njid] = 0;
00376                N_k = M2M->M2Nne_M1n[j]; 
00377                for (k=0; k<N_k; ++k) { /* for each neighbor */
00378                   nk = M2M->M2ne_M1n[j][k];
00379                   nkid = nk + i*nrow; /* index of nj's kth neighbor's data into far_data, COLUMN MAJOR BABY */ 
00380                   dt[njid] += far_data[nkid]* M2M->M2we_M1n[j][k];
00381                }
00382             }
00383          }
00384       } else {
00385          SUMA_LH("Using immediate neighbor, COLUMN MAJOR interpolation");
00386          k = 0;   /* just the closest neighbor  */
00387          for (i=0; i<ncol; ++i) { /* for each column */
00388             for (j=0; j<M2M->M1Nn; ++j) { /* for each node on M1 */
00389                nj = M2M->M1n[j]; 
00390                njid = j+i*M2M->M1Nn; /* index of nj's ith column entry into dt, COLUMN MAJOR BABY */
00391                dt[njid] = 0;
00392                N_k = M2M->M2Nne_M1n[j]; 
00393                if (M2M->M2Nne_M1n[j]) { /* Some nodes have no neighbors! */
00394                   /* k = 0, set above, just the closest neighbor  */
00395                   nk = M2M->M2ne_M1n[j][k];
00396                   nkid = nk + i*nrow; /* index of nj's kth neighbor's data into far_data, COLUMN MAJOR BABY */ 
00397                   dt[njid] += far_data[nkid];
00398                }
00399             }
00400          }
00401       }
00402    } else {
00403       SUMA_SL_Err("Bad order option");
00404       SUMA_free(dt); dt = NULL; SUMA_RETURN(dt); 
00405    }
00406 
00407    SUMA_RETURN(dt); 
00408 }