home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / xview / segal / threshold.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-12-24  |  11.4 KB  |  524 lines
  1. /*
  2.  *    threshold.c
  3.  *
  4.  *    Bryan Skene, LBL July 1991
  5.  *    part of 'segal'
  6.  *
  7.  *    Thresholding a grayscale image is pretty easy.  Thresholding a color
  8.  *    image involves a lot more user-input and processing.  To threshold a
  9.  *    color image, we define a closed surface in RGB space, call it S,
  10.  *    which contains the colors represented in some region of interest.
  11.  *    Any voxel with a color in S satisfies the threshold; others don't.  One
  12.  *    alternative approach is to build just some list L of colors in some
  13.  *    representative region of interest (cropped box perhaps) and then test
  14.  *    each voxel's color for membership in L.  I will implement both methods
  15.  *    and see which seems to be most useful.
  16.  *
  17.  *    BUT!  Before I do any of this, I will see what results I get by
  18.  *    simply thresholding on the MONO(r,g,b) value (= intensity).
  19.  */
  20.  
  21. #include "common.h"
  22.  
  23. #define PANELS 32
  24.  
  25. int hsum, hmu;
  26.  
  27. int histo[256];
  28. int num_panels;
  29. byte panel[PANELS];
  30. byte map[256];
  31.  
  32. /***********************************************/
  33. void
  34. refresh_histogram()
  35. {
  36.     void generate_histogram();
  37.     void draw_histogram();
  38.     void draw_thresh_bounds();
  39.     void generate_histo_stats();
  40.  
  41.     generate_histogram();
  42.     draw_histogram();
  43.     draw_thresh_bounds();
  44.     generate_histo_stats();
  45. }
  46.  
  47. /***********************************************/
  48. void
  49. generate_histogram()
  50. {
  51.     void traverse_region();
  52.     void add_to_histogram();
  53.     void scale_histogram();
  54.  
  55.     int i;
  56.  
  57.     for(i = 0; i < 256; i++) {
  58.         histo[i] = 0;
  59.     }
  60.  
  61.     traverse_region(threshold.roi, 0, 0, win[WIN_PAINT].f,
  62.         add_to_histogram);
  63.  
  64.     scale_histogram();
  65. }
  66.  
  67. /***********************************************/
  68. void
  69. add_to_histogram(x, y, f)
  70. int x, y, f;
  71. {
  72.     if(segal.color) {
  73.         if(threshold.plane == VAL_RGB)
  74.             histo[MONO(cbuf[RP][f][y][x], cbuf[GP][f][y][x],
  75.                 cbuf[BP][f][y][x])]++;
  76.         else histo[cbuf[threshold.plane][f][y][x]]++;
  77.     }
  78.     else histo[ibuf[f][y][x]]++;
  79. }
  80.  
  81. /***********************************************/
  82. void
  83. scale_histogram()
  84. {
  85.     char foo[14];
  86.     int i, min, max, max_i, old_max, delta;
  87. #define LEAVE_OUT_ZERO
  88.  
  89. #ifdef LEAVE_OUT_ZERO
  90. #define BEG_VALUE 1
  91. #else
  92. #define BEG_VALUE 0
  93. #endif
  94.  
  95.     /* sum pre-scaled histogram */
  96.     hsum = 0;
  97.     for(i = BEG_VALUE; i < 256; i++) hsum += histo[i];
  98.  
  99.     hmu = hsum / PANELS;
  100.  
  101.     /* everything greater than 4*hmu becomes 4*hmu */
  102.     for(i = BEG_VALUE; i < 256; i++)
  103.         if(histo[i] > 4*hmu) histo[i] = 4 * hmu;
  104.  
  105.     /* scale histogram */
  106.     max_i = 0;
  107.     min = max = old_max = histo[0];
  108.     for(i = BEG_VALUE; i < 256; i++) {
  109.         if(histo[i] < min) min = histo[i];
  110.         if(histo[i] > max) {
  111.             old_max = max;
  112.             max_i = i;
  113.             max = histo[i];
  114.         }
  115.     }
  116.  
  117.     if(max - old_max > old_max) {
  118.         histo[max_i] = old_max;
  119.         max = old_max;
  120.     }
  121.  
  122.     /* label the histogram with the "scaled" max value */
  123.     sprintf(foo, "%d", max);
  124.     xv_set(Threshold_pop_threshold->msg_max_histo,
  125.         PANEL_LABEL_STRING, foo,
  126.         NULL);
  127.  
  128.     delta = max - min;
  129.  
  130.     for(i = 0; i < 256; i++)
  131.         histo[i] = (int) ((float) 256 * (float) histo[i] /
  132.             (float) delta);
  133.  
  134.     /* sum post-scaled histogram */
  135.     hsum = 0;
  136.     for(i = BEG_VALUE; i < 256; i++) hsum += histo[i];
  137.  
  138.     hmu = hsum / PANELS;
  139. }
  140.  
  141. /***********************************************/
  142. void
  143. draw_histogram()
  144. {
  145.     unsigned long standout();
  146.  
  147.     int i;
  148.  
  149.     XClearArea(display, threshold.xid, 0, 0, 0, 0, FALSE);
  150.     XSetForeground(display, gc, standout(TURQUOISE));
  151.     for(i = 0; i < 256; i++) {
  152.         XDrawLine(display, threshold.xid, gc,
  153.             i, 255, i, 255 - histo[i]);
  154.     }
  155.  
  156.     xv_set(Threshold_pop_threshold->msg_mu,
  157.         XV_SHOW, FALSE,
  158.         NULL);
  159. }
  160.  
  161. /***********************************************/
  162. void
  163. erase_thresh_bounds()
  164. {
  165.     XSetForeground(display, gc, standout(CWHITE));
  166.     XDrawLine(display, threshold.xid, gc,
  167.         threshold.min, 0, threshold.min, 255 - histo[threshold.min]);
  168.     XDrawLine(display, threshold.xid, gc,
  169.         threshold.max, 0, threshold.max, 255 - histo[threshold.max]);
  170.  
  171.     XSetForeground(display, gc, standout(TURQUOISE));
  172.     XDrawLine(display, threshold.xid, gc,
  173.         threshold.min, 255 - histo[threshold.min], threshold.min, 255);
  174.     XDrawLine(display, threshold.xid, gc,
  175.         threshold.max, 255 - histo[threshold.max], threshold.max, 255);
  176. }
  177.  
  178. /***********************************************/
  179. void
  180. draw_thresh_bounds()
  181. {
  182.     char foo[20];
  183.     int min_color, max_color;
  184.  
  185.     sprintf(foo, "Max: %d", threshold.max);
  186.     xv_set(Threshold_pop_threshold->msg_max,
  187.         PANEL_LABEL_STRING, foo,
  188.         NULL);
  189.     sprintf(foo, "Min: %d", threshold.min);
  190.     xv_set(Threshold_pop_threshold->msg_min,
  191.         PANEL_LABEL_STRING, foo,
  192.         NULL);
  193.  
  194.     if(threshold.min == threshold.max) {
  195.         min_color = max_color = CBLACK;
  196.     }
  197.     else {
  198.         min_color = RED;
  199.         max_color = GREEN;
  200.     }
  201.  
  202.     XSetForeground(display, gc, standout(min_color));
  203.     XDrawLine(display, threshold.xid, gc,
  204.         threshold.min, 0, threshold.min, 255 - histo[threshold.min]);
  205.     XSetForeground(display, gc, standout(max_color));
  206.     XDrawLine(display, threshold.xid, gc,
  207.         threshold.max, 0, threshold.max, 255 - histo[threshold.max]);
  208.  
  209.     XSetForeground(display, gc, standout(PURPLE));
  210.     XDrawLine(display, threshold.xid, gc,
  211.         threshold.min, 255 - histo[threshold.min], threshold.min, 255);
  212.     XDrawLine(display, threshold.xid, gc,
  213.         threshold.max, 255 - histo[threshold.max], threshold.max, 255);
  214. }
  215.  
  216. /***********************************************/
  217. void
  218. generate_histo_stats()
  219. {
  220. /* it may be useful to replace BEG_VALUE w/ lower thresh and 256 w/ upper */
  221. #define LEAVE_OUT_ZERO
  222.  
  223. #ifdef LEAVE_OUT_ZERO
  224. #define BEG_VALUE 1
  225. #else
  226. #define BEG_VALUE 0
  227. #endif
  228.  
  229.     int i, j, sum, pan_w, sx0, sx1, dx0;
  230.  
  231.     pan_w = 256 / PANELS;
  232.  
  233.     /* build panel[] = array of panel's rightmost slice x-coords */
  234.     num_panels = 0;
  235.     sum = 0;
  236.     i = BEG_VALUE;
  237.     while(i < 256 && num_panels < PANELS) {
  238.         sum += histo[i];
  239.         if(sum > hmu) {
  240.             if(sum-hmu >= histo[i]/2
  241.             && i > 1
  242.             && ((num_panels > 0 && panel[num_panels - 1] != i - 1)
  243.               || num_panels == 0)) {
  244.                 /* give this slice to the next panel */
  245.                 panel[num_panels] = i - 1;
  246.             }
  247.             else {
  248.                 /* give this slice to this panel */
  249.                 panel[num_panels] = i;
  250.                 i++;
  251.             }
  252.             vprint"panel[%2d] = %d\n",
  253.                 num_panels, panel[num_panels]);
  254.  
  255.             num_panels++;
  256.             sum = 0;
  257.         }
  258.         else i++;
  259.     }
  260.  
  261.     /* make mappings for panels */
  262.     for(i = 0; i < num_panels; i++) {
  263.  
  264.         /* source */
  265.         if(i > 0) sx0 = panel[i - 1];
  266.         else sx0 = 0;
  267.  
  268.         sx1 = panel[i];
  269.  
  270.         /* destination */
  271.         dx0 = i * pan_w;
  272.  
  273.         for(j = sx0; j < sx1; j++)
  274.             map[j] = dx0 + (int)(((float)(j - sx0) / (float)(sx1 - sx0)) * (float)pan_w);
  275.     }
  276. }
  277.  
  278. /***********************************************/
  279. void
  280. draw_histo_stats()
  281. {
  282.     int i, y0;
  283.     char foo[20];
  284.  
  285.     /* draw mu */
  286.     XSetForeground(display, gc, standout(CBLACK));
  287.     XDrawLine(display, threshold.xid, gc,
  288.         0, 255 - hmu, 255, 255 - hmu);
  289.  
  290.     y0 = xv_get(Threshold_pop_threshold->canvas, XV_Y, NULL);
  291.     sprintf(foo, "Mu = %d\n", hmu);
  292.     xv_set(Threshold_pop_threshold->msg_mu,
  293.         XV_SHOW, TRUE,
  294.         XV_Y, y0 + 255 - hmu,
  295.         PANEL_LABEL_STRING, foo,
  296.         NULL);
  297.  
  298.     /* draw panels */
  299.     XSetForeground(display, gc, standout(ORANGE));
  300.     for(i = 0; i < num_panels; i++)
  301.         XDrawLine(display, threshold.xid, gc,
  302.             panel[i], 255, panel[i], 255 - histo[panel[i]]);
  303. }
  304.  
  305. /***********************************************/
  306. void
  307. histoeq()
  308. {
  309.     void save_image_undo();
  310.     void eq_i_pixel();
  311.     void refresh_histogram();
  312.     void redisplay_all();
  313.  
  314.     if(threshold.roi == R2d_WHOLE
  315.     || threshold.roi == R2d_CROP
  316.     || threshold.roi == R2d_PT_LIST)
  317.         save_image_undo(WIN_PAINT);
  318.  
  319.     /* map the image to it's new values */
  320.     traverse_region(threshold.roi, 0, 0, win[WIN_PAINT].f, eq_i_pixel);
  321.     img.changed_frame = TRUE;
  322.     refresh_histogram();
  323.     redisplay_all();
  324. }
  325.  
  326. /***********************************************/
  327. void
  328. eq_i_pixel(x, y, f)
  329. int x, y, f;
  330. {
  331.     if(segal.color) {
  332.         if(threshold.plane == VAL_RGB) {
  333.             cbuf[RP][f][y][x] += (11./32.)
  334.                 * (map[cbuf[RP][f][y][x]] - cbuf[RP][f][y][x]);
  335.             cbuf[GP][f][y][x] += (16./32.)
  336.                 * (map[cbuf[GP][f][y][x]] - cbuf[RP][f][y][x]);
  337.             cbuf[BP][f][y][x] += (5./32.)
  338.                 * (map[cbuf[BP][f][y][x]] - cbuf[RP][f][y][x]);
  339.         }
  340.     else cbuf[threshold.plane][f][y][x]
  341.         = map[cbuf[threshold.plane][f][y][x]];
  342.     }
  343.     else ibuf[f][y][x] = map[ibuf[f][y][x]];
  344. }
  345.  
  346. /***********************************************/
  347. void
  348. threshold_mask()
  349. {
  350.     void traverse_region();
  351.     void thresh_m_pixel();
  352.     void redisplay_all();
  353.  
  354.     traverse_region(threshold.roi, 0, 0, win[WIN_PAINT].f, thresh_m_pixel);
  355.     m[segal.e_m].changed_frame = TRUE;
  356.     redisplay_all();
  357. }
  358.  
  359. /***********************************************/
  360. void
  361. thresh_m_pixel(x, y, f)
  362. int x, y, f;
  363. {
  364.     int pval;
  365.  
  366.     if(segal.color) {
  367.         if(threshold.plane == VAL_RGB)
  368.             pval = MONO(cbuf[RP][f][y][x], cbuf[GP][f][y][x],
  369.                 cbuf[BP][f][y][x]);
  370.         else pval = cbuf[threshold.plane][f][y][x];
  371.     }
  372.     else pval = ibuf[f][y][x];
  373.  
  374.     if(pval >= threshold.min
  375.     && pval <= threshold.max)
  376.         switch(threshold.mask_effect) {
  377.         case THRESH_OVERWRITE :
  378.         case THRESH_ADD_TO :
  379.             TURN_BIT_ON(mbuf[f][y][x], m[segal.e_m].bit_key)
  380.             break;
  381.         case THRESH_REMOVE_FROM :
  382.             TURN_BIT_OFF(mbuf[f][y][x], m[segal.e_m].bit_key)
  383.             break;
  384.         default :
  385.             break;
  386.         }
  387.     else switch(threshold.mask_effect) {
  388.         case THRESH_OVERWRITE :
  389.             TURN_BIT_OFF(mbuf[f][y][x], m[segal.e_m].bit_key)
  390.             break;
  391.         case THRESH_REMOVE_FROM :
  392.         case THRESH_ADD_TO :
  393.         default :
  394.             break;
  395.     }
  396. }
  397.  
  398. /***********************************************/
  399. void
  400. threshold_image()
  401. {
  402.     void save_image_undo();
  403.     void traverse_region();
  404.     void thresh_i_pixel();
  405.     void refresh_histogram();
  406.     void redisplay_all();
  407.  
  408.     if(threshold.roi == R2d_WHOLE
  409.     || threshold.roi == R2d_CROP
  410.     || threshold.roi == R2d_PT_LIST)
  411.         save_image_undo(WIN_PAINT);
  412.  
  413.     traverse_region(threshold.roi, 0, 0, win[WIN_PAINT].f, thresh_i_pixel);
  414.     img.changed_frame = TRUE;
  415.     refresh_histogram();
  416.     redisplay_all();
  417. }
  418.  
  419. /***********************************************/
  420. void
  421. thresh_i_pixel(x, y, f)
  422. int x, y, f;
  423. {
  424. #define BIND_PVAL(p,d) {pval=p+d; RANGE(pval,0,255) p=pval;}
  425.     int pval;
  426.  
  427.     if(segal.color) {
  428.         if(threshold.plane == VAL_RGB)
  429.             pval = MONO(cbuf[RP][f][y][x], cbuf[GP][f][y][x],
  430.                 cbuf[BP][f][y][x]);
  431.         else pval = cbuf[threshold.plane][f][y][x];
  432.     }
  433.     else pval = ibuf[f][y][x];
  434.  
  435.     if(pval >= threshold.min
  436.     && pval <= threshold.max) {
  437.         if(segal.color) {
  438.             if(threshold.plane == VAL_RGB) {
  439.                 BIND_PVAL(cbuf[RP][f][y][x],
  440.                     threshold.image_effect)
  441.                 BIND_PVAL(cbuf[GP][f][y][x],
  442.                     threshold.image_effect)
  443.                 BIND_PVAL(cbuf[BP][f][y][x],
  444.                     threshold.image_effect)
  445.             }
  446.             else BIND_PVAL(cbuf[threshold.plane][f][y][x],
  447.                 threshold.image_effect)
  448.         }
  449.         else BIND_PVAL(ibuf[f][y][x], threshold.image_effect)
  450.     }
  451. #undef BIND_PVAL
  452. }
  453.  
  454. /***********************************************/
  455. void
  456. thresh_event(event)
  457. Event *event;
  458. {
  459. /* mouse events:
  460.  *    left button  : change min
  461.  *    right button : change max
  462.  */
  463.     void erase_thresh_bounds();
  464.     void draw_thresh_bounds();
  465.  
  466.     int x;
  467.     static int which_button;
  468.  
  469.     switch(event_id(event)) {
  470.     case MS_LEFT :
  471.         if(event_is_down(event)) {
  472.             which_button = MS_LEFT;
  473.             x = event_x(event);
  474.             if(x <= threshold.max
  475.             && x >= 0) {
  476.                 erase_thresh_bounds();
  477.                 threshold.min = x;
  478.                 draw_thresh_bounds();
  479.             }
  480.         }
  481.         break;
  482.     case MS_RIGHT :
  483.         if(event_is_down(event)) {
  484.             which_button = MS_RIGHT;
  485.             x = event_x(event);
  486.             if(x <= 255
  487.             && x >= threshold.min) {
  488.                 erase_thresh_bounds();
  489.                 threshold.max = x;
  490.                 draw_thresh_bounds();
  491.             }
  492.         }
  493.         break;
  494.     case MS_MIDDLE :
  495.         break;
  496.     case LOC_DRAG :
  497.         if(event_is_down(event)) switch(which_button) {
  498.         case MS_LEFT :
  499.             x = event_x(event);
  500.             if(x <= threshold.max
  501.             && x >= 0) {
  502.                 erase_thresh_bounds();
  503.                 threshold.min = x;
  504.                 draw_thresh_bounds();
  505.             }
  506.             break;
  507.         case MS_RIGHT :
  508.             x = event_x(event);
  509.             if(x <= 255
  510.             && x >= threshold.min) {
  511.                 erase_thresh_bounds();
  512.                 threshold.max = x;
  513.                 draw_thresh_bounds();
  514.             }
  515.             break;
  516.         default :
  517.             break;
  518.         }
  519.         break;
  520.     default :
  521.         break;
  522.     }
  523. }
  524.