#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#ifdef _WINDOWS
#define drand48()	((double) rand()/(double) RAND_MAX)
#endif

int Xsize, Ysize, Maxval ;

unsigned char *Image = NULL ;
unsigned char *Normal = NULL ;
unsigned char *Edge = NULL ;
float *Brightness = NULL ;
float *EdgeX ;
float *EdgeY ;
float Gamma = 4.0 ;
float Length = 6.0 ;
float Dither = 0.0 ;

#define SQR(x)  ((x)*(x))
#define MIN(a,b) ((a<b)?(a):(b))
#define MAX(a,b) ((a>b)?(a):(b))

/*
 * return the area defined by the intersection of the wide line from 
 * ex0, ey0 to ex1, ey1, of width "width", with the pixel going from 
 * x,y to x+1, y+1
 */

float
area(float x, float y, float ex0, float ey0, float ex1, float ey1, float width)
{
    float dx, dy ;
    float s0, s1, t0, t1 ;
    float tmp, l ;

    ex0 -= x ; ex1 -= x ;
    ey0 -= y ; ey1 -= y ;

    dx = ex1 - ex0 ;
    dy = ey1 - ey0 ;

    l = sqrt(SQR(dx)+SQR(dy)) ;

    if (fabs(dx) > 0.0001) {
	s0 = -ex0 / dx ;
	s1 = (1.0-ex0) / dx ;
	if (s0 > s1) {
	    tmp = s0 ; s0 = s1 ; s1 = tmp ;
	}
    } else {
	if (ex0 < 0.0 || ex0 > 1.0)
	    return 0.0 ;
	s0 = -1e6 ;
	s1 = -1e6 ;
    }

    if (fabs(dy) > 0.0001) {
	t0 = -ey0 / dy ;
	t1 = (1.0-ey0) / dy ;
	if (t0 > t1) {
	    tmp = t0 ; t0 = t1 ; t1 = tmp ;
	}
    } else {
	if (ey0 < 0.0 || ey0 > 1.0)
	    return 0.0 ;
	t0 = -1e6 ;
	t1 = 1e6 ;
    }

    t0 = MAX(s0, t0) ;
    t1 = MIN(s1, t1) ;

    if (t0 < 0.0) t0 = 0.0 ;
    if (t1 > 1.0) t1 = 1.0 ;
	
    if (t0 > t1)
	return 0.0 ;
    else
	return width * l * (t1 - t0) ;
}

/* find the two dimensional distance to the line segment */

float
linedist(float x, float y, float ex0, float ey0, float ex1, float ey1)
{
    float dx, dy ;
    float d0, d1, dm ;
    float xn, yn, t ;

    dx = ex1 - ex0 ; dy = ey1 - ey0 ;

    d0 = sqrt(SQR(x-ex0)+SQR(y-ey0)) ;
    d1 = sqrt(SQR(x-ex1)+SQR(y-ey1)) ;

    t = (dx * (x-ex0) + dy * (y-ey0)) /
            (SQR(dx) + SQR(dy)) ;

    xn = ex0 + t * dx ;
    yn = ey0 + t * dy ;

    dm = sqrt(SQR(x-xn)+SQR(y-yn)) ;

    return MIN(MIN(d0, d1), dm) ;
}

int
INDEX(int x, int y)
{
    if (x < 0) x = 0 ;
    if (x >= Xsize) x = Xsize-1 ;
    if (y < 0) y = 0 ;
    if (y >= Ysize) y = Ysize-1 ;
    return (y*Xsize+x) ;
}

struct t_bound {
    int minx, miny ;
    int maxx, maxy ;
} ;

void 
bound(struct t_bound *b, float x0, float y0, float x1, float y1, float width)
{
    b->minx = floor(x0-width) ;
    b->miny = floor(y0-width) ;
    b->maxx = ceil(x0+width) ;
    b->maxy = ceil(y0+width) ;

    if (floor(x1-width) < b->minx)
	b->minx = floor(x1-width) ;
    if (floor(y1-width) < b->miny)
	b->miny = floor(y1-width) ;
    if (ceil(x1+width) > b->maxx)
	b->maxx = floor(x1+width) ;
    if (ceil(y1+width) > b->maxy)
	b->maxy = floor(y1+width) ;

    if (b->minx < 0) b->minx = 0 ;
    if (b->maxx >= Xsize) b->maxx = Xsize-1 ;
    if (b->miny < 0) b->miny = 0 ;
    if (b->maxy >= Ysize) b->maxy = Ysize-1 ;
}

void
ink(char *ifname, char *efname, char *nfname, float width)
{
    FILE *fp ;
    int nxsize, nysize ;
    int i, j, k, x, y, failures ;
    float fx, fy ;
    float dx, dy, a ;
    unsigned char *ip ;
    float *bp, t ;
    float x0, y0, x1, y1 ;
    int nx, ny ;
    unsigned char val, newval ;
    struct t_bound b ;

    /* first, suck in the image file */

    if ((fp=fopen(ifname, "r")) == NULL) {
	perror(ifname) ;
	return ;
    }

    if (fscanf(fp, "P6\n%d %d\n%d\n", &Xsize, &Ysize, &Maxval) != 3) {
	fprintf(stderr, "%s: invalid format\n", ifname) ;
	fclose(fp) ;
	return ;
    }

    Image = (unsigned char *) malloc(Xsize*Ysize*3) ;
    fread(Image, sizeof(unsigned char)*Xsize*Ysize*3, 1, fp) ;
    fclose(fp) ;

    /* if we have a normal file, do the same for that */

    if (nfname) {
	if ((fp=fopen(nfname, "r")) == NULL) {
	    perror(ifname) ;
	    return ;
	}   

	if (fscanf(fp, "P6\n%d %d\n%d\n", &nxsize, &nysize, &Maxval) != 3) {
	    fprintf(stderr, "%s: invalid format\n", nfname) ;
	    fclose(fp) ;
	    return ;
	}   

	Normal = (unsigned char *) malloc(Xsize*Ysize*3) ;
	fread(Normal, Xsize*Ysize*3*sizeof(unsigned char), 1, fp) ;
	fclose(fp) ;
    }

    /* if we have a edge file, do the same for that */

    if (efname) {
	if ((fp=fopen(efname, "r")) == NULL) {
	    perror(efname) ;
	    return ;
	}   

	if (fscanf(fp, "P5\n%d %d\n%d\n", &nxsize, &nysize, &Maxval) != 3) {
	    fprintf(stderr, "%s: invalid format\n", efname) ;
	    fclose(fp) ;
	    return ;
	}   

	Edge = (unsigned char *) malloc(Xsize*Ysize) ;
	fread(Edge, Xsize*Ysize*sizeof(unsigned char), 1, fp) ;
	fclose(fp) ;
    }

    /* convert the real image to a luminance */

    Brightness = (float *) calloc(Xsize*Ysize, sizeof(float)) ;
    bp = Brightness ; ip = Image ;

    for (i=0; i<Xsize*Ysize; i++) {
	*bp++ = (0.30 * ip[0] + 0.59 * ip[1] + 0.11 * ip[2])/255.0 ;
	ip += 3 ;
    }

    /* now, invert the luminance, and apply some default gamma */

#define GAMMA_VALUE 4.0
    for (i=0, bp=Brightness; i<Xsize*Ysize; i++, bp++)
	*bp = pow(1.0-*bp, Gamma) ;

    /* we are ready to start, go ahead and print the postscript header
     * that we need to invert the image 
     */

    fp = fopen("www.pgm", "w") ;
    fprintf(fp, "P5\n%d %d\n%d\n", Xsize, Ysize, 255) ;
    for (i=0, bp=Brightness; i<Xsize*Ysize; i++, bp++)
	fputc(255.0 * *bp, fp) ;
    fclose(fp) ;

    printf("%%!\n") ;
    printf("%% Postscript Pen and Ink Drawing generated by ppmtopen\n") ;
    printf("%% written by Mark VandeWettering <markv@pixar.com>\n") ;
    printf("0 %d translate\n", Ysize) ;
    printf("1 -1 scale\n") ;
    printf("%f setlinewidth\n", width) ;
    printf("newpath\n") ;

#define NSTROKES        (1000000)

    failures = 0 ;
    for (i=0; i<NSTROKES && failures < NSTROKES ; i++) {

        fx = (drand48() * (Xsize-1)) ;
        fy = (drand48() * (Ysize-1)) ;

        x = (int) fx ;
        y = (int) fy ;

        if (drand48() > Brightness[INDEX(x,y)]) {
            failures ++ ;
            continue ;
        }

	if (Normal) {
	    dy = ((float) Normal[3*INDEX(x,y)] - 128.0)/128.0 ;
	    dx = ((float) Normal[3*INDEX(x,y)+1] - 128.0)/128.0 ;
	} else {
	    dx = drand48()*2.0-1.0 ;
	    dy = drand48()*2.0-1.0 ;
	}

	dx *= Length/2.0 ;
	dy *= Length/2.0 ;

#define DELTA_T 	(0.01)
	if (Edge) {
	    val = Edge[INDEX(x,y)] ;
	    for (t=DELTA_T; t<1.0; t+=DELTA_T) {
		nx = (int) (fx - t * dx) ;
		ny = (int) (fy - t * dy) ;
		newval = Edge[INDEX(nx,ny)] ;
		if (newval < val) {
		   x0 = fx - t * dx ;
		   y0 = fy - t * dy ;
		   goto skip_it_1;
		}
		val = newval ;
	    }
	    x0 = fx - dx ;
	    y0 = fy - dy ;
skip_it_1:
	    val = Edge[INDEX(x,y)] ;
	    for (t=DELTA_T; t<1.0; t+=DELTA_T) {
		nx = (int) (fx + t * dx) ;
		ny = (int) (fy + t * dy) ;
		newval = Edge[INDEX(nx,ny)] ;
		if (newval < val) {
		   x1 = fx + t * dx ;
		   y1 = fy + t * dy ;
		   goto skip_it_2;
		}
		val = newval ;
	    }
	    x1 = fx + dx ;
	    y1 = fy + dy ;
skip_it_2:
	    val = Edge[INDEX(x,y)] ;
	} else {
	    x0 = fx - dx + Dither * (drand48()-0.5) ;
	    x1 = fx + dx + Dither * (drand48()-0.5) ;
	    y0 = fy - dy + Dither * (drand48()-0.5) ;
	    y1 = fy + dy + Dither * (drand48()-0.5) ;
	}

        printf("%f %f moveto\n", x0, y0) ;
        printf("%f %f lineto\n", x1, y1) ;

	bound(&b, x0, y0, x1, y1, width) ;

        for (j=b.miny; j<b.maxy; j++) {
            for (k=b.minx; k<b.maxx ; k++) {
		a = area((float) k, (float) j,
			 fx-dx, fy-dy, fx+dx, fy+dy, width) ;
		Brightness[INDEX(k,j)] -= a ;
		if (Brightness[INDEX(k,j)]<0.0)
			Brightness[INDEX(k,j)] = 0.0 ;
            }
        }
     }

     fp = fopen("xxx.pgm", "w") ;
     fprintf(fp, "P5\n%d %d\n%d\n", Xsize, Ysize, 255) ;
     for (i=0, bp=Brightness; i<Xsize*Ysize; i++, bp++)
	fputc(255.0 * *bp, fp) ;
     fclose(fp) ;
     printf("stroke\n") ;
     printf("showpage\n") ;
}

int
main(int argc, char *argv[])
{
    extern int optind ;
    extern char *optarg ;
    int c ;
    char *normname = NULL ;
    char *edgename = NULL ;
    float width = 0.1 ;

    while ((c = getopt(argc, argv, "e:d:g:n:w:l:")) != EOF) {
	switch (c) {
	case 'l':
	    Length = atof(optarg) ;
	    break ;
	case 'e':
	    edgename = optarg ;
	    break ;
        case 'd':
	    Dither = atof(optarg) ;
	    break ;
	case 'n':
	    normname = optarg ;
	    break ;
	case 'g':
	    Gamma = atof(optarg) ;
	    break ;
	case 'w':
	    width = atof(optarg) ;
	    break ;
	default:
	    break ;
	}
    }

    for (; optind<argc; optind++)
	ink(argv[optind], edgename, normname, width) ;

    printf("quit\n") ;
    return(0);
}