[Commit] tess bentley.5c,NONE,1.1

Carl Worth commit at keithp.com
Sat Jul 3 00:41:24 PDT 2004


Committed by: cworth

Update of /local/src/CVS/tess
In directory home.keithp.com:/tmp/cvs-serv22121

Added Files:
	bentley.5c 
Log Message:
First working iterative Bentley-Ottman implementation.

--- NEW FILE: bentley.5c ---
autoload Sortlist
autoload Sort
autoload SVG

typedef struct {
    int x;
    int y;
} Point;

typedef struct {
    rational x;
    rational y;
} RationalPoint;

bool point_greater (Point a, Point b)
{
    if (a.y != b.y)
	return a.y > b.y;
    else
	return a.x > b.x;
}

typedef Edge;

typedef union {
    *Edge elt;
    void  nil;
} EdgePtr;

typedef struct {
    EdgePtr    prev, next;
    Point top;
    Point middle;
    Point bottom;
} Edge;

typedef enum {
    Start, Intersect, Stop
} EventType;

typedef Event;

typedef union {
    *Event	elt;
    void	nil;
} EventPtr;

typedef struct {
    EventPtr    prev, next;
    EventType type;
    EdgePtr e1;
    EdgePtr e2;
    Point point;
} Event;

bool slope_greater (Edge a, Edge b)
{
    int adx = a.bottom.x - a.top.x;
    int bdx = b.bottom.x - b.top.x;
    
    int ady = a.bottom.y - a.top.y;
    int bdy = b.bottom.y - b.top.y;
    
    return adx * bdy > ady * bdx;
}

bool event_greater (Event a, Event b)
{
    if (point_greater (a.point, b.point))
	return true;
    if (point_greater (b.point, a.point))
	return false;
    if (a.type != b.type)
	return (a.type == EventType.Start);

    /* Slope comparison */
    return slope_greater (*a.e1.elt, *b.e1.elt);
}

exception parallel (void);

RationalPoint intersect (Edge a, Edge b) {

    int det (int a, int b, int c, int d) {
	return a * d - b * c;
    }

    int dx1 = a.top.x - a.bottom.x;
    int dy1 = a.top.y - a.bottom.y;
	
    int dx2 = b.top.x - b.bottom.x;
    int dy2 = b.top.y - b.bottom.y;

    int den_det = det (dx1, dy1, dx2, dy2);

    if (den_det == 0)
	raise parallel ();

    int a_det = det (a.top.x, a.top.y,
		     a.bottom.x, a.bottom.y);
    int b_det = det (b.top.x, b.top.y,
		     b.bottom.x, b.bottom.y);
	    
    return (RationalPoint) {
	x = det (a_det, dx1,
		 b_det, dx2) / den_det,
	y = det (a_det, dy1,
		 b_det, dy2) / den_det
    };
}


bool event_ptr_greater (EventPtr a, EventPtr b)
{
    return event_greater (*a.elt, *b.elt);
}

Edge[*]
bentley_ottman (Edge[*] segments, &int intersection_count) {

    intersection_count = 0;
    EventPtr event_queue = Sortlist::List.nil;
    Edge[...] result = {};

    EventPtr initialize_event_queue (Edge[*] segments) {
	Event[...] events = {};

	for (int i=0; i < dim(segments); i++) {

	    if (point_greater (segments[i].top, segments[i].bottom)) {
		Point tmp;
		tmp = segments[i].top;
		segments[i].top = segments[i].bottom;
		segments[i].bottom = tmp;
	    }
	    segments[i].middle = segments[i].top;

	    events[i*2] = (Event) {
		type = EventType.Start,
		e1 = (EdgePtr.elt) &segments[i],
		point = segments[i].top,
		prev = Sortlist::List.nil,
		next = Sortlist::List.nil
	    };
	    events[i*2+1] = (Event) {
		type = EventType.Stop,
		e1 = (EdgePtr.elt) &segments[i],
		point = segments[i].bottom,
		prev = Sortlist::List.nil,
		next = Sortlist::List.nil
	    };
	}
	
	Sort::qsort (&events, event_greater);

	for (int i=0; i < dim(events) - 1; i++) {
	    events[i].next = (EventPtr.elt) &events[i+1];
	    events[i+1].prev = (EventPtr.elt) &events[i];
	}

	return (EventPtr.elt) &events[0];
    }

    event_queue = initialize_event_queue (segments);

    int current_y = event_queue.elt->point.y;

    bool edge_greater (EdgePtr a, EdgePtr b) {

	rational edge_x_for_y (EdgePtr e, int y) {
	    int dx = e.elt->bottom.x - e.elt->top.x;
	    int dy = e.elt->bottom.y - e.elt->top.y;

	    if (dy == 0)
		return e.elt->top.x;

	    return e.elt->top.x + (y - e.elt->top.y) * dx / dy;
	}

	rational ax = edge_x_for_y (a, current_y);
	rational bx = edge_x_for_y (b, current_y);

	if (ax != bx)
	    return ax > bx;

	/* Slope comparison */
	return slope_greater (*a.elt, *b.elt);
    }
/*
    bool edge_greater (EdgePtr a, EdgePtr b) {
	int adx = a.elt->bottom.x - a.elt->top.x;
	int bdx = b.elt->bottom.x - b.elt->top.x;

	int ady = a.elt->bottom.y - a.elt->top.y;
	int bdy = b.elt->bottom.y - b.elt->top.y;

	int ax = bdy*(adx*(current_y - a.elt->top.y) + a.elt->top.x * ady);
	int bx = ady*(bdx*(current_y - b.elt->top.y) + b.elt->top.x * bdy);

	return ax > bx;
    }
*/

    EdgePtr sweep_line = Sortlist::List.nil;

    void insert_event_if_intersect (EdgePtr a, EdgePtr b) {

	if (a == Sortlist::List.nil || b == Sortlist::List.nil)
	    return;

	RationalPoint intersection;
	try {
	    intersection = intersect (*a.elt, *b.elt);
	} catch parallel () {
	    return;
	}

	if (intersection.y > current_y
	    && intersection.y < a.elt->bottom.y 
	    && intersection.y < b.elt->bottom.y
	    )
	{
	    Sortlist::insert (&event_queue,
			      (EventPtr.elt) &(Event) {
		                  type = EventType.Intersect,
				      e1 = a,
				      e2 = b,
				      point = {
					  x = floor (intersection.x + 0.5),
					  y = floor (intersection.y + 0.5)
				      }
			      },
			      event_ptr_greater);
	}
    }

    static int frame_count = 0;

    while (event_queue != Sortlist::List.nil) {

	file svg_file;
	int text_y;

	void fprint_edge (file f, Edge e) {
	    fprintf (f, "(%2d, %2d)-(%2d, %2d)",
		    e.top.x, e.top.y,
		    e.bottom.x, e.bottom.y);
	}

	void print_edge_svg (Edge e) {
	    fprintf (svg_file, "<line x1='%d' y1='%d' x2='%d' y2='%d' />\n",
		    e.top.x, e.top.y,
		    e.bottom.x, e.bottom.y);
	}

	void fprint_event_queue (file f) {
	    EventPtr e = event_queue;
	    while (e != Sortlist::List.nil) {
		union switch (e.elt->type) {
		case Start:
		    fprintf (f, "Start: ");
		    break;
		case Stop:
		    fprintf (f, "Stop: ");
		    break;
		case Intersect:
		    fprintf (f, "Intersect: ");
		    break;
		}
		fprintf (f, "%g:\t", e.elt->point);
		fprint_edge (f, *e.elt->e1.elt);
		if (e.elt->type == EventType.Intersect) {
		    fprintf (f, " X ");
		    fprint_edge (f, *e.elt->e2.elt);
		}
		fprintf (f, "\n");
		e = e.elt->next;
	    }
	}

	void print_event_queue_svg () {
	    EventPtr e = event_queue;
	    while (e != Sortlist::List.nil) {
		fprintf (svg_file, "<text stroke='none' fill='black' x='0' y='%d'>\n", text_y);
		text_y += 14;
		union switch (e.elt->type) {
		case Start:
		    fprintf (svg_file, "Start: ");
		    break;
		case Stop:
		    fprintf (svg_file, "Stop: ");
		    break;
		case Intersect:
		    fprintf (svg_file, "Intersect: ");
		    break;
		}
		fprintf (svg_file, "%g:\t", e.elt->point);
		fprint_edge (svg_file, *e.elt->e1.elt);
		if (e.elt->type == EventType.Intersect) {
		    fprintf (svg_file, " X ");
		    fprint_edge (svg_file, *e.elt->e2.elt);
		}
		fprintf (svg_file, "\n");
		e = e.elt->next;
		fprintf (svg_file, "</text>\n");
	    }
	}

	void fprint_sweep_line (file f) {
	    EdgePtr e = sweep_line;
	    while (e != Sortlist::List.nil) {
		fprint_edge (f, *e.elt);
		e = e.elt->next;
		if (e != Sortlist::List.nil)
		    fprintf (f, ", ");
	    }
	    fprintf (f, "\n");
	}

	EventPtr event = event_queue;
	current_y = event.elt->point.y;

	void print_sweep_line_svg () {
	    fprintf (svg_file, "<?xml version='1.0' standalone='no'?>\n");
	    fprintf (svg_file, "<!DOCTYPE svg PUBLIC '-//W3C//DTD SVG 20001102//EN'\n");
	    fprintf (svg_file, "  'http://www.w3.org/TR/2000/CR-SVG-20001102/DTD/svg-20001102.dtd'>\n");
	    fprintf (svg_file, "<svg transform='translate(70,70)' fill='none' stroke='black' font='mono' font-size='12'>\n");

	    fprintf (svg_file, "<g transform='scale(30,30)' stroke-width='.05'>\n");
	    for (int i=0; i < dim(segments); i++) {
		print_edge_svg (segments[i]);
	    }
	    fprintf (svg_file, "</g>\n");

	    fprintf (svg_file, "<text stroke='none' fill='blue' x='0' y='%d'>\n", text_y);
	    text_y += 14;
	    EdgePtr e = sweep_line;
	    while (e != Sortlist::List.nil) {
		fprint_edge (svg_file, *e.elt);
		if (e != Sortlist::List.nil)
		    fprintf (svg_file, ", ");
		e = e.elt->next;
	    }
	    fprintf (svg_file, "</text>\n");

	    fprintf (svg_file, "<g transform='scale(30,30)' stroke='blue' stroke-width='.05'>\n");
	    EdgePtr e = sweep_line;
	    int edge_cnt = 0;
	    while (e != Sortlist::List.nil) {
		print_edge_svg (*e.elt);
		fprintf (svg_file, "<text font-size='.1' fill='blue' stroke='none' x='%d' y='%d'>%d</text>\n",
			 e.elt->top.x + 20, e.elt->top.y + 10, edge_cnt++);
		e = e.elt->next;
	    }
	    fprintf (svg_file, "<line x1='-100' y1='%d' x2='1000' y2='%d' stroke-width='.1' stroke='red' />\n",
		     current_y, current_y);
	    fprintf (svg_file, "<text font-size='1' fill='red' stroke='none' x='530' y='%d'>%d</text>\n", current_y - 2, current_y);

	    fprintf (svg_file, "<circle cx='%d' cy='%d' r='.15' stroke='none' fill='green'/>\n", event.elt->point.x, event.elt->point.y);
	    fprintf (svg_file, "</g>\n");

	    print_event_queue_svg ();

	    fprintf (svg_file, "</svg>\n");
	}

	string filename;
	filename = sprintf ("bentley_%03d.svg", frame_count++);
	svg_file = File::open (filename, "w");
	text_y = 200;
	print_sweep_line_svg ();
	File::close (svg_file);

	event_queue = event_queue.elt->next;
	if (event_queue != Sortlist::List.nil)
	    event_queue.elt->prev = Sortlist::List.nil;

	union switch (event.elt->type) {
	case Start:
	    EdgePtr edge, left, right;

	    edge = event.elt->e1;
	    Sortlist::insert (&sweep_line, edge, edge_greater);
	    left = edge.elt->prev;
	    right = edge.elt->next;
	    insert_event_if_intersect (left, edge);
	    insert_event_if_intersect (edge, right);
	    break;
	case Stop:
	    EdgePtr edge, left, right;

	    result[dim(result)] = (Edge) {
		top = event.elt->e1.elt->middle,
		bottom = event.elt->point
	    };

	    edge = event.elt->e1;
	    left = edge.elt->prev;
	    right = edge.elt->next;
	    Sortlist::delete (&sweep_line, edge);
	    insert_event_if_intersect (left, right);
	    break;
	case Intersect:
	    EdgePtr e1, e2, left, right;

	    intersection_count++;

	    result[dim(result)] = (Edge) {
		top = event.elt->e1.elt->middle,
		bottom = event.elt->point
	    };
	    result[dim(result)] = (Edge) {
		top = event.elt->e2.elt->middle,
		bottom = event.elt->point
	    };
	    event.elt->e1.elt->middle = event.elt->point;
	    event.elt->e2.elt->middle = event.elt->point;

	    Sortlist::swap (&sweep_line, event.elt->e1);

	    /* after the swap e2 is left of e1 */
	    e2 = event.elt->e2;
	    left = e2.elt->prev;
	    insert_event_if_intersect (left, e2);

	    e1 = event.elt->e1;
	    right = e1.elt->next;
	    insert_event_if_intersect (e1, right);
	    break;
	}
    }

    printf ("Found %d intersections\n", intersection_count);

    return result;
}

Edge[*] bentley_ottman_iterative (Edge[*] edges)
{
    Edge[*] old_edges;

    int count;
    do {
	old_edges = edges;
	edges = bentley_ottman(old_edges, &count);
    } while (count > 0);

    return edges;
}

/* Cool, degnerate case from Hobby's paper. Require's 3 passes until
 * we implement tolerance squares. */
Edge[*] hobby_edges = { { top = { x =   0, y =   0}, bottom = { x =   9, y =   5}},
			{ top = { x =   0, y =   0}, bottom = { x =  13, y =   6}},
			{ top = { x =   9, y =   5}, bottom = { x =  -1, y =  -2}}
};

bentley_ottman_iterative (hobby_edges);




More information about the Commit mailing list