home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / getobj3d / getobj3d.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-09-17  |  12.2 KB  |  459 lines
  1.  
  2. /*
  3.  *  getobj3d.c - 3d version of getobj
  4.  *
  5.  *  written 6/90 by Brian Tierney, Lawrence Berkeley Laboratory
  6.  *
  7.    Usage: getobj3d [-t NN][-b NN][-m NN][-g N][-c NN NN NN][-s]
  8.               [-a][-o][-v][-f fname] < inseq > outseq
  9.  
  10.   Options: [-t NN]  Surface threshold value (default = 50)
  11.            [-b NN]  background level (default = 0)
  12.            [-m NN]  minimum size object to keep (default = 10)
  13.            [-g N]  type of connectivity bridge: default = 0)
  14.                   ( 0 = none, 1 = 2D, 2 = 3D weak, 3 = 3D strong )
  15.            [-s]  computes stats on each object, creating file 'getobj.stat'
  16.            [-a]  find all objects in the image
  17.            [-o]  output a binary mask of selected objects.
  18.            [-v]  verbose mode
  19.            [-c NN NN NN] frame,row,col location of object
  20.            [-f fname ] file with list of seed values
  21.  */
  22.  
  23. #include "getobj.h"
  24.  
  25. /****************************************************************/
  26. main(argc, argv)
  27.     int       argc;
  28.     char    **argv;
  29. {
  30.     int       i, num_seed_vals;
  31.     int      *sx_list, *sy_list, *sz_list;
  32.     struct header hd;        /* hips header */
  33.  
  34.     void      clear_background_grid(), parse_args(), init_grid();
  35.     void      clear_background_image(), locate_surfaces(), show_slices();
  36.     u_char ***alloc_3d_byte_array();
  37.     u_short ***alloc_3d_short_array();
  38.     u_int  ***alloc_3d_int_array();
  39.  
  40.     Progname = strsave(*argv);
  41.     parse_args(argc, argv);
  42.  
  43.     read_header(&hd);
  44.     if (hd.pixel_format != PFSHORT && hd.pixel_format != PFBYTE
  45.     && hd.pixel_format != PFINT)
  46.     perr(HE_MSG,"image must be in int, short or byte format");
  47.  
  48.     ncol = hd.ocols;        /* x axis */
  49.     nrow = hd.orows;        /* y axis */
  50.     nframe = hd.num_frame;    /* z axis */
  51.     nvoxels = ncol * nrow * nframe;  /* total number of voxels */
  52.  
  53.     pix_format = hd.pixel_format;    /* byte, short , or int */
  54.  
  55.     if (nframe == 1) {
  56.     fprintf(stderr, "Error: only 1 frame, use getobj instead. \n\n");
  57.     exit(-1);
  58.     }
  59.     if (stats) {
  60.     if ((df = fopen("getobj.stats", "w")) == NULL) {
  61.         fprintf(stderr, "\n error opening stats output file \n\n");
  62.         exit(-1);
  63.     }
  64.     }
  65.     if (val_file) {
  66.     if ((vf = fopen(valfile, "r")) == NULL) {
  67.         fprintf(stderr, "\n error opening seed value file \n\n");
  68.         exit(-1);
  69.     }
  70.     num_seed_vals = get_value(vf);
  71.     } else
  72.     num_seed_vals = 1;
  73.  
  74.     if ((sx_list = Calloc(num_seed_vals, int)) == NULL)
  75.     perror("calloc error");
  76.     if ((sy_list = Calloc(num_seed_vals, int)) == NULL)
  77.     perror("calloc error");
  78.     if ((sz_list = Calloc(num_seed_vals, int)) == NULL)
  79.     perror("calloc error");
  80.  
  81.     if (val_file) {        /* read seed location file */
  82.     for (i = 0; i < num_seed_vals; i++) {
  83.         sx_list[i] = get_value(vf);
  84.         sy_list[i] = get_value(vf);
  85.         sz_list[i] = get_value(vf);
  86.     }
  87.     } else {
  88.     sx_list[0] = isx;
  89.     sy_list[0] = isy;
  90.     sz_list[0] = isz;
  91.     }
  92.  
  93.     for (i = 0; i < num_seed_vals; i++)    /* check if valid values */
  94.     if (sx_list[i] > ncol || sy_list[i] > nrow || sz_list[i] > nframe) {
  95.         fprintf(stderr, "\n Invalid seed value!(%d,%d,%d) \n\n",
  96.             sx_list[i], sy_list[i], sz_list[i]);
  97.         exit(-1);
  98.     }
  99.     fprintf(stderr, "\n rows: %d,  cols: %d,  frames: %d \n", nrow, ncol, nframe);
  100.  
  101.     if (verbose)
  102.     fprintf(stderr, "\n reading image... \n");
  103.  
  104.     grid = alloc_3d_byte_array(nframe, nrow, ncol);
  105.     if (hd.pixel_format == PFSHORT) {
  106.     sht_image = alloc_3d_short_array(nframe, nrow, ncol);
  107.     read_3d_short_array(stdin, sht_image, nframe, nrow, ncol);
  108.     } else if (hd.pixel_format == PFINT) {
  109.     int_image = alloc_3d_int_array(nframe, nrow, ncol);
  110.     read_3d_int_array(stdin, int_image, nframe, nrow, ncol);
  111.     } else {
  112.     image = alloc_3d_byte_array(nframe, nrow, ncol);
  113.     read_3d_byte_array(stdin, image, nframe, nrow, ncol);
  114.     }
  115.  
  116. /* #define DEBUG */
  117. #ifdef DEBUG
  118.     show_slices();
  119. #endif
  120.  
  121.     if (verbose)
  122.     fprintf(stderr, " initializing grid... \n");
  123.  
  124.     init_grid(0);
  125.  
  126.     if (verbose)
  127.     fprintf(stderr, " locating objects... \n");
  128.  
  129.     stack_size = identify_objects(tval);    /* label objects '1' */
  130.     alloc_stack(stack_size);    /* stack shouldn't be larger than the number
  131.                  * of pixels found by identify_objects */
  132.  
  133.     locate_surfaces();        /* label edges of object '2' */
  134.  
  135.     if (find_all) {
  136.     i = 1;
  137.     while (get_object(i, 0, 0, 0) >= 0)    /* label grid with '3' */
  138.         i++;
  139.     } else {
  140.     for (i = 0; i < num_seed_vals; i++) {    /* process each object */
  141.         (void) get_object(i + 1, sx_list[i], sy_list[i], sz_list[i]);
  142.     }
  143.     }
  144.  
  145.     if (verbose)
  146.     fprintf(stderr, "\n zeroing non-objects... ");
  147.  
  148.     clear_background_grid();    /* grid = 0 if grid != 3 */
  149.  
  150.     clear_background_image(bg);    /* set image to bg if grid != 3 */
  151.  
  152.     if (verbose)
  153.     fprintf(stderr, "\n writing image... ");
  154.  
  155.     if (output_binary_mask) {
  156.     init_header(&hd, "", "", nframe, "", nrow, ncol, PFBYTE,1, "");
  157.     update_header(&hd, argc, argv);
  158.     write_header(&hd);
  159.     write_3d_byte_array(stdout, grid, nframe, nrow, ncol);
  160.     } else {
  161.     update_header(&hd, argc, argv);
  162.     write_header(&hd);
  163.     if (pix_format == PFSHORT)
  164.         write_3d_short_array(stdout, sht_image, nrow, ncol, nframe);
  165.     else if (pix_format == PFINT)
  166.         write_3d_int_array(stdout, int_image, nframe, nrow, ncol);
  167.     else
  168.         write_3d_byte_array(stdout, image, nframe, nrow, ncol);
  169.     }
  170.  
  171.     if (pix_format == PFSHORT)
  172.     free_3d_short_array(sht_image);
  173.     else if (pix_format == PFINT)
  174.     free_3d_int_array((int) int_image);
  175.     else
  176.     free_3d_byte_array(image);
  177.     free_3d_byte_array(grid);
  178.     free((char *) stack);
  179.  
  180.     fprintf(stderr, "\n\n finished. \n\n");
  181.     return (0);
  182. }
  183.  
  184. /************************************************************/
  185. int
  186. get_object(num, sx, sy, sz)    /* locates and labels objects */
  187.     int       num;
  188.     int       sx, sy, sz;
  189. {
  190.     int       rval;
  191.     int       object_fill(), get_seed(), get_seed_guess();
  192.     void      reset_grid(), check_object_size();
  193.  
  194.     if (sx == 0 && sy == 0 && sz == 0)
  195.     rval = get_seed_guess(&sx, &sy, &sz);
  196.     else
  197.     rval = get_seed(&sx, &sy, &sz);
  198.  
  199.     if (rval < 0) {        /* didn't find an object */
  200.     if (find_all) {
  201.         if (verbose)
  202.         fprintf(stderr, "\n No more objects. \n");
  203.     } else
  204.         fprintf(stderr, "\n Object not found. \n");
  205.     return (-1);        /* give up */
  206.     }
  207.     pix_total = 0.;
  208.     pcnt = 0L;
  209.     sp = 0;            /* initialize stack pointer */
  210.     push(-1, -1, -1);        /* null stack */
  211.  
  212.     check_object_size(sx, sy, sz);
  213.     if (pcnt < min_size) {
  214.     if (verbose)
  215.         fprintf(stderr, "   skipping object %d: location %d,%d,%d; size of %d pixels \n",
  216.             num, sx, sy, sz, pcnt);
  217.     return (0);
  218.     } else {
  219.     reset_grid();
  220.     }
  221.  
  222.     pcnt = object_fill(sx, sy, sz);    /* grid = 3 if desired object */
  223.  
  224.     if (verbose)
  225.     fprintf(stderr, " Object: %d, location: (%d,%d,%d);  %d pixels \n",
  226.         num, sx, sy, sz, pcnt);
  227.  
  228.     if (stats) {
  229.     fprintf(stderr, "\n STATS: object %d at (%d,%d,%d) is %d pixels.",
  230.         num, sx, sy, sz, pcnt);
  231.     fprintf(df, "\n STATS: object at (%d,%d,%d) is %d pixels.",
  232.         sx, sy, sz, pcnt);
  233.     fprintf(stderr, "\n STATS: Total pixel intensity: %e,   log10: %f ",
  234.         pix_total, log10((double) pix_total));
  235.     fprintf(df, "\n STATS: Total pixel intensity:  %e,   log10: %f ",
  236.         pix_total, log10((double) pix_total));
  237.     fprintf(stderr, "\n STATS: average pixel value of this object is %.1f\n",
  238.         (double) pix_total / (double) pcnt);
  239.     fprintf(df, "\n STATS: average pixel value of this object is %.1f\n",
  240.         (double) pix_total / (double) pcnt);
  241.     }
  242.     return (0);
  243. }
  244.  
  245. /***************************************************************/
  246. void
  247. check_object_size(c, r, f)    /* recursive routine to check if an object is
  248.                  * big enough */
  249.     int       c, r, f;
  250.  /*
  251.   * Note: this routine doesnt count edge pixels, only pixels labeled '1'. It
  252.   * would require extra memory to store the 1 and the 2's (edges), and just
  253.   * counting the 1's is a good enough to get an idea of the objects size.
  254.   */
  255. {
  256.  
  257.     grid[f][r][c] = 9;
  258.     push(c, r, f);        /* add location to the stack */
  259.  
  260.     pcnt++;
  261.     if (pcnt >= min_size)
  262.     return;
  263.  
  264.     if (c < (ncol - 1))
  265.     if (grid[f][r][c + 1] == 1)
  266.         check_object_size(c + 1, r, f);
  267.  
  268.     if (c > 0)
  269.     if (grid[f][r][c - 1] == 1)
  270.         check_object_size(c - 1, r, f);
  271.  
  272.     if (r < (nrow - 1))
  273.     if (grid[f][r + 1][c] == 1)
  274.         check_object_size(c, r + 1, f);
  275.  
  276.     if (r > 0)
  277.     if (grid[f][r - 1][c] == 1)
  278.         check_object_size(c, r - 1, f);
  279.  
  280.     if (f < (nframe - 1))
  281.     if (grid[f + 1][r][c] == 1)
  282.         check_object_size(c, r, f + 1);
  283.  
  284.     if (f > 0)
  285.     if (grid[f - 1][r][c] == 1)
  286.         check_object_size(c, r, f - 1);
  287.  
  288.     return;
  289. }
  290.  
  291. /********************************************************/
  292. void
  293. reset_grid()
  294. {                /* restore grid to original state before
  295.                  * check_object_size routine */
  296.     int       c, r, f;
  297.  
  298.     pop(&c, &r, &f);
  299.     while (c >= 0) {
  300.     grid[f][r][c] = 1;
  301.     pop(&c, &r, &f);
  302.     }
  303. }
  304.  
  305. /***************************************************************/
  306. void
  307. sum_pixel(c, r, f)        /* used of stats flag is set */
  308.     int       c, r, f;
  309. {
  310.     if (pix_format == PFSHORT)
  311.     pix_total += (double) sht_image[f][r][c];
  312.     else if (pix_format == PFINT)
  313.     pix_total += (double) int_image[f][r][c];
  314.     else
  315.     pix_total += (double) image[f][r][c];
  316. }
  317.  
  318. /***************************************************************/
  319. void
  320. show_slices()
  321. {                /* for debugging  */
  322.     register int c, r, f;
  323.  
  324.     for (f = 0; f < nframe; f++) {
  325.     fprintf(stderr, "\n frame #: %d \n", f);
  326.     for (r = 0; r < nrow; r++) {
  327.         for (c = 0; c < ncol; c++)
  328.         fprintf(stderr, "%3d", image[f][r][c]);
  329.         fprintf(stderr, "\n");
  330.     }
  331.     }
  332. }
  333.  
  334. /******************************************************/
  335. static int
  336. get_value(vfile)        /* read seed values from file */
  337.     FILE     *vfile;
  338. {
  339.     int       val, tmp;
  340.  
  341.     if ((tmp = fscanf(vfile, "%d\n", &val)) == EOF) {
  342.     fprintf(stderr, "Out of values in seed file!\n\n");
  343.     exit(-1);
  344.     } else if (tmp != 1) {
  345.     fprintf(stderr, " Bad seed value file!\n\n");
  346.     exit(-1); 
  347.     }
  348.     return (val);
  349. }
  350.  
  351. /*************************************************/
  352. void
  353. parse_args(argc, argv)
  354.     int       argc;
  355.     char     *argv[];
  356. {
  357.     int       iv;
  358.     void      usageterm();
  359.  
  360.     /* set defaults */
  361.     tval = 50;
  362.     min_size = 10;
  363.     bg = isx = isy = isz = stats = verbose = find_all = 0;
  364.     output_binary_mask = bridges = 0;
  365.     valfile = NULL;
  366.  
  367.     /* Interpret options  */
  368.     while (--argc > 0 && (*++argv)[0] == '-') {
  369.     char     *s;
  370.     for (s = argv[0] + 1; *s; s++)
  371.         switch (*s) {
  372.         case 't':
  373.         if (argc < 2)
  374.             usageterm();
  375.         sscanf(*++argv, "%d", &iv);
  376.         tval = iv;
  377.         argc--;
  378.         break;
  379.         case 'b':
  380.         if (argc < 2)
  381.             usageterm();
  382.         sscanf(*++argv, "%d", &iv);
  383.         bg = iv;
  384.         argc--;
  385.         break;
  386.         case 'g':
  387.         if (argc < 2)
  388.             usageterm();
  389.         sscanf(*++argv, "%d", &iv);
  390.         bridges = iv;
  391.         argc--;
  392.         break;
  393.         case 'm':
  394.         if (argc < 2)
  395.             usageterm();
  396.         sscanf(*++argv, "%d", &iv);
  397.         min_size = iv;
  398.         argc--;
  399.         break;
  400.         case 's':
  401.         stats = 1;
  402.         break;
  403.         case 'a':
  404.         find_all = 1;
  405.         fprintf(stderr, " Looking for all objects \n");
  406.         break;
  407.         case 'o':
  408.         output_binary_mask = 1;
  409.         break;
  410.         case 'v':
  411.         verbose = 1;
  412.         break;
  413.         case 'c':
  414.         if (argc < 4)
  415.             usageterm();
  416.         sscanf(*++argv, "%d", &iv);
  417.         isx = iv;
  418.         sscanf(*++argv, "%d", &iv);
  419.         isy = iv;
  420.         sscanf(*++argv, "%d", &iv);
  421.         isz = iv;
  422.         argc -= 3;
  423.         break;
  424.         case 'f':
  425.         if (argc < 2)
  426.             usageterm();
  427.         valfile = *++argv;
  428.         val_file++;
  429.         argc--;
  430.         break;
  431.         case 'h':
  432.         usageterm();
  433.         break;
  434.         default:
  435.         usageterm();
  436.         break;
  437.         }
  438.     }                /* while */
  439. }
  440.  
  441. /******************************************************/
  442. void
  443. usageterm()
  444. {
  445.     fprintf(stderr, "Usage: getobj3d [-t NN][-b NN][-m NN][-g N][-c NN NN NN][-s][-a][-o][-v][-f fname] < inseq > outseq \n ");
  446.     fprintf(stderr, " Options: [-t NN]  Surface threshold value (default = 50)\n");
  447.     fprintf(stderr, "           [-b NN]  background level (default = 0)\n");
  448.     fprintf(stderr, "           [-m NN]  minimum size object to keep (default = 10)\n");
  449.     fprintf(stderr, "           [-g N]  type of connectivity bridge: default = 0) \n");
  450.     fprintf(stderr, "                  ( 0 = none, 1 = 2D, 2 = 3D weak, 3 = 3D strong ) \n");
  451.     fprintf(stderr, "           [-s]  computes stats on each object, creating file 'getobj.stat' \n");
  452.     fprintf(stderr, "           [-a]  find all objects in the image \n");
  453.     fprintf(stderr, "           [-o]  output a binary mask of selected objects. \n");
  454.     fprintf(stderr, "           [-v]  verbose mode \n");
  455.     fprintf(stderr, "           [-c NN NN NN] col,row,frame location of object \n");
  456.     fprintf(stderr, "           [-f fname ] file with list of seed values \n\n");
  457.     exit(0);
  458. }
  459.