1 /* pgmtopbm.c - read a portable graymap and write a portable bitmap
2 **
3 ** Copyright (C) 1989 by Jef Poskanzer.
4 **
5 ** Permission to use, copy, modify, and distribute this software and its
6 ** documentation for any purpose and without fee is hereby granted, provided
7 ** that the above copyright notice appear in all copies and that both that
8 ** copyright notice and this permission notice appear in supporting
9 ** documentation. This software is provided "as is" without express or
10 ** implied warranty.
11 */
12
13 #include <assert.h>
14
15 #include "pm_c_util.h"
16 #include "shhopt.h"
17 #include "pgm.h"
18 #include "dithers.h"
19 #include "mallocvar.h"
20
21 enum halftone {QT_FS, QT_THRESH, QT_DITHER8, QT_CLUSTER, QT_HILBERT};
22
23
24 struct cmdlineInfo {
25 /* All the information the user supplied in the command line,
26 in a form easy for the program to use.
27 */
28 const char * inputFilespec;
29 enum halftone halftone;
30 unsigned int clumpSize;
31 unsigned int clusterRadius;
32 /* Defined only for halftone == QT_CLUSTER */
33 float threshval;
34 };
35
36
37
38
39 static void
parseCommandLine(int argc,char ** argv,struct cmdlineInfo * cmdlineP)40 parseCommandLine(int argc, char ** argv,
41 struct cmdlineInfo *cmdlineP) {
42 /*----------------------------------------------------------------------------
43 Note that the file spec array we return is stored in the storage that
44 was passed to us as the argv array.
45 -----------------------------------------------------------------------------*/
46 optEntry *option_def;
47 /* Instructions to pm_optParseOptions3 on how to parse our options.
48 */
49 optStruct3 opt;
50
51 unsigned int option_def_index;
52 unsigned int floydOpt, hilbertOpt, thresholdOpt, dither8Opt,
53 cluster3Opt, cluster4Opt, cluster8Opt;
54 unsigned int valueSpec, clumpSpec;
55
56 MALLOCARRAY_NOFAIL(option_def, 100);
57
58 option_def_index = 0; /* incremented by OPTENTRY */
59 OPTENT3(0, "floyd", OPT_FLAG, NULL, &floydOpt, 0);
60 OPTENT3(0, "fs", OPT_FLAG, NULL, &floydOpt, 0);
61 OPTENT3(0, "threshold", OPT_FLAG, NULL, &thresholdOpt, 0);
62 OPTENT3(0, "hilbert", OPT_FLAG, NULL, &hilbertOpt, 0);
63 OPTENT3(0, "dither8", OPT_FLAG, NULL, &dither8Opt, 0);
64 OPTENT3(0, "d8", OPT_FLAG, NULL, &dither8Opt, 0);
65 OPTENT3(0, "cluster3", OPT_FLAG, NULL, &cluster3Opt, 0);
66 OPTENT3(0, "c3", OPT_FLAG, NULL, &cluster3Opt, 0);
67 OPTENT3(0, "cluster4", OPT_FLAG, NULL, &cluster4Opt, 0);
68 OPTENT3(0, "c4", OPT_FLAG, NULL, &cluster4Opt, 0);
69 OPTENT3(0, "cluster8", OPT_FLAG, NULL, &cluster8Opt, 0);
70 OPTENT3(0, "c8", OPT_FLAG, NULL, &cluster8Opt, 0);
71 OPTENT3(0, "value", OPT_FLOAT, &cmdlineP->threshval,
72 &valueSpec, 0);
73 OPTENT3(0, "clump", OPT_UINT, &cmdlineP->clumpSize,
74 &clumpSpec, 0);
75
76 opt.opt_table = option_def;
77 opt.short_allowed = FALSE; /* We have no short (old-fashioned) options */
78 opt.allowNegNum = FALSE; /* We may have parms that are negative numbers */
79
80 pm_optParseOptions3(&argc, argv, opt, sizeof(opt), 0);
81 /* Uses and sets argc, argv, and some of *cmdlineP and others. */
82
83 if (floydOpt + thresholdOpt + hilbertOpt + dither8Opt +
84 cluster3Opt + cluster4Opt + cluster8Opt == 0)
85 cmdlineP->halftone = QT_FS;
86 else if (floydOpt + thresholdOpt + dither8Opt +
87 cluster3Opt + cluster4Opt + cluster8Opt > 1)
88 pm_error("No cannot specify more than one halftoning type");
89 else {
90 if (floydOpt)
91 cmdlineP->halftone = QT_FS;
92 else if (thresholdOpt)
93 cmdlineP->halftone = QT_THRESH;
94 else if (hilbertOpt)
95 cmdlineP->halftone = QT_HILBERT;
96 else if (dither8Opt)
97 cmdlineP->halftone = QT_DITHER8;
98 else if (cluster3Opt) {
99 cmdlineP->halftone = QT_CLUSTER;
100 cmdlineP->clusterRadius = 3;
101 } else if (cluster4Opt) {
102 cmdlineP->halftone = QT_CLUSTER;
103 cmdlineP->clusterRadius = 4;
104 } else if (cluster8Opt) {
105 cmdlineP->halftone = QT_CLUSTER;
106 cmdlineP->clusterRadius = 8;
107 } else
108 pm_error("INTERNAL ERROR. No halftone option");
109 }
110
111 if (!valueSpec)
112 cmdlineP->threshval = 0.5;
113 else {
114 if (cmdlineP->threshval < 0.0)
115 pm_error("-value cannot be negative. You specified %f",
116 cmdlineP->threshval);
117 if (cmdlineP->threshval > 1.0)
118 pm_error("-value cannot be greater than one. You specified %f",
119 cmdlineP->threshval);
120 }
121
122 if (!clumpSpec)
123 cmdlineP->clumpSize = 5;
124 else {
125 if (cmdlineP->clumpSize < 2)
126 pm_error("-clump must be at least 2. You specified %u",
127 cmdlineP->clumpSize);
128 }
129
130 if (argc-1 > 1)
131 pm_error("Too many arguments (%d). There is at most one "
132 "non-option argument: the file name",
133 argc-1);
134 else if (argc-1 == 1)
135 cmdlineP->inputFilespec = argv[1];
136 else
137 cmdlineP->inputFilespec = "-";
138 }
139
140
141
142 /* Hilbert curve tracer */
143
144 #define MAXORD 18
145
146 static int hil_order,hil_ord;
147 static int hil_turn;
148 static int hil_dx,hil_dy;
149 static int hil_x,hil_y;
150 static int hil_stage[MAXORD];
151 static int hil_width,hil_height;
152
153 static void
init_hilbert(int const w,int const h)154 init_hilbert(int const w,
155 int const h) {
156 /*----------------------------------------------------------------------------
157 Initialize the Hilbert curve tracer
158 -----------------------------------------------------------------------------*/
159 int big,ber;
160 hil_width = w;
161 hil_height = h;
162 big = w > h ? w : h;
163 for (ber = 2, hil_order = 1; ber < big; ber <<= 1, hil_order++);
164 if (hil_order > MAXORD)
165 pm_error("Sorry, hilbert order is too large");
166 hil_ord = hil_order;
167 hil_order--;
168 }
169
170
171
172 static int
hilbert(int * const px,int * const py)173 hilbert(int * const px, int * const py) {
174 /*----------------------------------------------------------------------------
175 Return non-zero if got another point
176 -----------------------------------------------------------------------------*/
177 int temp;
178 if (hil_ord > hil_order) {
179 /* have to do first point */
180
181 hil_ord--;
182 hil_stage[hil_ord] = 0;
183 hil_turn = -1;
184 hil_dy = 1;
185 hil_dx = hil_x = hil_y = 0;
186 *px = *py = 0;
187 return 1;
188 }
189
190 /* Operate the state machine */
191 for(;;) {
192 switch (hil_stage[hil_ord]) {
193 case 0:
194 hil_turn = -hil_turn;
195 temp = hil_dy;
196 hil_dy = -hil_turn * hil_dx;
197 hil_dx = hil_turn * temp;
198 if (hil_ord > 0) {
199 hil_stage[hil_ord] = 1;
200 hil_ord--;
201 hil_stage[hil_ord]=0;
202 continue;
203 }
204 case 1:
205 hil_x += hil_dx;
206 hil_y += hil_dy;
207 if (hil_x < hil_width && hil_y < hil_height) {
208 hil_stage[hil_ord] = 2;
209 *px = hil_x;
210 *py = hil_y;
211 return 1;
212 }
213 case 2:
214 hil_turn = -hil_turn;
215 temp = hil_dy;
216 hil_dy = -hil_turn * hil_dx;
217 hil_dx = hil_turn * temp;
218 if (hil_ord > 0) {
219 /* recurse */
220
221 hil_stage[hil_ord] = 3;
222 hil_ord--;
223 hil_stage[hil_ord]=0;
224 continue;
225 }
226 case 3:
227 hil_x += hil_dx;
228 hil_y += hil_dy;
229 if (hil_x < hil_width && hil_y < hil_height) {
230 hil_stage[hil_ord] = 4;
231 *px = hil_x;
232 *py = hil_y;
233 return 1;
234 }
235 case 4:
236 if (hil_ord > 0) {
237 /* recurse */
238 hil_stage[hil_ord] = 5;
239 hil_ord--;
240 hil_stage[hil_ord]=0;
241 continue;
242 }
243 case 5:
244 temp = hil_dy;
245 hil_dy = -hil_turn * hil_dx;
246 hil_dx = hil_turn * temp;
247 hil_turn = -hil_turn;
248 hil_x += hil_dx;
249 hil_y += hil_dy;
250 if (hil_x < hil_width && hil_y < hil_height) {
251 hil_stage[hil_ord] = 6;
252 *px = hil_x;
253 *py = hil_y;
254 return 1;
255 }
256 case 6:
257 if (hil_ord > 0) {
258 /* recurse */
259 hil_stage[hil_ord] = 7;
260 hil_ord--;
261 hil_stage[hil_ord]=0;
262 continue;
263 }
264 case 7:
265 temp = hil_dy;
266 hil_dy = -hil_turn * hil_dx;
267 hil_dx = hil_turn * temp;
268 hil_turn = -hil_turn;
269 /* Return from a recursion */
270 if (hil_ord < hil_order)
271 hil_ord++;
272 else
273 return 0;
274 }
275 }
276 }
277
278
279
doHilbert(FILE * const ifP,unsigned int const clump_size)280 static void doHilbert(FILE * const ifP,
281 unsigned int const clump_size) {
282 /*----------------------------------------------------------------------------
283 Use hilbert space filling curve dithering
284 -----------------------------------------------------------------------------*/
285 /*
286 * This is taken from the article "Digital Halftoning with
287 * Space Filling Curves" by Luiz Velho, proceedings of
288 * SIGRAPH '91, page 81.
289 *
290 * This is not a terribly efficient or quick version of
291 * this algorithm, but it seems to work. - Graeme Gill.
292 * graeme@labtam.labtam.OZ.AU
293 *
294 */
295
296 int cols, rows;
297 gray maxval;
298 gray **grays;
299 bit **bits;
300 int end;
301 int *x,*y;
302 int sum = 0;
303
304 grays = pgm_readpgm(ifP, &cols,&rows, &maxval);
305 bits = pbm_allocarray(cols,rows);
306
307 MALLOCARRAY(x, clump_size);
308 MALLOCARRAY(y, clump_size);
309 if (x == NULL || y == NULL)
310 pm_error("out of memory");
311 init_hilbert(cols,rows);
312
313 end = clump_size;
314 while (end == clump_size) {
315 int i;
316 /* compute the next clust co-ordinates along hilbert path */
317 for (i = 0; i < end; i++) {
318 if (hilbert(&x[i],&y[i])==0)
319 end = i; /* we reached the end */
320 }
321 /* sum levels */
322 for (i = 0; i < end; i++)
323 sum += grays[y[i]][x[i]];
324 /* dither half and half along path */
325 for (i = 0; i < end; i++) {
326 if (sum >= maxval) {
327 bits[y[i]][x[i]] = PBM_WHITE;
328 sum -= maxval;
329 } else
330 bits[y[i]][x[i]] = PBM_BLACK;
331 }
332 }
333 pbm_writepbm(stdout, bits, cols, rows, 0);
334 }
335
336
337
338 struct converter {
339 void (*convertRow)(struct converter * const converterP,
340 unsigned int const row,
341 gray grayrow[],
342 bit bitrow[]);
343 void (*destroy)(struct converter * const converterP);
344 unsigned int cols;
345 gray maxval;
346 void * stateP;
347 };
348
349
350
351 /*=============================================================================
352 Converter: fs
353 =============================================================================*/
354
355 unsigned int const fs_scale = 1024;
356 unsigned int const half_fs_scale = 512;
357
358 struct fsState {
359 long* thiserr;
360 long* nexterr;
361 bool fs_forward;
362 long threshval; /* Threshold gray value, scaled by FS_SCALE */
363 };
364
365
366 static void
fsConvertRow(struct converter * const converterP,unsigned int const row,gray grayrow[],bit bitrow[])367 fsConvertRow(struct converter * const converterP,
368 unsigned int const row,
369 gray grayrow[],
370 bit bitrow[]) {
371
372 struct fsState * const stateP = converterP->stateP;
373
374 long * const thiserr = stateP->thiserr;
375 long * const nexterr = stateP->nexterr;
376
377 bit* bP;
378 gray* gP;
379
380 unsigned int limitcol;
381 unsigned int col;
382
383 for (col = 0; col < converterP->cols + 2; ++col)
384 nexterr[col] = 0;
385 if (stateP->fs_forward) {
386 col = 0;
387 limitcol = converterP->cols;
388 gP = grayrow;
389 bP = bitrow;
390 } else {
391 col = converterP->cols - 1;
392 limitcol = -1;
393 gP = &(grayrow[col]);
394 bP = &(bitrow[col]);
395 }
396 do {
397 long sum;
398 sum = ((long) *gP * fs_scale) / converterP->maxval +
399 thiserr[col + 1];
400 if (sum >= stateP->threshval) {
401 *bP = PBM_WHITE;
402 sum = sum - stateP->threshval - half_fs_scale;
403 } else
404 *bP = PBM_BLACK;
405
406 if (stateP->fs_forward) {
407 thiserr[col + 2] += (sum * 7) / 16;
408 nexterr[col ] += (sum * 3) / 16;
409 nexterr[col + 1] += (sum * 5) / 16;
410 nexterr[col + 2] += (sum ) / 16;
411
412 ++col;
413 ++gP;
414 ++bP;
415 } else {
416 thiserr[col ] += (sum * 7) / 16;
417 nexterr[col + 2] += (sum * 3) / 16;
418 nexterr[col + 1] += (sum * 5) / 16;
419 nexterr[col ] += (sum ) / 16;
420
421 --col;
422 --gP;
423 --bP;
424 }
425 } while (col != limitcol);
426
427 stateP->thiserr = nexterr;
428 stateP->nexterr = thiserr;
429 stateP->fs_forward = ! stateP->fs_forward;
430 }
431
432
433
434 static void
fsDestroy(struct converter * const converterP)435 fsDestroy(struct converter * const converterP) {
436 free(converterP->stateP);
437 }
438
439
440
441 static struct converter
createFsConverter(unsigned int const cols,gray const maxval,float const threshFraction)442 createFsConverter(unsigned int const cols,
443 gray const maxval,
444 float const threshFraction) {
445
446 struct fsState * stateP;
447 struct converter converter;
448
449 MALLOCVAR_NOFAIL(stateP);
450
451 /* Initialize Floyd-Steinberg error vectors. */
452 MALLOCARRAY_NOFAIL(stateP->thiserr, cols + 2);
453 MALLOCARRAY_NOFAIL(stateP->nexterr, cols + 2);
454 srand(pm_randseed());
455
456 {
457 /* (random errors in [-fs_scale/8 .. fs_scale/8]) */
458 unsigned int col;
459 for (col = 0; col < cols + 2; ++col)
460 stateP->thiserr[col] =
461 (long)(rand() % fs_scale - half_fs_scale) / 4;
462 }
463
464 stateP->fs_forward = TRUE;
465 stateP->threshval = threshFraction * fs_scale;
466 converter.stateP = stateP;
467 converter.cols = cols;
468 converter.maxval = maxval;
469 converter.convertRow = &fsConvertRow;
470 converter.destroy = &fsDestroy;
471
472 return converter;
473 }
474
475
476
477 /*=============================================================================
478 Converter: thresh
479 =============================================================================*/
480
481 struct threshState {
482 gray threshval;
483 };
484
485
486 static void
threshConvertRow(struct converter * const converterP,unsigned int const row,gray grayrow[],bit bitrow[])487 threshConvertRow(struct converter * const converterP,
488 unsigned int const row,
489 gray grayrow[],
490 bit bitrow[]) {
491
492 struct threshState * const stateP = converterP->stateP;
493
494 unsigned int col;
495 for (col = 0; col < converterP->cols; ++col)
496 if (grayrow[col] >= stateP->threshval)
497 bitrow[col] = PBM_WHITE;
498 else
499 bitrow[col] = PBM_BLACK;
500 }
501
502
503
504 static void
threshDestroy(struct converter * const converterP)505 threshDestroy(struct converter * const converterP) {
506 free(converterP->stateP);
507 }
508
509
510
511 static struct converter
createThreshConverter(unsigned int const cols,gray const maxval,float const threshFraction)512 createThreshConverter(unsigned int const cols,
513 gray const maxval,
514 float const threshFraction) {
515
516 struct threshState * stateP;
517 struct converter converter;
518
519 MALLOCVAR_NOFAIL(stateP);
520
521 converter.cols = cols;
522 converter.maxval = maxval;
523 converter.convertRow = &threshConvertRow;
524 converter.destroy = &threshDestroy;
525
526 stateP->threshval = ROUNDU(maxval * threshFraction);
527 converter.stateP = stateP;
528
529 return converter;
530 }
531
532
533
534 /*=============================================================================
535 Converter: dither8
536 =============================================================================*/
537
538 struct dither8State {
539 int dither8[16][16];
540 };
541
542
543
544 static void
dither8ConvertRow(struct converter * const converterP,unsigned int const row,gray grayrow[],bit bitrow[])545 dither8ConvertRow(struct converter * const converterP,
546 unsigned int const row,
547 gray grayrow[],
548 bit bitrow[]) {
549
550 struct dither8State * const stateP = converterP->stateP;
551
552 unsigned int col;
553
554 for (col = 0; col < converterP->cols; ++col)
555 if (grayrow[col] > stateP->dither8[row % 16][col % 16])
556 bitrow[col] = PBM_WHITE;
557 else
558 bitrow[col] = PBM_BLACK;
559 }
560
561
562
563 static void
dither8Destroy(struct converter * const converterP)564 dither8Destroy(struct converter * const converterP) {
565
566 struct dither8State * const stateP = converterP->stateP;
567
568 free(stateP);
569 }
570
571
572
573 static struct converter
createDither8Converter(unsigned int const cols,gray const maxval)574 createDither8Converter(unsigned int const cols,
575 gray const maxval) {
576
577 struct converter converter;
578 struct dither8State * stateP;
579
580 unsigned int row;
581
582 MALLOCVAR_NOFAIL(stateP);
583
584 converter.cols = cols;
585 converter.convertRow = &dither8ConvertRow;
586 converter.destroy = dither8Destroy;
587 converter.stateP = stateP;
588 converter.maxval = maxval;
589
590 /* Scale dither matrix. */
591 for (row = 0; row < 16; ++row) {
592 unsigned int col;
593 for (col = 0; col < 16; ++col)
594 stateP->dither8[row][col] = dither8[row][col] * maxval / 256;
595 }
596 return converter;
597 }
598
599
600
601 /*=============================================================================
602 Converter: cluster
603 =============================================================================*/
604
605 struct clusterState {
606 unsigned int radius;
607 int ** clusterMatrix;
608 };
609
610
611
612 static void
clusterConvertRow(struct converter * const converterP,unsigned int const row,gray grayrow[],bit bitrow[])613 clusterConvertRow(struct converter * const converterP,
614 unsigned int const row,
615 gray grayrow[],
616 bit bitrow[]) {
617
618 struct clusterState * const stateP = converterP->stateP;
619 unsigned int const diameter = 2 * stateP->radius;
620
621 unsigned int col;
622
623 for (col = 0; col < converterP->cols; ++col)
624 if (grayrow[col] >
625 stateP->clusterMatrix[row % diameter][col % diameter])
626 bitrow[col] = PBM_WHITE;
627 else
628 bitrow[col] = PBM_BLACK;
629 }
630
631
632
633 static void
clusterDestroy(struct converter * const converterP)634 clusterDestroy(struct converter * const converterP) {
635
636 struct clusterState * const stateP = converterP->stateP;
637 unsigned int const diameter = 2 * stateP->radius;
638
639 unsigned int row;
640
641 for (row = 0; row < diameter; ++row)
642 free(stateP->clusterMatrix[row]);
643
644 free(stateP->clusterMatrix);
645
646 free(stateP);
647 }
648
649
650
651 static struct converter
createClusterConverter(unsigned int const radius,unsigned int const cols,gray const maxval)652 createClusterConverter(unsigned int const radius,
653 unsigned int const cols,
654 gray const maxval) {
655
656 int const clusterNormalizer = radius * radius * 2;
657 unsigned int const diameter = 2 * radius;
658
659 struct converter converter;
660 struct clusterState * stateP;
661 unsigned int row;
662
663 converter.cols = cols;
664 converter.maxval = maxval;
665 converter.convertRow = &clusterConvertRow;
666 converter.destroy = &clusterDestroy;
667
668 MALLOCVAR_NOFAIL(stateP);
669
670 stateP->radius = radius;
671
672 MALLOCARRAY_NOFAIL(stateP->clusterMatrix, diameter);
673 for (row = 0; row < diameter; ++row) {
674 unsigned int col;
675
676 MALLOCARRAY_NOFAIL(stateP->clusterMatrix[row], diameter);
677
678 for (col = 0; col < diameter; ++col) {
679 int val;
680 switch (radius) {
681 case 3: val = cluster3[row][col]; break;
682 case 4: val = cluster4[row][col]; break;
683 case 8: val = cluster8[row][col]; break;
684 default:
685 pm_error("INTERNAL ERROR: invalid radius");
686 }
687 stateP->clusterMatrix[row][col] = val * maxval / clusterNormalizer;
688 }
689 }
690
691 converter.stateP = stateP;
692
693 return converter;
694 }
695
696
697
698 int
main(int argc,char * argv[])699 main(int argc, char *argv[]) {
700
701 struct cmdlineInfo cmdline;
702 FILE* ifP;
703 gray* grayrow;
704 bit* bitrow;
705
706 pgm_init(&argc, argv);
707
708 parseCommandLine(argc, argv, &cmdline);
709
710 ifP = pm_openr(cmdline.inputFilespec);
711
712 if (cmdline.halftone == QT_HILBERT)
713 doHilbert(ifP, cmdline.clumpSize);
714 else {
715 struct converter converter;
716 int cols, rows;
717 gray maxval;
718 int format;
719 int row;
720
721 pgm_readpgminit(ifP, &cols, &rows, &maxval, &format);
722
723 pbm_writepbminit(stdout, cols, rows, 0);
724
725 switch (cmdline.halftone) {
726 case QT_FS:
727 converter = createFsConverter(cols, maxval, cmdline.threshval);
728 break;
729 case QT_THRESH:
730 converter = createThreshConverter(cols, maxval, cmdline.threshval);
731 break;
732 case QT_DITHER8:
733 converter = createDither8Converter(cols, maxval);
734 break;
735 case QT_CLUSTER:
736 converter =
737 createClusterConverter(cmdline.clusterRadius, cols, maxval);
738 break;
739 case QT_HILBERT:
740 pm_error("INTERNAL ERROR: halftone is QT_HILBERT where it "
741 "shouldn't be.");
742 break;
743 }
744
745 grayrow = pgm_allocrow(cols);
746 bitrow = pbm_allocrow(cols);
747
748 for (row = 0; row < rows; ++row) {
749 pgm_readpgmrow(ifP, grayrow, cols, maxval, format);
750
751 converter.convertRow(&converter, row, grayrow, bitrow);
752
753 pbm_writepbmrow(stdout, bitrow, cols, 0);
754 }
755 pbm_freerow(bitrow);
756 pgm_freerow(grayrow);
757
758 if (converter.destroy)
759 converter.destroy(&converter);
760 }
761
762 pm_close(ifP);
763
764 return 0;
765 }
766