home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / netpbma.zip / pgm / pgmtopbm.c < prev    next >
C/C++ Source or Header  |  1994-02-18  |  11KB  |  478 lines
  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 "pgm.h"
  14. #include "dithers.h"
  15.  
  16. static void init_hilbert ARGS((int w, int h));
  17. static int hilbert ARGS((int *px, int *py));
  18.  
  19.  
  20. int
  21. main( argc, argv )
  22.     int argc;
  23.     char* argv[];
  24. {
  25.     FILE* ifp;
  26.     gray* grayrow;
  27.     register gray* gP;
  28.     bit* bitrow;
  29.     register bit* bP;
  30.     int argn, rows, cols, format, row, col, limitcol;
  31.     float fthreshval;
  32.     int clump_size;
  33.     gray maxval;
  34.     char* usage = "[-floyd|-fs  | -hilbert | -threshold | -dither8|-d8 |\n     -cluster3|-c3|-cluster4|-c4|-cluster8|-c8] [-value <val>] [-clump <size>] [pgmfile]";
  35.     int halftone;
  36. #define QT_FS 1
  37. #define QT_THRESH 2
  38. #define QT_DITHER8 3
  39. #define QT_CLUSTER3 4
  40. #define QT_CLUSTER4 5
  41. #define QT_CLUSTER8 6
  42. #define QT_HILBERT 7
  43.     long threshval, sum;
  44.     long* thiserr;
  45.     long* nexterr;
  46.     long* temperr;
  47. #define FS_SCALE 1024
  48. #define HALF_FS_SCALE 512
  49.     int fs_direction;
  50.  
  51.     pgm_init( &argc, argv );
  52.  
  53.     argn = 1;
  54.     halftone = QT_FS;    /* default quantization is Floyd-Steinberg */
  55.     fthreshval = 0.5;
  56.     clump_size = 5;    /* default hilbert curve clump size */
  57.  
  58.     while ( argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0' )
  59.     {
  60.     if ( pm_keymatch( argv[argn], "-fs", 2 ) ||
  61.         pm_keymatch( argv[argn], "-floyd", 2 ) )
  62.         halftone = QT_FS;
  63.     else if ( pm_keymatch( argv[argn], "-hilbert", 2 ))
  64.         halftone = QT_HILBERT;
  65.     else if ( pm_keymatch( argv[argn], "-threshold", 2 ) )
  66.         halftone = QT_THRESH;
  67.     else if ( pm_keymatch( argv[argn], "-dither8", 2 ) ||
  68.           pm_keymatch( argv[argn], "-d8", 3 ) )
  69.         halftone = QT_DITHER8;
  70.     else if ( pm_keymatch( argv[argn], "-cluster3", 9 ) ||
  71.           pm_keymatch( argv[argn], "-c3", 3 ) )
  72.         halftone = QT_CLUSTER3;
  73.     else if ( pm_keymatch( argv[argn], "-cluster4", 9 ) ||
  74.               pm_keymatch( argv[argn], "-c4", 3 ) )
  75.         halftone = QT_CLUSTER4;
  76.     else if ( pm_keymatch( argv[argn], "-cluster8", 9 ) ||
  77.               pm_keymatch( argv[argn], "-c8", 3 ) )
  78.         halftone = QT_CLUSTER8;
  79.     else if ( pm_keymatch( argv[argn], "-value", 2 ) )
  80.     {
  81.         ++argn;
  82.         if ( argn == argc || sscanf( argv[argn], "%f", &fthreshval ) != 1 ||
  83.         fthreshval < 0.0 || fthreshval > 1.0 )
  84.         pm_usage( usage );
  85.     }
  86.     else if ( pm_keymatch( argv[argn], "-clump", 2 ) )
  87.     {
  88.         ++argn;
  89.         if ( argn == argc || sscanf( argv[argn], "%d", &clump_size ) != 1 ||
  90.         clump_size < 2)
  91.         pm_usage( usage );
  92.     }
  93.     else
  94.         pm_usage( usage );
  95.     ++argn;
  96.     }
  97.  
  98.     if ( argn != argc )
  99.     {
  100.     ifp = pm_openr( argv[argn] );
  101.     ++argn;
  102.     }
  103.     else
  104.     ifp = stdin;
  105.  
  106.     if ( argn != argc )
  107.     pm_usage( usage );
  108.  
  109.     if ( halftone != QT_HILBERT )
  110.     {
  111.     pgm_readpgminit( ifp, &cols, &rows, &maxval, &format );
  112.     grayrow = pgm_allocrow( cols );
  113.  
  114.     pbm_writepbminit( stdout, cols, rows, 0 );
  115.     bitrow = pbm_allocrow( cols );
  116.  
  117.     /* Initialize. */
  118.     switch ( halftone )
  119.     {
  120.     case QT_FS:
  121.         /* Initialize Floyd-Steinberg error vectors. */
  122.         thiserr = (long*) pm_allocrow( cols + 2, sizeof(long) );
  123.         nexterr = (long*) pm_allocrow( cols + 2, sizeof(long) );
  124.         srandom( (int) ( time( 0 ) ^ getpid( ) ) );
  125.         for ( col = 0; col < cols + 2; ++col )
  126.         thiserr[col] = ( random( ) % FS_SCALE - HALF_FS_SCALE ) / 4;
  127.         /* (random errors in [-FS_SCALE/8 .. FS_SCALE/8]) */
  128.         fs_direction = 1;
  129.         threshval = fthreshval * FS_SCALE;
  130.         break;
  131.  
  132.     case QT_HILBERT:
  133.         break;
  134.  
  135.     case QT_THRESH:
  136.         threshval = fthreshval * maxval + 0.999999;
  137.         break;
  138.  
  139.     case QT_DITHER8:
  140.         /* Scale dither matrix. */
  141.         for ( row = 0; row < 16; ++row )
  142.         for ( col = 0; col < 16; ++col )
  143.             dither8[row][col] = dither8[row][col] * ( maxval + 1 ) / 256;
  144.         break;
  145.  
  146.     case QT_CLUSTER3:
  147.         /* Scale order-3 clustered dither matrix. */
  148.         for ( row = 0; row < 6; ++row )
  149.         for ( col = 0; col < 6; ++col )
  150.             cluster3[row][col] = cluster3[row][col] * ( maxval + 1 ) / 18;
  151.         break;
  152.  
  153.     case QT_CLUSTER4:
  154.         /* Scale order-4 clustered dither matrix. */
  155.         for ( row = 0; row < 8; ++row )
  156.         for ( col = 0; col < 8; ++col )
  157.             cluster4[row][col] = cluster4[row][col] * ( maxval + 1 ) / 32;
  158.         break;
  159.  
  160.     case QT_CLUSTER8:
  161.         /* Scale order-8 clustered dither matrix. */
  162.         for ( row = 0; row < 16; ++row )
  163.         for ( col = 0; col < 16; ++col )
  164.             cluster8[row][col] = cluster8[row][col] * ( maxval + 1 ) / 128;
  165.         break;
  166.  
  167.     default:
  168.         pm_error( "can't happen" );
  169.         exit( 1 );
  170.     }
  171.  
  172.     for ( row = 0; row < rows; ++row )
  173.     {
  174.         pgm_readpgmrow( ifp, grayrow, cols, maxval, format );
  175.     
  176.         switch ( halftone )
  177.         {
  178.         case QT_FS:
  179.         for ( col = 0; col < cols + 2; ++col )
  180.             nexterr[col] = 0;
  181.         if ( fs_direction )
  182.         {
  183.             col = 0;
  184.             limitcol = cols;
  185.             gP = grayrow;
  186.             bP = bitrow;
  187.         }
  188.         else
  189.         {
  190.             col = cols - 1;
  191.             limitcol = -1;
  192.             gP = &(grayrow[col]);
  193.             bP = &(bitrow[col]);
  194.         }
  195.         do
  196.         {
  197.             sum = ( (long) *gP * FS_SCALE ) / maxval + thiserr[col + 1];
  198.             if ( sum >= threshval )
  199.             {
  200.             *bP = PBM_WHITE;
  201.             sum = sum - threshval - HALF_FS_SCALE;
  202.             }
  203.             else
  204.             *bP = PBM_BLACK;
  205.     
  206.             if ( fs_direction )
  207.             {
  208.             thiserr[col + 2] += ( sum * 7 ) / 16;
  209.             nexterr[col    ] += ( sum * 3 ) / 16;
  210.             nexterr[col + 1] += ( sum * 5 ) / 16;
  211.             nexterr[col + 2] += ( sum     ) / 16;
  212.     
  213.             ++col;
  214.             ++gP;
  215.             ++bP;
  216.             }
  217.             else
  218.             {
  219.             thiserr[col    ] += ( sum * 7 ) / 16;
  220.             nexterr[col + 2] += ( sum * 3 ) / 16;
  221.             nexterr[col + 1] += ( sum * 5 ) / 16;
  222.             nexterr[col    ] += ( sum     ) / 16;
  223.     
  224.             --col;
  225.             --gP;
  226.             --bP;
  227.             }
  228.         }
  229.         while ( col != limitcol );
  230.         temperr = thiserr;
  231.         thiserr = nexterr;
  232.         nexterr = temperr;
  233.         fs_direction = ! fs_direction;
  234.         break;
  235.     
  236.         case QT_THRESH:
  237.         for ( col = 0, gP = grayrow, bP = bitrow; col < cols; ++col, ++gP, ++bP )
  238.             if ( *gP >= threshval )
  239.             *bP = PBM_WHITE;
  240.             else
  241.             *bP = PBM_BLACK;
  242.         break;
  243.     
  244.         case QT_DITHER8:
  245.         for ( col = 0, gP = grayrow, bP = bitrow; col < cols; ++col, ++gP, ++bP )
  246.             if ( *gP >= dither8[row % 16][col % 16] )
  247.             *bP = PBM_WHITE;
  248.             else
  249.             *bP = PBM_BLACK;
  250.         break;
  251.     
  252.         case QT_CLUSTER3:
  253.         for ( col = 0, gP = grayrow, bP = bitrow; col < cols; ++col, ++gP, ++bP )
  254.             if ( *gP >= cluster3[row % 6][col % 6] )
  255.             *bP = PBM_WHITE;
  256.             else
  257.             *bP = PBM_BLACK;
  258.         break;
  259.     
  260.         case QT_CLUSTER4:
  261.         for ( col = 0, gP = grayrow, bP = bitrow; col < cols; ++col, ++gP, ++bP )
  262.             if ( *gP >= cluster4[row % 8][col % 8] )
  263.             *bP = PBM_WHITE;
  264.             else
  265.             *bP = PBM_BLACK;
  266.         break;
  267.     
  268.         case QT_CLUSTER8:
  269.         for ( col = 0, gP = grayrow, bP = bitrow; col < cols; ++col, ++gP, ++bP )
  270.             if ( *gP >= cluster8[row % 16][col % 16] )
  271.             *bP = PBM_WHITE;
  272.             else
  273.             *bP = PBM_BLACK;
  274.         break;
  275.     
  276.         default:
  277.         pm_error( "can't happen" );
  278.         exit( 1 );
  279.         }
  280.     
  281.         pbm_writepbmrow( stdout, bitrow, cols, 0 );
  282.     }
  283.     }
  284.     else    /* else use hilbert space filling curve dithering */
  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.         gray **grays;
  297.         bit **bits;
  298.     int xx,yy;
  299.     int end;
  300.     int *x,*y;
  301.     int sum = 0;
  302.  
  303.     grays = pgm_readpgm( ifp, &cols,&rows, &maxval );
  304.     bits = pbm_allocarray(cols,rows);
  305.  
  306.     x = (int *) malloc( sizeof(int) * clump_size);
  307.     y = (int *) malloc( sizeof(int) * clump_size);
  308.     if (!x || !y)
  309.         pm_error( "Run out of memory" );
  310.     init_hilbert(cols,rows);
  311.  
  312.     end = clump_size;
  313.     while (end == clump_size)
  314.     {
  315.         int i;
  316.         /* compute the next clust co-ordinates along hilbert path */
  317.         for (i = 0; i < end; i++)
  318.         {
  319.         if (hilbert(&x[i],&y[i])==0)
  320.             end = i;    /* we reached the end */
  321.         }
  322.         /* sum levels */
  323.         for (i = 0; i < end; i++)
  324.         sum += grays[y[i]][x[i]];
  325.         /* dither half and half along path */
  326.         for (i = 0; i < end; i++)
  327.         {
  328.         if (sum >= maxval)
  329.         {
  330.             bits[y[i]][x[i]] = PBM_WHITE;
  331.             sum -= maxval;
  332.         }
  333.         else
  334.             bits[y[i]][x[i]] = PBM_BLACK;
  335.         }
  336.     }
  337.     pbm_writepbm( stdout, bits, cols, rows, 0 );
  338.     }
  339.  
  340.     pm_close( ifp );
  341.  
  342.     exit( 0 );
  343. }
  344.  
  345.  
  346. /* Hilbert curve tracer */
  347.  
  348. #define MAXORD 18
  349.  
  350. static int hil_order,hil_ord;
  351. static int hil_turn;
  352. static int hil_dx,hil_dy;
  353. static int hil_x,hil_y;
  354. static int hil_stage[MAXORD];
  355. static int hil_width,hil_height;
  356.  
  357. /* Initialise the Hilbert curve tracer */
  358. static void init_hilbert(w,h)
  359.     int w,h;
  360. {
  361.     int big,ber;
  362.     hil_width = w;
  363.     hil_height = h;
  364.     big = w > h ? w : h;
  365.     for(ber = 2, hil_order = 1; ber < big; ber <<= 1, hil_order++);
  366.     if (hil_order > MAXORD)
  367.     pm_error( "Sorry, hilbert order is too large" );
  368.     hil_ord = hil_order;
  369.     hil_order--;
  370. }
  371.  
  372. /* Return non-zero if got another point */
  373. static int hilbert(px,py)
  374.     int *px,*py;
  375. {
  376.     int temp;
  377.     if (hil_ord > hil_order)    /* have to do first point */
  378.     {
  379.     hil_ord--;
  380.     hil_stage[hil_ord] = 0;
  381.     hil_turn = -1;
  382.     hil_dy = 1;
  383.     hil_dx = hil_x = hil_y = 0;
  384.     *px = *py = 0;
  385.     return 1;
  386.     }
  387.     for(;;) /* Operate the state machine */
  388.     {
  389.     switch (hil_stage[hil_ord])
  390.     {
  391.     case 0:
  392.         hil_turn = -hil_turn;
  393.         temp = hil_dy;
  394.         hil_dy = -hil_turn * hil_dx;
  395.         hil_dx = hil_turn * temp;
  396.         if (hil_ord > 0) 
  397.         {
  398.         hil_stage[hil_ord] = 1;
  399.         hil_ord--;
  400.         hil_stage[hil_ord]=0;
  401.         continue;
  402.         }
  403.     case 1:
  404.         hil_x += hil_dx;
  405.         hil_y += hil_dy;
  406.         if (hil_x < hil_width && hil_y < hil_height)
  407.         {
  408.         hil_stage[hil_ord] = 2;
  409.         *px = hil_x;
  410.         *py = hil_y;
  411.         return 1;
  412.         }
  413.     case 2:
  414.         hil_turn = -hil_turn;
  415.         temp = hil_dy;
  416.         hil_dy = -hil_turn * hil_dx;
  417.         hil_dx = hil_turn * temp;
  418.         if (hil_ord > 0) /* recurse */
  419.         {
  420.         hil_stage[hil_ord] = 3;
  421.         hil_ord--;
  422.         hil_stage[hil_ord]=0;
  423.         continue;
  424.         }
  425.     case 3:
  426.         hil_x += hil_dx;
  427.         hil_y += hil_dy;
  428.         if (hil_x < hil_width && hil_y < hil_height)
  429.         {
  430.         hil_stage[hil_ord] = 4;
  431.         *px = hil_x;
  432.         *py = hil_y;
  433.         return 1;
  434.         }
  435.     case 4:
  436.         if (hil_ord > 0) /* recurse */
  437.         {
  438.         hil_stage[hil_ord] = 5;
  439.         hil_ord--;
  440.         hil_stage[hil_ord]=0;
  441.         continue;
  442.         }
  443.     case 5:
  444.         temp = hil_dy;
  445.         hil_dy = -hil_turn * hil_dx;
  446.         hil_dx = hil_turn * temp;
  447.         hil_turn = -hil_turn;
  448.         hil_x += hil_dx;
  449.         hil_y += hil_dy;
  450.         if (hil_x < hil_width && hil_y < hil_height)
  451.         {
  452.         hil_stage[hil_ord] = 6;
  453.         *px = hil_x;
  454.         *py = hil_y;
  455.         return 1;
  456.         }
  457.     case 6:
  458.         if (hil_ord > 0) /* recurse */
  459.         {
  460.         hil_stage[hil_ord] = 7;
  461.         hil_ord--;
  462.         hil_stage[hil_ord]=0;
  463.         continue;
  464.         }
  465.     case 7:
  466.         temp = hil_dy;
  467.         hil_dy = -hil_turn * hil_dx;
  468.         hil_dx = hil_turn * temp;
  469.         hil_turn = -hil_turn;
  470.         /* Return from a recursion */
  471.         if (hil_ord < hil_order)
  472.         hil_ord++;
  473.         else
  474.         return 0;
  475.     }
  476.     }
  477. }
  478.