#include #include #include #include #include "colourspace.h" static double intensity_to_voltage_srgb( double val ) { /* Handle invalid values before doing a gamma transform. */ if( val < 0.0 ) return 0.0; if( val > 1.0 ) return 1.0; /* sRGB uses an exponent of 2.4 in its nonstandard curve, but it * still advertises a display gamma of 2.2. */ if( val <= 0.0031308 ) { return 12.92 * val; } return ( ( 1.055 * pow( val, 1.0 / 2.4 ) ) - 0.055 ); } static double voltage_to_intensity_srgb( double val ) { /* Handle invalid values before doing a gamma transform. */ if( val < 0.0 ) return 0.0; if( val > 1.0 ) return 1.0; if( val <= 0.04045 ) { return val / 12.92; } return pow( ( val + 0.055 ) / 1.055, 2.4 ); } /** * Rec-709 specifies a gamma of 2.2 and the following curve with linear * segment. */ static double intensity_to_voltage_rec709( double val ) { /* Handle invalid values before doing a gamma transform. */ if( val < 0.0 ) return 0.0; if( val > 1.0 ) return 1.0; if( val <= 0.018 ) { return 4.5 * val; } return ( ( 1.099 * pow( val, 1.0 / 2.2 ) ) - 0.099 ); } static double voltage_to_intensity_rec709( double val ) { /* Handle invalid values before doing a gamma transform. */ if( val < 0.0 ) return 0.0; if( val > 1.0 ) return 1.0; if( val <= 0.081 ) { return val / 4.5; } return pow( ( val + 0.099 ) / 1.099, 2.2 ); } void linear_to_nonlinear_srgb( double r, double g, double b, double *rp, double *gp, double *bp ) { *rp = intensity_to_voltage_srgb( r ); *gp = intensity_to_voltage_srgb( g ); *bp = intensity_to_voltage_srgb( b ); } void nonlinear_to_linear_srgb( double rp, double gp, double bp, double *r, double *g, double *b ) { *r = voltage_to_intensity_srgb( rp ); *g = voltage_to_intensity_srgb( gp ); *b = voltage_to_intensity_srgb( bp ); } void linear_to_nonlinear_rec709( double r, double g, double b, double *rp, double *gp, double *bp ) { *rp = intensity_to_voltage_rec709( r ); *gp = intensity_to_voltage_rec709( g ); *bp = intensity_to_voltage_rec709( b ); } void nonlinear_to_linear_rec709( double rp, double gp, double bp, double *r, double *g, double *b ) { *r = voltage_to_intensity_rec709( rp ); *g = voltage_to_intensity_rec709( gp ); *b = voltage_to_intensity_rec709( bp ); } void xyz_to_rgb_rec709( double x, double y, double z, double *r, double *g, double *b ) { *r = ( 3.240479 * x ) + ( -1.537150 * y ) + ( -0.498535 * z ); *g = ( -0.969256 * x ) + ( 1.875992 * y ) + ( 0.041556 * z ); *b = ( 0.055648 * x ) + ( -0.204043 * y ) + ( 1.057311 * z ); } void rgb_to_xyz_rec709( double r, double g, double b, double *x, double *y, double *z ) { *x = ( 0.412453 * r ) + ( 0.357580 * g ) + ( 0.180423 * b ); *y = ( 0.212671 * r ) + ( 0.715160 * g ) + ( 0.072169 * b ); *z = ( 0.019334 * r ) + ( 0.119193 * g ) + ( 0.950227 * b ); } void xyz_to_srgb( double x, double y, double z, double *r, double *g, double *b ) { *r = ( 3.2406 * x ) + ( -1.5372 * y ) + ( -0.4986 * z ); *g = ( -0.9689 * x ) + ( 1.8758 * y ) + ( 0.0415 * z ); *b = ( 0.0557 * x ) + ( -0.2040 * y ) + ( 1.0570 * z ); } void srgb_to_xyz( double r, double g, double b, double *x, double *y, double *z ) { *x = ( 0.4124 * r ) + ( 0.3576 * g ) + ( 0.1805 * b ); *y = ( 0.2126 * r ) + ( 0.7152 * g ) + ( 0.0722 * b ); *z = ( 0.0193 * r ) + ( 0.1192 * g ) + ( 0.9505 * b ); } void xyz_to_fakelab( double x, double y, double z, double *lstar, double *astar, double *bstar ) { if( lstar ) { *lstar = pow( y, 1.0 / 3.0 ); if( *lstar < 0.0 ) *lstar = 0.0; } if( astar ) *astar = pow( x, 1.0 / 3.0 ); if( bstar ) *bstar = pow( z, 1.0 / 3.0 ); } void fakelab_to_xyz( double lstar, double astar, double bstar, double *x, double *y, double *z ) { if( y ) *y = lstar * lstar * lstar; if( x ) *x = astar * astar * astar; if( z ) *z = bstar * bstar * bstar; } void rgb_to_fakelab_rec709( double r, double g, double b, double *lstar, double *astar, double *bstar ) { double x, y, z; rgb_to_xyz_rec709( r, g, b, &x, &y, &z ); xyz_to_fakelab( x, y, z, lstar, astar, bstar ); } void fakelab_to_rgb_rec709( double lstar, double astar, double bstar, double *r, double *g, double *b ) { double x, y, z; fakelab_to_xyz( lstar, astar, bstar, &x, &y, &z ); xyz_to_rgb_rec709( x, y, z, r, g, b ); } const double uprimen = 4.0 / ( 1.0 + 15.0 + 3.0 ); const double vprimen = 9.0 / ( 1.0 + 15.0 + 3.0 ); void xyz_to_luv( double x, double y, double z, double *lstar, double *ustar, double *vstar ) { double uprime = ( 4.0 * x ) / ( x + ( 15.0 * y ) + ( 3.0 * z ) ); double vprime = ( 9.0 * y ) / ( x + ( 15.0 * y ) + ( 3.0 * z ) ); *lstar = ( 116.0 * pow( y, 1.0 / 3.0 ) ) - 16.0; if( *lstar < 0.0 ) *lstar = 0.0; *ustar = 13.0 * (*lstar) * ( uprime - uprimen ); *vstar = 13.0 * (*lstar) * ( vprime - vprimen ); } void luv_to_xyz( double lstar, double ustar, double vstar, double *x, double *y, double *z ) { double uprime = ( ustar / ( 13.0 * lstar ) ) + uprimen; double vprime = ( vstar / ( 13.0 * lstar ) ) + vprimen; double ytmp = ( lstar + 16.0 ) / 116.0; double newy = ( ytmp * ytmp * ytmp ); uvy_to_xyz( uprime, vprime, newy, x, y, z ); } void xyz_to_lupvp( double x, double y, double z, double *lstar, double *uprime, double *vprime ) { *uprime = ( 4.0 * x ) / ( x + ( 15.0 * y ) + ( 3.0 * z ) ); *vprime = ( 9.0 * y ) / ( x + ( 15.0 * y ) + ( 3.0 * z ) ); *lstar = ( 116.0 * pow( y, 1.0 / 3.0 ) ) - 16.0; if( *lstar < 0.0 ) *lstar = 0.0; } void lupvp_to_xyz( double lstar, double uprime, double vprime, double *x, double *y, double *z ) { double ytmp = ( lstar + 16.0 ) / 116.0; double newy = ( ytmp * ytmp * ytmp ); uvy_to_xyz( uprime, vprime, newy, x, y, z ); } void xyy_to_xyz( double x, double y, double Y, double *rx, double *ry, double *rz ) { *ry = Y; *rx = ( x / y ) * Y; *rz = ( ( 1.0 - x - y ) / y ) * Y; /* *ry = Y; *rx = ( x / y ) * Y; *rz = ( *rx / x ) - *rx - Y; */ } void get_natural_Y( double x, double y, double Y, double *ry ) { double ix, iy, iz, t; xyy_to_xyz( x, y, Y, &ix, &iy, &iz ); t = - ( ( -1.0 ) / ( ix + iy + iz ) ); *ry = iy * t; } void uvy_to_xyz( double up, double vp, double Y, double *rx, double *ry, double *rz ) { *ry = Y; *rx = ( 9.0 * Y * up ) / ( 4.0 * vp ); *rz = ( ( ( 4.0 * *rx ) / up ) - *rx - ( 15.0 * Y ) ) / 3.0; } int intersect_two_line_segments( double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, double *x, double *y ) { double minx1, maxx1, miny1, maxy1; double minx2, maxx2, miny2, maxy2; if( x1 == x2 && x4 == x3 ) return 0; if( x2 == x1 ) { double s2 = ( y4 - y3 ) / ( x4 - x3 ); double b2 = y3 - ( s2 * x3 ); *x = x1; *y = ( s2 * (*x) ) + b2; } else if( x4 == x3 ) { return intersect_two_line_segments( x3, y3, x4, y4, x1, y1, x2, y2, x, y ); } else { double s1 = ( y2 - y1 ) / ( x2 - x1 ); double s2 = ( y4 - y3 ) / ( x4 - x3 ); double b1 = y1 - ( s1 * x1 ); double b2 = y3 - ( s2 * x3 ); if( s1 == s2 ) return 0; *x = ( b2 - b1 ) / ( s1 - s2 ); *y = ( s1 * (*x) ) + b1; } minx1 = x1; if( x2 < minx1 ) minx1 = x2; maxx1 = x1; if( x2 > maxx1 ) maxx1 = x2; miny1 = y1; if( y2 < miny1 ) miny1 = y2; maxy1 = y1; if( y2 > maxy1 ) maxy1 = y2; minx2 = x3; if( x4 < minx2 ) minx2 = x4; maxx2 = x3; if( x4 > maxx2 ) maxx2 = x4; miny2 = y3; if( y4 < miny2 ) miny2 = y4; maxy2 = y3; if( y4 > maxy2 ) maxy2 = y4; if( *x >= minx1 && *x <= maxx1 && *y >= miny1 && *y <= maxy1 ) if( *x >= minx2 && *x <= maxx2 && *y >= miny2 && *y <= maxy2 ) return 1; return 0; } void closest_in_gamut_xy( double px, double py, double wx, double wy, double rx, double ry, double gx, double gy, double bx, double by, double *cx, double *cy ) { if( intersect_two_line_segments( px, py, wx, wy, rx, ry, gx, gy, cx, cy ) ) return; if( intersect_two_line_segments( px, py, wx, wy, bx, by, rx, ry, cx, cy ) ) return; if( intersect_two_line_segments( px, py, wx, wy, gx, gy, bx, by, cx, cy ) ) return; *cx = 1.0; *cy = 1.0; // assert( 0 ); } void closest_in_gamut_srgb( double r, double g, double b, double wr, double wg, double wb, double *ri, double *gi, double *bi ) { *ri = r; *gi = g; *bi = b; if( *ri < 0.0 ) *ri = 0.0; if( *gi < 0.0 ) *gi = 0.0; if( *bi < 0.0 ) *bi = 0.0; return; if( *ri >= 0.0 && *gi >= 0.0 && *bi >= 0.0 ) return; // wr = wg = wb = 1.0; for(;;) { // XXX wtf?! // double vr = wr - *ri; // double vg = wg - *gi; // double vb = wb - *bi; // double scale; if( *ri >= 0.0 && *gi >= 0.0 && *bi >= 0.0 ) return; if( *ri < *gi && *ri < *bi ) { *gi += -*ri; *bi += -*ri; *ri = 0.0; } else if( *gi < *ri && *gi < *bi ) { *bi += -*gi; *ri += -*gi; *gi = 0.0; } else { *ri += -*bi; *gi += -*bi; *bi = 0.0; } /* if( *ri < *gi && *ri < *bi ) { scale = ( -(*ri) ) / vr; } else if( *gi < *ri && *gi < *bi ) { scale = ( -(*gi) ) / vg; } else { scale = ( -(*bi) ) / vb; } *ri = *ri + ( scale * vr ); *gi = *gi + ( scale * vg ); *bi = *bi + ( scale * vb ); */ } } /* chromaticities: S170M: x y g 0.310 0.595 b 0.155 0.070 r 0.630 0.340 D65 0.3127 0.3290 BT.470 System M: g 0.21 0.71 b 0.14 0.08 r 0.67 0.33 C 0.310 0.316 BT.470 System B,G: g 0.29 0.60 b 0.15 0.06 r 0.64 0.33 D65 0.313 0.329 Poynton on apple/sony monitors: http://www.inforamp.net/~poynton/notes/short_subjects/comp.graphics/color_temperature_for_Apple x y tol persistance R .625 .340 +-.030 1 ms G .280 .595 +-.030 40 us B .155 .070 +-.016 30 us W .283 .298 +-.030 (9300 K +8 MPCD) Name Red Green Blue White Point Source xr yr xg yg xb yb xw yw Short-Persistence 0.61 0.35 0.29 0.59 0.15 0.063 N/A N/A [Foley96, p. 583] Long-Persistence 0.62 0.33 0.21 0.685 0.15 0.063 N/A N/A [Foley96, p. 583] NTSC 0.67 0.33 0.21 0.71 0.14 0.08 Illuminant C [Walker98] CCIR 601-1 is the old NTSC standard. Now called Rec 601-1. EBU 0.64 0.33 0.30 0.60 0.15 0.06 Illuminant D65 [Walker98] Dell (all monitors except 21" Mitsubishi p/n 65532) 0.625 0.340 0.275 0.605 0.150 0.065 9300 K [Dell, E-mail, 12 Jan 99] SMPTE 0.630 0.340 0.310 0.595 0.155 0.070 Illuminant D65 [Walker98] P22 phosphor in NEC Multisync C400 0.610 0.350 0.307 0.595 0.150 0.065 0.280 0.315 [NEC98] 9300 K Gamma = 2.2 P22 phosphor in KDS VS19 0.625 0.340 0.285 0.605 0.150 0.065 0.281 0.311 High Brightness LEDs 0.700 0.300 0.170 0.700 0.130 0.075 0.310 0.320 Nichia Corporation Illuminant C 6774 0.31006 0.31616 0.3104 0.3191 [Wyszecki82, p 139] [Agoston87, p. 103] */ struct matrix_s { double d[ 3 ][ 3 ]; }; double matrix_determinant( matrix_t *m ) { return m->d[0][0]*m->d[1][1]*m->d[2][2]+ m->d[0][1]*m->d[1][2]*m->d[2][0]+ m->d[0][2]*m->d[1][0]*m->d[2][1]- m->d[0][0]*m->d[1][2]*m->d[2][1]- m->d[0][1]*m->d[1][0]*m->d[2][2]- m->d[0][2]*m->d[1][1]*m->d[2][0]; } void matrix_scale( matrix_t *r, matrix_t *a, double b ) { int x, y; for( x = 0; x < 3; x++ ) for( y = 0; y < 3; y++ ) r->d[x][y] = a->d[x][y] * b; } void matrix_transpose( matrix_t *r, matrix_t *a ) { int x, y; for( x = 0; x < 3; x++ ) for( y = 0; y < 3; y++ ) r->d[x][y] = a->d[y][x]; } void matrix_print( matrix_t *m ) { int x, y; for( x = 0; x < 3; x++ ) { for( y = 0; y < 3; y++ ) { fprintf( stderr, "%f ", m->d[x][y] ); } fprintf( stderr, "\n"); } } void matrix_multiply( matrix_t *r, matrix_t *a, matrix_t *b ) { int x, y; for( x = 0; x < 3; x++ ) for( y = 0; y < 3; y++ ) r->d[x][y] = (a->d[x][0] * b->d[0][y]) + (a->d[x][1] * b->d[1][y]) + (a->d[x][2] * b->d[2][y]); } int matrix_invert( matrix_t *r, matrix_t *m ) { /* Inverse matrix is cofactor matrix divided by the determinant */ double det; det = matrix_determinant( m ); if( !det ) return 0; /* matrix of cofactors */ r->d[0][0] = m->d[1][1]*m->d[2][2]-m->d[1][2]*m->d[2][1]; r->d[0][1] = m->d[2][0]*m->d[1][2]-m->d[1][0]*m->d[2][2]; r->d[0][2] = m->d[1][0]*m->d[2][1]-m->d[2][0]*m->d[1][1]; r->d[1][0] = m->d[2][1]*m->d[0][2]-m->d[0][1]*m->d[2][2]; r->d[1][1] = m->d[0][0]*m->d[2][2]-m->d[0][2]*m->d[2][0]; r->d[1][2] = m->d[2][0]*m->d[0][1]-m->d[0][0]*m->d[2][1]; r->d[2][0] = m->d[0][1]*m->d[1][2]-m->d[1][1]*m->d[0][2]; r->d[2][1] = m->d[1][0]*m->d[0][2]-m->d[0][0]*m->d[1][2]; r->d[2][2] = m->d[0][0]*m->d[1][1]-m->d[1][0]*m->d[0][1]; /* Multiply it by 1/det */ det = 1.0 / det; matrix_scale( r, r, det ); return 1; } /** * Using the method described in SMPTE RP-177. */ matrix_t *build_rgb_xyz_matrix_from_chromaticities( double rx, double ry, double gx, double gy, double bx, double by, double wx, double wy ) { matrix_t *final; matrix_t P, Ptrans, Pinv, C, NPM, NPMinv, NPMtrans; double wX, wY, wZ; double rz, gz, bz; double Cr, Cg, Cb; /** * We need the XYZ coordinates for our white point first. We use a * Y value of 1 so that our RGB value for white will be at (1,1,1). */ xyy_to_xyz( wx, wy, 1.0, &wX, &wY, &wZ ); fprintf( stderr, "White point is:\n" "Xw: %f, Yw: %f, Zw: %f\n", wX, wY, wZ ); /** * Construct the matrix of primaries. */ rz = 1.0 - rx - ry; gz = 1.0 - gx - gy; bz = 1.0 - bx - by; P.d[0][0] = rx; P.d[0][1] = gx; P.d[0][2] = bx; P.d[1][0] = ry; P.d[1][1] = gy; P.d[1][2] = by; P.d[2][0] = rz; P.d[2][1] = gz; P.d[2][2] = bz; if( !matrix_invert( &Ptrans, &P ) ) { fprintf( stderr, "colourspace: Can't invert colour primaries matrix!\n" ); return 0; } matrix_transpose( &Pinv, &Ptrans ); fprintf( stderr, "P matrix is:\n" ); matrix_print( &P ); fprintf( stderr, "P inverse is:\n" ); matrix_print( &Pinv ); Cr = (Pinv.d[0][0] * wX) + (Pinv.d[0][1] * wY) + (Pinv.d[0][2] * wZ); Cg = (Pinv.d[1][0] * wX) + (Pinv.d[1][1] * wY) + (Pinv.d[1][2] * wZ); Cb = (Pinv.d[2][0] * wX) + (Pinv.d[2][1] * wY) + (Pinv.d[2][2] * wZ); fprintf( stderr, "C values are:\n" "Cr: %f, Cg: %f, Cb: %f\n", Cr, Cg, Cb ); C.d[0][0] = Cr; C.d[0][1] = 0.0; C.d[0][2] = 0.0; C.d[1][0] = 0.0; C.d[1][1] = Cg; C.d[1][2] = 0.0; C.d[2][0] = 0.0; C.d[2][1] = 0.0; C.d[1][2] = Cb; matrix_multiply( &NPM, &Ptrans, &C ); NPM.d[0][0] = rx * Cr; NPM.d[0][1] = gx * Cg; NPM.d[0][2] = bx * Cb; NPM.d[1][0] = ry * Cr; NPM.d[1][1] = gy * Cg; NPM.d[1][2] = by * Cb; NPM.d[2][0] = rz * Cr; NPM.d[2][1] = gz * Cg; NPM.d[2][2] = bz * Cb; fprintf( stderr, "matrix is:\n" ); matrix_print( &NPM ); if( !matrix_invert( &NPMtrans, &NPM ) ) { fprintf( stderr, "colourspace: Can't invert conversion matrix!\n" ); return 0; } matrix_transpose( &NPMinv, &NPMtrans ); fprintf( stderr, "inv matrix is:\n" ); matrix_print( &NPMinv ); final = (matrix_t *) malloc( sizeof( matrix_t ) ); *final = NPM; return final; } /** * 601 numbers: * * Y' = 0.299*R' + 0.587*G' + 0.114*B' (in 0.0 to 1.0) * Cb = -0.169*R' - 0.331*G' + 0.500*B' (in -0.5 to +0.5) * Cr = 0.500*R' - 0.419*G' - 0.081*B' (in -0.5 to +0.5) * * S170M numbers: * Y' = 0.299*R' + 0.587*G' + 0.114*B' (in 0.0 to 1.0) * R-Y' = 0.701*R' - 0.587*G' - 0.114*B' * B-Y' = -0.299*R' - 0.587*G' + 0.886*B' */ /** * numbers from MPEG2 spec: * output Y' = ( 219 * Y' ) + 16 * output Cb = ( 224 * Pb ) + 128 * output Cr = ( 224 * Pr ) + 128 * * Rec-709: * Y' = 0.7154*G' + 0.0721*B' + 0.2125*R' * Pb = -0.386 *G' + 0.500 *B' - 0.115 *R' * Pr = -0.454 *G' - 0.046 *B' + 0.500 *R' * * 'FCC' (corresponds to 470-M?): * Y' = 0.59 *G' + 0.11 *B' + 0.30 *R' * Pb = -0.331*G' + 0.500*B' - 0.169*R' * Pr = -0.421*G' - 0.079*B' + 0.500*R' * * 470-2 System B,G: * Y' = 0.587*G' + 0.114*B' + 0.299*R' * Pb = -0.331*G' + 0.500*B' - 0.169*R' * Pr = -0.419*G' - 0.081*B' + 0.500*R' * * SMPTE 170M: * Y' = 0.587*G' + 0.114*B' + 0.299*R' * Pb = -0.331*G' + 0.500*B' - 0.169*R' * Pr = -0.419*G' - 0.081*B' + 0.500*R' * * SMPTE 240M: * Y' = 0.701*G' + 0.087*B' + 0.212*R' * Pb = -0.384*G' + 0.500*B' - 0.116*R' * Pr = -0.445*G' - 0.055*B' + 0.500*R' */ void matrix_rgb_to_xyz( matrix_t *conv, double r, double g, double b, double *x, double *y, double *z ) { *x = (conv->d[0][0] * r) + (conv->d[0][1] * g) + (conv->d[0][2] * b); *y = (conv->d[1][0] * r) + (conv->d[1][1] * g) + (conv->d[1][2] * b); *z = (conv->d[2][0] * r) + (conv->d[2][1] * g) + (conv->d[2][2] * b); } /** * Converts gamma-corrected RGB in [0,1] to Y'Pbpr in [0,1]/[-0.5,0.5]. */ void mpeg2_rgb_to_ypbpr_rec709( double r, double g, double b, double *yp, double *pb, double *pr ) { *yp = (0.7154 * g) + (0.0721 * b) + (0.2125 * r); *pb = (-0.386 * g) + (0.500 * b) - (0.115 * r); *pr = (-0.454 * g) - (0.046 * b) + (0.500 * r); } void mpeg2_ypbpr_to_rgb_rec709( double yp, double pb, double pr, double *r, double *g, double *b ) { *g = (0.99981 * yp) - (0.18723 * pb) - (0.46798 * pr); *b = (1.0019 * yp) + (1.8556 * pb) + (0.0010051 * pr); *r = (1.00000 * yp) + (0.00071500 * pb) + (1.5752 * pr); } void mpeg2_rgb_to_ypbpr_fcc( double r, double g, double b, double *yp, double *pb, double *pr ) { /** * MPEG2 gives this formula with less precision. I'm going to * instead use the numbers from Rec.601. * *yp = (0.59 * g) + (0.11 * b) + (0.30 * r); */ *yp = ( 0.587 * g) + (0.114 * b) + (0.299 * r); *pb = (-0.331 * g) + (0.500 * b) - (0.169 * r); *pr = (-0.421 * g) - (0.079 * b) + (0.500 * r); } void mpeg2_ypbpr_to_rgb_fcc( double yp, double pb, double pr, double *r, double *g, double *b ) { /** * MPEG2 gives this formula with less precision. I'm going to * instead use the numbers from Rec.601. * *yp = (0.59 * g) + (0.11 * b) + (0.30 * r); */ *g = (1.0 * yp) - (0.34068 * pb) - (0.71315 * pr); *b = (1.0 * yp) + (1.7722 * pb) + (0.0009881 * pr); *r = (1.0 * yp) - (0.0068498 * pb) + (1.3997 * pr); } /** * POS bicubic resampler. */ double bicubic( double x, double y, double *image, int width, int height, int linestride ) { double output; int i, j, nx, ny; nx = (int) floor( x ); ny = (int) floor( y ); output = 0.0; for( i = 0; i < 4; i++ ) { int cy = ny + i - 1; for( j = 0; j < 4; j++ ) { int cx = nx + j - 1; double tx, hx; double ty, hy; double curpixel; if( cx < 0 || cy < 0 || cx >= width || cy >= height ) { continue; } /* Calculate delta x and delta y. */ tx = fabs( x - (double) cx ); ty = fabs( y - (double) cy ); /* Calculate the x contribution. */ if( tx < 1.0 ) { hx = 1.0 - ( 2.0 * tx * tx ) + ( tx * tx * tx ); } else { hx = 4.0 - ( 8.0 * tx ) + ( 5.0 * tx * tx ) - ( tx * tx * tx ); } /* Calculate the y contribution. */ if( ty < 1.0 ) { hy = 1.0 - ( 2.0 * ty * ty ) + ( ty * ty * ty ); } else { hy = 4.0 - ( 8.0 * ty ) + ( 5.0 * ty * ty ) - ( ty * ty * ty ); } /* Add in the contribution from this pixel. */ curpixel = image[ (cy * linestride) + cx ]; output += curpixel * ( hx * hy ); } } return output; } #if 1 void resample_chroma_420_to_444( double *c420, double *c444, int width, int height ) { double poffh = 1.0 / 3.0; double poffw = 0.5; int i, j; for( i = 0; i < height; i++ ) { for( j = 0; j < width; j++ ) { /* we do four samples per input sample. */ double curposw = (double) j; double curposh = (double) i; c444[ 0 ] = bicubic( curposw, curposh - poffh, c420, width, height, width ); c444[ 1 ] = bicubic( curposw + poffw, curposh - poffh, c420, width, height, width ); c444[ (width*2) ] = bicubic( curposw, curposh + poffh, c420, width, height, width ); c444[ (width*2)+1 ] = bicubic( curposw + poffw, curposh + poffh, c420, width, height, width ); c444 += 2; } c444 += width*2; } } #else void resample_chroma_420_to_444( double *c420, double *c444, int width, int height ) { int i, j; for( i = 0; i < height; i++ ) { for( j = 0; j < width; j++ ) { c444[ 0 ] = *c420; c444[ 1 ] = *c420; c444[ (width*2) ] = *c420; c444[ (width*2)+1 ] = *c420; c444 += 2; c420++; } c444 += width*2; } } #endif