1 #include <stdbool.h>
2 #include <math.h>
3 #include <stdio.h>
4 #include <string.h>
5 #include <stdlib.h>
6 
7 #include "netpbm/pam.h"
8 #include "netpbm/pm_system.h"
9 #include "netpbm/pm_gamma.h"
10 #include "netpbm/nstring.h"
11 #include "netpbm/ppm.h"
12 
13 #include "shhopt.h"
14 #include "mallocvar.h"
15 
16 /* ----------------------------- Type aliases ------------------------------ */
17 
18 typedef unsigned char uchar;
19 typedef unsigned int  uint;
20 typedef struct   pam  pam;
21 
22 typedef struct {
23 /*----------------------------------------------------------------------------
24   An RGB triple, in linear intensity or linear brightness; user's choice.
25 -----------------------------------------------------------------------------*/
26     double _[3];
27 } Rgb;
28 
29 typedef struct {
30 /*----------------------------------------------------------------------------
31   A quadratic polynomial
32 -----------------------------------------------------------------------------*/
33     double coeff[3];
34 } Polynomial;
35 
36 typedef struct {
37 /*----------------------------------------------------------------------------
38   A set of source or target sample values, in some plane.
39 
40   These are either intensity-linear or brightness-linear; user's choice.
41 
42   There could be two or three values; user must know which.
43 -----------------------------------------------------------------------------*/
44     double _[3];
45 } SampleSet;
46 
47 /* ------------------------- Parse transformations ------------------------- */
48 
49 typedef struct {
50 /*----------------------------------------------------------------------------
51   A mapping of one source color to one target color, encoded in linear RGB
52 -----------------------------------------------------------------------------*/
53     tuplen from;
54     tuplen to;
55 } Trans;
56 
57 typedef struct {
58     const char * from;  /* color specifications  */
59     const char * to;    /*   as they appear on commandline */
60     unsigned int  hasFrom;      /* "from" part is present */
61     unsigned int  hasTo;        /* "to" part is present */
62     char  nameFromS[3]; /* short option name */
63     char  nameToS  [3];
64     char  nameFromL[6]; /* long option name */
65     char  nameToL  [6];
66 } TransArg;
67 
68 typedef struct {
69     TransArg _[3];
70 } TransArgSet;
71 
72 typedef struct {
73     unsigned int n;
74         /* Number of elements in 't', 2 for linear transformation; 3 for
75            quadratic.
76         */
77     Trans t[3];
78 } TransSet;
79 
80 typedef struct {
81     unsigned int linear;
82     unsigned int fitbrightness;
83     TransSet     xlats; /* color mappings (-from1, -to1, etc.) */
84     const char * inputFileName;  /* the input file name, "-" for stdin     */
85 } CmdlineInfo;
86 
87 
88 
89 static void
optAddTrans(optEntry * const option_def,unsigned int * const option_def_indexP,TransArg * const xP,char const index)90 optAddTrans (optEntry *     const option_def,
91              unsigned int * const option_def_indexP,
92              TransArg *     const xP,
93              char           const index) {
94 
95     char indexc;
96     uint option_def_index;
97 
98     option_def_index = *option_def_indexP;
99 
100     indexc = '0' + index;
101 
102     STRSCPY(xP->nameFromL, "from "); xP->nameFromL[4] = indexc;
103     STRSCPY(xP->nameToL,   "to "  ); xP->nameToL  [2] = indexc;
104     STRSCPY(xP->nameFromS, "f "   ); xP->nameFromS[1] = indexc;
105     STRSCPY(xP->nameToS,   "t "   ); xP->nameToS  [1] = indexc;
106 
107     OPTENT3(0, xP->nameFromL, OPT_STRING, &xP->from, &xP->hasFrom, 0);
108     OPTENT3(0, xP->nameFromS, OPT_STRING, &xP->from, &xP->hasFrom, 0);
109     OPTENT3(0, xP->nameToL,   OPT_STRING, &xP->to,   &xP->hasTo,   0);
110     OPTENT3(0, xP->nameToS,   OPT_STRING, &xP->to,   &xP->hasTo,   0);
111 
112     *option_def_indexP = option_def_index;
113 }
114 
115 
116 
117 static void
parseColor(const char * const text,tuplen * const colorP)118 parseColor(const char * const text,
119            tuplen *     const colorP) {
120 /*----------------------------------------------------------------------------
121   Parses color secification in <text>, converts it into linear RGB,
122   and stores the result in <colorP>.
123 -----------------------------------------------------------------------------*/
124     const char * const lastsc = strrchr(text, ':');
125 
126     const char * colorname;
127     double mul;
128     tuplen unmultipliedColor;
129     tuplen color;
130 
131     if (lastsc) {
132         /* Specification contains a colon.  It might be the colon that
133            introduces the optional multiplier, or it might just be the colon
134            after the type specifier, e.g. "rgbi:...".
135         */
136 
137         if (strstr(text, "rgb") == text && strchr(text, ':') == lastsc) {
138             /* The only colon present is the one on the type specifier.
139                So there is no multiplier.
140             */
141             mul = 1.0;
142             colorname = pm_strdup(text);
143         } else {
144             /* There is a multiplier (possibly invalid, though). */
145             const char * const mulstart = lastsc + 1;
146 
147             char * endP;
148             char colorbuf[50];
149 
150             errno = 0;
151             mul = strtod(mulstart, &endP);
152             if (errno != 0 || endP == mulstart)
153                 pm_error("Invalid sample multiplier: '%s'", mulstart);
154 
155             strncpy(colorbuf, text, lastsc - text);
156             colorbuf[lastsc - text] = '\0';
157             colorname = pm_strdup(colorbuf);
158         }
159     } else {
160         mul = 1.0;
161         colorname = pm_strdup(text);
162     }
163 
164     unmultipliedColor = pnm_parsecolorn(colorname);
165 
166     pm_strfree(colorname);
167 
168     MALLOCARRAY_NOFAIL(color, 3);
169 
170     {
171         /* Linearize and apply multiplier */
172         unsigned int i;
173         for (i = 0; i < 3; ++i)
174             color[i] = pm_ungamma709(unmultipliedColor[i]) * mul;
175     }
176     free(unmultipliedColor);
177 
178     *colorP = color;
179 }
180 
181 
182 
183 static void
parseTran(TransArg const transArg,Trans * const rP)184 parseTran (TransArg const transArg,
185            Trans *  const rP) {
186 
187     parseColor(transArg.from, &rP->from);
188     parseColor(transArg.to,   &rP->to);
189 }
190 
191 
192 
193 static void
calcTrans(TransArgSet const transArgs,TransSet * const transP)194 calcTrans(TransArgSet   const transArgs,
195           TransSet *    const transP) {
196 /*----------------------------------------------------------------------------
197    Interpret transformation option (-from1, etc.) values 'transArg'
198    as transformations, *transP.
199 -----------------------------------------------------------------------------*/
200     unsigned int xi;
201 
202     for (transP->n = 0, xi = 0; xi < 3; ++xi) {
203         const TransArg * const xformP = &transArgs._[xi];
204 
205         if (xformP->hasFrom || xformP->hasTo) {
206             if (!xformP->hasFrom || !xformP->hasTo)
207                 pm_error("Mapping %u incompletely specified - "
208                          "you specified -fromN or -toN but not the other",
209                     xi + 1);
210             parseTran(*xformP, &transP->t[transP->n++]);
211         }
212     }
213     if (transP->n < 2)
214         pm_error("You must specify at least two mappings with "
215                  "-from1, -to1, etc.  You specified %u", transP->n);
216 }
217 
218 
219 
220 static void
parseCommandLine(int argc,const char ** argv,CmdlineInfo * const cmdlineP)221 parseCommandLine(int                 argc,
222                  const char **       argv,
223                  CmdlineInfo * const cmdlineP) {
224 /*----------------------------------------------------------------------------
225    Parse program command line described in Unix standard form by argc
226    and argv.  Return the information in the options as *cmdlineP.
227 
228    If command line is internally inconsistent (invalid options, etc.),
229    issue error message to stderr and abort program.
230 
231    Note that the strings we return are stored in the storage that
232    was passed to us as the argv array.  We also trash *argv.
233 -----------------------------------------------------------------------------*/
234     optEntry * option_def;
235         /* Instructions to pm_optParseOptions3 on how to parse our options.
236          */
237     optStruct3 opt;
238 
239     unsigned int option_def_index;
240 
241     TransArgSet xlations; /* color mapping as read from command line */
242 
243     MALLOCARRAY_NOFAIL(option_def, 100);
244 
245     option_def_index = 0;  /* incremented by OPTENT3 */
246 
247     OPTENT3(0, "fitbrightness",          OPT_FLAG, NULL,
248             &cmdlineP->fitbrightness, 0);
249     OPTENT3(0, "linear",                 OPT_FLAG, NULL,
250             &cmdlineP->linear,        0);
251 
252     {
253         unsigned int i;
254         for (i = 0; i < 3; ++i)
255             optAddTrans(option_def, &option_def_index,
256                         &xlations._[i], i + 1);
257     }
258 
259     opt.opt_table     = option_def;
260     opt.short_allowed = 0;
261     opt.allowNegNum   = 0;
262 
263     pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);
264 
265     if (cmdlineP->linear && cmdlineP->fitbrightness) {
266         pm_error("You cannot use -linear and -fitbrightness together");
267         /* Note: It actually makes sense to use them together; we're just not
268            willing to put the effort into something it's unlikely anyone will
269            want.
270         */
271     }
272 
273     calcTrans(xlations, &cmdlineP->xlats);
274 
275     if (argc-1 < 1)
276         cmdlineP->inputFileName = "-";
277     else {
278         cmdlineP->inputFileName = argv[1];
279 
280         if (argc-1 > 1)
281             pm_error("Too many arguments.  "
282                      "The only possible non-option argument "
283                      "is the input file name");
284     }
285 
286     free(option_def);
287 }
288 
289 
290 
291 static void
errResolve(void)292 errResolve(void) {
293     pm_error( "Cannot resolve the transformations");
294 }
295 
296 
297 
298 static double
sqr(double const x)299 sqr(double const x) {
300     return x * x;
301 }
302 
303 
304 
305 static void
solveOnePlane(SampleSet const f,SampleSet const t,unsigned int const n,Polynomial * const solutionP)306 solveOnePlane(SampleSet    const f,
307               SampleSet    const t,
308               unsigned int const n,
309               Polynomial * const solutionP) {
310 /*----------------------------------------------------------------------------
311   Find the transformation that maps f[i] to t[i] for 0 <= i < n.
312 -----------------------------------------------------------------------------*/
313     double const eps = 0.00001;
314 
315     double a, b, c;
316 
317     /* I have decided against generic methods of solving systems of linear
318        equations in favour of simple explicit formulas, with no memory
319        allocation and tedious matrix processing.
320     */
321 
322     switch (n) {
323     case 3: {
324         double const aDenom =
325             sqr( f._[0] ) * ( f._[1] - f._[2] ) -
326             sqr( f._[2] ) * ( f._[1] - f._[0] ) -
327             sqr( f._[1] ) * ( f._[0] - f._[2] );
328 
329         if (fabs(aDenom) < eps)
330             errResolve();
331 
332         a = (t._[1] * (f._[2] - f._[0]) - t._[0] * (f._[2] - f._[1]) -
333              t._[2] * (f._[1] - f._[0]))
334             / aDenom;
335     } break;
336     case 2:
337         a = 0.0;
338         break;
339     default:
340         a = 0.0; /* to avoid a warning that <a> "may be uninitialized". */
341         pm_error("INTERNAL ERROR: solve(): impossible value of n: %u", n);
342     }
343 
344     {
345         double const bDenom = f._[1] - f._[0];
346 
347         if (fabs(bDenom) < eps)
348             errResolve();
349 
350         b = (t._[1] - t._[0] + a * (sqr(f._[0]) - sqr(f._[1]))) / bDenom;
351     }
352 
353     c = -a * sqr(f._[0]) - b * f._[0] + t._[0];
354 
355     solutionP->coeff[0] = a; solutionP->coeff[1] = b; solutionP->coeff[2] = c;
356 }
357 
358 
359 
360 static void
chanData(TransSet const ta,bool const fittingBrightness,unsigned int const plane,SampleSet * const fromP,SampleSet * const toP)361 chanData(TransSet     const ta,
362          bool         const fittingBrightness,
363          unsigned int const plane,
364          SampleSet *  const fromP,
365          SampleSet *  const toP) {
366 /*----------------------------------------------------------------------------
367   Collate transformations from 'ta' for plane 'plane'.
368 -----------------------------------------------------------------------------*/
369     unsigned int i;
370 
371     for (i = 0; i < ta.n; ++i) {
372         if (fittingBrightness) { /* working with gamma-compressed values */
373             fromP->_[i] = pm_gamma709(ta.t[i].from[plane]);
374             toP->  _[i] = pm_gamma709(ta.t[i].to  [plane]);
375         } else { /* working in linear RGB */
376             fromP->_[i] = ta.t[i].from[plane];
377             toP->  _[i] = ta.t[i].to  [plane];
378         }
379     }
380 }
381 
382 
383 
384 typedef struct {
385     Polynomial _[3];  /* One per plane */
386 } Solution;
387 
388 
389 
390 static void
solveFmCmdlineOpts(CmdlineInfo const cmdline,unsigned int const depth,Solution * const solutionP)391 solveFmCmdlineOpts(CmdlineInfo  const cmdline,
392                    unsigned int const depth,
393                    Solution *   const solutionP) {
394 /*----------------------------------------------------------------------------
395    Compute the function that will transform the tuples, based on what the user
396    requested ('cmdline').
397 
398    The function takes intensity-linear tuples for the normal levels function,
399    or brightness-linear for the brightness approximation levels function.
400 
401    The transformed image has 'depth' planes.
402 -----------------------------------------------------------------------------*/
403     unsigned int plane;
404     SampleSet from, to;
405     /* This initialization to bypass the "may be uninitialized" warning: */
406     to  ._[0] = 0; to.  _[1] = 0; to  ._[2] = 0;
407     from._[0] = 1; from._[1] = 0; from._[2] = 0;
408 
409     for (plane = 0; plane < depth; ++plane) {
410 
411         chanData(cmdline.xlats, cmdline.fitbrightness, plane, &from, &to);
412         solveOnePlane(from, to, cmdline.xlats.n, &solutionP->_[plane]);
413     }
414 }
415 
416 
417 
418 static samplen
xformedSample(samplen const value,Polynomial const polynomial)419 xformedSample(samplen    const value,
420               Polynomial const polynomial) {
421 /*----------------------------------------------------------------------------
422   'sample' transformed by 'polynomial'.
423 -----------------------------------------------------------------------------*/
424     double const res =
425         (polynomial.coeff[0] * value + polynomial.coeff[1]) * value +
426         polynomial.coeff[2];
427 
428     return MAX(0.0f, MIN(1.0f, res));
429 }
430 
431 
432 
433 static void
pamlevels(CmdlineInfo const cmdline)434 pamlevels(CmdlineInfo const cmdline) {
435 
436     unsigned int row;
437     pam      inPam, outPam;
438     Solution solution;
439     tuplen * tuplerown;
440     FILE   * ifP;
441 
442     ifP = pm_openr(cmdline.inputFileName);
443 
444     pnm_readpaminit(ifP, &inPam, PAM_STRUCT_SIZE(tuple_type));
445 
446     outPam = inPam;
447     outPam.file = stdout;
448 
449     solveFmCmdlineOpts(cmdline, inPam.depth, &solution);
450 
451     tuplerown = pnm_allocpamrown(&inPam);
452 
453     pnm_writepaminit(&outPam);
454 
455     for (row = 0; row < inPam.height; ++row) {
456         unsigned int col;
457 
458         pnm_readpamrown(&inPam, tuplerown);
459 
460         if (!cmdline.linear && !cmdline.fitbrightness)
461             pnm_ungammarown(&inPam, tuplerown);
462 
463         for (col = 0; col < inPam.width; ++col) {
464             unsigned int plane;
465 
466             for (plane = 0; plane < inPam.depth; ++plane) {
467                 tuplerown[col][plane] =
468                     xformedSample(tuplerown[col][plane], solution._[plane]);
469             }
470         }
471         if (!cmdline.linear && !cmdline.fitbrightness)
472             pnm_gammarown(&inPam, tuplerown);
473 
474         pnm_writepamrown(&outPam, tuplerown);
475     }
476     pnm_freepamrown(tuplerown);
477     pm_close(ifP);
478 }
479 
480 
481 
482 static void
freeCmdLineInfo(CmdlineInfo cmdline)483 freeCmdLineInfo(CmdlineInfo cmdline) {
484 /*----------------------------------------------------------------------------
485   Free any memory that has been dynamically allocated in <cmdline>.
486 -----------------------------------------------------------------------------*/
487     TransSet * const xxP = &cmdline.xlats;
488 
489     uint x;
490 
491     for (x = 0; x < xxP->n; ++x) {
492         free(xxP->t[x].from);
493         free(xxP->t[x].to);
494     }
495 }
496 
497 
498 
main(int argc,const char * argv[])499 int main(int    argc, const char * argv[]) {
500 
501     CmdlineInfo cmdline;
502 
503     pm_proginit(&argc, argv);
504 
505     parseCommandLine(argc, argv, &cmdline);
506 
507     pamlevels(cmdline);
508 
509     freeCmdLineInfo(cmdline);
510 
511     return 0;
512 }
513 
514 
515 
516