/* pgmtopbm.c - read a portable graymap and write a portable bitmap
**
** Copyright (C) 1989 by 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 <assert.h>
#include "pgm.h"
#include "dithers.h"
#include "mallocvar.h"

static void init_hilbert ARGS((int w, int h));
static int hilbert ARGS((int *px, int *py));


int
main( argc, argv )
    int argc;
    char* argv[];
{
    FILE* ifp;
    gray* grayrow;
    register gray* gP;
    bit* bitrow;
    register bit* bP;
    int argn, rows, cols, format, row, col, limitcol;
    float fthreshval;
    int clump_size;
    gray maxval;
    const char* const usage = "[-floyd|-fs  | -hilbert | -threshold | -dither8|-d8 |\n     -cluster3|-c3|-cluster4|-c4|-cluster8|-c8] [-value <val>] [-clump <size>] [pgmfile]";
    int halftone;
#define QT_FS 1
#define QT_THRESH 2
#define QT_DITHER8 3
#define QT_CLUSTER3 4
#define QT_CLUSTER4 5
#define QT_CLUSTER8 6
#define QT_HILBERT 7
    long threshval, sum;
    long* thiserr;
    long* nexterr;
    long* temperr;
#define FS_SCALE 1024
#define HALF_FS_SCALE 512
    int fs_direction;

    pgm_init( &argc, argv );

    argn = 1;
    halftone = QT_FS;   /* default quantization is Floyd-Steinberg */
    fthreshval = 0.5;
    clump_size = 5; /* default hilbert curve clump size */

    while ( argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0' )
    {
        if ( pm_keymatch( argv[argn], "-fs", 2 ) ||
             pm_keymatch( argv[argn], "-floyd", 2 ) )
            halftone = QT_FS;
        else if ( pm_keymatch( argv[argn], "-hilbert", 2 ))
            halftone = QT_HILBERT;
        else if ( pm_keymatch( argv[argn], "-threshold", 2 ) )
            halftone = QT_THRESH;
        else if ( pm_keymatch( argv[argn], "-dither8", 2 ) ||
                  pm_keymatch( argv[argn], "-d8", 3 ) )
            halftone = QT_DITHER8;
        else if ( pm_keymatch( argv[argn], "-cluster3", 9 ) ||
                  pm_keymatch( argv[argn], "-c3", 3 ) )
            halftone = QT_CLUSTER3;
        else if ( pm_keymatch( argv[argn], "-cluster4", 9 ) ||
                  pm_keymatch( argv[argn], "-c4", 3 ) )
            halftone = QT_CLUSTER4;
        else if ( pm_keymatch( argv[argn], "-cluster8", 9 ) ||
                  pm_keymatch( argv[argn], "-c8", 3 ) )
            halftone = QT_CLUSTER8;
        else if ( pm_keymatch( argv[argn], "-value", 2 ) )
        {
            /* Before 2001-06-08, this strtod() was a sscanf("%f").  Arjen
               Bax found on his CYGWIN system this causes a comma to be
               used as a decimal separator.  But on 2001-09.03, Charles
               Wilson found this not to be the case.  
            */
            char *ep;
            ++argn;
            if ( argn == argc
                 || (fthreshval = strtod(argv[argn], &ep), *ep != 0)
                 || fthreshval < 0.0 || fthreshval > 1.0 )
            {
                pm_usage( usage );
            }
        }
        else if ( pm_keymatch( argv[argn], "-clump", 2 ) )
        {
            ++argn;
            if ( argn == argc || 
                 sscanf( argv[argn], "%d", &clump_size ) != 1 ||
                 clump_size < 2)
                pm_usage( usage );
        }
        else
            pm_usage( usage );
        ++argn;
    }

    if ( argn != argc )
    {
        ifp = pm_openr( argv[argn] );
        ++argn;
    }
    else
        ifp = stdin;

    if ( argn != argc )
        pm_usage( usage );

    if ( halftone != QT_HILBERT )
    {
        pgm_readpgminit( ifp, &cols, &rows, &maxval, &format );
        grayrow = pgm_allocrow( cols );

        pbm_writepbminit( stdout, cols, rows, 0 );
        bitrow = pbm_allocrow( cols );

        /* Initialize. */
        switch ( halftone )
        {
        case QT_FS:
            /* Initialize Floyd-Steinberg error vectors. */
            thiserr = (long*) pm_allocrow( cols + 2, sizeof(long) );
            nexterr = (long*) pm_allocrow( cols + 2, sizeof(long) );
            srand( (int) ( time( 0 ) ^ getpid( ) ) );
            for ( col = 0; col < cols + 2; ++col )
                thiserr[col] = ( rand( ) % FS_SCALE - HALF_FS_SCALE ) / 4;
            /* (random errors in [-FS_SCALE/8 .. FS_SCALE/8]) */
            fs_direction = 1;
            threshval = fthreshval * FS_SCALE;
            break;

        case QT_HILBERT:
            break;

        case QT_THRESH:
            threshval = fthreshval * maxval + 0.999999;
            break;

        case QT_DITHER8:
            /* Scale dither matrix. */
            for ( row = 0; row < 16; ++row )
                for ( col = 0; col < 16; ++col )
                    dither8[row][col] = 
                        dither8[row][col] * maxval / 256;
            break;

        case QT_CLUSTER3:
            /* Scale order-3 clustered dither matrix. */
            for ( row = 0; row < 6; ++row )
                for ( col = 0; col < 6; ++col )
                    cluster3[row][col] = 
                        cluster3[row][col] * maxval / 18;
            break;

        case QT_CLUSTER4:
            /* Scale order-4 clustered dither matrix. */
            for ( row = 0; row < 8; ++row )
                for ( col = 0; col < 8; ++col )
                    cluster4[row][col] = 
                        cluster4[row][col] * maxval / 32;
            break;

        case QT_CLUSTER8:
            /* Scale order-8 clustered dither matrix. */
            for ( row = 0; row < 16; ++row )
                for ( col = 0; col < 16; ++col )
                    cluster8[row][col] = 
                        cluster8[row][col] * maxval / 128;
            break;

        default:
            pm_error( "can't happen" );
            exit( 1 );
        }

        for ( row = 0; row < rows; ++row )
        {
            pgm_readpgmrow( ifp, grayrow, cols, maxval, format );
    
            switch ( halftone )
            {
            case QT_FS:
                for ( col = 0; col < cols + 2; ++col )
                    nexterr[col] = 0;
                if ( fs_direction )
                {
                    col = 0;
                    limitcol = cols;
                    gP = grayrow;
                    bP = bitrow;
                }
                else
                {
                    col = cols - 1;
                    limitcol = -1;
                    gP = &(grayrow[col]);
                    bP = &(bitrow[col]);
                }
                do
                {
                    sum = 
                        ( (long) *gP * FS_SCALE ) / maxval + thiserr[col + 1];
                    if ( sum >= threshval )
                    {
                        *bP = PBM_WHITE;
                        sum = sum - threshval - HALF_FS_SCALE;
                    }
                    else
                        *bP = PBM_BLACK;
    
                    if ( fs_direction )
                    {
                        thiserr[col + 2] += ( sum * 7 ) / 16;
                        nexterr[col    ] += ( sum * 3 ) / 16;
                        nexterr[col + 1] += ( sum * 5 ) / 16;
                        nexterr[col + 2] += ( sum     ) / 16;
    
                        ++col;
                        ++gP;
                        ++bP;
                    }
                    else
                    {
                        thiserr[col    ] += ( sum * 7 ) / 16;
                        nexterr[col + 2] += ( sum * 3 ) / 16;
                        nexterr[col + 1] += ( sum * 5 ) / 16;
                        nexterr[col    ] += ( sum     ) / 16;
    
                        --col;
                        --gP;
                        --bP;
                    }
                }
                while ( col != limitcol );
                temperr = thiserr;
                thiserr = nexterr;
                nexterr = temperr;
                fs_direction = ! fs_direction;
                break;
    
            case QT_THRESH: {
                unsigned int col;
                for ( col = 0; col < cols; ++col )
                    if ( grayrow[col] >= threshval )
                        bitrow[col] = PBM_WHITE;
                    else
                        bitrow[col] = PBM_BLACK;
            }
            break;
    
            case QT_DITHER8: {
                unsigned int col;
                for ( col = 0; col < cols; ++col )
                    if ( grayrow[col] > dither8[row % 16][col % 16] )
                        bitrow[col] = PBM_WHITE;
                    else
                        bitrow[col] = PBM_BLACK;
            }
            break;
    
            case QT_CLUSTER3: {
                unsigned int col;
                for ( col = 0; col < cols; ++col )
                    if ( grayrow[col] > cluster3[row % 6][col % 6] )
                        bitrow[col] = PBM_WHITE;
                    else
                        bitrow[col] = PBM_BLACK;
            }
            break;
    
            case QT_CLUSTER4: {
                unsigned int col;
                for ( col = 0; col < cols; ++col )
                    if ( grayrow[col] > cluster4[row % 8][col % 8] )
                        bitrow[col] = PBM_WHITE;
                    else
                        bitrow[col] = PBM_BLACK;
            }
            break;
                
            case QT_CLUSTER8: {
                unsigned int col;
                for ( col = 0; col < cols; ++col )
                    if ( grayrow[col] > cluster8[row % 16][col % 16] )
                        bitrow[col] = PBM_WHITE;
                    else
                        bitrow[col] = PBM_BLACK;
            }
            break;
    
            default:
                pm_error( "Internal error: impossible value of 'halftone'" );
            }
            pbm_writepbmrow( stdout, bitrow, cols, 0 );
        }
    }
    else    /* else use hilbert space filling curve dithering */
        /*
         * This is taken from the article "Digital Halftoning with
         * Space Filling Curves" by Luiz Velho, proceedings of
         * SIGRAPH '91, page 81.
         *
         * This is not a terribly efficient or quick version of
         * this algorithm, but it seems to work. - Graeme Gill.
         * graeme@labtam.labtam.OZ.AU
         *
         */
    {
        gray **grays;
        bit **bits;
        int end;
        int *x,*y;
        int sum = 0;

        grays = pgm_readpgm( ifp, &cols,&rows, &maxval );
        bits = pbm_allocarray(cols,rows);

        MALLOCARRAY(x, clump_size);
        MALLOCARRAY(y, clump_size);
        if (x == NULL  || y == NULL)
            pm_error( "Run out of memory" );
        init_hilbert(cols,rows);

        end = clump_size;
        while (end == clump_size)
        {
            int i;
            /* compute the next clust co-ordinates along hilbert path */
            for (i = 0; i < end; i++)
            {
                if (hilbert(&x[i],&y[i])==0)
                    end = i;    /* we reached the end */
            }
            /* sum levels */
            for (i = 0; i < end; i++)
                sum += grays[y[i]][x[i]];
            /* dither half and half along path */
            for (i = 0; i < end; i++)
            {
                if (sum >= maxval)
                {
                    bits[y[i]][x[i]] = PBM_WHITE;
                    sum -= maxval;
                }
                else
                    bits[y[i]][x[i]] = PBM_BLACK;
            }
        }
        pbm_writepbm( stdout, bits, cols, rows, 0 );
    }

    pm_close( ifp );

    exit( 0 );
}


/* Hilbert curve tracer */

#define MAXORD 18

static int hil_order,hil_ord;
static int hil_turn;
static int hil_dx,hil_dy;
static int hil_x,hil_y;
static int hil_stage[MAXORD];
static int hil_width,hil_height;

/* Initialise the Hilbert curve tracer */
static void init_hilbert(w,h)
    int w,h;
{
    int big,ber;
    hil_width = w;
    hil_height = h;
    big = w > h ? w : h;
    for(ber = 2, hil_order = 1; ber < big; ber <<= 1, hil_order++);
    if (hil_order > MAXORD)
        pm_error( "Sorry, hilbert order is too large" );
    hil_ord = hil_order;
    hil_order--;
}

/* Return non-zero if got another point */
static int hilbert(px,py)
    int *px,*py;
{
    int temp;
    if (hil_ord > hil_order)    /* have to do first point */
    {
        hil_ord--;
        hil_stage[hil_ord] = 0;
        hil_turn = -1;
        hil_dy = 1;
        hil_dx = hil_x = hil_y = 0;
        *px = *py = 0;
        return 1;
    }
    for(;;) /* Operate the state machine */
    {
        switch (hil_stage[hil_ord])
        {
        case 0:
            hil_turn = -hil_turn;
            temp = hil_dy;
            hil_dy = -hil_turn * hil_dx;
            hil_dx = hil_turn * temp;
            if (hil_ord > 0) 
            {
                hil_stage[hil_ord] = 1;
                hil_ord--;
                hil_stage[hil_ord]=0;
                continue;
            }
        case 1:
            hil_x += hil_dx;
            hil_y += hil_dy;
            if (hil_x < hil_width && hil_y < hil_height)
            {
                hil_stage[hil_ord] = 2;
                *px = hil_x;
                *py = hil_y;
                return 1;
            }
        case 2:
            hil_turn = -hil_turn;
            temp = hil_dy;
            hil_dy = -hil_turn * hil_dx;
            hil_dx = hil_turn * temp;
            if (hil_ord > 0) /* recurse */
            {
                hil_stage[hil_ord] = 3;
                hil_ord--;
                hil_stage[hil_ord]=0;
                continue;
            }
        case 3:
            hil_x += hil_dx;
            hil_y += hil_dy;
            if (hil_x < hil_width && hil_y < hil_height)
            {
                hil_stage[hil_ord] = 4;
                *px = hil_x;
                *py = hil_y;
                return 1;
            }
        case 4:
            if (hil_ord > 0) /* recurse */
            {
                hil_stage[hil_ord] = 5;
                hil_ord--;
                hil_stage[hil_ord]=0;
                continue;
            }
        case 5:
            temp = hil_dy;
            hil_dy = -hil_turn * hil_dx;
            hil_dx = hil_turn * temp;
            hil_turn = -hil_turn;
            hil_x += hil_dx;
            hil_y += hil_dy;
            if (hil_x < hil_width && hil_y < hil_height)
            {
                hil_stage[hil_ord] = 6;
                *px = hil_x;
                *py = hil_y;
                return 1;
            }
        case 6:
            if (hil_ord > 0) /* recurse */
            {
                hil_stage[hil_ord] = 7;
                hil_ord--;
                hil_stage[hil_ord]=0;
                continue;
            }
        case 7:
            temp = hil_dy;
            hil_dy = -hil_turn * hil_dx;
            hil_dx = hil_turn * temp;
            hil_turn = -hil_turn;
            /* Return from a recursion */
            if (hil_ord < hil_order)
                hil_ord++;
            else
                return 0;
        }
    }
}