/* pnmgamma.c - perform gamma correction on a PNM image ** ** Copyright (C) 1991 by Bill Davidson and Jef Poskanzer. ** ** Permission to use, copy, modify, and distribute this software and its ** documentation for any purpose and without fee is hereby granted, provided ** that the above copyright notice appear in all copies and that both that ** copyright notice and this permission notice appear in supporting ** documentation. This software is provided "as is" without express or ** implied warranty. */ #include #include #include #include "pm_c_util.h" #include "mallocvar.h" #include "shhopt.h" #include "pnm.h" enum transferFunction { XF_EXP, XF_EXP_INVERSE, XF_BT709RAMP, XF_BT709RAMP_INVERSE, XF_SRGBRAMP, XF_SRGBRAMP_INVERSE, XF_BT709_TO_SRGB, XF_SRGB_TO_BT709 }; struct cmdlineInfo { /* All the information the user supplied in the command line, in a form easy for the program to use. */ const char * filespec; /* '-' if stdin */ enum transferFunction transferFunction; float rgamma, ggamma, bgamma; unsigned int maxval; unsigned int makeNewMaxval; }; static void interpretOldArguments(int const argc, char ** const argv, float const defaultGamma, struct cmdlineInfo * const cmdlineP) { /* Use the old syntax wherein the gamma values come from arguments. If there is one argument, it's a gamma value for all three components. If 3 arguments, it's separate gamma values. If 2, it's a single gamma value plus a file name. If 4, it's separate gamma values plus a file name. */ if (argc-1 == 0) { cmdlineP->rgamma = defaultGamma; cmdlineP->ggamma = defaultGamma; cmdlineP->bgamma = defaultGamma; cmdlineP->filespec = "-"; } else if (argc-1 == 1) { cmdlineP->rgamma = atof(argv[1]); cmdlineP->ggamma = atof(argv[1]); cmdlineP->bgamma = atof(argv[1]); cmdlineP->filespec = "-"; } else if (argc-1 == 2) { cmdlineP->rgamma = atof(argv[1]); cmdlineP->ggamma = atof(argv[1]); cmdlineP->bgamma = atof(argv[1]); cmdlineP->filespec = argv[2]; } else if (argc-1 == 3) { cmdlineP->rgamma = atof(argv[1]); cmdlineP->ggamma = atof(argv[2]); cmdlineP->bgamma = atof(argv[3]); cmdlineP->filespec = "-"; } else if (argc-1 == 4) { cmdlineP->rgamma = atof(argv[1]); cmdlineP->ggamma = atof(argv[2]); cmdlineP->bgamma = atof(argv[3]); cmdlineP->filespec = argv[4]; } else pm_error("Wrong number of arguments. " "You may have 0, 1, or 3 gamma values " "plus zero or one filename"); if (cmdlineP->rgamma <= 0.0 || cmdlineP->ggamma <= 0.0 || cmdlineP->bgamma <= 0.0 ) pm_error("Invalid gamma value. Must be positive floating point " "number."); } static void getGammaFromOpts(struct cmdlineInfo * const cmdlineP, bool const gammaSpec, float const gammaOpt, bool const rgammaSpec, bool const ggammaSpec, bool const bgammaSpec, float const defaultGamma) { if (gammaSpec) if (gammaOpt < 0.0) pm_error("Invalid gamma value %f. Must be positive.", gammaOpt); if (rgammaSpec) { if (cmdlineP->rgamma < 0.0) pm_error("Invalid red gamma value %f. Must be positive.", cmdlineP->rgamma); } else { if (gammaSpec) cmdlineP->rgamma = gammaOpt; else cmdlineP->rgamma = defaultGamma; } if (ggammaSpec) { if (cmdlineP->ggamma < 0.0) pm_error("Invalid green gamma value %f. Must be positive.", cmdlineP->ggamma); } else { if (gammaSpec) cmdlineP->ggamma = gammaOpt; else cmdlineP->ggamma = defaultGamma; } if (bgammaSpec) { if (cmdlineP->bgamma < 0.0) pm_error("Invalid blue gamma value %f. Must be positive.", cmdlineP->bgamma); } else { if (gammaSpec) cmdlineP->bgamma = gammaOpt; else cmdlineP->bgamma = defaultGamma; } } static void parseCommandLine(int argc, char ** argv, struct cmdlineInfo * const cmdlineP) { optEntry *option_def; /* Instructions to pm_optParseOptions3 on how to parse our options. */ optStruct3 opt; unsigned int bt709ramp, srgbramp, ungamma, bt709tosrgb, srgbtobt709; unsigned int bt709tolinear, lineartobt709; unsigned int gammaSpec, rgammaSpec, ggammaSpec, bgammaSpec; float gammaOpt; float defaultGamma; unsigned int option_def_index; MALLOCARRAY_NOFAIL(option_def, 100); option_def_index = 0; /* incremented by OPTENT3 */ OPTENT3(0, "ungamma", OPT_FLAG, NULL, &ungamma, 0); OPTENT3(0, "bt709tolinear", OPT_FLAG, NULL, &bt709tolinear, 0); OPTENT3(0, "lineartobt709", OPT_FLAG, NULL, &lineartobt709, 0); OPTENT3(0, "bt709ramp", OPT_FLAG, NULL, &bt709ramp, 0); OPTENT3(0, "cieramp", OPT_FLAG, NULL, &bt709ramp, 0); OPTENT3(0, "srgbramp", OPT_FLAG, NULL, &srgbramp, 0); OPTENT3(0, "bt709tosrgb", OPT_FLAG, NULL, &bt709tosrgb, 0); OPTENT3(0, "srgbtobt709", OPT_FLAG, NULL, &srgbtobt709, 0); OPTENT3(0, "maxval", OPT_UINT, &cmdlineP->maxval, &cmdlineP->makeNewMaxval, 0); OPTENT3(0, "gamma", OPT_FLOAT, &gammaOpt, &gammaSpec, 0); OPTENT3(0, "rgamma", OPT_FLOAT, &cmdlineP->rgamma, &rgammaSpec, 0); OPTENT3(0, "ggamma", OPT_FLOAT, &cmdlineP->ggamma, &ggammaSpec, 0); OPTENT3(0, "bgamma", OPT_FLOAT, &cmdlineP->bgamma, &bgammaSpec, 0); opt.opt_table = option_def; opt.short_allowed = FALSE; /* We have no short (old-fashioned) options */ opt.allowNegNum = TRUE; pm_optParseOptions3(&argc, argv, opt, sizeof(opt), 0); /* Uses and sets argc, argv, and some of *cmdline_p and others. */ if (bt709tolinear + lineartobt709 + bt709ramp + srgbramp + bt709tosrgb + srgbtobt709 > 1) pm_error("You may specify only one function option"); else { if (bt709tolinear) { if (ungamma) pm_error("You cannot specify -ungamma with -bt709tolinear"); else cmdlineP->transferFunction = XF_BT709RAMP_INVERSE; } else if (lineartobt709) { if (ungamma) pm_error("You cannot specify -ungamma with -lineartobt709"); else cmdlineP->transferFunction = XF_BT709RAMP; } else if (bt709tosrgb) { if (ungamma) pm_error("You cannot specify -ungamma with -bt709tosrgb"); else cmdlineP->transferFunction = XF_BT709_TO_SRGB; } else if (srgbtobt709) { if (ungamma) pm_error("You cannot specify -ungamma with -srgbtobt709"); else cmdlineP->transferFunction = XF_SRGB_TO_BT709; } else if (bt709ramp) { if (ungamma) cmdlineP->transferFunction = XF_BT709RAMP_INVERSE; else cmdlineP->transferFunction = XF_BT709RAMP; } else if (srgbramp) { if (ungamma) cmdlineP->transferFunction = XF_SRGBRAMP_INVERSE; else cmdlineP->transferFunction = XF_SRGBRAMP; } else { if (ungamma) cmdlineP->transferFunction = XF_EXP_INVERSE; else cmdlineP->transferFunction = XF_EXP; } } if (cmdlineP->makeNewMaxval) { if (cmdlineP->maxval > PNM_OVERALLMAXVAL) pm_error("Largest possible maxval is %u. You specified %u", PNM_OVERALLMAXVAL, cmdlineP->maxval); } switch (cmdlineP->transferFunction) { case XF_BT709RAMP: case XF_BT709RAMP_INVERSE: case XF_SRGB_TO_BT709: defaultGamma = 1.0/0.45; break; case XF_SRGBRAMP: case XF_SRGBRAMP_INVERSE: case XF_BT709_TO_SRGB: /* The whole function is often approximated with exponent 2.2 and no linear piece. We do the linear piece, so we use the real exponent of 2.4. */ defaultGamma = 2.4; break; case XF_EXP: case XF_EXP_INVERSE: defaultGamma = 2.2; break; } if (bt709tolinear || lineartobt709 || bt709tosrgb || srgbtobt709) { /* Use the new syntax wherein the gamma values come from options, not arguments. So if there's an argument, it's a file name. */ getGammaFromOpts(cmdlineP, gammaSpec, gammaOpt, rgammaSpec, ggammaSpec, bgammaSpec, defaultGamma); if (argc-1 < 1) cmdlineP->filespec = "-"; else { cmdlineP->filespec = argv[1]; if (argc-1 > 1) pm_error("Too many arguments (%u). With this function, there " "is at most one argument: the file name", argc-1); } } else { if (gammaSpec || rgammaSpec || ggammaSpec || bgammaSpec) pm_error("With this function, you specify the gamma values in " "arguments, not with the -gamma, etc."); interpretOldArguments(argc, argv, defaultGamma, cmdlineP); } } static void buildPowGamma(xelval table[], xelval const maxval, xelval const newMaxval, double const gamma) { /*---------------------------------------------------------------------------- Build a gamma table of size maxval+1 for the given gamma value. This function depends on pow(3m). If you don't have it, you can simulate it with '#define pow(x,y) exp((y)*log(x))' provided that you have the exponential function exp(3m) and the natural logarithm function log(3m). I can't believe I actually remembered my log identities. -----------------------------------------------------------------------------*/ xelval i; double const oneOverGamma = 1.0 / gamma; for (i = 0 ; i <= maxval; ++i) { double const normalized = ((double) i) / maxval; /* Xel sample value normalized to 0..1 */ double const v = pow(normalized, oneOverGamma); table[i] = MIN((xelval)(v * newMaxval + 0.5), newMaxval); /* denormalize, round and clip */ } } static void buildBt709Gamma(xelval table[], xelval const maxval, xelval const newMaxval, double const gamma) { /*---------------------------------------------------------------------------- Build a gamma table of size maxval+1 for the ITU Recommendation BT.709 gamma transfer function. 'gamma' must be 1/0.45 for true Rec. 709. -----------------------------------------------------------------------------*/ double const oneOverGamma = 1.0 / gamma; xelval i; /* This transfer function is linear for sample values 0 .. maxval*.018 and an exponential for larger sample values. The exponential is slightly stretched and translated, though, unlike the popular pure exponential gamma transfer function. */ xelval const linearCutoff = (xelval) (maxval * 0.018 + 0.5); double const linearExpansion = (1.099 * pow(0.018, oneOverGamma) - 0.099) / 0.018; double const maxvalScaler = (double)newMaxval/maxval; for (i = 0; i <= linearCutoff; ++i) table[i] = i * linearExpansion * maxvalScaler + 0.5; for (; i <= maxval; ++i) { double const normalized = ((double) i) / maxval; /* Xel sample value normalized to 0..1 */ double const v = 1.099 * pow(normalized, oneOverGamma) - 0.099; table[i] = MIN((xelval)(v * newMaxval + 0.5), newMaxval); /* denormalize, round, and clip */ } } static void buildBt709GammaInverse(xelval table[], xelval const maxval, xelval const newMaxval, double const gamma) { /*---------------------------------------------------------------------------- Build a gamma table of size maxval+1 for the Inverse of the ITU Rec. BT.709 gamma transfer function. 'gamma' must be 1/0.45 for true Rec. 709. -----------------------------------------------------------------------------*/ double const oneOverGamma = 1.0 / gamma; xelval i; /* This transfer function is linear for sample values 0 .. maxval*.018 and an exponential for larger sample values. The exponential is slightly stretched and translated, though, unlike the popular pure exponential gamma transfer function. */ xelval const linearCutoff = (xelval) (maxval * 0.018 + 0.5); double const linearCompression = 0.018 / (1.099 * pow(0.018, oneOverGamma) - 0.099); double const maxvalScaler = (double)newMaxval/maxval; for (i = 0; i <= linearCutoff / linearCompression; ++i) table[i] = i * linearCompression * maxvalScaler + 0.5; for (; i <= maxval; ++i) { double const normalized = ((double) i) / maxval; /* Xel sample value normalized to 0..1 */ double const v = pow((normalized + 0.099) / 1.099, gamma); table[i] = MIN((xelval)(v * newMaxval + 0.5), newMaxval); /* denormalize, round, and clip */ } } static void buildSrgbGamma(xelval table[], xelval const maxval, xelval const newMaxval, double const gamma) { /*---------------------------------------------------------------------------- Build a gamma table of size maxval+1 for the IEC SRGB gamma transfer function (Standard IEC 61966-2-1). 'gamma' must be 2.4 for true SRGB -----------------------------------------------------------------------------*/ double const oneOverGamma = 1.0 / gamma; xelval i; /* This transfer function is linear for sample values 0 .. maxval*.040405 and an exponential for larger sample values. The exponential is slightly stretched and translated, though, unlike the popular pure exponential gamma transfer function. */ xelval const linearCutoff = (xelval) maxval * 0.0031308 + 0.5; double const linearExpansion = (1.055 * pow(0.0031308, oneOverGamma) - 0.055) / 0.0031308; double const maxvalScaler = (double)newMaxval/maxval; for (i = 0; i <= linearCutoff; ++i) table[i] = i * linearExpansion * maxvalScaler + 0.5; for (; i <= maxval; ++i) { double const normalized = ((double) i) / maxval; /* Xel sample value normalized to 0..1 */ double const v = 1.055 * pow(normalized, oneOverGamma) - 0.055; table[i] = MIN((xelval)(v * newMaxval + 0.5), newMaxval); /* denormalize, round, and clip */ } } static void buildSrgbGammaInverse(xelval table[], xelval const maxval, xelval const newMaxval, double const gamma) { /*---------------------------------------------------------------------------- Build a gamma table of size maxval+1 for the Inverse of the IEC SRGB gamma transfer function (Standard IEC 61966-2-1). 'gamma' must be 2.4 for true SRGB -----------------------------------------------------------------------------*/ double const oneOverGamma = 1.0 / gamma; xelval i; /* This transfer function is linear for sample values 0 .. maxval*.040405 and an exponential for larger sample values. The exponential is slightly stretched and translated, though, unlike the popular pure exponential gamma transfer function. */ xelval const linearCutoff = (xelval) maxval * 0.0031308 + 0.5; double const linearCompression = 0.0031308 / (1.055 * pow(0.0031308, oneOverGamma) - 0.055); double const maxvalScaler = (double)newMaxval/maxval; for (i = 0; i <= linearCutoff / linearCompression; ++i) table[i] = i * linearCompression * maxvalScaler + 0.5; for (; i <= maxval; ++i) { double const normalized = ((double) i) / maxval; /* Xel sample value normalized to 0..1 */ double const v = pow((normalized + 0.055) / 1.055, gamma); table[i] = MIN((xelval)(v * newMaxval + 0.5), newMaxval); /* denormalize, round, and clip */ } } static void buildBt709ToSrgbGamma(xelval table[], xelval const maxval, xelval const newMaxval, double const gammaSrgb) { /*---------------------------------------------------------------------------- Build a gamma table of size maxval+1 for the combination of the inverse of ITU Rec BT.709 and the forward SRGB gamma transfer functions. I.e. this converts from Rec 709 to SRGB. 'gammaSrgb' must be 2.4 for true SRGB. -----------------------------------------------------------------------------*/ double const oneOverGamma709 = 0.45; double const gamma709 = 1.0 / oneOverGamma709; double const oneOverGammaSrgb = 1.0 / gammaSrgb; double const normalizer = 1.0 / maxval; /* This transfer function is linear for sample values 0 .. maxval*.018 and an exponential for larger sample values. The exponential is slightly stretched and translated, though, unlike the popular pure exponential gamma transfer function. */ xelval const linearCutoff709 = (xelval) (maxval * 0.018 + 0.5); double const linearCompression709 = 0.018 / (1.099 * pow(0.018, oneOverGamma709) - 0.099); double const linearCutoffSrgb = 0.0031308; double const linearExpansionSrgb = (1.055 * pow(0.0031308, oneOverGammaSrgb) - 0.055) / 0.0031308; xelval i; for (i = 0; i <= maxval; ++i) { double const normalized = i * normalizer; /* Xel sample value normalized to 0..1 */ double radiance; double srgb; if (i < linearCutoff709 / linearCompression709) radiance = normalized * linearCompression709; else radiance = pow((normalized + 0.099) / 1.099, gamma709); assert(radiance <= 1.0); if (radiance < linearCutoffSrgb * normalizer) srgb = radiance * linearExpansionSrgb; else srgb = 1.055 * pow(radiance, oneOverGammaSrgb) - 0.055; assert(srgb <= 1.0); table[i] = srgb * newMaxval + 0.5; } } static void buildSrgbToBt709Gamma(xelval table[], xelval const maxval, xelval const newMaxval, double const gamma709) { /*---------------------------------------------------------------------------- Build a gamma table of size maxval+1 for the combination of the inverse of SRGB and the forward ITU Rec BT.709 gamma transfer functions. I.e. this converts from SRGB to Rec 709. 'gamma709' must be 1/0.45 for true Rec. 709. -----------------------------------------------------------------------------*/ double const oneOverGamma709 = 1.0 / gamma709; double const gammaSrgb = 2.4; double const oneOverGammaSrgb = 1.0 / gammaSrgb; double const normalizer = 1.0 / maxval; /* This transfer function is linear for sample values 0 .. maxval*.040405 and an exponential for larger sample values. The exponential is slightly stretched and translated, though, unlike the popular pure exponential gamma transfer function. */ xelval const linearCutoffSrgb = (xelval) maxval * 0.0031308 + 0.5; double const linearCompressionSrgb = 0.0031308 / (1.055 * pow(0.0031308, oneOverGammaSrgb) - 0.055); xelval const linearCutoff709 = (xelval) (maxval * 0.018 + 0.5); double const linearExpansion709 = (1.099 * pow(0.018, oneOverGamma709) - 0.099) / 0.018; xelval i; for (i = 0; i <= maxval; ++i) { double const normalized = i * normalizer; /* Xel sample value normalized to 0..1 */ double radiance; double bt709; if (i < linearCutoffSrgb / linearCompressionSrgb) radiance = normalized * linearCompressionSrgb; else radiance = pow((normalized + 0.055) / 1.055, gammaSrgb); assert(radiance <= 1.0); if (radiance < linearCutoff709 * normalizer) bt709 = radiance * linearExpansion709; else bt709 = 1.099 * pow(radiance, oneOverGamma709) - 0.099; assert(bt709 <= 1.0); table[i] = bt709 * newMaxval + 0.5; } } static void createGammaTables(enum transferFunction const transferFunction, xelval const maxval, xelval const newMaxval, double const rgamma, double const ggamma, double const bgamma, xelval ** const rtableP, xelval ** const gtableP, xelval ** const btableP) { /* Allocate space for the tables. */ MALLOCARRAY(*rtableP, maxval+1); MALLOCARRAY(*gtableP, maxval+1); MALLOCARRAY(*btableP, maxval+1); if (*rtableP == NULL || *gtableP == NULL || *btableP == NULL) pm_error("Can't get memory to make gamma transfer tables"); /* Build the gamma correction tables. */ switch (transferFunction) { case XF_BT709RAMP: { buildBt709Gamma(*rtableP, maxval, newMaxval, rgamma); buildBt709Gamma(*gtableP, maxval, newMaxval, ggamma); buildBt709Gamma(*btableP, maxval, newMaxval, bgamma); } break; case XF_BT709RAMP_INVERSE: { buildBt709GammaInverse(*rtableP, maxval, newMaxval, rgamma); buildBt709GammaInverse(*gtableP, maxval, newMaxval, ggamma); buildBt709GammaInverse(*btableP, maxval, newMaxval, bgamma); } break; case XF_SRGBRAMP: { buildSrgbGamma(*rtableP, maxval, newMaxval, rgamma); buildSrgbGamma(*gtableP, maxval, newMaxval, ggamma); buildSrgbGamma(*btableP, maxval, newMaxval, bgamma); } break; case XF_SRGBRAMP_INVERSE: { buildSrgbGammaInverse(*rtableP, maxval, newMaxval, rgamma); buildSrgbGammaInverse(*gtableP, maxval, newMaxval, ggamma); buildSrgbGammaInverse(*btableP, maxval, newMaxval, bgamma); } break; case XF_EXP: { buildPowGamma(*rtableP, maxval, newMaxval, rgamma); buildPowGamma(*gtableP, maxval, newMaxval, ggamma); buildPowGamma(*btableP, maxval, newMaxval, bgamma); } break; case XF_EXP_INVERSE: { buildPowGamma(*rtableP, maxval, newMaxval, 1.0/rgamma); buildPowGamma(*gtableP, maxval, newMaxval, 1.0/ggamma); buildPowGamma(*btableP, maxval, newMaxval, 1.0/bgamma); } break; case XF_BT709_TO_SRGB: { buildBt709ToSrgbGamma(*rtableP, maxval, newMaxval, rgamma); buildBt709ToSrgbGamma(*gtableP, maxval, newMaxval, ggamma); buildBt709ToSrgbGamma(*btableP, maxval, newMaxval, bgamma); } break; case XF_SRGB_TO_BT709: { buildSrgbToBt709Gamma(*rtableP, maxval, newMaxval, rgamma); buildSrgbToBt709Gamma(*gtableP, maxval, newMaxval, ggamma); buildSrgbToBt709Gamma(*btableP, maxval, newMaxval, bgamma); } break; } } static void convertRaster(FILE * const ifP, FILE * const ofP, int const cols, int const rows, xelval const maxval, int const format, xelval const outputMaxval, int const outputFormat, xelval * const rtable, xelval * const gtable, xelval * const btable) { xel * xelrow; unsigned int row; xelrow = pnm_allocrow(cols); for (row = 0; row < rows; ++row) { pnm_readpnmrow(ifP, xelrow, cols, maxval, format); pnm_promoteformatrow(xelrow, cols, maxval, format, maxval, outputFormat); switch (PNM_FORMAT_TYPE(outputFormat)) { case PPM_TYPE: { unsigned int col; for (col = 0; col < cols; ++col) { xelval const r = PPM_GETR(xelrow[col]); xelval const g = PPM_GETG(xelrow[col]); xelval const b = PPM_GETB(xelrow[col]); PPM_ASSIGN(xelrow[col], rtable[r], gtable[g], btable[b]); } } break; case PGM_TYPE: { unsigned int col; for (col = 0; col < cols; ++col) { xelval const xel = PNM_GET1(xelrow[col]); PNM_ASSIGN1(xelrow[col], gtable[xel]); } } break; default: pm_error("Internal error. Impossible format type"); } pnm_writepnmrow(ofP, xelrow, cols, outputMaxval, outputFormat, 0); } pnm_freerow(xelrow); } int main(int argc, char *argv[]) { struct cmdlineInfo cmdline; FILE * ifP; xelval maxval; int rows, cols, format; xelval outputMaxval; int outputFormat; xelval * rtable; xelval * gtable; xelval * btable; pnm_init(&argc, argv); parseCommandLine(argc, argv, &cmdline); ifP = pm_openr(cmdline.filespec); pnm_readpnminit(ifP, &cols, &rows, &maxval, &format); if (PNM_FORMAT_TYPE(format) == PPM_TYPE) outputFormat = PPM_TYPE; else if (cmdline.rgamma != cmdline.ggamma || cmdline.ggamma != cmdline.bgamma) outputFormat = PPM_TYPE; else outputFormat = PGM_TYPE; if (PNM_FORMAT_TYPE(format) != outputFormat) { if (outputFormat == PPM_TYPE) pm_message("Promoting to PPM"); if (outputFormat == PGM_TYPE) pm_message("Promoting to PGM"); } outputMaxval = cmdline.makeNewMaxval ? cmdline.maxval : maxval; createGammaTables(cmdline.transferFunction, maxval, outputMaxval, cmdline.rgamma, cmdline.ggamma, cmdline.bgamma, &rtable, >able, &btable); pnm_writepnminit(stdout, cols, rows, outputMaxval, outputFormat, 0); convertRaster(ifP, stdout, cols, rows, maxval, format, outputMaxval, outputFormat, rtable, gtable, btable); pm_close(ifP); pm_close(stdout); return 0; }