/*
======================================================================
hot.c

Ernie Wright  7 Jul 98
MSVC 4.0

Calculate whether an RGB pixel value will exceed certain broadcast
signal limits after encoding.

Based on David Martindale and Alan W. Paeth, "Television Color
Encoding and 'Hot' Broadcast Colors", in GRAPHICS GEMS II, James Arvo
ed., Academic Press, 1991, pp. 147-158.

   init_hot()     initialize
   hot()          determine whether a pixel is hot
====================================================================== */

#include <math.h>

#define SCALE       8192         /* scale factor: do floats with int math */
#define MAXPIXEL    255          /* white value */
#define GAMMA_NTSC  2.2
#define GAMMA_PAL   2.8

#define gc( x )          pow( x, gam )
#define pix_decode( v )  (( double ) v / MAXPIXEL )


static double code_matrix_ntsc[ 3 ][ 3 ] = {    // RGB to YIQ
    0.2989,  0.5866,  0.1144,
    0.5959, -0.2741, -0.3218,
    0.2113, -0.5227,  0.3113,
};

static double code_matrix_pal[ 3 ][ 3 ] = {     //  RGB to YUV
    0.2989,  0.5866,  0.1144,
   -0.1473, -0.2891,  0.4364,
    0.6149, -0.5145, -0.1004,
};

static int    tab[ 3 ][ 3 ][ MAXPIXEL + 1 ];
static int    ichroma_lim2;
static int    icompos_lim;
static double gam;


/*
======================================================================
init_hot()

Initialize the lookup table and parameters.

INPUTS
   encoding       0 = NTSC, 1 = PAL
   ped            pedestal, normally 7.5 for NTSC and 0.0 for PAL
   com            composite limit, IRE units (e.g. 110)
   chr            chroma limit, IRE units (e.g., 50)

RESULTS
   Initializes the lookup table and parameters used during the hot()
   calculation.

Composite and chroma limits need to be converted into the units used
internally for YIQ (YUV).  The conversion depends on the pedestal,
since as Y goes from 0 to 1, the signal goes from the pedestal level
to 100 IRE.  Chroma is always scaled to remain consistent with Y.
====================================================================== */

void init_hot( int encoding, double ped, double com, double chr )
{
   double f;
   int pv, i, j;

   if ( encoding == 0 )
      gam = 1.0 / GAMMA_NTSC;
   else
      gam = 1.0 / GAMMA_PAL;

   for ( pv = 0; pv <= MAXPIXEL; pv++ ) {
      f = SCALE * gc( pix_decode( pv ));

      for ( i = 0; i < 3; i++ )
         for ( j = 0; j < 3; j++ )
            if ( encoding == 0 ) tab[ i ][ j ][ pv ] =
               ( int )( f * code_matrix_ntsc[ i ][ j ] + 0.5 );
            else tab[ i ][ j ][ pv ] =
               ( int )( f * code_matrix_pal[ i ][ j ] + 0.5 );
   }

   icompos_lim  = ( int )((( com - ped ) / ( 100.0 - ped )) * SCALE + 0.5 );
   ichroma_lim2 = ( int )(( chr / ( 100.0 - ped )) * SCALE + 0.5 );
   ichroma_lim2 *= ichroma_lim2;
}


/*
======================================================================
hot()

Check to see if the chrominance vector is too long or the composite
waveform amplitude is too large.

Pixel decoding, gamma correction, and matrix multiplication are all
done by lookup table.

i and q are the two chrominance components.  For NTSC, they are I and
Q.  For PAL, i is U (scaled B-Y) and q is V (scaled R-Y).  Since we
only care about the length of the chroma vector, not its angle, we
don't care which is which.

Chrominance is too large if

   sqrt( i^2 + q^2 ) > chroma_lim.

The composite signal amplitude is too large if

   y + sqrt( i^2 + q^2 ) > compos_lim.

We avoid doing the sqrt by checking

   i^2 + q^2 > chroma_lim^2

and

   y + sqrt( i^2 + q^2 ) > compos_lim
   sqrt( i^2 + q^2 )     > compos_lim - y
   i^2 + q^2             > ( compos_lim - y )^2

Note that this fails when y > compos_lim, a possibility I don't think
Martindale and Paeth contemplated.  But we can test for that case
directly.
====================================================================== */

int hot( int r, int g, int b )
{
   int y, i, q, y2, c2;

   y = tab[ 0 ][ 0 ][ r ] + tab[ 0 ][ 1 ][ g ] + tab[ 0 ][ 2 ][ b ];
   i = tab[ 1 ][ 0 ][ r ] + tab[ 1 ][ 1 ][ g ] + tab[ 1 ][ 2 ][ b ];
   q = tab[ 2 ][ 0 ][ r ] + tab[ 2 ][ 1 ][ g ] + tab[ 2 ][ 2 ][ b ];

   if ( y > icompos_lim ) return 1;

   c2 = i * i + q * q;
   if ( c2 > ichroma_lim2 ) return 1;

   y2 = icompos_lim - y;
   y2 *= y2;
   return ( c2 > y2 );
}