home *** CD-ROM | disk | FTP | other *** search
/ Photo CD Demo 1 / Demo.bin / gems / gemsii / hot.c < prev    next >
C/C++ Source or Header  |  1991-04-21  |  11KB  |  396 lines
  1. /*
  2.  * hot.c - Scan an image for pixels with RGB values that will give
  3.  *    "unsafe" values of chrominance signal or composite signal
  4.  *    amplitude when encoded into an NTSC or PAL colour signal.
  5.  *    (This happens for certain high-intensity high-saturation colours
  6.  *    that are rare in real scenes, but can easily be present
  7.  *    in synthetic images.)
  8.  *
  9.  *     Such pixels can be flagged so the user may then choose other
  10.  *    colours.  Or, the offending pixels can be made "safe"
  11.  *    in a manner that preserves hue.
  12.  *
  13.  *    There are two reasonable ways to make a pixel "safe":
  14.  *    We can reduce its intensity (luminance) while leaving
  15.  *    hue and saturation the same.  Or, we can reduce saturation
  16.  *    while leaving hue and luminance the same.  A #define selects
  17.  *    which strategy to use.
  18.  *
  19.  * Note to the user: You must add your own read_pixel() and write_pixel()
  20.  *    routines.  You may have to modify pix_decode() and pix_encode().
  21.  *    MAXPIX, WID, and HGT are likely to need modification.
  22.  */
  23.  
  24. /*
  25.  * Originally written as "ikNTSC.c" by Alan Wm Paeth,
  26.  *    University of Waterloo, August, 1985
  27.  * Updated by Dave Martindale, Imax Systems Corp., December 1990
  28.  */
  29.  
  30. /*
  31.  * Compile-time options.
  32.  *
  33.  * Define either NTSC or PAL as 1 to select the colour system.
  34.  * Define the other one as zero, or leave it undefined.
  35.  *
  36.  * Define FLAG_HOT as 1 if you want "hot" pixels set to black
  37.  *    to identify them.  Otherwise they will be made safe.
  38.  *
  39.  * Define REDUCE_SAT as 1 if you want hot pixels to be repaired by
  40.  *    reducing their saturation.  By default, luminance is reduced.
  41.  *
  42.  * CHROMA_LIM is the limit (in IRE units) of the overall
  43.  *    chrominance amplitude; it should be 50 or perhaps
  44.  *    very slightly higher.
  45.  * 
  46.  * COMPOS_LIM is the maximum amplitude (in IRE units) allowed for
  47.  *    the composite signal.  A value of 100 is the maximum
  48.  *    monochrome white, and is always safe.  120 is the absolute
  49.  *    limit for NTSC broadcasting, since the transmitter's carrier
  50.  *    goes to zero with 120 IRE input signal.  Generally, 110
  51.  *    is a good compromise - it allows somewhat brighter colours
  52.  *    than 100, while staying safely away from the hard limit.
  53.  */
  54.  
  55. #define    NTSC        1
  56. #define    PAL         0
  57. #define    FLAG_HOT    0
  58. #define    REDUCE_SAT  0
  59.  
  60.  
  61. #define    CHROMA_LIM      50.0        /* chroma amplitude limit */
  62. #define    COMPOS_LIM      110.0        /* max IRE amplitude */
  63.  
  64. #if NTSC
  65. /*
  66.  * RGB to YIQ encoding matrix.
  67.  */
  68. double code_matrix[3][3] = {
  69.      0.2989,     0.5866,     0.1144,
  70.      0.5959,    -0.2741,    -0.3218,
  71.      0.2113,    -0.5227,     0.3113,
  72. };
  73.  
  74. #define    PEDESTAL    7.5        /* 7.5 IRE black pedestal */
  75. #define    GAMMA        2.2
  76. #endif NTSC
  77.  
  78. #if PAL
  79. /*
  80.  * RGB to YUV encoding matrix.
  81.  */
  82. double code_matrix[3][3] = {
  83.      0.2989,     0.5866,     0.1144,
  84.     -0.1473,    -0.2891,     0.4364,
  85.      0.6149,    -0.5145,    -0.1004,
  86. };
  87.  
  88. #define    PEDESTAL    0.0        /* no pedestal in PAL */
  89. #define    GAMMA        2.8
  90. #endif PAL
  91.  
  92.  
  93. #define SCALE    8192            /* scale factor: do floats with int math */
  94. #define MAXPIX     255            /* white value */
  95. #define WID     1024           /* FB dimensions */
  96. #define HGT     768
  97.  
  98. typedef struct {
  99.     unsigned char    r, g, b;
  100. } Pixel;
  101.  
  102. int    tab[3][3][MAXPIX+1];    /* multiply lookup table */
  103. double    chroma_lim;             /* chroma limit */
  104. double    compos_lim;             /* composite amplitude limit */
  105. long    ichroma_lim2;           /* chroma limit squared (scaled integer) */
  106. int    icompos_lim;            /* composite amplitude limit (scaled integer) */
  107.  
  108. double    pix_decode(), gc(), inv_gc();
  109. int    pix_encode(), hot();
  110.  
  111. main()
  112. {
  113.     Pixel    p;
  114.     int    row, col;
  115.  
  116.     build_tab();
  117.  
  118.     for (col=0; col<WID; col++) {
  119.         for(row=0; row<HGT; row++) {
  120.             read_pixel(row, col, &p);
  121.             if (hot(&p)) {
  122.                 write_pixel(row, col, &p);
  123.             }
  124.         }
  125.     }
  126. }
  127.  
  128. /*
  129.  * build_tab: Build multiply lookup table.
  130.  *
  131.  * For each possible pixel value, decode value into floating-point
  132.  * intensity.  Then do gamma correction required by the video
  133.  * standard.  Scale the result by our fixed-point scale factor.
  134.  * Then calculate 9 lookup table entries for this pixel value.
  135.  *
  136.  * We also calculate floating-point and scaled integer versions
  137.  * of our limits here.  This prevents evaluating expressions every pixel
  138.  * when the compiler is too stupid to evaluate constant-valued
  139.  * floating-point expressions at compile time.
  140.  *
  141.  * For convenience, the limits are #defined using IRE units.
  142.  * We must convert them here into the units in which YIQ
  143.  * are measured.  The conversion from IRE to internal units
  144.  * depends on the pedestal level in use, since as Y goes from
  145.  * 0 to 1, the signal goes from the pedestal level to 100 IRE.
  146.  * Chroma is always scaled to remain consistent with Y.
  147.  */
  148.  
  149. build_tab()
  150. {
  151.     register double    f;
  152.     register int    pv;
  153.  
  154.     for (pv = 0; pv <= MAXPIX; pv++) {
  155.         f = SCALE * gc(pix_decode(pv));
  156.         tab[0][0][pv] = (int)(f * code_matrix[0][0] + 0.5);
  157.         tab[0][1][pv] = (int)(f * code_matrix[0][1] + 0.5);
  158.         tab[0][2][pv] = (int)(f * code_matrix[0][2] + 0.5);
  159.         tab[1][0][pv] = (int)(f * code_matrix[1][0] + 0.5);
  160.         tab[1][1][pv] = (int)(f * code_matrix[1][1] + 0.5);
  161.         tab[1][2][pv] = (int)(f * code_matrix[1][2] + 0.5);
  162.         tab[2][0][pv] = (int)(f * code_matrix[2][0] + 0.5);
  163.         tab[2][1][pv] = (int)(f * code_matrix[2][1] + 0.5);
  164.         tab[2][2][pv] = (int)(f * code_matrix[2][2] + 0.5);
  165.     }
  166.  
  167.     chroma_lim = (double)CHROMA_LIM / (100.0 - PEDESTAL);
  168.     compos_lim = ((double)COMPOS_LIM - PEDESTAL) / (100.0 - PEDESTAL);
  169.  
  170.     ichroma_lim2 = (int)(chroma_lim * SCALE + 0.5);
  171.     ichroma_lim2 *= ichroma_lim2;
  172.     icompos_lim = (int)(compos_lim * SCALE + 0.5);
  173. }
  174.  
  175. int
  176. hot(p)
  177. Pixel    *p;
  178. {
  179.     register int    r, g, b;
  180.     register int    y, i, q;
  181.     register long    y2, c2;
  182.     double        pr, pg, pb;
  183. #if REDUCE_SAT
  184.     double        py;
  185. #endif
  186.     register double    fy, fc, t, scale;
  187. #if !FLAG_HOT
  188.     static int    prev_r = 0, prev_g = 0, prev_b = 0;
  189.     static int    new_r, new_g, new_b;
  190. #endif
  191.     extern double    pow(), hypot();
  192.  
  193.     r = p->r;
  194.     g = p->g;
  195.     b = p->b;
  196.  
  197.     /*
  198.      * Pixel decoding, gamma correction, and matrix multiplication
  199.      * all done by lookup table.
  200.      *
  201.      * "i" and "q" are the two chrominance components;
  202.      * they are I and Q for NTSC.
  203.      * For PAL, "i" is U (scaled B-Y) and "q" is V (scaled R-Y).
  204.      * Since we only care about the length of the chroma vector,
  205.      * not its angle, we don't care which is which.
  206.      */
  207.     y = tab[0][0][r] + tab[0][1][g] + tab[0][2][b];
  208.     i = tab[1][0][r] + tab[1][1][g] + tab[1][2][b];
  209.     q = tab[2][0][r] + tab[2][1][g] + tab[2][2][b];
  210.  
  211.     /*
  212.      * Check to see if the chrominance vector is too long or the
  213.      * composite waveform amplitude is too large.
  214.      *
  215.      * Chrominance is too large if
  216.      *
  217.      *    sqrt(i^2, q^2)  >  chroma_lim.
  218.      *
  219.      * The composite signal amplitude is too large if
  220.      *
  221.      *    y + sqrt(i^2, q^2)  >  compos_lim.
  222.      *
  223.      * We avoid doing the sqrt by checking
  224.      *
  225.      *    i^2 + q^2  >  chroma_lim^2
  226.      * and
  227.      *    y + sqrt(i^2 + q^2)  >  compos_lim
  228.      *    sqrt(i^2 + q^2)  >  compos_lim - y
  229.      *    i^2 + q^2  >  (compos_lim - y)^2
  230.      *
  231.      */
  232.  
  233.     c2 = (long)i * i + (long)q * q;
  234.     y2 = (long)icompos_lim - y;
  235.     y2 *= y2;
  236.     if (c2 <= ichroma_lim2 && c2 <= y2)    /* no problems */
  237.         return 0;
  238.  
  239.     /*
  240.      * Pixel is hot, choose desired (compilation time controlled) strategy
  241.      */
  242. #if FLAG_HOT
  243.     /*
  244.      * Set the hot pixel to black to identify it.
  245.      */
  246.     p->r = p->g = p->b = 0;
  247. #else FLAG_HOT
  248.     /*
  249.      * Optimization: cache the last-computed hot pixel.
  250.      */
  251.     if (r == prev_r && g == prev_g && b == prev_b) {
  252.         p->r = new_r;
  253.         p->g = new_g;
  254.         p->b = new_b;
  255.         return 1;
  256.     }
  257.     prev_r = r;
  258.     prev_g = g;
  259.     prev_b = b;
  260.  
  261.     /*
  262.      * Get Y and chroma amplitudes in floating point.
  263.      *
  264.      * If your C library doesn't have hypot(), just use
  265.      * hypot(a,b) = sqrt(a*a, b*b);
  266.      *
  267.      * Then extract linear (un-gamma-corrected) floating-point
  268.      * pixel RGB values.
  269.      */
  270.     fy = (double)y / SCALE;
  271.     fc = hypot((double)i / SCALE, (double)q / SCALE);
  272.  
  273.     pr = pix_decode(r);
  274.     pg = pix_decode(g);
  275.     pb = pix_decode(b);
  276.  
  277.     /*
  278.      * Reducing overall pixel intensity by scaling
  279.      * R, G, and B reduces Y, I, and Q by the same factor.
  280.      * This changes luminance but not saturation, since saturation
  281.      * is determined by the chroma/luminance ratio.
  282.      *
  283.      * On the other hand, by linearly interpolating between the
  284.      * original pixel value and a grey pixel with the same
  285.      * luminance (R=G=B=Y), we change saturation without
  286.      * affecting luminance.
  287.      */
  288.  
  289. #if !REDUCE_SAT
  290.     /*
  291.      * Calculate a scale factor that will bring the pixel
  292.      * within both chroma and composite limits, if we scale
  293.      * luminance and chroma simultaneously.
  294.      *
  295.      * The calculated chrominance reduction applies to the
  296.      * gamma-corrected RGB values that are the input to
  297.      * the RGB-to-YIQ operation.  Multiplying the
  298.      * original un-gamma-corrected pixel values by
  299.      * the scaling factor raised to the "gamma" power
  300.      * is equivalent, and avoids calling gc() and inv_gc()
  301.      * three times each.
  302.      */
  303.     scale = chroma_lim / fc;
  304.     t = compos_lim / (fy + fc);
  305.     if (t < scale)
  306.         scale = t;
  307.     scale = pow(scale, GAMMA);
  308.  
  309.     r = pix_encode(scale * pr);
  310.     g = pix_encode(scale * pg);
  311.     b = pix_encode(scale * pb);
  312. #else REDUCE_SAT
  313.     /*
  314.      * Calculate a scale factor that will bring the pixel
  315.      * within both chroma and composite limits, if we scale
  316.      * chroma while leaving luminance unchanged.
  317.      *
  318.      * We have to interpolate gamma-corrected RGB values,
  319.      * so we must convert from linear to gamma-corrected
  320.      * before interpolation and then back to linear afterwards.
  321.      */
  322.     scale = chroma_lim / fc;
  323.     t = (compos_lim - fy) / fc;
  324.     if (t < scale)
  325.         scale = t;
  326.  
  327.     pr = gc(pr);
  328.     pg = gc(pg);
  329.     pb = gc(pb);
  330.     py = pr * code_matrix[0][0] + pg * code_matrix[0][1]
  331.         + pb * code_matrix[0][2];
  332.     r = pix_encode(inv_gc(py + scale * (pr - py)));
  333.     g = pix_encode(inv_gc(py + scale * (pg - py)));
  334.     b = pix_encode(inv_gc(py + scale * (pb - py)));
  335. #endif REDUCE_SAT
  336.  
  337.     p->r = new_r = r;
  338.     p->g = new_g = g;
  339.     p->b = new_b = b;
  340. #endif FLAG_HOT
  341.     return 1;
  342. }
  343.  
  344. /*
  345.  * gc: apply the gamma correction specified for this video standard.
  346.  * inv_gc: inverse function of gc.
  347.  *
  348.  * These are generally just a call to pow(), but be careful!
  349.  * Future standards may use more complex functions.
  350.  * (e.g. SMPTE 240M's "electro-optic transfer characteristic").
  351.  */
  352.  
  353. double
  354. gc(x)
  355. double    x;
  356. {
  357.     extern double    pow();
  358.  
  359.     return pow(x, 1.0 / GAMMA);
  360. }
  361.  
  362. double
  363. inv_gc(x)
  364. double    x;
  365. {
  366.     extern double    pow();
  367.  
  368.     return pow(x, GAMMA);
  369. }
  370.  
  371. /*
  372.  * pix_decode: decode an integer pixel value into a floating-point
  373.  *    intensity in the range [0, 1].
  374.  *
  375.  * pix_encode: encode a floating-point intensity into an integer
  376.  *    pixel value.
  377.  *
  378.  * The code given here assumes simple linear encoding; you must change
  379.  * these routines if you use a different pixel encoding technique.
  380.  */
  381.  
  382. double
  383. pix_decode(v)
  384. int    v;
  385. {
  386.     return (double)v / MAXPIX;
  387. }
  388.  
  389. int
  390. pix_encode(v)
  391. double    v;
  392. {
  393.     return (int)(v * MAXPIX + 0.5);
  394. }
  395.  
  396.