00001 /*Program to compute the minimum distance metric between two surfaces */
00002 
00003 /* Header FILES */
00004    
00005 #include "SUMA_suma.h"
00006 
00007 SUMA_SurfaceViewer *SUMAg_cSV; /*!< Global pointer to current Surface Viewer structure*/
00008 SUMA_SurfaceViewer *SUMAg_SVv; /*!< Global pointer to the vector containing the various Surface Viewer Structures 
00009                                     SUMAg_SVv contains SUMA_MAX_SURF_VIEWERS structures */
00010 int SUMAg_N_SVv = 0; /*!< Number of SVs realized by X */
00011 SUMA_DO *SUMAg_DOv;     /*!< Global pointer to Displayable Object structure vector*/
00012 int SUMAg_N_DOv = 0; /*!< Number of DOs stored in DOv */
00013 SUMA_CommonFields *SUMAg_CF; /*!< Global pointer to structure containing info common to all viewers */
00014 
00015 /* function declarations */
00016 float maximum(int N, float *inarray);
00017 float minimum(int N, float *inarray);
00018 void cmp_surf_usage ();
00019 
00020 /*
00021 File : SUMA_compare_surfaces_direction.c
00022 authors : Shruti Japee and Ziad Saad
00023 Date : Thu August 1 10:24:01 EST 2002
00024 Purpose : Two determine the min distance between nodes on one surface to faces on another.
00025 */
00026 
00027 int main (int argc,char *argv[])
00028 {
00029   static char FuncName[]={"SUMA_compare_surfaces"}; 
00030   SUMA_SurfSpecFile Spec;
00031   char *specfilename = NULL;
00032   int kar, id;
00033   SUMA_Boolean brk;
00034   SUMA_Boolean SurfIn = NOPE;
00035 
00036   SUMA_SurfaceObject *Surf1 = NULL, *Surf2=NULL;
00037   char *Surf1_FileName=NULL;
00038   char *Surf2_FileName=NULL;
00039   char *Vol1Parent_FileName=NULL;
00040   char *Vol2Parent_FileName=NULL;
00041   
00042   /* for SureFit surface */
00043   SUMA_SFname *Surf1_SFName=NULL, *Surf2_SFName=NULL;
00044    
00045   /* other variables */
00046   int trouble;
00047   int i,j,k;
00048   int num_nodes1;
00049   int num_nodes2;
00050   int onenode = -1, istart = -1, istop = -1, FailedDistance = -1; 
00051   float P0[3];
00052   float delta_t; 
00053   float P1[3];
00054   float P2[3];
00055   float N0[3];
00056   float maxdistance, mindistance;
00057   float *distance, distanceneg, distancepos;
00058   float Points[2][3];
00059   SUMA_COLOR_MAP *MyColMap;
00060   SUMA_SCALE_TO_MAP_OPT *MyOpt;
00061   SUMA_COLOR_SCALED_VECT * MySV;
00062   SUMA_MT_INTERSECT_TRIANGLE *triangle = NULL;
00063   SUMA_SurfaceObject *SO1, *SO2;
00064   struct timeval tt; 
00065   FILE *colorfile;
00066   FILE *distancefile;
00067   FILE *segfile1, *segfile2, *segfile3,*trianglesfile;
00068   char colorfilename[1000];
00069   char distancefilename[1000];
00070   char segmentfilename1[1000],segmentfilename2[1000],segmentfilename3[1000],trianglesfilename[1000];
00071   char *tag1 = NULL;
00072   char *tag2 = NULL;
00073   char *state1 = NULL;
00074   char *state2 = NULL;
00075   char *hemi = NULL;
00076   float B_dim[3];
00077   char *fout = NULL, *fname=NULL;
00078   SUMA_Boolean KeepMTI = YUP, Partial = NOPE, SkipConsistent = NOPE;
00079 
00080   SUMA_mainENTRY;
00081   SUMA_STANDALONE_INIT;
00082   
00083   if (argc < 7) {
00084     cmp_surf_usage();
00085     exit (1);
00086   }
00087   
00088   /* read in the surfaces */
00089   kar = 1;
00090   brk = NOPE;
00091   SurfIn = NOPE;
00092    Partial = NOPE;
00093    SkipConsistent = NOPE;
00094   while (kar < argc) {
00095     /* loop accross command line options */
00096     if ((strcmp(argv[kar], "-h") == 0) || (strcmp(argv[kar], "-help") == 0)) {
00097       cmp_surf_usage ();
00098       exit (1);
00099     }
00100    
00101    if (!brk && (strcmp(argv[kar], "-nocons") == 0)) {
00102       SkipConsistent = YUP;
00103       brk = YUP;
00104    }
00105            
00106     if (!brk && (strcmp(argv[kar], "-sv1")) == 0) {
00107       kar ++;
00108       if (kar >= argc) {
00109         fprintf (SUMA_STDERR, "need argument after -sv1");
00110         exit (1);
00111       }
00112       Vol1Parent_FileName = argv[kar];
00113       brk = YUP;
00114     }
00115     if (!brk && (strcmp(argv[kar], "-sv2")) == 0) {
00116       kar ++;
00117       if (kar >= argc) {
00118         fprintf (SUMA_STDERR, "need argument after -sv2");
00119         exit (1);
00120       }
00121       Vol2Parent_FileName = argv[kar];
00122       brk = YUP;
00123     }
00124     if (!brk && (strcmp(argv[kar], "-prefix")) == 0) {
00125       kar ++;
00126       if (kar >= argc) {
00127         fprintf (SUMA_STDERR, "need argument after -prefix");
00128         exit (1);
00129       }
00130       fout = argv[kar];
00131       brk = YUP;
00132     }
00133     if (!brk && (strcmp(argv[kar], "-hemi")) == 0) {
00134       kar ++;
00135       if (kar >= argc) {
00136         fprintf (SUMA_STDERR, "need argument after -hemi");
00137         exit (1);
00138          }
00139       hemi = argv[kar];
00140       brk = YUP;
00141     }
00142     if (!brk && (strcmp(argv[kar], "-spec")) == 0) {
00143       kar ++;
00144       if (kar >= argc) {
00145         fprintf (SUMA_STDERR, "need argument after -spec ");
00146         exit (1);
00147       }
00148       specfilename = argv[kar];
00149       brk = YUP;
00150     }
00151     
00152     if (!brk && (strcmp(argv[kar], "-onenode")) == 0) {
00153       if (Partial) {
00154          fprintf (SUMA_STDERR, "-onenode is incompatible with -noderange");
00155               exit (1);
00156       }
00157       kar ++;
00158       if (kar >= argc) {
00159         fprintf (SUMA_STDERR, "need argument after -onenode");
00160         exit (1);
00161            }
00162       istart = atoi(argv[kar]);
00163       istop = istart;
00164       Partial = YUP;
00165       brk = YUP;
00166     }
00167     
00168     if (!brk && (strcmp(argv[kar], "-noderange")) == 0) {
00169       if (Partial) {
00170          fprintf (SUMA_STDERR, "-noderange is incompatible with -onenode");
00171               exit (1);
00172       }
00173       kar ++;
00174       if (kar+1 >= argc) {
00175         fprintf (SUMA_STDERR, "need 2 arguments after -noderange");
00176         exit (1);
00177            }
00178       istart = atoi(argv[kar]); kar ++;
00179       istop = atoi(argv[kar]);
00180       Partial = YUP;
00181       brk = YUP;
00182     }
00183     
00184     
00185     if (!brk) {
00186       fprintf (SUMA_STDERR,"Error %s: Option %s not understood. Try -help for usage\n", FuncName, argv[kar]);
00187       exit (1);
00188     } 
00189     else {      
00190       brk = NOPE;
00191       kar ++;
00192     }
00193   }/* loop across command line options */
00194   
00195 
00196   /* allocate for the surface objects */
00197   Surf1 = (SUMA_SurfaceObject *) SUMA_malloc(sizeof(SUMA_SurfaceObject));       
00198   SO1 = (SUMA_SurfaceObject *) SUMA_malloc(sizeof(SUMA_SurfaceObject)); 
00199 
00200         
00201   if (specfilename == NULL) {
00202     fprintf (SUMA_STDERR,"Error %s: No spec filename specified.\n", FuncName);
00203     exit(1);
00204   }
00205   if (!SUMA_Read_SpecFile (specfilename, &Spec)) {
00206     fprintf(SUMA_STDERR,"Error %s: Error in SUMA_Read_SpecFile\n", FuncName);
00207     exit(1);
00208   }     
00209 
00210   /**** loading the first surface *****/
00211   if (SUMA_iswordin(Spec.SurfaceType[0], "FreeSurfer") == 1) {
00212     Surf1_FileName = Spec.SurfaceFile[0];
00213     Surf1 = SUMA_Load_Surface_Object(Surf1_FileName, SUMA_FREE_SURFER, SUMA_ASCII, Vol1Parent_FileName);
00214     tag1 =  "FS";
00215   }
00216   else 
00217     if (SUMA_iswordin(Spec.SurfaceType[0], "SureFit") == 1) {
00218       Surf1_SFName = SUMA_malloc(sizeof(SUMA_SFname));
00219       strcpy(Surf1_SFName->name_coord,Spec.CoordFile[0]);
00220       strcpy(Surf1_SFName->name_topo, Spec.TopoFile[0]);
00221       strcpy(Surf1_SFName->name_param, Spec.SureFitVolParam[0]);
00222       Surf1 = SUMA_Load_Surface_Object(Surf1_SFName, SUMA_SUREFIT, SUMA_ASCII,Vol1Parent_FileName);
00223       tag1 =  "SF";
00224     }
00225   state1 = Spec.State[0];
00226   
00227   /*** loading the second surface ****/
00228   fprintf(SUMA_STDERR,"loading the next surface \n");
00229   if (SUMA_iswordin(Spec.SurfaceType[1], "FreeSurfer") == 1) {
00230     Surf2_FileName = Spec.SurfaceFile[1];
00231     Surf2 = SUMA_Load_Surface_Object(Surf2_FileName, SUMA_FREE_SURFER, SUMA_ASCII, Vol2Parent_FileName);
00232     tag2 =  "FS";
00233   }
00234   else 
00235     if (SUMA_iswordin(Spec.SurfaceType[1], "SureFit") == 1) {
00236       Surf2_SFName = SUMA_malloc(sizeof(SUMA_SFname));
00237       strcpy(Surf2_SFName->name_coord, Spec.CoordFile[1]);
00238       strcpy(Surf2_SFName->name_topo, Spec.TopoFile[1]);
00239       strcpy(Surf2_SFName->name_param, Spec.SureFitVolParam[1]);
00240       Surf2 = SUMA_Load_Surface_Object(Surf2_SFName, SUMA_SUREFIT, SUMA_ASCII,Vol2Parent_FileName);
00241       tag2 =  "SF";
00242     }
00243   state2 = Spec.State[1];
00244 
00245   /* check on the output filename */
00246   if (fout == NULL) { /* form default name */
00247     fname = (char *)SUMA_malloc (
00248             ( strlen(hemi) + strlen(tag1) + strlen(state1) + strlen(tag2) + strlen(state2) + 20 ) *
00249             sizeof(char)); 
00250     sprintf(fname, "%s_%s_%s_%s_%s",hemi,tag1,state1,tag2,state2);
00251 
00252   } else {
00253      fname = SUMA_copy_string(fout); 
00254   }
00255   
00256     sprintf (colorfilename, "%s.col", fname);
00257     sprintf (distancefilename, "%s.dist", fname);
00258    
00259   /* section that will create the files that contain 6 values corresponding to two points to draw segments in SUMA */
00260 
00261   sprintf(segmentfilename1, "%s.allsegs.txt",fname);
00262   sprintf(segmentfilename2, "%s.longsegs.txt",fname);
00263   sprintf(segmentfilename3, "%s.badsegs.txt",fname); 
00264   sprintf(trianglesfilename, "%s.triangles.txt",fname); 
00265 
00266  if((segfile1 = fopen(segmentfilename1, "w"))==NULL) {
00267     fprintf(SUMA_STDERR, "Could not open segment file 1.\n");
00268     exit(1);
00269   }
00270  if((segfile2 = fopen(segmentfilename2, "w"))==NULL) {
00271     fprintf(SUMA_STDERR, "Could not open segment file 2.\n");
00272     exit(1);
00273   }
00274  if((segfile3 = fopen(segmentfilename3, "w"))==NULL) {
00275     fprintf(SUMA_STDERR, "Could not open segment file 3.\n");
00276     exit(1);
00277   }
00278 
00279  if((trianglesfile = fopen(trianglesfilename, "w"))==NULL) {
00280     fprintf(SUMA_STDERR, "Could not open triangles file.\n");
00281     exit(1);
00282   }
00283 
00284 
00285   if (SUMA_filexists(colorfilename) || SUMA_filexists(distancefilename)) {
00286     fprintf (SUMA_STDERR,"Error %s: One or both of output files %s, %s exists.\nWill not overwrite.\n", \
00287              FuncName, distancefilename, colorfilename);
00288     exit(1);
00289   }
00290          
00291   /***********************************************************************************************/
00292 
00293   SO1 = Surf1;
00294   SO2 = Surf2;
00295   
00296 
00297   if (SkipConsistent) {
00298    fprintf (SUMA_STDERR,"Skipping consistency check.\n");
00299   } else {
00300      if (!SO1->EL) SO1->EL = SUMA_Make_Edge_List (SO1->FaceSetList, SO1->N_FaceSet, SO1->N_Node,SO1->NodeList, SO1->idcode_str); 
00301      if (SUMA_MakeConsistent (SO1->FaceSetList, SO1->N_FaceSet, SO1->EL, 1, &trouble) == YUP)
00302        fprintf(SUMA_STDERR,"faces are consistent\n");
00303      else
00304        fprintf(SUMA_STDERR,"faces are not consistent\n");
00305   }
00306  
00307   if (!SO2->EL) SO2->EL = SUMA_Make_Edge_List (SO2->FaceSetList, SO2->N_FaceSet, SO2->N_Node, SO2->NodeList, SO2->idcode_str); 
00308   if (SUMA_MakeConsistent (SO2->FaceSetList, SO2->N_FaceSet, SO2->EL, 1, &trouble) == YUP)
00309     fprintf(SUMA_STDERR,"faces are consistent\n");
00310   else
00311     fprintf(SUMA_STDERR,"faces are not consistent\n");
00312    
00313 
00314    num_nodes1 = SO1->N_Node;
00315    num_nodes2 = SO2->N_Node;
00316    
00317    fprintf(SUMA_STDERR, "Number of nodes in surface 1: %d \n", num_nodes1);
00318    fprintf(SUMA_STDERR, "Number of nodes in surface 2: %d \n", num_nodes2);
00319    fprintf(SUMA_STDERR, "Number of faces in surface 1: %d \n", SO1->N_FaceSet);
00320    fprintf(SUMA_STDERR, "Number of faces in surface 2: %d \n", SO2->N_FaceSet);
00321    if (!SO1->NodeNormList) SUMA_RECOMPUTE_NORMALS(SO1);
00322    if (!SO2->NodeNormList) SUMA_RECOMPUTE_NORMALS(SO2);
00323    
00324    
00325    /* add some noise to surface 2 */
00326    /*  for (i = 0; i < SO2->N_Node; i ++)
00327      {
00328        id = SO2->NodeDim*i;
00329        SO2->NodeList[id] += 0.0002;
00330        SO2->NodeList[id+1] += 0.0003;
00331        SO2->NodeList[id+2] -= 0.0004;
00332      }
00333      
00334    */
00335     
00336   /* Take each node in the SO1-> Nodelist and its corresponding SO1->NodeNormList.  This is the normalized normal vector to the node on the first surface.
00337      So each node is P0.  P1 is computed as some point along the normal vector to that node.  Lets say P0 is P1 + 20 mm along the normal. 
00338      Now feed P0 and P1 into the intersect triangle routine and feed the node and face list of surface 2 */
00339   
00340   distance = SUMA_malloc(num_nodes1*sizeof(float));
00341   /* *****YOU SHOULD ALLOCATE FOR triangle that is done in the function****** triangle = SUMA_malloc(sizeof(SUMA_MT_INTERSECT_TRIANGLE)); */
00342   /* for each node on the first surface do the following */
00343   SUMA_etime (&tt, 0);
00344   
00345    if (!Partial){
00346       istart = 0;
00347       istop = SO1->N_Node-1;
00348    } else {
00349       if (istart > istop) {
00350          fprintf (SUMA_STDERR,"Error %s: starting node %d > stopping node %d\n", 
00351            FuncName, istart, istop);
00352          exit(1);
00353       }
00354       if (istart < 0 || istop > SO1->N_Node-1) {
00355          fprintf (SUMA_STDERR,"Error %s: starting node %d is either < 0 or stopping node > %d (N_Node -1)\n", 
00356            FuncName, onenode, SO1->N_Node-1);
00357          exit(1);
00358       }
00359    }
00360   
00361   FailedDistance = 0; 
00362   for (i = istart; i <= istop; i++) {
00363     id = SO1->NodeDim * i;
00364     P0[0] = SO1->NodeList[id];
00365     P0[1] = SO1->NodeList[id+1];
00366     P0[2] = SO1->NodeList[id+2];
00367 
00368     N0[0] = SO1->NodeNormList[id];
00369     N0[1] = SO1->NodeNormList[id+1];
00370     N0[2] = SO1->NodeNormList[id+2];
00371 
00372    SUMA_POINT_AT_DISTANCE(N0, P0, 100, Points);
00373    P1[0] = Points[0][0];
00374    P1[1] = Points[0][1];
00375    P1[2] = Points[0][2];
00376    P2[0] = Points[1][0];
00377    P2[1] = Points[1][1];
00378    P2[2] = Points[1][2];
00379 
00380    /* now determine the distance along normal */
00381    triangle = SUMA_MT_intersect_triangle(P0,P1, SO2->NodeList, SO2->N_Node, SO2->FaceSetList, SO2->N_FaceSet, triangle);
00382    /* fprintf(SUMA_STDERR,"number of hits for node %d : %d\n", i,triangle->N_hits); */ 
00383    if (triangle->N_hits ==0) {
00384    fprintf(SUMA_STDERR, "Could not find hit for node %d in either direction.\n", i);
00385    fprintf(segfile3,"%f %f %f %f %f %f\n",P0[0],P0[1],P0[2],P1[0],P1[1],P1[2]);
00386    distance[i] = 0.0;
00387    }
00388    else {
00389    fprintf(trianglesfile,"distance for surf 1 node %d:\n",i);
00390    for (k = 0; k < triangle->N_el; k++) {
00391    if (triangle->isHit[k] == YUP)
00392    fprintf(trianglesfile, "hit %d: %f (%f, %f)\n",k,triangle->t[k], triangle->u[k], triangle->v[k]);
00393    }
00394    /* distance[i] = sqrtf(pow(triangle->P[0]-P0[0],2)+pow(triangle->P[1]-P0[1],2)+pow(triangle->P[2]-P0[2],2)); */
00395    distance[i] = triangle->t[triangle->ifacemin];
00396    fprintf(segfile2,"%f %f %f %f %f %f\n",P0[0],P0[1],P0[2],P1[0],P1[1],P1[2]);
00397    }
00398 
00399    if (!KeepMTI) triangle = SUMA_Free_MT_intersect_triangle(triangle); 
00400     
00401 
00402    fprintf(segfile1,"%f %f %f %f %f %f\n",P0[0],P0[1],P0[2],P1[0],P1[1],P1[2]);
00403 
00404     if (!(i%100)) {
00405       delta_t = SUMA_etime(&tt, 1);
00406       fprintf (SUMA_STDERR, " [%d]/[%d] %.2f/100%% completed. Dt = %.2f min done of %.2f min total\r" ,  i, num_nodes1, (float)i / num_nodes1 * 100, delta_t/60, delta_t/i * num_nodes1/60);
00407     }
00408    
00409    if (Partial) {
00410       /* output distance to screen */
00411       fprintf(SUMA_STDERR, "\nDistance at node %d is %f\n", i, distance[i]);
00412    } 
00413         
00414   }
00415   
00416 
00417  
00418   
00419      /* write out the distance file */
00420      if((distancefile = fopen(distancefilename, "w"))==NULL) {
00421        fprintf(SUMA_STDERR, "Could not open file distance.txt.\n");
00422        exit(1);
00423      }
00424      else {  
00425        for (i=0; i < num_nodes1; ++i) {
00426          fprintf (distancefile,"%d\t%f\n", i, distance[i]);
00427        }
00428        fclose (distancefile);
00429      }
00430 
00431      /* output this distance as a color file */
00432      MyColMap = SUMA_GetStandardMap(SUMA_CMAP_MATLAB_DEF_BYR64);
00433      MyOpt = SUMA_ScaleToMapOptInit();
00434      MySV = SUMA_Create_ColorScaledVect(num_nodes1);
00435      mindistance = minimum(num_nodes1, distance);
00436      maxdistance = maximum(num_nodes1, distance);
00437      SUMA_ScaleToMap(distance,num_nodes1,mindistance, maxdistance, MyColMap,MyOpt,MySV);
00438 
00439 
00440      /* write out the distance color file */
00441      if((colorfile = fopen(colorfilename, "w"))==NULL) {
00442        fprintf(SUMA_STDERR, "Could not open file distance.col.\n");
00443        exit(1);
00444      }
00445      else {
00446        for (i=0; i < num_nodes1; ++i) {
00447          fprintf (colorfile,"%d\t%f\t%f\t%f\n", i, MySV->cV[3*i  ], MySV->cV[3*i+1], MySV->cV[3*i+2]);
00448        }
00449        fclose (colorfile);
00450      }
00451 
00452   if (fname) SUMA_free(fname);
00453   
00454   if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);  
00455   SUMA_RETURN (1);
00456 }
00457 
00458 /*************************** FUNCTION DEFINITIONS **********************************/
00459 
00460 void cmp_surf_usage ()
00461 {
00462   static char FuncName[]={"cmp_surf_usage"};
00463   char * s = NULL;
00464   s = SUMA_help_basics();
00465   printf ("\n"
00466           "   Usage:    CompareSurfaces \n"
00467           "             -spec <Spec file>\n"
00468           "             -hemi <L or R>\n"
00469           "             -sv1 <volparentaligned1.BRIK>\n"
00470           "             -sv2 <volparentaligned2.BRIK> \n"
00471           "             [-prefix <fileprefix>]\n"
00472           "\n"
00473           "   NOTE: This program is now superseded by SurfToSurf\n"
00474           "\n"
00475           "   This program calculates the distance, at each node in Surface 1 (S1) to Surface 2 (S2)\n"
00476           "   The distances are computed along the local surface normal at each node in S1.\n"
00477           "   S1 and S2 are the first and second surfaces encountered in the spec file, respectively.\n"
00478           "\n"
00479           "   -spec <Spec file>: File containing surface specification. This file is typically \n"
00480           "                      generated by @SUMA_Make_Spec_FS (for FreeSurfer surfaces) or \n"
00481           "                      @SUMA_Make_Spec_SF (for SureFit surfaces).\n"
00482           "   -hemi <left or right>: specify the hemisphere being processed \n"
00483           "   -sv1 <volume parent BRIK>:volume parent BRIK for first surface \n"
00484           "   -sv2 <volume parent BRIK>:volume parent BRIK for second surface \n"
00485           "\n"
00486           "Optional parameters:\n"
00487           "   [-prefix <fileprefix>]: Prefix for distance and node color output files.\n"
00488           "                           Existing file will not be overwritten.\n"
00489           "   [-onenode <index>]: output results for node index only. \n"
00490           "                       This option is for debugging.\n"
00491           "   [-noderange <istart> <istop>]: output results from node istart to node istop only. \n"
00492           "                                  This option is for debugging.\n"
00493           "   NOTE: -noderange and -onenode are mutually exclusive\n"
00494           "   [-nocons]: Skip mesh orientation consistency check.\n"
00495           "              This speeds up the start time so it is useful\n"
00496           "              for debugging runs.\n"
00497           "\n"
00498           "%s"
00499           "\n   For more help: http://afni.nimh.nih.gov/ssc/ziad/SUMA/SUMA_doc.htm\n"
00500           "\n"
00501           "\n   If you can't get help here, please get help somewhere.\n", s);
00502   SUMA_free(s); s = NULL;
00503   s = SUMA_New_Additions(0, 1); printf("%s\n", s);SUMA_free(s); s = NULL;
00504   printf ("\n    Shruti Japee LBC/NIMH/NIH shruti@codon.nih.gov Ziad S. Saad SSSC/NIMH/NIH ziad@nih.gov \n\n");
00505   /*
00506     printf ("\n   [-dev]: This option will give access to options that are not well polished for consumption.\n");
00507     printf ("\n           \n");
00508   */
00509   exit (0);
00510 }
00511 
00512 
00513 
00514 float minimum (int N, float *inarray)
00515 {
00516   float min = inarray[0];
00517   int i;
00518   for (i = 1; i < N; i++)
00519     {
00520       if (inarray[i] < min)
00521         min = inarray[i];
00522     }
00523   return min;
00524 }
00525 /* determine the max in a vector of floats */
00526 float maximum (int N, float *inarray)
00527 {
00528   float max = inarray[0];
00529   int i;
00530   for (i = 1; i < N; i++)
00531     {
00532       if (inarray[i] > max)
00533         max = inarray[i];
00534     }
00535   return max;
00536 }