00001 /********************************************************************************
00002 3dEdge                                                                                  
00003 Justin Knoll 5/23/00                                                                    
00004 
00005 Generates an edge-detected mask from the input 3d+time dataset.                          
00006 Needs boundary checking for non-sagittal datasets.                                      
00007 
00008 ********************************************************************************/
00009 
00010 #include "mrilib.h"
00011 #include <stdio.h>
00012 #include <stdlib.h>
00013 #include <string.h>
00014 #include <math.h>
00015 
00016 /*******************************************************************************/
00017 void errorExit(char *message) { 
00018         
00019         fprintf(stderr, "\n\nError in 3dEdge:\n%s\n\nTry 3dEdge -help\n",message);
00020         exit(1);
00021 }
00022 
00023 /*******************************************************************************/
00024 void helpMessage()
00025 /*******************************************************************************/
00026 {
00027   printf(
00028          "Generates an edge-detected mask from the input dataset.\n"
00029          "Usage: 3dEdge [options] dataset\n"
00030          "Where the options are:\n"
00031          "  -prefix name     The prefix for the filename of the edge-masked\n"
00032          "                   dataset.\n"
00033          "  -mask name       The filename for the separately output mask.\n"
00034          "  -rms float       The floating point value (0 <= rms <= 1) for RMS\n"
00035          "                   threshhold. Voxels in the time-averaged volume with\n"
00036          "                   >= this fraction of the peak RMS for that volume\n"
00037          "                   are provisionally part of the mask, though other\n"
00038          "                   tests may subsequently remove them.\n"
00039          "  -neighbor float  The floating point value for the neighbor\n"
00040          "                   threshhold. Voxels in the RMS-derived mask are\n"
00041          "                   removed if less than this fraction of their nearest\n"
00042          "                   neighbors are included in the RMS-derived mask,\n"
00043          "                   else they are added.\n"
00044   );
00045 }
00046 
00047 /*******************************************************************************/
00048 int edgeCount(int l, int nxx, int nyy, int nvox)
00049 /*
00050   computes the number of edges the voxel has - e.g. a voxel in the corner of a sub-brick has three
00051   edges. this is necessary to apply the neighborCount to boundary voxels. 
00052 */
00053 {
00054   int count=0;
00055   if (l < (nxx*nyy)) count++;
00056   if (l < (nxx*nyy + 1)) count++;
00057 
00058   if (l > (nvox - (nxx*nyy))) count++;
00059   if (l > (nvox - (nxx*nyy) - 1)) count++;
00060 
00061 
00062   if (l % (nxx*nyy) < nxx) count++;
00063   if (l % (nxx*nyy) < (nxx + 1)) count++;
00064 
00065   if (l % (nxx*nyy) > (nxx * (nyy - 1))) count++;
00066   if (l % (nxx*nyy) > (nxx * (nyy - 1) - 1)) count++;
00067 
00068 
00069   if (((l % (nxx*nyy)) %  nxx) == 0) count++;
00070   if (((l % (nxx*nyy)) %  nxx) == 1) count++;
00071 
00072   if (((l % (nxx*nyy)) %  nxx) == (nxx - 1)) count++;
00073   if (((l % (nxx*nyy)) %  nxx) == (nxx - 2)) count++;
00074  
00075   return (count*25) - (count*count)-1;
00076 }
00077 
00078 /*******************************************************************************/
00079 float neighborCount(unsigned char *mask, int location, int nxx, int nyy, int nvox)
00080 {
00081   int c=0; /* number of neighbors */
00082   int i, j, k; /*axes */
00083   int tmpLocation;
00084   float percent;  
00085   for (i=-2 ; i <= 2 ; i++)
00086     for (j=-2 ; j <= 2 ; j++)
00087       for (k=-2 ; k <= 2 ; k++)
00088 /* 
00089    the if condition below is the bitwise OR of the three variables. It is true if at least one variable is not zero.
00090    if all three are zero, we are at the voxel itself, not a neighbor.
00091 */
00092         if (i | j | k)
00093           {
00094             tmpLocation = location;
00095             tmpLocation += i;
00096             tmpLocation += (j * nxx);
00097             if (location < (nxx*nyy)) 
00098               {
00099                 if (k < 0) continue;
00100                   tmpLocation += (k * nxx * nyy);
00101               }
00102             else if (location < ((nxx*nyy) + 1))
00103               {
00104                 if (k < -1) continue;
00105                   tmpLocation += (k * nxx * nyy);
00106               }
00107             else if (location > (nvox - (nxx*nyy)))
00108               {
00109                 if (k > 0) continue;
00110                   tmpLocation += (k * nxx * nyy);
00111               }
00112             else if (location > (nvox - (nxx*nyy) - 1))
00113               {
00114                 if (k > 1) continue;
00115                   tmpLocation += (k * nxx * nyy);
00116               }
00117             if (mask[tmpLocation]) 
00118               c++;
00119           }
00120   percent = ( (float) c / (124 - edgeCount(location, nxx, nyy, nvox)));
00121   return percent;
00122 }
00123 
00124 /*******************************************************************************/
00125 void recurseSelect(unsigned char *mask, int nxx, int nyy, int nvox, int i)
00126 {
00127   float x;
00128   float y;
00129   float z;
00130 
00131   if (mask[i] == 1)
00132     {
00133       mask[i] = 2;
00134         {
00135           recurseSelect(mask, nxx, nyy, nvox, (i + 1));
00136           recurseSelect(mask, nxx, nyy, nvox, (i - 1));
00137       
00138           recurseSelect(mask, nxx, nyy, nvox, (i + nxx));
00139           recurseSelect(mask, nxx, nyy, nvox, (i - nxx));
00140       
00141           recurseSelect(mask, nxx, nyy, nvox, (i + (nxx * nyy)));
00142           recurseSelect(mask, nxx, nyy, nvox, (i - (nxx * nyy)));
00143         }
00144     }
00145   else
00146     return;
00147 }
00148 
00149 /*******************************************************************************/
00150 int centerOfMass(unsigned char *mask, int nxx, int nyy, int nvox)
00151 {
00152   int i;
00153   float xMass=0;
00154   float yMass=0;
00155   float zMass=0;
00156   int mass=0;
00157 
00158   for(i=0 ; i < nvox ; i++)
00159     {
00160       if (mask[i])
00161         {
00162           xMass += (i % (nxx * nyy) % nxx);
00163           yMass += (i % (nxx * nyy) / nxx);
00164           zMass += (i / (nxx * nyy));
00165           mass++;
00166         }
00167     }
00168   xMass /= mass;
00169   yMass /= mass;
00170   zMass /= mass;
00171   return (int) ((int) (xMass+0.49) + ((int) (yMass+0.49) * nxx) + ((int) (zMass+0.49) * (nxx * nyy))); 
00172 }
00173 
00174 /*******************************************************************************/
00175 unsigned char* edgeDetect(THD_3dim_dataset *inputDataset, float rmsThresh, float neighbors)
00176 {
00177   int nvox, ntimes;
00178   int nxx, nyy;
00179   int i,j;
00180   unsigned char *inputTimeStep_b;  
00181   short *inputTimeStep_s;  
00182   float *inputTimeStep_f;
00183   unsigned char *mask;
00184   double *sums;
00185   double rms=0;
00186   unsigned char *newMask; /* temporary mask used to prevent order-dependent errors in neighbor counting. */
00187   THD_3dim_dataset *outputDataset=NULL;
00188 
00189   DSET_load(inputDataset);
00190   
00191   ntimes = DSET_NUM_TIMES(inputDataset);
00192   nvox = DSET_NVOX(inputDataset);
00193 
00194   nxx=inputDataset->daxes->nxx;
00195   nyy=inputDataset->daxes->nyy;
00196  
00197   sums = (double *) malloc(nvox * sizeof(double));
00198   mask = (unsigned char *) malloc(nvox * sizeof(unsigned char));
00199   newMask = (unsigned char *) malloc(nvox * sizeof(unsigned char));
00200 
00201 /*
00202   the for loop below switches to handle different subbrick types (float, byte, short), and stores the
00203   square of the values at each voxel in sums[]. 
00204 */
00205   for (j=0 ; j < ntimes ; j++)
00206     {
00207       switch (DSET_BRICK_TYPE(inputDataset, j))
00208         {
00209         case MRI_byte:
00210           {
00211             inputTimeStep_b = (unsigned char *) DSET_ARRAY(inputDataset, j);
00212             for (i=0 ; i < nvox ; i++)
00213               {  
00214                 sums[i] += inputTimeStep_b[i]; 
00215               }
00216             break;
00217           }
00218         case MRI_short:
00219           {
00220             inputTimeStep_s = (short *) DSET_ARRAY(inputDataset, j);
00221             for (i=0 ; i < nvox ; i++)
00222               {  
00223                 sums[i] += inputTimeStep_s[i]; 
00224               }
00225             break;
00226           }
00227         case MRI_float:
00228           {
00229             inputTimeStep_f = (float *) DSET_ARRAY(inputDataset, j);
00230             for (i=0 ; i < nvox ; i++)
00231               {  
00232                 sums[i] += inputTimeStep_f[i];
00233               }
00234             break;
00235           }
00236         default:
00237           errorExit("Dataset sub-brick is of unrecognized type.");
00238         }
00239 
00240      /* calculate average of voxels over time */   
00241     } 
00242   for (i=0 ; i < nvox ; i++)
00243     {
00244       sums[i] /= ntimes;
00245       rms += sums[i]*sums[i];
00246     }
00247   rms /= nvox;
00248   rms = sqrt(rms);
00249   
00250 
00251   for (i=0 ; i < nvox ; i++)
00252     {
00253       mask[i] = (sums[i] >= (rms * rmsThresh));
00254     }
00255 /* fill holes and remove outlayers based on number of nearest neighbors */
00256    for (i = 0 ; i < nvox ; i++)  
00257      if (neighborCount(mask, i, nxx, nyy, nvox) < neighbors) newMask[i]=0;
00258      else
00259        newMask[i]=mask[i];
00260    for (i = 0 ; i < nvox ; i++) 
00261      mask[i] = newMask[i];
00262    for (i = 0 ; i < nvox ; i++)  
00263       if (neighborCount(mask, i, nxx, nyy, nvox) > neighbors) newMask[i]=1;
00264    for (i = 0 ; i < nvox ; i++) 
00265      mask[i] = newMask[i];
00266 
00267    i = centerOfMass(mask, nxx, nyy, nvox);
00268 /*   printf("i=%d\n", i);*/
00269    recurseSelect(mask, nxx, nyy, nvox, i);
00270 
00271    for (i = 0 ; i < nvox ; i++) 
00272      mask[i] = mask[i]==2;
00273    
00274    free(sums);
00275    free(newMask);
00276    return mask;
00277 }
00278 
00279 /*******************************************************************************/
00280 int applyMask(THD_3dim_dataset *inputDataset, unsigned char *mask, char *prefix, int argc, char *argv[])
00281 /* argc and argv are needed for to put the command-line args in the history. */
00282 {
00283   int ierr;
00284   int i,j;
00285   int ntimes, nvox;
00286   double* timeStep;
00287   THD_3dim_dataset *outputDataset;
00288 
00289   unsigned char *inputTimeStep_b;  
00290   short *inputTimeStep_s;  
00291   float *inputTimeStep_f;
00292 
00293   unsigned char *outputTimeStep_b;  
00294   short *outputTimeStep_s;  
00295   float *outputTimeStep_f;
00296 
00297   ntimes = DSET_NUM_TIMES(inputDataset);
00298   nvox = DSET_NVOX(inputDataset);
00299   outputDataset = EDIT_empty_copy(inputDataset);
00300 
00301   ierr = EDIT_dset_items(outputDataset,
00302                            ADN_type, HEAD_ANAT_TYPE,     /* dataset type */
00303                            ADN_func_type, FUNC_FIM_TYPE, /* dataset type */
00304                            ADN_prefix, prefix,
00305                          ADN_none);
00306 
00307   for (j=0 ; j < ntimes ; j++)
00308     {
00309       switch (DSET_BRICK_TYPE(inputDataset, j))
00310         {
00311         case MRI_byte:
00312           {
00313             inputTimeStep_b = (unsigned char *) DSET_ARRAY(inputDataset, j);
00314             outputTimeStep_b = (unsigned char *) malloc(sizeof(unsigned char) * nvox);
00315             for (i=0 ; i < nvox ; i++)
00316               {  
00317                 if (mask[i])
00318                   outputTimeStep_b[i] = inputTimeStep_b[i];
00319                 else
00320                   outputTimeStep_b[i] = 0;
00321               }
00322             EDIT_substitute_brick(outputDataset, j, DSET_BRICK_TYPE(inputDataset, j), outputTimeStep_b);
00323             break;
00324           }
00325         case MRI_short:
00326           {
00327             inputTimeStep_s = (short *) DSET_ARRAY(inputDataset, j);
00328             outputTimeStep_s = (short *) malloc(sizeof(short) * nvox);
00329             for (i=0 ; i < nvox ; i++)
00330               {  
00331                 if (mask[i])
00332                   outputTimeStep_s[i] = inputTimeStep_s[i];
00333                 else
00334                   outputTimeStep_s[i] = 0;
00335               }
00336             EDIT_substitute_brick(outputDataset, j, DSET_BRICK_TYPE(inputDataset, j), outputTimeStep_s);
00337             break;
00338           }
00339         case MRI_float:
00340           {
00341             inputTimeStep_f = (float *) DSET_ARRAY(inputDataset, j);
00342             outputTimeStep_f = (float *) malloc(sizeof(inputTimeStep_f));
00343             for (i=0 ; i < nvox ; i++)
00344               {  
00345                 if (mask[i])
00346                   outputTimeStep_f[i] = inputTimeStep_f[i];
00347                 else
00348                   outputTimeStep_f[i] = 0;
00349               }
00350             EDIT_substitute_brick(outputDataset, j, DSET_BRICK_TYPE(inputDataset, j), outputTimeStep_f);
00351             break;
00352           }
00353         default:
00354           errorExit("Dataset sub-brick is of unrecognized type.");
00355         }
00356     } 
00357     tross_Copy_History( inputDataset, outputDataset );
00358     tross_Make_History( "3dedge" , argc, argv, outputDataset) ;
00359     DSET_write(outputDataset);  
00360 /*    DSET_unload(inputDataset);
00361     DSET_unload(outputDataset); */
00362 }
00363 
00364 void saveMask(THD_3dim_dataset *inputDataset, unsigned char *mask, char *prefix)
00365 {
00366   int ierr;
00367   THD_3dim_dataset *outputDataset;
00368   
00369   outputDataset = EDIT_empty_copy(inputDataset);
00370   ierr = EDIT_dset_items(outputDataset,
00371                            ADN_prefix, prefix,      /* prefix */
00372                            ADN_datum_all, MRI_byte,      /* datum type */
00373                            ADN_nvals, 1,                 /* number of 3d sub-bricks */
00374                            ADN_type, HEAD_FUNC_TYPE,     /* dataset type */
00375                            ADN_func_type, FUNC_FIM_TYPE, /* dataset type */
00376                            ADN_ntt, 1,                   /* number of time steps */
00377                          ADN_none);
00378   EDIT_substitute_brick(outputDataset, 0, MRI_byte, mask);
00379   DSET_write(outputDataset);
00380 }
00381 
00382 
00383 THD_3dim_dataset* openDataset(char *name)
00384 {
00385   THD_3dim_dataset *newDataset = THD_open_dataset(name);
00386 
00387 /* check for cases which indicate errors or are not supported. */
00388   if (newDataset == NULL)
00389     errorExit("Cannot open new dataset!");
00390   if (DSET_BRICK_TYPE(newDataset, 0) == MRI_complex)
00391     errorExit("I can't deal with complex datasets.");
00392   if (DSET_BRICK_FACTOR(newDataset, 0))
00393     errorExit("I can't deal with datasets containing scaling factors");
00394 /* else, make it so   */
00395   return newDataset; 
00396 }
00397 
00398 /*******************************************************************************/
00399 int main(int argc, char *argv[])
00400 {
00401   int narg=1;
00402   int stringLength;
00403   float rms=1.0;
00404   float neighbors=0.5;
00405   char *prefix = "masked";
00406   char *maskName = "mask";
00407   unsigned char *mask;
00408   THD_3dim_dataset *inputDataset=NULL;
00409   
00410   if (argc < 2 || (strcmp(argv[1], "-help") == 0))
00411     {
00412       helpMessage();
00413       exit(1);
00414     }
00415 
00416   while((narg < argc) && (argv[narg][0] == '-'))
00417     {
00418       if (strcmp(argv[narg], "-prefix") == 0)
00419         {
00420           narg++;
00421           stringLength = strlen(argv[narg]);
00422           prefix = (char *) malloc((stringLength + 1) * sizeof(char));
00423           strcpy(prefix, argv[narg++]);
00424 /*        printf("prefix: %s\n", prefix); */
00425           continue;
00426         }
00427       if (strcmp(argv[narg], "-mask") == 0)
00428         {
00429           narg++;
00430           stringLength = strlen(argv[narg]);
00431           maskName = (char *) malloc((stringLength + 1) * sizeof(char));
00432           strcpy(maskName, argv[narg++]);
00433 /*        printf("mask: %s\n", maskName); */
00434           continue;
00435         }
00436       if (strcmp(argv[narg], "-rms") == 0)
00437         {
00438           narg++;
00439           rms = (float)atof(argv[narg++]);
00440 /*        printf("rms: %f\n", rms);*/
00441           continue;
00442         }
00443       if (strcmp(argv[narg], "-neighbor") == 0)
00444         {
00445           narg++;
00446           neighbors = (float)atof(argv[narg++]);
00447 /*        printf("neighbors: %f\n", neighbors);*/
00448           continue;
00449         }
00450       else
00451         {
00452           printf("*** unrecognized option, %s\n", argv[narg]);
00453           narg++;
00454         }
00455     }
00456   inputDataset = openDataset(argv[narg]);
00457   mask = edgeDetect(inputDataset, rms, neighbors);
00458   if (strcmp(maskName, "mask") != 0)
00459         saveMask(inputDataset, mask, maskName);
00460   applyMask(inputDataset, mask, prefix, argc, argv);
00461   return 0;
00462 }
00463 
00464 
00465 
00466 
00467 
00468 
00469 
00470 
00471 
00472 
00473 
00474 
00475 
00476 
00477 
00478 
00479 
00480 
00481 
00482 
00483 
00484