/* unmscheme.c -- Lance R. Williams, April 15, 2006. Copyright 2006 University of New Mexico. All rights reserved. Permission to copy and modify this software and its documen- tation only for internal use in your organization is hereby granted, provided that this notice is retained thereon and on all copies. UNM makes no representations as to the sui- tability and operability of this software for any purpose. It is provided "as is" without express or implied warranty. UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT- NESS. IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT- SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER- FORMANCE OF THIS SOFTWARE. No other rights, including, for example, the right to redis- tribute this software and its documentation or the right to prepare derivative works, are granted unless specifically provided in a separate license agreement. Copyright 2006, University of New Mexico. All rights reserved. */ /* gcc unmscheme.c -o unmscheme -lm -O3 -Wall */ #define GL 1 /* gcc unmscheme.c -o unmscheme -lm -lglut -lGL -lGLU -L/x11/R6/lib -O3 -Wall */ #include #include #include #include #include #include #include #include #include #include #ifdef GL #include #include #endif #include "complex.c" #include "svd.c" #include "four1.c" typedef struct sexpr { int type; union { int i; double x; char c; char *text; fcomplex z; struct symbol *symbol; struct pair *pair; struct closure *closure; struct primitive *primitive; struct sexpr *sexpr; struct plumber *plumber; struct image *image; struct complex_image *complex_image; struct line *line; struct sketch *sketch; struct vector *vector; struct continuation *continuation; FILE *file; } u; } sexpr; typedef struct node { sexpr *sp; struct node *next; } node; typedef struct pair { sexpr *car; sexpr *cdr; } pair; typedef struct closure { int n; sexpr *lookup; sexpr *compiled; sexpr *env; } closure; typedef struct symbol { char *name; closure *macro; } symbol; typedef struct primitive { int arity; char *name; sexpr *(*func)(); } primitive; typedef struct vector { int length; sexpr **vector; } vector; typedef struct continuation { int s; sexpr **args; sexpr **cv_stack; sexpr ***pc_stack; sexpr **lu_stack; sexpr **e_stack; int *r_stack; } continuation; typedef struct plumber { double x; double y; double heading; } plumber; #include "unmscheme.h" typedef struct image { int rows; int cols; float *data; } image; typedef struct complex_image { int rows; int cols; fcomplex *data; } complex_image; typedef struct line { float x[2]; float y[2]; float theta; float contrast; float length; float coverage; sexpr *links[2]; int replaced; int merged; int age; } line; typedef struct sketch { int rows; int cols; int bixel_rows; int bixel_cols; float bixel_dx; float bixel_dy; sexpr **grid; int count; sexpr *ls; } sketch; #define PI 3.141592653589 #define TWOPI 6.283185307179 #define HALFPI 1.570796326794 #define S30 0.5 #define C30 0.86602540378 #define C45 0.70710678118 #define max(a,b) ((a) > (b) ? (a) : (b)); #define min(a,b) ((a) < (b) ? (a) : (b)); #define mod(x,n) ((x) >= (n) ? (x)-(n) : ((x) < 0 ? (n)+(x) : (x))) #define opposite(x) ((x == 1) ? 0 : 1) #define even(x) ((x) % 2 ? 0 : 1) #define odd(x) ((x) % 2 ? 1 : 0) #define print2stdout(x) display_or_print((x),PRINT,TOPLEVEL,stdout); #define STRLEN 256 #define NIL 0 #define NUMBER 1 #define TOKEN 2 #define CHARACTER 3 #define STRING 4 #define VIRGIN 5 #define CLOSURE 6 #define PRIMITIVE 7 #define BOOLEAN 8 #define VECTOR 9 #define INPUT_PORT 10 #define OUTPUT_PORT 11 #define EOF_OBJECT 12 #define GRAPHIC 13 #define PLUMBER 14 #define IMAGE 15 #define COLOR_IMAGE 16 #define COMPLEX_IMAGE 17 #define LINE 18 #define SKETCH 19 #define CONTINUATION 20 #define COMPLEX 21 #define UNDEFINED 22 #define SYMBOL 100 #define PAIR 200 #define COPIED 300 #define HALT 401 #define LOCAL_REFER 402 #define GLOBAL_REFER 403 #define REFER 404 #define CONSTANT 405 #define CLOSE 406 #define RETURN 407 #define TEST 408 #define LOCAL_ASSIGN 409 #define GLOBAL_ASSIGN 410 #define ASSIGN 411 #define DEFINE 412 #define APPLY 413 #define FUNCALL 414 #define FRAME 415 #define ARGUMENT 416 #define CONTI 417 #define NUATE 418 #define THROW 419 #define CATCH 420 #define null(sp) ((sp)->type == NIL) #define pair(sp) ((sp)->type == PAIR) #define number(sp) ((sp)->type == NUMBER) #define token(sp) ((sp)->type == TOKEN) #define character(sp) ((sp)->type == CHARACTER) #define string(sp) ((sp)->type == STRING) #define symbol(sp) ((sp)->type == SYMBOL) #define virgin(sp) ((sp)->type == VIRGIN) #define closure(sp) ((sp)->type == CLOSURE) #define primitive(sp) ((sp)->type == PRIMITIVE) #define boolean(sp) ((sp)->type == BOOLEAN) #define vector(sp) ((sp)->type == VECTOR) #define graphic(sp) ((sp)->type == GRAPHIC) #define plumber(sp) ((sp)->type == PLUMBER) #define image(sp) ((sp)->type == IMAGE) #define color_image(sp) ((sp)->type == COLOR_IMAGE) #define complex_image(sp) ((sp)->type == COMPLEX_IMAGE) #define line(sp) ((sp)->type == LINE) #define sketch(sp) ((sp)->type == SKETCH) #define input_port(sp) ((sp)->type == INPUT_PORT) #define output_port(sp) ((sp)->type == OUTPUT_PORT) #define eof_object(sp) ((sp)->type == EOF_OBJECT) #define komplex(sp) ((sp)->type == COMPLEX) #define undefined(sp) ((sp)->type == UNDEFINED) #define copied(sp) ((sp)->type == COPIED) #define bytecode(sp) (((sp)->type >= HALT) && ((sp)->type <= NUATE)) #define self_evaluating(sp) ((sp)->type < PAIR) #define HASHSIZE 10001 #define SEXPRS 80000000 #define STACKSIZE 10000 static sexpr *global_env_vars; static sexpr *global_env_vals; static sexpr *global_vars; static sexpr *global_vals; static sexpr *macros; static sexpr *nil; static sexpr *eof_object; static sexpr *undefined; static sexpr *true; static sexpr *false; static sexpr *current_input_port; static sexpr *current_output_port; static sexpr *symbol$nil; static sexpr *symbol$straight; static sexpr *symbol$spot; static sexpr *symbol$transparent; static sexpr *symbol$bend; static sexpr *symbol$adjoin; static sexpr *symbol$adorn; static sexpr *symbol$text; static sexpr *graphic$default_color; static sexpr *graphic$nil; static sexpr *bytecode$halt; static sexpr *bytecode$local_refer; static sexpr *bytecode$global_refer; static sexpr *bytecode$refer; static sexpr *bytecode$constant; static sexpr *bytecode$close; static sexpr *bytecode$return; static sexpr *bytecode$test; static sexpr *bytecode$local_assign; static sexpr *bytecode$global_assign; static sexpr *bytecode$assign; static sexpr *bytecode$define; static sexpr *bytecode$funcall; static sexpr *bytecode$apply; static sexpr *bytecode$frame; static sexpr *bytecode$argument; static sexpr *bytecode$conti; static sexpr *bytecode$nuate; static sexpr *bytecode$throw; static sexpr *bytecode$catch; static sexpr *symbol$lambda; static sexpr *symbol$begin; static sexpr *symbol$apply; static sexpr *symbol$define; static sexpr *symbol$define_macro; static sexpr *symbol$if; static sexpr *symbol$callcc; static sexpr *symbol$callec; static sexpr *primitive$append; static sexpr *primitive$cons; static sexpr *symbol$let; static sexpr *symbol$let_star; static sexpr *symbol$letrec; static sexpr *symbol$set_bang; static sexpr *symbol$quote; static sexpr *symbol$quasiquote; static sexpr *symbol$unquote; static sexpr *symbol$unquote_splicing; static sexpr *symbol$hyphen; static sexpr *token$left; static sexpr *token$right; static sexpr *token$dot; static sexpr *token$quote; static sexpr *token$quasiquote; static sexpr *token$unquote; static sexpr *token$unquote_splicing; static sexpr *token$pound_sign; static node *hashtab[HASHSIZE]; static int gargc; static char **gargv; static int gensyms; static int gc; static sexpr *smalloc; static sexpr *smalloc0[SEXPRS]; static sexpr *smalloc1[SEXPRS]; static int loads; static jmp_buf esc; #ifdef GL static sexpr *graphic_sp; static sexpr *image_sp; static sexpr *sketch_sp; #endif static double image_display_scale = 1.0; int main(int argc, char **argv) { int i; sexpr *input, *result; gargc = argc; gargv = argv; char c; gc = 0; smalloc = (sexpr *) smalloc0; for (i = 0; i < HASHSIZE; i++) hashtab[i] = NULL; gensyms = 0; define_scheme_constants(); make_global_env(); load(str2exp("compiler-basics.scm")); load(str2exp("plumbing.scm")); load(str2exp("boldt.scm")); printf("Welcome to UNM Scheme 2.2 Copyright (c) 2006 The University of New Mexico\n"); setjmp(esc); while (loads < argc) load(str2exp(argv[loads++])); printf("> "); while ((c = getc(stdin)) != EOF) { ungetc(c,stdin); input = parse(stdin,'\n'); if (!null(input)) { result = eval(car(input)); while (!null(input = cdr(input))) { print(result); printf("> "); result = eval(car(input)); } print(result); garbage_collect(); printf("> "); } } exit(0); } void print_value_escape(char *msg, struct sexpr *sp) { printf("%s",msg); if (undefined(sp)) printf("\n"); else print(sp); longjmp(esc,1); } int isfalse(sexpr *sp) { return (boolean(sp) && (sp->u.i == 0)); } sexpr *color(double r, double g, double b) { return cons(num2exp(r),cons(num2exp(g),cons(num2exp(b),nil))); } sexpr *make_bytecode(char *name, int i) { sexpr *sp = malloc(sizeof(sexpr)); sp->type = i; sp->u.symbol = malloc(sizeof(symbol)); sp->u.symbol->name = malloc((strlen(name) + 1)*sizeof(char)); strcpy(sp->u.symbol->name,name); return sp; } void define_scheme_constants() { /* Scheme constants */ loads = 0; nil = malloc(sizeof(sexpr)); nil->type = NIL; eof_object = malloc(sizeof(sexpr)); eof_object->type = EOF_OBJECT; current_input_port = malloc(sizeof(sexpr)); current_input_port->type = INPUT_PORT; current_input_port->u.file = stdin; current_output_port = malloc(sizeof(sexpr)); current_output_port->type = OUTPUT_PORT; current_output_port->u.file = stdout; undefined = malloc(sizeof(sexpr)); undefined->type = UNDEFINED; true = malloc(sizeof(sexpr)); true->type = BOOLEAN; true->u.i = 1; false = malloc(sizeof(sexpr)); false->type = BOOLEAN; false->u.i = 0; global_vars = make_vector(1000); global_vars->u.vector->length = 0; global_vals = make_vector(1000); global_vals->u.vector->length = 0; global_env_vars = cons(global_vars,nil); global_env_vals = cons(global_vals,nil); macros=cons(make_vector(0),nil); /* Plumbing Graphics constants */ symbol$nil = str2symbol("nil"); symbol$straight = str2symbol("straight"); symbol$spot = str2symbol("spot"); symbol$transparent = str2symbol("transparent"); symbol$bend = str2symbol("bend"); symbol$adjoin = str2symbol("adjoin"); symbol$adorn = str2symbol("adorn"); symbol$text = str2symbol("text"); graphic$default_color = color(255.0,255.0,255.0); graphic$nil = malloc(sizeof(sexpr)); graphic$nil->type = GRAPHIC; graphic$nil->u.sexpr = nil; /* Constants needed for compilation */ symbol$quote = str2symbol("quote"); symbol$quasiquote = str2symbol("quasiquote"); symbol$unquote = str2symbol("unquote"); symbol$unquote_splicing = str2symbol("unquote-splicing"); primitive$cons = make_primitive(cons,2,"cons"); primitive$append = make_primitive(append,2,"scheme:append"); symbol$lambda = str2symbol("lambda"); symbol$begin = str2symbol("scheme:begin"); symbol$apply = str2symbol("apply"); symbol$define = str2symbol("scheme:define"); symbol$define_macro = str2symbol("define-macro"); symbol$if = str2symbol("if"); symbol$callcc = str2symbol("call/cc"); symbol$callec = str2symbol("call/ec"); symbol$let = str2symbol("let"); symbol$let_star = str2symbol("let*"); symbol$letrec = str2symbol("letrec"); symbol$set_bang = str2symbol("set!"); /* Scanner constants */ token$left = char2token('('); token$right = char2token(')'); token$dot = char2token('.'); token$quote = char2token('\''); token$quasiquote = char2token('`'); token$unquote = char2token(','); token$unquote_splicing = char2token('@'); token$pound_sign = char2token('#'); symbol$hyphen = str2symbol("-"); /* Virtual machine constants */ bytecode$halt = make_bytecode("halt",HALT); bytecode$local_refer = make_bytecode("local-refer",LOCAL_REFER); bytecode$global_refer = make_bytecode("global-refer",GLOBAL_REFER); bytecode$refer = make_bytecode("refer",REFER); bytecode$constant = make_bytecode("constant",CONSTANT); bytecode$close = make_bytecode("close",CLOSE); bytecode$return = make_bytecode("return",RETURN); bytecode$test = make_bytecode("test",TEST); bytecode$local_assign = make_bytecode("local-assign",LOCAL_ASSIGN); bytecode$global_assign = make_bytecode("global-assign",GLOBAL_ASSIGN); bytecode$assign = make_bytecode("assign",ASSIGN); bytecode$define = make_bytecode("define",DEFINE); bytecode$apply = make_bytecode("apply",APPLY); bytecode$funcall = make_bytecode("funcall",FUNCALL); bytecode$frame = make_bytecode("frame",FRAME); bytecode$argument = make_bytecode("argument",ARGUMENT); bytecode$conti = make_bytecode("conti",CONTI); bytecode$nuate = make_bytecode("nuate",NUATE); bytecode$throw = make_bytecode("throw",THROW); bytecode$catch = make_bytecode("catch",CATCH); } sexpr *load(sexpr *filename) { FILE *file; sexpr *input, *result; if (!string(filename)) { printf("load: Illegal filename.\n"); longjmp(esc,1); } file = fopen(filename->u.text,"r"); if (file) { input = parse(file,EOF); if (!pair(input)) longjmp(esc,1); result = nil; while (!null(input)) { result = eval(car(input)); input = cdr(input); } fclose(file); return result; } else { printf("load: File not found.\n"); longjmp(esc,1); } } void eatline(FILE *file, char c) { while (c != EOF && c != '\n') c=getc(file); return; } sexpr *open_input_file(sexpr *filename) { sexpr *sp = smalloc++; sp->type = INPUT_PORT; if (!string(filename)) { printf("open-input-file: Illegal filename.\n"); longjmp(esc,1); } sp->u.file = fopen(filename->u.text,"r"); if (sp->u.file) return sp; else { printf("open-input-file: Attempt to open input file failed.\n"); longjmp(esc,1); } } sexpr *open_output_file(sexpr *filename) { sexpr *sp = smalloc++; sp->type = OUTPUT_PORT; if (!string(filename)) { printf("open-output-file: Illegal filename.\n"); longjmp(esc,1); } sp->u.file = fopen(filename->u.text,"w"); if (sp->u.file) return sp; printf("open-output-file: Attempt to open output file failed.\n"); longjmp(esc,1); } sexpr *close_input_port(sexpr *port) { if (!input_port(port)) { print_value_escape("close-input-port: Argument is not an input-port: ",port); } fclose(port->u.file); return undefined; } sexpr *close_output_port(sexpr *port) { if (!output_port(port)) { print_value_escape("close-output-port: Argument is not an output-port: ",port); } fclose(port->u.file); return undefined; } sexpr *read_char(sexpr *port) { char c; if (!input_port(port)) { print_value_escape("read-char: Argument is not an input-port: ",port); } else { c = getc(port->u.file); if (c != EOF) return char2exp(c); else return eof_object; } } sexpr *write_char(sexpr *c, sexpr *port) { if (!output_port(port)) { print_value_escape("write-char: Argument is not an output-port: ",port); } putc(c->u.c,port->u.file); return undefined; } sexpr *peek_char(sexpr *port) { char c; if (!input_port(port)) { print_value_escape("peek-char: Argument is not an input-port: ",port); } c = getc(port->u.file); ungetc(c,port->u.file); if (c != EOF) return char2exp(c); else return eof_object; } int isspecialer(char c) { return(c=='~' || c=='!' || c=='#' || c=='$' || c=='%' || c=='^' || c=='&' || c=='*' || c=='-' || c=='_' || c=='+' || c=='=' || c==':' || c=='<' || c=='>' || c=='/' || c=='?' || c == '@'); } int isletter(char c) { return((c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z')); } int islegal1st(char c) { return(isletter(c) || isspecialer(c)); } int islegal2nd(char c) { return(isletter(c) || isspecialer(c) || isdigit(c)); } int isdigit_or_pt(char c) { return(isdigit(c) || c=='.'); } sexpr *lex_error(FILE *file, char c) { printf("Lexical error reading character '%c'\n",c); eatline(file, c); longjmp(esc,1); return nil; } int eatit(int (*func)(char), char *name, FILE *file, char stop) { char c; int m = 0; c = getc(file); while (func(c) && (c != stop) && (c != EOF)) { if (m >= STRLEN) lex_error(file, c); name[m++] = c; c = getc(file); } name[m] = '\0'; ungetc(c,file); return 1; } double eat_number(FILE *file, char stop) { char c; int m = 0; double x; int state = 1; char *name = malloc(STRLEN*sizeof(char)); c = getc(file); while ((c != stop) && (c != EOF)) { switch(state) { case 1: if (isdigit(c)) state = 2; else if (c == '.') state = 3; else { printf("Malformed number.\n"); longjmp(esc,1); } break; case 2: if (isdigit(c)) state = 2; else if (c == '.') state = 3; else { ungetc(c, file); name[m] = '\0'; sscanf(name,"%lg",&x); free(name); return x; } break; case 3: if (isdigit(c)) state = 3; else { ungetc(c, file); name[m] = '\0'; sscanf(name,"%lg",&x); free(name); return x; } break; } name[m++] = c; c = getc(file); } if (state == 2 || state == 3) { name[m] = '\0'; ungetc(c,file); sscanf(name,"%lg",&x); free(name); return x; } printf("Malformed number.\n"); longjmp(esc,1); } sexpr *cons(sexpr *car, sexpr *cdr) { sexpr *sp; sp = smalloc++; sp->type = PAIR; sp->u.pair=(pair *) malloc(sizeof(pair)); sp->u.pair->car=car; sp->u.pair->cdr=cdr; return sp; } sexpr *car(sexpr *sp) { if (pair(sp)) return sp->u.pair->car; print_value_escape("Attempt to compute car of non-pair: ",sp); } sexpr *cdr(sexpr *sp) { if (pair(sp)) return sp->u.pair->cdr; print_value_escape("Attempt to compute cdr of non-pair: ",sp); } sexpr *caar(sexpr *sp) { return car(car(sp)); } sexpr *cadr(sexpr *sp) { return car(cdr(sp)); } sexpr *cdar(sexpr *sp) { return cdr(car(sp)); } sexpr *cddr(sexpr *sp) { return cdr(cdr(sp)); } sexpr *caddr(sexpr *sp) { return car(cdr(cdr(sp))); } sexpr *cadddr(sexpr *sp) { return car(cdr(cdr(cdr(sp)))); } double remainder(double x, double y) { double ratio = x/y; return y*(ratio-floor(ratio)); } double deg2rad(double deg) { return (remainder(deg,360)/180.0)*PI; } double angle_difference(double rad1, double rad2) { return remainder(rad1 - rad2,TWOPI); } #ifdef GL void display_graphic() { sexpr *p = make_plumber(0.0,-0.66666666667,HALFPI); glClear(GL_COLOR_BUFFER_BIT); draw_plumber(p,graphic$default_color); p=draw_graphic(p,graphic_sp); draw_plumber(p,graphic$default_color); glFlush(); } void display_image() { glClear(GL_COLOR_BUFFER_BIT); glPixelZoom(image_display_scale,-image_display_scale); glRasterPos2i(-1,1); glDrawPixels(image_sp->u.image->cols,image_sp->u.image->rows,GL_LUMINANCE,GL_FLOAT,image_sp->u.image->data); glFlush(); } void display_color_image() { glClear(GL_COLOR_BUFFER_BIT); glPixelZoom(image_display_scale,-image_display_scale); glRasterPos2i(-1,1); glDrawPixels(image_sp->u.image->cols,image_sp->u.image->rows,GL_RGB,GL_FLOAT,image_sp->u.image->data); glFlush(); } void draw_line(sexpr *lp, float sx, float sy) { glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); glEnable(GL_BLEND); glEnable(GL_LINE_SMOOTH); glBegin(GL_LINES); glVertex2f(lp->u.line->x[0]/sx-1.0,1.0-lp->u.line->y[0]/sy); glVertex2f(lp->u.line->x[1]/sx-1.0,1.0-lp->u.line->y[1]/sy); glEnd(); } #define RGB 0 #define HOT 1 void draw_lines(sexpr *ls0, float sx, float sy, int color, float r, float g, float b) { sexpr *first, *ls; float contrast, avg, delta, total, stddev; int n; ls = ls0; total = 0.0; n = 0; while (!null(ls)) { total += car(ls)->u.line->contrast; n++; ls = ls->u.pair->cdr; } avg = total/n; ls = ls0; total = 0.0; while (!null(ls)) { delta = car(ls)->u.line->contrast - avg; total += delta*delta; ls = ls->u.pair->cdr; } stddev = avg + sqrt(total/n); ls = ls0; while (!null(ls)) { first = car(ls); if (color == HOT) { contrast = first->u.line->contrast; r = rhot(0.0,stddev,contrast); g = ghot(0.0,stddev,contrast); b = bhot(0.0,stddev,contrast); } glColor3f(r,g,b); draw_line(first,sx,sy); ls = ls->u.pair->cdr; } } void draw_link(float x1, float y1, float x2, float y2, float sx, float sy) { float cos_theta, sin_theta, cos_phi, sin_phi, d, xc, yc, t; float r = 0.5; float rad1, rad2; d = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)); cos_theta = d/(r*2.0); sin_theta = sin(acos(cos_theta)); cos_phi = (x2-x1)/d; sin_phi = (y2-y1)/d; xc = x1 + r*(cos_theta*cos_phi - sin_theta*sin_phi); yc = y1 + r*(sin_theta*cos_phi + cos_theta*sin_phi); rad1 = atan2(y1-yc,x1-xc); rad2 = atan2(y2-yc,x2-xc); glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); glEnable(GL_BLEND); glEnable(GL_LINE_SMOOTH); if (angle_difference((double) rad1,(double) rad2) > angle_difference((double) rad2,(double) rad1)) { glBegin(GL_LINE_STRIP); for (t = rad2; angle_difference((double) t,(double) rad1) > 0.05; t = t + 0.05) glVertex2f((xc+cos(t)*r)/sx-1.0,1.0-(yc+sin(t)*r)/sy); glEnd(); } else { glBegin(GL_LINE_STRIP); for (t = rad1; angle_difference((double) t,(double) rad2) > 0.05; t = t + 0.05) glVertex2f((xc+cos(t)*r)/sx-1.0,1.0-(yc+sin(t)*r)/sy); glEnd(); } } void draw_links(sexpr *ls, float sx, float sy) { line *ln1, *ln2; sexpr *lns; int end1, end2; float x1, y1, x2, y2; glColor3f(1.0,1.0,1.0); if (!null(ls)) { for (end1 = 0; end1 <= 1; end1++) { ln1 = car(ls)->u.line; lns = ln1->links[end1]; while (!null(lns)) { ln2 = car(lns)->u.line; end2 = opposite(end1); x1 = ln1->x[end1]; y1 = ln1->y[end1]; x2 = ln2->x[end2]; y2 = ln2->y[end2]; draw_link(x1,y1,x2,y2,sx,sy); lns = cdr(lns); } } draw_links(cdr(ls),sx,sy); } } void display_sketch() { sketch *lp = sketch_sp->u.sketch; glClear(GL_COLOR_BUFFER_BIT); glLineWidth(1.0); draw_lines(lp->ls,lp->cols/2.0,lp->rows/2.0,HOT,0.0,0.0,0.0); glFlush(); } void display_link_graph() { sketch *lp = sketch_sp->u.sketch; glClear(GL_COLOR_BUFFER_BIT); glLineWidth(1.0); draw_lines(lp->ls,lp->cols/2.0,lp->rows/2.0,RGB,0.0,0.0,1.0); draw_links(lp->ls,lp->cols/2.0,lp->rows/2.0); glFlush(); } sexpr *show_link_graph(sexpr *sp) { pid_t pid; switch (pid = fork()) { case -1: printf("fork failed"); break; case 0: glutInit(&gargc,gargv); glutInitWindowSize(sp->u.sketch->cols,sp->u.sketch->rows); glutInitWindowPosition(0,0); glutCreateWindow("link graph"); sketch_sp = sp; glutDisplayFunc(display_link_graph); glutMainLoop(); exit(0); break; } return undefined; } #endif float bhot(float xmin, float xmax, float x) { x = (x - xmin)/(xmax - xmin); if (x < 0.666666667) return 0.0; else return (x - 0.666666667)*3.0; } float ghot(float xmin, float xmax, float x) { x = (x - xmin)/(xmax - xmin); if (x < 0.333333333) return 0.0; else if (x < 0.666666667) return (x - 0.333333333)*3.0; else return 1.0; } float rhot(float xmin, float xmax, float x) { x = (x - xmin)/(xmax - xmin); if (x < 0.333333333) return x*3.0; else return 1.0; } #define DISPLAY 0 #define PRINT 1 #define EMBEDDED 0 #define TOPLEVEL 1 sexpr *display_or_print(sexpr *sp, int print_flag, int toplevel_flag, FILE *file) { #ifdef GL pid_t pid; #endif char text[STRLEN]; int r, c; float x0, y0, x1, y1; switch(sp->type) { case NIL: fprintf(file,"()"); break; case PAIR: pair_print(sp,print_flag,file); break; case NUMBER: fprintf(file,"%g",sp->u.x); break; case TOKEN: if (eq(sp,token$unquote_splicing)) fprintf(file,",@"); else fprintf(file,"%c",sp->u.c); break; case CHARACTER: if (print_flag) { switch (sp->u.c) { case '\n': fprintf(file,"#\\newline"); break; case ' ': fprintf(file,"#\\space"); break; case '\t': fprintf(file,"#\\tab"); break; default: fprintf(file,"#\\%c",sp->u.c); } } else fprintf(file,"%c",sp->u.c); break; case STRING: if (print_flag) fprintf(file,"\"%s\"",sp->u.text); else fprintf(file,"%s",sp->u.text); break; case SYMBOL: fprintf(file,"%s",sp->u.symbol->name); break; case VIRGIN: case CLOSURE: fprintf(file,"#"); break; case PRIMITIVE: fprintf(file,"#",sp->u.primitive->name); break; case BOOLEAN: fprintf(file,"%s",sp->u.i ? "#t" : "#f"); break; case VECTOR: vector_print(sp,print_flag,file); break; case INPUT_PORT: fprintf(file,"#"); break; case OUTPUT_PORT: fprintf(file,"#"); break; case EOF_OBJECT: fprintf(file,"#"); break; case GRAPHIC: sprintf(text,"#",gtype(sp)->u.symbol->name); if (toplevel_flag) { #ifdef GL switch (pid = fork()) { case -1: printf("fork failed"); break; case 0: glutInit(&gargc,gargv); glutInitWindowSize(700,700); glutCreateWindow(text); graphic_sp = sp; glutDisplayFunc(display_graphic); glutMainLoop(); exit(0); default: if (print_flag) fprintf(file,text); return undefined; } #else if (print_flag) fprintf(file,text); return undefined; #endif } break; case PLUMBER: x0 = sp->u.plumber->x; y0 = sp->u.plumber->y; x1 = sp->u.plumber->heading; fprintf(file,"#u.image->rows; c = sp->u.image->cols; sprintf(text,"#",r,c); if (toplevel_flag) { #ifdef GL switch (pid = fork()) { case -1: printf("fork failed"); break; case 0: glutInit(&gargc,gargv); glutInitWindowSize(c*image_display_scale,r*image_display_scale); glutCreateWindow(text); image_sp = image_normalize(sp); glutDisplayFunc(display_image); glutMainLoop(); exit(0); default: if (print_flag) fprintf(file,text); return undefined; } #else if (print_flag) fprintf(file,text); return undefined; #endif } break; case COMPLEX_IMAGE: r = sp->u.complex_image->rows; c = sp->u.complex_image->cols; sprintf(text,"#",r,c); if (toplevel_flag) { #ifdef GL switch (pid = fork()) { case -1: printf("fork failed"); break; case 0: glutInit(&gargc,gargv); glutInitWindowSize(c*image_display_scale,r*image_display_scale); glutCreateWindow(text); image_sp = color_image_normalize(complex2color(sp)); glutDisplayFunc(display_color_image); glutMainLoop(); exit(0); default: if (print_flag) fprintf(file,text); return undefined; } #else if (print_flag) fprintf(file,text); return undefined; #endif } else fprintf(file,text); break; case COLOR_IMAGE: r = sp->u.image->rows; c = sp->u.image->cols; sprintf(text,"#",r,c); if (toplevel_flag) { #ifdef GL switch (pid = fork()) { case -1: printf("fork failed"); break; case 0: glutInit(&gargc,gargv); glutInitWindowSize(c*image_display_scale,r*image_display_scale); glutCreateWindow(text); image_sp = color_image_normalize(sp); glutDisplayFunc(display_color_image); glutMainLoop(); exit(0); default: if (print_flag) fprintf(file,text); return undefined; } #else if (print_flag) fprintf(file,text); return undefined; #endif } break; case LINE: x0 = sp->u.line->x[0]; y0 = sp->u.line->y[0]; x1 = sp->u.line->x[1]; y1 = sp->u.line->y[1]; fprintf(file,"#",x0,y0,x1,y1); break; case SKETCH: r = sp->u.sketch->rows; c = sp->u.sketch->cols; sprintf(text,"#",r,c,sp->u.sketch->count); if (toplevel_flag) { #ifdef GL switch (pid = fork()) { case -1: printf("fork failed"); break; case 0: glutInit(&gargc,gargv); glutInitWindowSize(c,r); glutInitWindowPosition(0,0); glutCreateWindow(text); sketch_sp = sp; glutDisplayFunc(display_sketch); glutMainLoop(); exit(0); default: if (print_flag) fprintf(file,text); return undefined; } #else if (print_flag) fprintf(file,text); return undefined; #endif } break; case CONTINUATION: fprintf(file,"#"); break; case UNDEFINED: fprintf(file,"#"); break; case COMPLEX: if (sp->u.z.i < 0) fprintf(file,"%g-%gi",(double) sp->u.z.r, (double) fabs(sp->u.z.i)); else fprintf(file,"%g+%gi",(double) sp->u.z.r, (double) sp->u.z.i); break; default: fprintf(file,"%s",sp->u.symbol->name); } return undefined; } void ReadString(FILE *ff, char *str) { char ch; int i; for(i=0, ch=fgetc(ff); ch!=EOF && ch!=10 && ch!=13; i++) { str[i]=ch; ch = fgetc(ff); } str[i]=0; } int ReadCommentedString(FILE *ff, char *str) { char trash[256]; int i; if (0==fscanf(ff,"%s",str)) return(3); /* EOF */ for (i=0; str[i] != '#' && str[i] != 0; i++); if (str[i]=='#') ReadString(ff,trash); str[i]=0; return 0; } int ReadNextNumber(FILE *ff, unsigned *x) { int res; char str[256]; for (res=ReadCommentedString(ff,str); !res && sscanf(str,"%u",x)!=1; res=ReadCommentedString(ff,str)); return res; } sexpr *read_image_helper(FILE *file, int image_type) { unsigned x = 0, y = 0, greylevels = 0; unsigned i, n, data; sexpr *sp; if (ReadNextNumber(file,&x) || ReadNextNumber(file,&y) || ReadNextNumber(file,&greylevels)) { printf("read-image: Bad image header.\n"); longjmp(esc,1); } if (image_type == IMAGE) { sp = make_image((int) y, (int) x); n = y*x; } else { sp = make_color_image((int) y, (int) x); n = y*x*3; } for (i = 0; i < n; i++) { if (ReadNextNumber(file,&data)) { printf("read-image: Bad data.\n"); longjmp(esc,1); } sp->u.image->data[i] = (float) data; } return sp; } sexpr *readpgm(sexpr *filename) { FILE *file; sexpr *sp; if (!string(filename)) { printf("read-image: Illegal filename.\n"); longjmp(esc,1); } file=fopen(filename->u.text,"r"); if (!file) { printf("read-image: File not found.\n"); longjmp(esc,1); } if (getc(file) != 'P' || getc(file) != '2') { printf("read-image: Image file must be ASCII .pgm format.\n"); longjmp(esc,1); } sp = read_image_helper(file,IMAGE); fclose(file); return sp; } sexpr *readppm(sexpr *filename) { FILE *file; sexpr *sp; if (!string(filename)) { printf("read-color-image: Illegal filename.\n"); longjmp(esc,1); } file=fopen(filename->u.text,"r"); if (!file) { printf("read-color-image: File not found.\n"); longjmp(esc,1); } if (getc(file) != 'P' || getc(file) != '3') { printf("read-color-image: Image file must be ASCII .ppm format.\n"); longjmp(esc,1); } sp = read_image_helper(file,COLOR_IMAGE); fclose(file); return sp; } sexpr *writepgm(sexpr *sp, sexpr *filename) { FILE *file; int i, j, n, grey_value, max_grey_value, rows, cols; if (!image(sp)) { printf("write-image: First argument must be an image.\n"); longjmp(esc,1); } rows = sp->u.image->rows; cols = sp->u.image->cols; n = rows*cols; if (!string(filename)) { printf("write-image: Illegal filename.\n"); longjmp(esc,1); } max_grey_value = 0; for (i = 0; i < n; i++) { grey_value = sp->u.image->data[i]; if (grey_value < 0) { print_value_escape("write-image: Negative grey value: ",num2exp(grey_value)); } if (grey_value > max_grey_value) max_grey_value = grey_value; } file = fopen(filename->u.text,"w"); fprintf(file,"P2\n"); fprintf(file,"# Creator: University of New Mexico Scheme 2.0\n"); fprintf(file,"%d %d\n",cols,rows); fprintf(file,"%d\n",max_grey_value); j = 1; for (i = 0; i < n; i++) { fprintf(file,"%5d ",(int) floor(sp->u.image->data[i] + 0.5)); if (j++ == 15) { j = 1; fprintf(file,"\n"); } } fprintf(file,"\n"); fclose(file); return undefined; } sexpr *writeppm(sexpr *sp, sexpr *filename) { FILE *file; int i, j, n, grey_value, max_grey_value, rows, cols; if (!color_image(sp)) { printf("write-color-image: First argument must be a color image.\n"); longjmp(esc,1); } rows = sp->u.image->rows; cols = sp->u.image->cols; n = rows*cols; if (!string(filename)) { printf("write-color-image: Illegal filename.\n"); longjmp(esc,1); } max_grey_value = 0; for (i = 0; i < 3*n; i++) { grey_value = sp->u.image->data[i]; if (grey_value < 0) { print_value_escape("write-color-image: Negative grey value: ",num2exp(grey_value)); } if (grey_value > max_grey_value) max_grey_value = grey_value; } file=fopen(filename->u.text,"w"); fprintf(file,"P3\n"); fprintf(file,"# Creator: University of New Mexico Scheme 2.0\n"); fprintf(file,"%d %d\n",cols,rows); fprintf(file,"%d\n",max_grey_value); j = 1; for (i = 0; i < n; i++) { fprintf(file,"%5d ",(int) floor(sp->u.image->data[i*3] + 0.5)); fprintf(file,"%5d ",(int) floor(sp->u.image->data[i*3+1] + 0.5)); fprintf(file,"%5d ",(int) floor(sp->u.image->data[i*3+2] + 0.5)); if (j++ == 5) { j = 1; fprintf(file,"\n"); } } fprintf(file,"\n"); fclose(file); return undefined; } sexpr *image_ref(sexpr *sp1, sexpr *sp2, sexpr *sp3) { int i, j, rows, cols; float x, y; float f00, f01, f10, f11, *data; fcomplex c00, c01, c10, c11, *cdata; if ((!image(sp1) && !complex_image(sp1)) || !number(sp2) || !number(sp3)) { printf("image-ref: Illegal argument.\n"); longjmp(esc,1); } rows = sp1->u.image->rows; cols = sp1->u.image->cols; i = (int) floor((float) sp2->u.x); j = (int) floor((float) sp3->u.x); if (i < 0 || i >= rows || j < 0 || j >= cols) { return num2exp(0.0); } y = sp2->u.x - (float) i; x = sp3->u.x - (float) j; if (x > 0 || y > 0) { switch (sp1->type) { case IMAGE: data = sp1->u.image->data; f00 = data[i*cols + j]; f01 = data[i*cols + j + 1]; f10 = data[(i+1)*cols + j]; f11 = data[(i+1)*cols + j + 1]; return num2exp((f01-f00)*x+(f10-f00)*y+(f11+f00-f10-f01)*x*y+f00); case COMPLEX_IMAGE: cdata = sp1->u.complex_image->data; c00 = cdata[i*cols + j]; c01 = cdata[i*cols + j + 1]; c10 = cdata[(i+1)*cols + j]; c11 = cdata[(i+1)*cols + j + 1]; return complex2exp(Cadd(RCmul(x,Csub(c01,c00)),Cadd(RCmul(y,Csub(c10,c00)),Cadd(RCmul(x*y,Cadd(c11,Csub(c00,Cadd(c10,c01)))),c00)))); default: print_value_escape("image-ref: Argument is not an image: ",sp1); } } else { switch (sp1->type) { case IMAGE: data = sp1->u.image->data; return num2exp(data[i*cols + j]); case COMPLEX_IMAGE: cdata = sp1->u.complex_image->data; return complex2exp(cdata[i*cols + j]); default: print_value_escape("image-ref: Argument is not an image: ",sp1); } } } sexpr *image_set(sexpr *sp1, sexpr *sp2, sexpr *sp3, sexpr *sp4) { int i, j, index; if ((!image(sp1) && !complex_image(sp1)) || !number(sp2) || !number(sp3)) { printf("image-set!: Illegal argument.\n"); longjmp(esc,1); } i = (int) floor((float) sp2->u.x); j = (int) floor((float) sp3->u.x); if (i < 0 || i >= sp1->u.image->rows) { printf("image-set!: Row index out of range: %d\n",i); longjmp(esc,1); } if (j < 0 || j >= sp1->u.image->cols) { printf("image-set!: Column index out of range: %d\n",j); longjmp(esc,1); } switch (sp1->type) { case IMAGE: if (number(sp4)) { index = i*sp1->u.image->cols + j; sp1->u.image->data[index] = (float) sp4->u.x; return undefined; } print_value_escape("image-set!: Argument is not a real number: ",sp4); case COMPLEX_IMAGE: switch (sp4->type) { case NUMBER: index = i*sp1->u.image->cols + j; sp1->u.complex_image->data[index] = Complex(sp4->u.x,0.0); return undefined; case COMPLEX: index = i*sp1->u.image->cols + j; sp1->u.complex_image->data[index] = sp4->u.z; return undefined; default: print_value_escape("image-set!: Argument is not a number: ",sp4); } default: print_value_escape("image-set!: Argument is not an image: ",sp1); } } int member(sexpr *item, sexpr *ls) { while (!null(ls)) { if (item == ls->u.pair->car) return 1; ls=ls->u.pair->cdr; } return 0; } void check_image_sizes(char *name, int row, int col, image *ip) { if (row != ip->rows && row != -1) { printf("%s: Row size %d not equal to %d.\n",name,row,ip->rows); longjmp(esc,1); } if (col != ip->cols && col != -1) { printf("%s: Column size %d not equal to %d.\n",name,col,ip->cols); longjmp(esc,1); } return; } sexpr *image2array(sexpr *sp1) { int i, j, rows, cols; int index; sexpr *sp2; vector *row, *col; if (!image(sp1) && !complex_image(sp1)) { printf("image->array: Argument is not an image.\n"); longjmp(esc,1); } rows = sp1->u.image->rows; cols = sp1->u.image->cols; sp2 = smalloc++; sp2->type = VECTOR; col = sp2->u.vector = malloc(sizeof(vector)); col->length = rows; col->vector = malloc(rows*sizeof(sexpr)); for (i = 0; i < rows; i++) { col->vector[i] = smalloc++; col->vector[i]->type = VECTOR; row = col->vector[i]->u.vector = malloc(sizeof(vector)); row->length = cols; row->vector = malloc(cols*sizeof(sexpr)); if (image(sp1)) { for (j = 0; j < cols; j++) { index = cols*i + j; row->vector[j] = num2exp(sp1->u.image->data[index]); } } else { for (j = 0; j < cols; j++) { index = cols*i + j; row->vector[j] = complex2exp(sp1->u.complex_image->data[index]); } } } return sp2; } sexpr *make_image(int rows, int cols) { sexpr *sp = smalloc++; sp->type = IMAGE; sp->u.image = (image *) malloc(sizeof(image)); sp->u.image->rows = rows; sp->u.image->cols = cols; sp->u.image->data = malloc(cols*rows*sizeof(float)); return sp; } sexpr *make_color_image(int rows, int cols) { sexpr *sp = smalloc++; sp->type = COLOR_IMAGE; sp->u.image = (image *) malloc(sizeof(image)); sp->u.image->rows = rows; sp->u.image->cols = cols; sp->u.image->data = malloc(3*cols*rows*sizeof(float)); return sp; } sexpr *make_complex_image(int rows, int cols) { sexpr *sp = smalloc++; sp->type = COMPLEX_IMAGE; sp->u.complex_image = (complex_image *) malloc(sizeof(complex_image)); sp->u.complex_image->rows = rows; sp->u.complex_image->cols = cols; sp->u.complex_image->data = malloc(cols*rows*sizeof(fcomplex)); return sp; } sexpr *make_covariance_matrix(sexpr *sp1) { int i, j, k, rows, cols, n, m; float x, y; float **images; float *data3; sexpr *sp2, *sp3; if (!vector(sp1)) { printf("make-covariance-matrix: Illegal argument.\n"); longjmp(esc,1); } /* find size */ m = sp1->u.vector->length; sp2 = sp1->u.vector->vector[0]; if (!image(sp2)) { printf("make-covariance-matrix: Illegal argument.\n"); longjmp(esc,1); } rows = sp2->u.image->rows; cols = sp2->u.image->cols; images = (float **) malloc(m*sizeof(float *)); images[0] = sp2->u.image->data; /* check sizes */ for (j = 1; j < m; j++) { sp2 = sp1->u.vector->vector[j]; if (!image(sp2)) { printf("make-covariance-matrix: Illegal argument.\n"); free(images); longjmp(esc,1); } check_image_sizes("make-covariance-matrix",rows,cols,sp2->u.image); images[j] = sp2->u.image->data; } /* make image */ sp3 = make_image(m,m); data3 = sp3->u.image->data; /* for every pixel */ n = rows*cols; for (i = 0; i < n; i++) { /* for every image */ for (j = 0; j < m; j++) { x = images[j][i]; for (k = 0; k < m; k++) { y = images[k][i]; data3[j*m + k] += x*y; } } } /* normalize */ for (j = 0; j < m; j++) { for (k = 0; k < m; k++) { data3[j*m + k] /= n; } } /* return matrix */ free(images); return sp3; } sexpr *test_svd(sexpr *sp) { int i, j, rows; sexpr *singular_vectors, *singular_values; double **A, *s; if (!image(sp)) { printf("svd: Argument is not an image.\n"); longjmp(esc,1); } rows = sp->u.image->rows; if (rows != sp->u.image->cols) { printf("svd: Argument is not a square image.\n"); longjmp(esc,1); } A = (double **) malloc(2*rows*sizeof(double *)); s = (double *) malloc(rows*sizeof(double)); for (i = 0; i < 2*rows; i++) A[i] = malloc(rows*sizeof(double)); for (i = 0; i < rows; i++) { for (j = 0; j < rows; j++) A[i][j] = (double) sp->u.image->data[i*rows+j]; } svd(A,s,rows); singular_vectors = make_image(rows,rows); singular_values = make_vector(rows); for (i = 0; i < rows; i++) { singular_values->u.vector->vector[i] = num2exp(s[i]); for (j = 0; j < rows; j++) { singular_vectors->u.image->data[i*rows+j] = (float) A[rows+i][j]; } } for (i = 0; i < 2*rows; i++) free(A[i]); free(A); free(s); return list2(singular_vectors,singular_values); } float proj_pt_on_line_rel_dist(sexpr *ln, int end, float x, float y) { float x0 = ln->u.line->x[end]; float y0 = ln->u.line->y[end]; float x1 = ln->u.line->x[opposite(end)]; float y1 = ln->u.line->y[opposite(end)]; float length = ln->u.line->length; float dx = x1-x0; float dy = y1-y0; return ((x-x0)*dx + (y-y0)*dy)/(length*length); } float pt_line_dist_sq(float x0, float y0, float x1, float y1, float x, float y) { float dx = x1-x0; float dy = y1-y0; float length_squared = dx*dx + dy*dy; float dot = ((x-x0)*dx + (y-y0)*dy)/length_squared; float u = x0 + dx*dot; float v = y0 + dy*dot; return (u-x)*(u-x) + (v-y)*(v-y); } int line_intersects_circle(float x0, float y0, float x1, float y1, float x, float y, float r) { float dx=x1-x0; float dy=y1-y0; float m, b, d, f, g, h; if (fabs(dx) >= fabs(dy) || dy == 0) { m = dy/dx; b = -m*(x0-x)+(y0-y); f = -2*m*b; g = 2*(m*m+1); h = f*f+2*g*r*r; if (h < 0) return 0; d = sqrt(h); return ((x1 >= (f+d)/g+x >= x0) && (x1 >= (f-d)/g+x >= x0)); } else { m = dx/dy; b = -m*(y0-y)+(x0-x); f = -2*m*b; g = 2*(m*m+1); h = f*f+2*g*r*r; if (h < 0) return 0; d = sqrt(h); return ((y1 >= (f+d)/g+y >= y0) && (y1 >= (f-d)/g+y >= y0)); } } sexpr *get_all_lns_end_in_circle(sketch *sketch, float x, float y, int end, float r) { int index; float x0, y0; float rr = r*r; float bixel_dx = sketch->bixel_dx; float bixel_dy = sketch->bixel_dy; float bixel_rows = (float) sketch->bixel_rows; float bixel_cols = (float) sketch->bixel_cols; int bixel_i0 = (int) floor((y-r)/bixel_dy); int bixel_j0 = (int) floor((x-r)/bixel_dx); int bixel_i1 = (int) ceil((y+r)/bixel_dy); int bixel_j1 = (int) ceil((x+r)/bixel_dx); int bixel_i, bixel_j; sexpr *lns_in_bixel, *line; sexpr *lns_in_circle = nil; bixel_i0 = max(0,bixel_i0); bixel_j0 = max(0,bixel_j0); bixel_i1 = min(bixel_i1,bixel_rows); bixel_j1 = min(bixel_j1,bixel_cols); for (bixel_i = bixel_i0; bixel_i < bixel_i1; bixel_i++) { for (bixel_j = bixel_j0; bixel_j < bixel_j1; bixel_j++) { index = bixel_i*bixel_cols+bixel_j; lns_in_bixel = sketch->grid[index]; while (!null(lns_in_bixel)) { line = lns_in_bixel->u.pair->car; x0 = line->u.line->x[end]; y0 = line->u.line->y[end]; if (((x0-x)*(x0-x) + (y0-y)*(y0-y) < rr) && !member(line,lns_in_circle)) lns_in_circle=cons(line,lns_in_circle); lns_in_bixel=lns_in_bixel->u.pair->cdr; } } } return lns_in_circle; } sexpr *get_all_lns_in_circle(sketch *sketch, float x, float y, float r) { int index; float x0, y0, x1, y1; float rr=r*r; float bixel_dx = sketch->bixel_dx; float bixel_dy = sketch->bixel_dy; float bixel_rows = (float) sketch->bixel_rows; float bixel_cols = (float) sketch->bixel_cols; int bixel_i0 = (int) floor((y-r)/bixel_dy); int bixel_j0 = (int) floor((x-r)/bixel_dx); int bixel_i1 = (int) ceil((y+r)/bixel_dy); int bixel_j1 = (int) ceil((x+r)/bixel_dx); int bixel_i, bixel_j; sexpr *lns_in_bixel, *line; sexpr *lns_in_circle = nil; bixel_i0 = max(0,bixel_i0); bixel_j0 = max(0,bixel_j0); bixel_i1 = min(bixel_i1,bixel_rows); bixel_j1 = min(bixel_j1,bixel_cols); for (bixel_i = bixel_i0; bixel_i < bixel_i1; bixel_i++) { for (bixel_j = bixel_j0; bixel_j < bixel_j1; bixel_j++) { index = bixel_i*bixel_cols+bixel_j; lns_in_bixel = sketch->grid[index]; while (!null(lns_in_bixel)) { line = lns_in_bixel->u.pair->car; x0 = line->u.line->x[0]; y0 = line->u.line->y[0]; x1 = line->u.line->x[1]; y1 = line->u.line->y[1]; if (((x0-x)*(x0-x) + (y0-y)*(y0-y) < rr || (x1-x)*(x1-x) + (y1-y)*(y1-y) < rr || line_intersects_circle(x0,y0,x1,y1,x,y,r)) && !member(line,lns_in_circle)) lns_in_circle=cons(line,lns_in_circle); lns_in_bixel=lns_in_bixel->u.pair->cdr; } } } return lns_in_circle; } sexpr *po_ln_filter_magnitude(sexpr *ln, sexpr *lns, float contrast_ratio_min, float contrast_ratio_max) { float contrast = ln->u.line->contrast; float contrast_ratio; sexpr *previous; sexpr *current = lns; previous = nil; while (!null(current)) { contrast_ratio = car(current)->u.line->contrast/contrast; if ((contrast_ratio < contrast_ratio_min) || (contrast_ratio > contrast_ratio_max)) { if (null(previous)) { current = lns = lns->u.pair->cdr; continue; } else { previous->u.pair->cdr = current = current->u.pair->cdr; continue; } } previous = current; current = current->u.pair->cdr; } return lns; } sexpr *po_ln_filter_direction(sexpr *ln, sexpr *lns, float delta_theta_max) { float theta0 = ln->u.line->theta; float theta1; float delta_theta; sexpr *previous; sexpr *current = lns; previous = nil; while (!null(current)) { theta1 = current->u.pair->car->u.line->theta; delta_theta = (float) fabs(angle_difference((double) theta0,(double) theta1)); delta_theta = (float) min(delta_theta, TWOPI-delta_theta); if (delta_theta > delta_theta_max) { if (null(previous)) { current = lns = lns->u.pair->cdr; continue; } else { previous->u.pair->cdr = current = current->u.pair->cdr; continue; } } previous = current; current = current->u.pair->cdr; } return lns; } sexpr *po_ln_filter_end_project(sexpr *ln, int end, sexpr *lns, float end_proj_rel_max) { sexpr *previous; sexpr *current = lns; sexpr *first; float d1, d2; previous = nil; while (!null(current)) { first = car(current); /* No need to divide by length! */ d1 = proj_pt_on_line_rel_dist(ln,end,first->u.line->x[opposite(end)],first->u.line->y[opposite(end)]); d2 = proj_pt_on_line_rel_dist(ln,end,first->u.line->x[end],first->u.line->y[end]); if (d1 > end_proj_rel_max || d1 <= d2) { if (null(previous)) { current = lns = lns->u.pair->cdr; continue; } else { previous->u.pair->cdr = current = current->u.pair->cdr; continue; } } previous = current; current = current->u.pair->cdr; } return lns; } sexpr *po_ln_filter_lateral_distance(sexpr *ln, sexpr *lns, float lateral_dist_abs_max_sq) { sexpr *previous; sexpr *current = lns; sexpr *first; float xm, ym; float x0 = ln->u.line->x[0]; float y0 = ln->u.line->y[0]; float x1 = ln->u.line->x[1]; float y1 = ln->u.line->y[1]; previous = nil; while (!null(current)) { first = car(current); xm = (first->u.line->x[0] + first->u.line->x[1])/2.0; ym = (first->u.line->y[0] + first->u.line->y[1])/2.0; if (pt_line_dist_sq(x0,y0,x1,y1,xm,ym) > lateral_dist_abs_max_sq) { if (null(previous)) { current = lns = lns->u.pair->cdr; continue; } else { previous->u.pair->cdr = current = current->u.pair->cdr; continue; } } previous = current; current = current->u.pair->cdr; } return lns; } sexpr *delete(sexpr *item, sexpr *ls) { sexpr *previous; sexpr *current = ls; previous = nil; while (!null(current)) { if (eq(item,current->u.pair->car)) { if (null(previous)) { current = ls = ls->u.pair->cdr; continue; } else { previous->u.pair->cdr = current = current->u.pair->cdr; continue; } } previous = current; current = current->u.pair->cdr; } return ls; } void po_initialize_grid_for_ln_ends(sexpr *sp1, int bixel_rows, int bixel_cols) { int index, bixels; int rows = sp1->u.sketch->rows; int cols = sp1->u.sketch->cols; sexpr *ls = sp1->u.sketch->ls; sexpr *line, **grid; float bixel_dx, bixel_dy; float x0, y0, x1, y1; bixels = bixel_rows*bixel_cols; bixel_dx = (float) cols/(float) bixel_cols; bixel_dy = (float) rows/(float) bixel_rows; grid = malloc(bixels*sizeof(sexpr)); for (index=0; index < bixels; index++) grid[index]=nil; while (!null(ls)) { line = ls->u.pair->car; x0 = line->u.line->x[0]; y0 = line->u.line->y[0]; x1 = line->u.line->x[1]; y1 = line->u.line->y[1]; index = (int) floor(y0/bixel_dy)*bixel_cols+floor(x0/bixel_dx); if (index >= 0 && index < bixels) grid[index]=cons(line,grid[index]); index = (int) floor(y1/bixel_dy)*bixel_cols+floor(x1/bixel_dx); if (index >= 0 && index < bixels) grid[index] = cons(line,grid[index]); ls = ls->u.pair->cdr; } sp1->u.sketch->bixel_rows = bixel_rows; sp1->u.sketch->bixel_cols = bixel_cols; sp1->u.sketch->bixel_dx = bixel_dx; sp1->u.sketch->bixel_dy = bixel_dy; sp1->u.sketch->grid = grid; } void po_initialize_grid_for_all_lns(sexpr *sp1, int bixel_rows, int bixel_cols) { int index, bixels; int rows = sp1->u.sketch->rows; int cols = sp1->u.sketch->cols; sexpr *ls = sp1->u.sketch->ls; sexpr *line, **grid; float dx, dy, bixel_dx, bixel_dy; float x, y, x0, y0, x1, y1; bixels = bixel_rows*bixel_cols; bixel_dx = (float) cols/(float) bixel_cols; bixel_dy = (float) rows/(float) bixel_rows; grid=malloc(bixels*sizeof(sexpr)); for (index=0; index < bixels; index++) grid[index]=nil; while (!null(ls)) { line = ls->u.pair->car; x0 = line->u.line->x[0]; y0 = line->u.line->y[0]; x1 = line->u.line->x[1]; y1 = line->u.line->y[1]; dy = y1-y0; dx = x1-x0; x = x0; y = y0; if (fabs(dx) > fabs(dy)) { if (x <= x1) { while (x <= x1) { index = (int) floor(y/bixel_dy)*bixel_cols+floor(x/bixel_dx); if (index >= 0 && index < bixels) grid[index]=cons(line,grid[index]); x += bixel_dx; y = y0 + (x-x0)*(dy/dx); } } else { while (x >= x1) { index = (int) floor(y/bixel_dy)*bixel_cols+floor(x/bixel_dx); if (index >= 0 && index < bixels) grid[index]=cons(line,grid[index]); x -= bixel_dx; y = y0 + (x-x0)*(dy/dx); } } } else { if (y <= y1) { while (y <= y1) { index = (int) floor(y/bixel_dy)*bixel_cols+floor(x/bixel_dx); if (index >= 0 && index < bixels) grid[index]=cons(line,grid[index]); y += bixel_dy; x = x0 + (y-y0)*(dx/dy); } } else { while (y >= y1) { index = (int) floor(y/bixel_dy)*bixel_cols+floor(x/bixel_dx); if (index >= 0 && index < bixels) grid[index]=cons(line,grid[index]); y -= bixel_dy; x = x0 + (y-y0)*(dx/dy); } } } index = (int) floor(y1/bixel_dy)*bixel_cols+floor(x1/bixel_dx); if (index >= 0 && index < bixels) grid[index]=cons(line,grid[index]); ls=ls->u.pair->cdr; } sp1->u.sketch->bixel_rows = bixel_rows; sp1->u.sketch->bixel_cols = bixel_cols; sp1->u.sketch->bixel_dx = bixel_dx; sp1->u.sketch->bixel_dy = bixel_dy; sp1->u.sketch->grid = grid; } sexpr *po_link(sexpr *sp1, sexpr *sp2) { sexpr *ls = sp1->u.sketch->ls; sexpr **link_parameters = sp2->u.vector->vector; int end; sexpr *ln, *lns; int rows, cols, get_lns_mode; float link_radius, contrast_ratio_min, contrast_ratio_max; float delta_theta_max, end_proj_rel_max, lateral_dist_abs_max_sq; float x, y; if (!sketch(sp1)) { printf("po-link!: First argument must be sketch.\n"); longjmp(esc,1); } rows = sp1->u.sketch->rows; cols = sp1->u.sketch->cols; if (!vector(sp2) || sp2->u.vector->length != 7) { printf("po-link!: Link-parameters must be length seven vector.\n"); longjmp(esc,1); } if (!boolean(link_parameters[0])) { printf("po-link!: Get-lines-mode parameter must be boolean.\n"); longjmp(esc,1); } get_lns_mode = link_parameters[0]->u.i; if (!number(link_parameters[1])) { printf("po-link!: Link-radius must be number.\n"); longjmp(esc,1); } link_radius = link_parameters[1]->u.x; if (!number(link_parameters[2])) { printf("po-link!: Contrast-ratio-minimum parameter must be number.\n"); longjmp(esc,1); } contrast_ratio_min = link_parameters[2]->u.x; if (!number(link_parameters[3])) { printf("po-link!: Contrast-ratio-maximum parameter must be number.\n"); longjmp(esc,1); } contrast_ratio_max = link_parameters[3]->u.x; if (!number(link_parameters[4])) { printf("po-link!: Delta-theta-maximum parameter must be number.\n"); longjmp(esc,1); } delta_theta_max = link_parameters[4]->u.x; if (!number(link_parameters[5])) { printf("po-link!: End-project-relative-maximum parameter must be number.\n"); longjmp(esc,1); } end_proj_rel_max = link_parameters[5]->u.x; if (!number(link_parameters[6])) { printf("po-link!: Lateral-distance-absolute-maximum-squared parameter must be number.\n"); longjmp(esc,1); } lateral_dist_abs_max_sq = link_parameters[6]->u.x; if (get_lns_mode == 1) po_initialize_grid_for_all_lns(sp1,rows/(link_radius*2.0),cols/(link_radius*2.0)); else po_initialize_grid_for_ln_ends(sp1,rows/(link_radius*2.0),cols/(link_radius*2.0)); while (!null(ls)) { ln = ls->u.pair->car; for (end = 0; end <= 1; end++) { x = ln->u.line->x[end]; y = ln->u.line->y[end]; if (get_lns_mode == 1) lns=get_all_lns_in_circle(sp1->u.sketch,x,y,link_radius); else lns=get_all_lns_end_in_circle(sp1->u.sketch,x,y,opposite(end),link_radius); /* delete yourself */ lns = delete(ln,lns); lns = po_ln_filter_magnitude(ln,lns,contrast_ratio_min,contrast_ratio_max); lns = po_ln_filter_direction(ln,lns,delta_theta_max); lns = po_ln_filter_end_project(ln,end,lns,end_proj_rel_max); lns = po_ln_filter_lateral_distance(ln,lns,lateral_dist_abs_max_sq); ln->u.line->links[end]=lns; } ls = ls->u.pair->cdr; } return sp1; } sexpr *make_sketch(sexpr *rows, sexpr *cols, sexpr *ls) { sexpr *result; sexpr *ln; if (!number(rows) || !number(cols)) { printf("make-sketch: First and second arguments must be numbers.\n"); longjmp(esc,1); } result = smalloc++; result->type = SKETCH; result->u.sketch = malloc(sizeof(sketch)); result->u.sketch->rows = rows->u.x; result->u.sketch->cols = cols->u.x; result->u.sketch->bixel_rows = 0; result->u.sketch->bixel_cols = 0; result->u.sketch->count = 0; result->u.sketch->grid = NULL; result->u.sketch->ls = nil; while (!null(ls)) { if (!pair(ls)) { printf("make-sketch: Third argument must be list.\n"); longjmp(esc,1); } ln = ls->u.pair->car; if (!line(ln)) { printf("make-sketch: Third argument must be list of lines.\n"); longjmp(esc,1); } result->u.sketch->ls = cons(ln,result->u.sketch->ls); result->u.sketch->count++; ls = ls->u.pair->cdr; } return result; } void add_new_line(float x0, float y0, float x1, float y1, float contrast, float coverage, sexpr *sp1) { sexpr *sp2; float dx = x1-x0; float dy = y1-y0; sp2 = smalloc++; sp2->type = LINE; sp2->u.line = malloc(sizeof(line)); sp2->u.line->x[0] = x0; sp2->u.line->y[0] = y0; sp2->u.line->x[1] = x1; sp2->u.line->y[1] = y1; sp2->u.line->theta = atan2(dy,dx); sp2->u.line->contrast = contrast; sp2->u.line->length = sqrt(dx*dx+dy*dy); sp2->u.line->coverage = coverage; sp2->u.line->age = 0; sp2->u.line->replaced = 0; sp2->u.line->merged = 0; sp2->u.line->links[0] = nil; sp2->u.line->links[1] = nil; sp1->u.sketch->count++; sp1->u.sketch->ls=cons(sp2,sp1->u.sketch->ls); } void add_old_line(sexpr *sp0, sexpr *sp1) { sp0->u.line->age++; sp0->u.line->links[0] = nil; sp0->u.line->links[1] = nil; sp1->u.sketch->count++; sp1->u.sketch->ls=cons(sp0,sp1->u.sketch->ls); } sexpr *po_merge(sexpr *sp1, sexpr *sp2) { sexpr *ls, *ln; if (!sketch(sp1) || !sketch(sp2) || sp1->u.sketch->rows != sp2->u.sketch->rows || sp1->u.sketch->cols != sp2->u.sketch->cols) { printf("po-merge!: Arguments must be sketches with equal dimensions.\n"); longjmp(esc,1); } ls = sp2->u.sketch->ls; while (!null(ls)) { ln = ls->u.pair->car; if (ln->u.line->replaced == 0 && ln->u.line->merged == 0 && !member(ln,sp1->u.sketch->ls)) { ln->u.line->merged = 1; sp1->u.sketch->count++; sp1->u.sketch->ls = cons(ln,sp1->u.sketch->ls); } ls = ls->u.pair->cdr; } return sp1; } sexpr *sketch_union(sexpr *sp1, sexpr *sp2) { sexpr *result, *ls1, *ls2, *ln; if (!sketch(sp1) || !sketch(sp2) || sp1->u.sketch->rows != sp2->u.sketch->rows || sp1->u.sketch->cols != sp2->u.sketch->cols) { printf("sketch-union: Arguments must be sketches with equal dimensions.\n"); longjmp(esc,1); } result = make_sketch(num2exp(sp1->u.sketch->rows), num2exp(sp1->u.sketch->cols),nil); ls1 = sp1->u.sketch->ls; ls2 = sp2->u.sketch->ls; while (!null(ls1)) { ln = ls1->u.pair->car; result->u.sketch->count++; result->u.sketch->ls = cons(ln,result->u.sketch->ls); ls1 = ls1->u.pair->cdr; } while (!null(ls2)) { ln = ls2->u.pair->car; if (!member(ln,result->u.sketch->ls)) { result->u.sketch->count++; result->u.sketch->ls = cons(ln,result->u.sketch->ls); } ls2 = ls2->u.pair->cdr; } return result; } sexpr *sketch_lines(sexpr *sp) { if (!sketch(sp)) { printf("sketch-lines: Argument is not a sketch.\n"); longjmp(esc,1); } return sp->u.sketch->ls; } sexpr *sketch_rows(sexpr *sp) { if (!sketch(sp)) { printf("sketch-rows: Argument is not a sketch.\n"); longjmp(esc,1); } return num2exp(sp->u.sketch->rows); } sexpr *sketch_cols(sexpr *sp) { if (!sketch(sp)) { printf("sketch-cols: Argument is not a sketch.\n"); longjmp(esc,1); } return num2exp(sp->u.sketch->cols); } sexpr *sketch_count(sexpr *sp) { if (!sketch(sp)) { printf("sketch-count: Argument is not a sketch.\n"); longjmp(esc,1); } return num2exp(sp->u.sketch->count); } sexpr *line_x0(sexpr *sp) { if (!line(sp)) { printf("line-x0: Argument is not a line.\n"); longjmp(esc,1); } return num2exp(sp->u.line->x[0]); } sexpr *line_y0(sexpr *sp) { if (!line(sp)) { printf("line-y0: Argument is not a line.\n"); longjmp(esc,1); } return num2exp(sp->u.line->y[0]); } sexpr *line_x1(sexpr *sp) { if (!line(sp)) { printf("line-x1: Argument is not a line.\n"); longjmp(esc,1); } return num2exp(sp->u.line->x[1]); } sexpr *line_y1(sexpr *sp) { if (!line(sp)) { printf("line-y1: Argument is not a line.\n"); longjmp(esc,1); } return num2exp(sp->u.line->y[1]); } sexpr *line_contrast(sexpr *sp) { if (!line(sp)) { printf("line-contrast: Argument is not a line.\n"); longjmp(esc,1); } return num2exp(sp->u.line->contrast); } sexpr *line_theta(sexpr *sp) { if (!line(sp)) { printf("line-theta: Argument is not a line.\n"); longjmp(esc,1); } return num2exp(sp->u.line->theta); } sexpr *line_length(sexpr *sp) { if (!line(sp)) { printf("line-length: Argument is not a line.\n"); longjmp(esc,1); } return num2exp(sp->u.line->length); } sexpr *line_links(sexpr *sp) { if (!line(sp)) { printf("line-links: Argument is not a line.\n"); longjmp(esc,1); } return list2(sp->u.line->links[0],sp->u.line->links[1]); } sexpr *set_line_contrast(sexpr *sp1, sexpr *sp2) { if (!line(sp1)) { printf("set-line-contrast!: First argument is not a line.\n"); longjmp(esc,1); } if (!number(sp2)) { printf("set-line-contrast!: Second argument is not a number.\n"); longjmp(esc,1); } sp1->u.line->contrast = sp2->u.x; return undefined; } sexpr *po_replace(sexpr *sp1, sexpr *sp2) { sexpr **replace_parameters = sp2->u.vector->vector; int i; float x, y, xm, ym, dx, dy; float x0, y0, x1, y1, u0, v0, u1, v1; float x00, y00, x01, y01, x10, y10, x11, y11; float u00, v00, u01, v01, u10, v10, u11, v11; float dot00, dot01, dot0, dot1, dot10, dot11; float cumulative_coverage, coverage; float scale, mse, best_mse; float best_avg_contrast = 0, best_cumulative_coverage = 0; float c, c0, c1, l, l0, l1, ll; float m, m0, m1, n0, n1, o0, o1; float k, k0, k1; sexpr *best_ln = nil, *best_ln0 = nil, *best_ln1 = nil; float best_u00 = 0, best_v00 = 0, best_u0 = 0, best_v0 = 0; float best_u1 = 0, best_v1 = 0, best_u11 = 0, best_v11 = 0; int new, copied, paths, path_histogram[10]; sexpr *ls, *ls0, *ls1; sexpr *ln, *ln0, *ln1; double **A, *s; sexpr *result; float replace_radius; float straightness; int coverage_filter_on; float min_cumulative_coverage; float min_coverage; int delta_abst_lev_max; int replace_straightest_path; float rr; if (!sketch(sp1)) { printf("po-replace!: First argument must be sketch.\n"); longjmp(esc,1); } result = make_sketch(num2exp(sp1->u.sketch->rows),num2exp(sp1->u.sketch->cols),nil); if (!vector(sp2) || sp2->u.vector->length != 7) { printf("po-replace!: Replace parameters must be length seven vector.\n"); longjmp(esc,1); } if (!number(replace_parameters[0])) { printf("po-replace!: Replace-radius must be number.\n"); longjmp(esc,1); } replace_radius = replace_parameters[0]->u.x; if (!number(replace_parameters[1])) { printf("po-replace!: Straightness-threshold must be number.\n"); longjmp(esc,1); } straightness = replace_parameters[1]->u.x; if (!boolean(replace_parameters[2])) { printf("po-replace!: Coverage-filter-on parameter must be boolean.\n"); longjmp(esc,1); } coverage_filter_on = replace_parameters[2]->u.i; if (!number(replace_parameters[3])) { printf("po-replace!: Minimum-cumulative-coverage parameter must be number.\n"); longjmp(esc,1); } min_cumulative_coverage = replace_parameters[3]->u.x; if (!number(replace_parameters[4])) { printf("po-replace!: Minimum-coverage-parameter must be number.\n"); longjmp(esc,1); } min_coverage = replace_parameters[4]->u.x; if (!number(replace_parameters[5])) { printf("po-replace!: Delta-abstraction-level-maximum parameter must be number.\n"); longjmp(esc,1); } delta_abst_lev_max = (int) replace_parameters[5]->u.x; if (!boolean(replace_parameters[6])) { printf("po-replace!: replace-straightest-path parameter must be boolean.\n"); longjmp(esc,1); } replace_straightest_path = replace_parameters[6]->u.i; rr = replace_radius*replace_radius; A = (double **) malloc(4*sizeof(double *)); s = (double *) malloc(2*sizeof(double)); new = copied = 0; for (i = 0; i < 4; i++) A[i] = malloc(2*sizeof(double)); for (i = 0; i <= 9; i++) path_histogram[i] = 0; /* Enumerate length three paths */ ls = sp1->u.sketch->ls; while (!null(ls)) { paths = 0; best_mse = FLT_MAX; ln = ls->u.pair->car; if (ln->u.line->replaced) { ls = ls->u.pair->cdr; continue; } x0 = ln->u.line->x[0]; y0 = ln->u.line->y[0]; x1 = ln->u.line->x[1]; y1 = ln->u.line->y[1]; x = (x0 + x1)/2.0; y = (y0 + y1)/2.0; ls0 = ln->u.line->links[0]; while (!null(ls0)) { ln0 = ls0->u.pair->car; if (ln0->u.line->replaced) { ls0 = ls0->u.pair->cdr; continue; } x00 = ln0->u.line->x[0]; y00 = ln0->u.line->y[0]; if (((x-x00)*(x-x00) + (y-y00)*(y-y00)) > rr) { ls0 = ls0->u.pair->cdr; continue; } x01 = ln0->u.line->x[1]; y01 = ln0->u.line->y[1]; ls1 = ln->u.line->links[1]; while (!null(ls1)) { ln1 = ls1->u.pair->car; if (ln1->u.line->replaced) { ls1 = ls1->u.pair->cdr; continue; } x11 = ln1->u.line->x[1]; y11 = ln1->u.line->y[1]; if (((x-x11)*(x-x11) + (y-y11)*(y-y11)) > rr) { ls1 = ls1->u.pair->cdr; continue; } x10 = ln1->u.line->x[0]; y10 = ln1->u.line->y[0]; xm = (x00 + x01 + x0 + x1 + x10 + x11)/6.0; ym = (y00 + y01 + y0 + y1 + y10 + y11)/6.0; A[0][0] = (x00-xm)*(x00-xm) + (x01-xm)*(x01-xm) + (x0-xm)*(x0-xm) + (x1-xm)*(x1-xm) + (x10-xm)*(x10-xm) + (x11-xm)*(x11-xm); A[1][1] = (y00-ym)*(y00-ym) + (y01-ym)*(y01-ym) + (y0-ym)*(y0-ym) + (y1-ym)*(y1-ym) + (y10-ym)*(y10-ym) + (y11-ym)*(y11-ym); A[0][1] = A[1][0] = (x00-xm)*(y00-ym) + (x01-xm)*(y01-ym) + (x0-xm)*(y0-ym) + (x1-xm)*(y1-ym) + (x10-xm)*(y10-ym) + (x11-xm)*(y11-ym); svd(A,s,2); l = A[2][0]*A[2][0] + A[2][1]*A[2][1]; dx = -A[2][0]/l; dy = A[2][1]/l; dot00 = dx*(x00-xm)+dy*(y00-ym); u00 = xm + dot00*dx; v00 = ym + dot00*dy; dot01 = dx*(x01-xm)+dy*(y01-ym); u01 = xm + dot01*dx; v01 = ym + dot01*dy; dot0 = dx*(x0-xm)+dy*(y0-ym); u0 = xm + dot0*dx; v0 = ym + dot0*dy; dot1 = dx*(x1-xm)+dy*(y1-ym); u1 = xm + dot1*dx; v1 = ym + dot1*dy; dot10 = dx*(x10-xm)+dy*(y10-ym); u10 = xm + dot10*dx; v10 = ym + dot10*dy; dot11 = dx*(x11-xm)+dy*(y11-ym); u11 = xm + dot11*dx; v11 = ym + dot11*dy; l0 = ln0->u.line->length; l = ln->u.line->length; l1 = ln1->u.line->length; scale = ((x00-x11)*(x00-x11) + (y00-y11)*(y00-y11))*4.0; mse = (pt_line_dist_sq(u00,v00,u11,v11,x00,y00)*l0 + pt_line_dist_sq(u00,v00,u11,v11,x01,y01)*l0 + pt_line_dist_sq(u00,v00,u11,v11,x0,y0)*l + pt_line_dist_sq(u00,v00,u11,v11,x1,y1)*l + pt_line_dist_sq(u00,v00,u11,v11,x10,y10)*l1 + pt_line_dist_sq(u00,v00,u11,v11,x11,y11)*l1)/scale; paths++; if (mse < best_mse && mse < straightness) { m0 = sqrt((u01-u00)*(u01-u00)+(v01-v00)*(v01-v00)); m = sqrt((u1-u0)*(u1-u0)+(v1-v0)*(v1-v0)); m1 = sqrt((u11-u10)*(u11-u10)+(v11-v10)*(v11-v10)); n0 = sqrt((u00-u0)*(u00-u0)+(v00-v0)*(v00-v0)); n1 = sqrt((u11-u1)*(u11-u1)+(v11-v1)*(v11-v1)); if (n0 < m0) o0 = m0 - n0; else o0 = 0.0; if (n1 < m1) o1 = m1 - n1; else o1 = 0.0; k0 = ln0->u.line->coverage; k = ln->u.line->coverage; k1 = ln1->u.line->coverage; ll = sqrt((u11-u00)*(u11-u00)+(v11-v00)*(v11-v00)); cumulative_coverage = (m0*k0 + m*k + m1*k1 - o0 - o1)/ll; coverage = (m0 + m + m1 - o0 - o1)/ll; if (coverage_filter_on == 0 || (coverage_filter_on == 1 && cumulative_coverage > min_cumulative_coverage && coverage > min_coverage)) { best_mse = mse; c0 = ln0->u.line->contrast; c = ln->u.line->contrast; c1 = ln1->u.line->contrast; best_avg_contrast = (c0*l0 + c*l + c1*l1)/(l0 + l + l1); best_cumulative_coverage = cumulative_coverage; best_u00 = u00; best_v00 = v00; best_u11 = u11; best_v11 = v11; best_ln0 = ln0; best_ln = ln; best_ln1 = ln1; if (!replace_straightest_path) goto break_break; } } ls1 = ls1->u.pair->cdr; } ls0 = ls0->u.pair->cdr; } break_break: if (best_mse < FLT_MAX) { new++; add_new_line(best_u00,best_v00,best_u11,best_v11,best_avg_contrast,best_cumulative_coverage,result); best_ln0->u.line->replaced = 1; best_ln->u.line->replaced = 1; best_ln1->u.line->replaced = 1; if (paths < 9) path_histogram[paths]++; else path_histogram[10]++; } ls = ls->u.pair->cdr; } /* printf("length three path histogram: \n"); for (i = 0; i < 10; i++) printf("%d(%d) ",i,path_histogram[i]); printf("\n"); */ for (i = 0; i <= 9; i++) path_histogram[i] = 0; /* Enumerate length two left-paths */ ls = sp1->u.sketch->ls; while (!null(ls)) { paths = 0; best_mse = FLT_MAX; ln = ls->u.pair->car; if (ln->u.line->replaced) { ls = ls->u.pair->cdr; continue; } x0 = ln->u.line->x[0]; y0 = ln->u.line->y[0]; x1 = ln->u.line->x[1]; y1 = ln->u.line->y[1]; x = (x0 + x1)/2.0; y = (y0 + y1)/2.0; ls0 = ln->u.line->links[0]; while(!null(ls0)) { ln0 = ls0->u.pair->car; if (ln0->u.line->replaced) { ls0 = ls0->u.pair->cdr; continue; } x00 = ln0->u.line->x[0]; y00 = ln0->u.line->y[0]; if (((x-x00)*(x-x00) + (y-y00)*(y-y00)) > rr) { ls0 = ls0->u.pair->cdr; continue; } x01 = ln0->u.line->x[1]; y01 = ln0->u.line->y[1]; xm = (x00 + x01 + x0 + x1)/4.0; ym = (y00 + y01 + y0 + y1)/4.0; A[0][0] = (x00-xm)*(x00-xm) + (x01-xm)*(x01-xm) + (x0-xm)*(x0-xm) + (x1-xm)*(x1-xm); A[1][1] = (y00-ym)*(y00-ym) + (y01-ym)*(y01-ym) + (y0-ym)*(y0-ym) + (y1-ym)*(y1-ym); A[0][1] = A[1][0] = (x00-xm)*(y00-ym) + (x01-xm)*(y01-ym) + (x0-xm)*(y0-ym) + (x1-xm)*(y1-ym); svd(A,s,2); l = A[2][0]*A[2][0] + A[2][1]*A[2][1]; dx = -A[2][0]/l; dy = A[2][1]/l; dot00 = dx*(x00-xm)+dy*(y00-ym); u00 = xm + dot00*dx; v00 = ym + dot00*dy; dot01 = dx*(x01-xm)+dy*(y01-ym); u01 = xm + dot01*dx; v01 = ym + dot01*dy; dot0 = dx*(x0-xm)+dy*(y0-ym); u0 = xm + dot0*dx; v0 = ym + dot0*dy; dot1 = dx*(x1-xm)+dy*(y1-ym); u1 = xm + dot1*dx; v1 = ym + dot1*dy; l0 = ln0->u.line->length; l = ln->u.line->length; scale = ((x00-x1)*(x00-x1) + (y00-y1)*(y00-y1))*2.0; mse = (pt_line_dist_sq(u00,v00,u1,v1,x00,y00)*l0 + pt_line_dist_sq(u00,v00,u1,v1,x01,y01)*l0 + pt_line_dist_sq(u00,v00,u1,v1,x0,y0)*l + pt_line_dist_sq(u00,v00,u1,v1,x1,y1)*l)/scale; paths++; if (mse < best_mse && mse < straightness) { m0 = sqrt((u01-u00)*(u01-u00)+(v01-v00)*(v01-v00)); m = sqrt((u1-u0)*(u1-u0)+(v1-v0)*(v1-v0)); n0 = sqrt((u00-u0)*(u00-u0)+(v00-v0)*(v00-v0)); if (n0 < m0) o0 = m0 - n0; else o0 = 0.0; k0 = ln0->u.line->coverage; k = ln->u.line->coverage; ll = sqrt((u1-u00)*(u1-u00)+(v1-v00)*(v1-v00)); cumulative_coverage = (m0*k0 + m*k - o0)/ll; coverage = (m0 + m - o0)/ll; if (coverage_filter_on == 0 || (coverage_filter_on == 1 && cumulative_coverage > min_cumulative_coverage && coverage > min_coverage)) { best_mse = mse; c0 = ln0->u.line->contrast; c = ln->u.line->contrast; best_avg_contrast = (c0*l0 + c*l)/(l0 + l); best_cumulative_coverage = cumulative_coverage; best_u00 = u00; best_v00 = v00; best_u1 = u1; best_v1 = v1; best_ln0 = ln0; best_ln = ln; if (!replace_straightest_path) break; } } ls0 = ls0->u.pair->cdr; } if (best_mse < FLT_MAX) { new++; add_new_line(best_u00,best_v00,best_u1,best_v1,best_avg_contrast,best_cumulative_coverage,result); best_ln0->u.line->replaced = 1; best_ln->u.line->replaced = 1; if (paths < 9) path_histogram[paths]++; else path_histogram[10]++; } ls = ls->u.pair->cdr; } /* Enumerate length two right-paths */ ls = sp1->u.sketch->ls; while (!null(ls)) { paths = 0; best_mse = FLT_MAX; ln = ls->u.pair->car; if (ln->u.line->replaced) { ls = ls->u.pair->cdr; continue; } x0 = ln->u.line->x[0]; y0 = ln->u.line->y[0]; x1 = ln->u.line->x[1]; y1 = ln->u.line->y[1]; x = (x0 + x1)/2.0; y = (y0 + y1)/2.0; ls1 = ln->u.line->links[1]; while (!null(ls1)) { ln1 = ls1->u.pair->car; if (ln1->u.line->replaced) { ls1 = ls1->u.pair->cdr; continue; } x11 = ln1->u.line->x[1]; y11 = ln1->u.line->y[1]; if (((x-x11)*(x-x11) + (y-y11)*(y-y11)) > rr) { ls1 = ls1->u.pair->cdr; continue; } x10 = ln1->u.line->x[0]; y10 = ln1->u.line->y[0]; xm = (x0 + x1 + x10 + x11)/4.0; ym = (y0 + y1 + y10 + y11)/4.0; A[0][0] = (x0-xm)*(x0-xm) + (x1-xm)*(x1-xm) + (x10-xm)*(x10-xm) + (x11-xm)*(x11-xm); A[1][1] = (y0-ym)*(y0-ym) + (y1-ym)*(y1-ym) + (y10-ym)*(y10-ym) + (y11-ym)*(y11-ym); A[0][1] = A[1][0] = (x0-xm)*(y0-ym) + (x1-xm)*(y1-ym) + (x10-xm)*(y10-ym) + (x11-xm)*(y11-ym); svd(A,s,2); l = A[2][0]*A[2][0] + A[2][1]*A[2][1]; dx = -A[2][0]/l; dy = A[2][1]/l; dot0 = dx*(x0-xm)+dy*(y0-ym); u0 = xm + dot0*dx; v0 = ym + dot0*dy; dot1 = dx*(x1-xm)+dy*(y1-ym); u1 = xm + dot1*dx; v1 = ym + dot1*dy; dot10 = dx*(x10-xm)+dy*(y10-ym); u10 = xm + dot10*dx; v10 = ym + dot10*dy; dot11 = dx*(x11-xm)+dy*(y11-ym); u11 = xm + dot11*dx; v11 = ym + dot11*dy; l = ln->u.line->length; l1 = ln1->u.line->length; scale = ((x0-x11)*(x0-x11) + (y0-y11)*(y0-y11))*2.0; mse = (pt_line_dist_sq(u0,v0,u11,v11,x0,y0)*l + pt_line_dist_sq(u0,v0,u11,v11,x1,y1)*l + pt_line_dist_sq(u0,v0,u11,v11,x10,y10)*l1 + pt_line_dist_sq(u0,v0,u11,v11,x11,y11)*l1)/scale; paths++; if (mse < best_mse && mse < straightness) { m = sqrt((u1-u0)*(u1-u0)+(v1-v0)*(v1-v0)); m1 = sqrt((u11-u10)*(u11-u10)+(v11-v10)*(v11-v10)); n1 = sqrt((u11-u1)*(u11-u1)+(v11-v1)*(v11-v1)); if (n1 < m1) o1 = m1 - n1; else o1 = 0.0; k = ln->u.line->coverage; k1 = ln1->u.line->coverage; ll = sqrt((u11-u0)*(u11-u0)+(v11-v0)*(v11-v0)); cumulative_coverage = (m*k + m1*k1 - o1)/ll; coverage = (m + m1 - o1)/ll; if (coverage_filter_on == 0 || (coverage_filter_on == 1 && cumulative_coverage > min_cumulative_coverage && coverage > min_coverage)) { best_mse = mse; c = ln->u.line->contrast; c1 = ln1->u.line->contrast; best_avg_contrast = (c*l + c1*l1)/(l + l1); best_cumulative_coverage = cumulative_coverage; best_u0 = u0; best_v0 = v0; best_u11 = u11; best_v11 = v11; best_ln = ln; best_ln1 = ln1; if (!replace_straightest_path) break; } } ls1 = ls1->u.pair->cdr; } if (best_mse < FLT_MAX) { new++; add_new_line(best_u0,best_v0,best_u11,best_v11,best_avg_contrast,best_cumulative_coverage,result); best_ln->u.line->replaced = 1; best_ln1->u.line->replaced = 1; if (paths < 9) path_histogram[paths]++; else path_histogram[10]++; } ls = ls->u.pair->cdr; } /* printf("length two path histogram: \n"); for (i = 0; i < 10; i++) printf("%d(%d) ",i,path_histogram[i]); printf("\n"); */ /* Copy non-replaced lines */ ls = sp1->u.sketch->ls; while (!null(ls)) { ln = ls->u.pair->car; if (ln->u.line->replaced == 0 && ln->u.line->age < delta_abst_lev_max) { add_old_line(ln,result); copied++; } ls = ls->u.pair->cdr; } /* printf("new: %d copied: %d total: %d\n",new,copied,new+copied); */ return result; } void add_zc0(float x0, float y0, float x1, float y1, float contrast, sexpr *sp1) { sexpr *sp2; float dx = x1-x0; float dy = y1-y0; sp2 = smalloc++; sp2->type = LINE; sp2->u.line = malloc(sizeof(line));; sp2->u.line->x[0] = x0; sp2->u.line->y[0] = y0; sp2->u.line->x[1] = x1; sp2->u.line->y[1] = y1; sp2->u.line->theta = atan2(dy,dx); sp2->u.line->contrast = contrast; sp2->u.line->length = sqrt(dx*dx + dy*dy); sp2->u.line->coverage = 1.0; sp2->u.line->age = 0; sp2->u.line->replaced = 0; sp2->u.line->merged = 0; sp2->u.line->links[0] = nil; sp2->u.line->links[1] = nil; sp1->u.sketch->count++; sp1->u.sketch->ls=cons(sp2,sp1->u.sketch->ls); } void add_zc1(float x, float y, float theta, float contrast, sexpr *sp1) { sexpr *sp2; float dx = cos(theta)/2.0; float dy = sin(theta)/2.0; sp2 = smalloc++; sp2->type = LINE; sp2->u.line = malloc(sizeof(line)); sp2->u.line->x[0] = x - dx; sp2->u.line->y[0] = y - dy; sp2->u.line->x[1] = x + dx; sp2->u.line->y[1] = y + dy; sp2->u.line->theta = theta; sp2->u.line->contrast = contrast; sp2->u.line->length = 1.0; sp2->u.line->coverage = 1.0; sp2->u.line->age = 0; sp2->u.line->replaced = 0; sp2->u.line->merged = 0; sp2->u.line->links[0] = nil; sp2->u.line->links[1] = nil; sp1->u.sketch->count++; sp1->u.sketch->ls=cons(sp2,sp1->u.sketch->ls); } sexpr *make_line(sexpr *x0, sexpr *y0, sexpr *x1, sexpr *y1, sexpr *contrast) { sexpr *result; if (!number(x0) || !number(y0) || !number(x1) || !number(y1) || !number(contrast)) { printf("make-line: Illegal argument.\n"); longjmp(esc,1); } result = smalloc++; result->type = LINE; result->u.sketch = malloc(sizeof(line)); result->u.line->x[0] = x0->u.x; result->u.line->y[0] = y0->u.x; result->u.line->x[1] = x1->u.x; result->u.line->y[1] = y1->u.x; result->u.line->theta = atan2(y1->u.x-y0->u.x,x1->u.x-x0->u.x); result->u.line->contrast = contrast->u.x; result->u.line->length = sqrt((x1->u.x-x0->u.x)*(x1->u.x-x0->u.x) + (y1->u.x-y0->u.x)*(y1->u.x-y0->u.x)); result->u.line->coverage = 1.0; result->u.line->age = 0; result->u.line->replaced = 0; result->u.line->links[0] = nil; result->u.line->links[1] = nil; return result; } #define interpolate(f00,f01,f10,f11,x,y) (f01-f00)*x+(f10-f00)*y+(f11+f00-f10-f01)*x*y+f00; sexpr *zerocrossings(sexpr *sp1, sexpr *sp2, sexpr *sp3, sexpr *sp4) { int i, j, rows, cols; unsigned code; float *lapl, *dx, *dy; float contrast; float lapl00, lapl01, lapl10, lapl11; float dx00, dx01, dx10, dx11; float dy00, dy01, dy10, dy11; float x0, y0, x1, y1, dxc, dyc, xc, yc; int method; sexpr *result; x0 = y0 = x1 = y1 = 0.0; if (!image(sp1) || !image(sp2) || !image(sp3)) { printf("zero-crossings: First three arguments must be images.\n"); longjmp(esc,1); } rows=sp1->u.image->rows; cols=sp1->u.image->cols; if (sp2->u.image->rows != rows || sp2->u.image->cols != cols || sp3->u.image->rows != rows || sp3->u.image->cols != cols) { printf("zero-crossings: Image argument size mismatch.\n"); longjmp(esc,1); } if (!number(sp4)) { printf("zero-crossings: Method parameter must be number.\n"); longjmp(esc,1); } method = sp4->u.x; if (method != 0 && method != 1 && method != 2) { printf("zero-crossings: Method parameter must equal 0, 1, or 2.\n"); longjmp(esc,1); } result = smalloc++; result->type = SKETCH; result->u.sketch = malloc(sizeof(sketch)); result->u.sketch->rows=rows; result->u.sketch->cols=cols; result->u.sketch->bixel_rows = 0; result->u.sketch->bixel_cols = 0; result->u.sketch->count=0; result->u.sketch->grid=NULL; result->u.sketch->ls=nil; lapl = sp1->u.image->data; dx = sp2->u.image->data; dy = sp3->u.image->data; for (i=0; i < rows-1; i++) { for (j=0; j < cols-1; j++) { /* 00 -- 01 */ /* | | */ /* 10 -- 11 */ lapl00 = lapl[i*cols + j]; lapl01 = lapl[i*cols + j + 1]; lapl10 = lapl[(i+1)*cols + j]; lapl11 = lapl[(i+1)*cols + j + 1]; code = 0; if (lapl00*lapl01*lapl10*lapl11 != 0) { if (lapl00 > 0) code += 1; if (lapl01 > 0) code += 2; if (lapl10 > 0) code += 4; if (lapl11 > 0) code += 8; } /* + to the right */ switch (code) { case 0: /* - - */ /* - - */ continue; case 1: /* + - */ /* - - */ x0 = j+lapl00/(lapl00-lapl01); y0 = i; x1 = j; y1 = i+lapl00/(lapl00-lapl10); break; case 2: /* - + */ /* - - */ x0 = j+1; y0 = i+lapl01/(lapl01-lapl11); x1 = j+lapl00/(lapl00-lapl01); y1 = i; break; case 3: /* + + */ /* - - */ x0 = j+1; y0 = i+lapl01/(lapl01-lapl11); x1 = j; y1 = i+lapl00/(lapl00-lapl10); break; case 4: /* - - */ /* + - */ x0 = j; y0 = i+lapl00/(lapl00-lapl10); x1 = j+lapl10/(lapl10-lapl11); y1 = i+1; break; case 5: /* + - */ /* + - */ x0 = j+lapl00/(lapl00-lapl01); y0 = i; x1 = j+lapl10/(lapl10-lapl11); y1 = i+1; break; case 6: /* - + */ /* + - */ break; case 7: /* + + */ /* + - */ x0 = j+1; y0 = i+lapl01/(lapl01-lapl11); x1 = j+lapl10/(lapl10-lapl11); y1 = i+1; break; case 8: /* - - */ /* - + */ x0 = j+lapl10/(lapl10-lapl11); y0 = i+1; x1 = j+1; y1 = i+lapl01/(lapl01-lapl11); break; case 9: /* + - */ /* - + */ continue; case 10: /* - + */ /* - + */ x0 = j+lapl10/(lapl10-lapl11); y0 = i+1; x1 = j+lapl00/(lapl00-lapl01); y1 = i; break; case 11: /* + + */ /* - + */ x0 = j+lapl10/(lapl10-lapl11); y0 = i+1; x1 = j; y1 = i+lapl00/(lapl00-lapl10); break; case 12: /* - - */ /* + + */ x0 = j; y0 = i+lapl00/(lapl00-lapl10); x1 = j+1; y1 = i+lapl01/(lapl01-lapl11); break; case 13: /* + - */ /* + + */ x0 = j+lapl00/(lapl00-lapl01); y0 = i; x1 = j+1; y1 = i+lapl01/(lapl01-lapl11); break; case 14: /* - + */ /* + + */ x0 = j; y0 = i+lapl00/(lapl00-lapl10); x1 = j+lapl00/(lapl00-lapl01); y1 = i; break; case 15: /* + + */ /* + + */ continue; default: continue; } xc=(x0+x1)/2-j; yc=(y0+y1)/2-i; dx00 = dx[i*cols + j]; dx01 = dx[i*cols + j + 1]; dx10 = dx[(i+1)*cols + j]; dx11 = dx[(i+1)*cols + j + 1]; dy00 = dy[i*cols + j]; dy01 = dy[i*cols + j + 1]; dy10 = dy[(i+1)*cols + j]; dy11 = dy[(i+1)*cols + j + 1]; dxc = interpolate(dx00,dx01,dx10,dx11,xc,yc); dyc = interpolate(dy00,dy01,dy10,dy11,xc,yc); contrast = sqrt(dxc*dxc + dyc*dyc); switch (method) { case 0: add_zc0(x0,y0,x1,y1,contrast,result); break; case 1: add_zc1((x0+x1)/2,(y0+y1)/2,atan2(-dxc,dyc),contrast,result); break; case 2: add_zc1(x0,y0,atan2(-dxc,dyc),contrast,result); } } } return result; } sexpr *filter_on_contrast(sexpr *sp1, sexpr *sp2, sexpr *sp3) { sexpr *result; sexpr *ls1, *ls, *first; float lower, upper, contrast; int count; if (!sketch(sp1) || !number(sp2) || !number(sp3)) { printf("filter-on-contrast: Illegal arguments.\n"); longjmp(esc,1); } lower = sp2->u.x; upper = sp3->u.x; result = smalloc++; result->type = SKETCH; result->u.sketch = malloc(sizeof(sketch)); result->u.sketch->rows = sp1->u.sketch->rows; result->u.sketch->cols = sp1->u.sketch->cols; result->u.sketch->bixel_rows = 0; result->u.sketch->bixel_cols = 0; result->u.sketch->grid = NULL; count = 0; ls = nil; ls1 = sp1->u.sketch->ls; while (!null(ls1)) { first = ls1->u.pair->car; contrast = first->u.line->contrast; if (contrast > lower && contrast < upper) { count++; ls = cons(first,ls); } ls1 = ls1->u.pair->cdr; } result->u.sketch->ls=ls; result->u.sketch->count=count; return result; } sexpr *filter_on_length(sexpr *sp1, sexpr *sp2, sexpr *sp3) { sexpr *result; sexpr *ls1, *ls, *first; float lower, upper, length; int count; if (!sketch(sp1) || !number(sp2) || !number(sp3)) { printf("filter-on-length: Illegal arguments.\n"); longjmp(esc,1); } lower = sp2->u.x; upper = sp3->u.x; result = smalloc++; result->type = SKETCH; result->u.sketch = malloc(sizeof(sketch)); result->u.sketch->rows = sp1->u.sketch->rows; result->u.sketch->cols = sp1->u.sketch->cols; result->u.sketch->bixel_rows = 0; result->u.sketch->bixel_cols = 0; result->u.sketch->grid = NULL; count = 0; ls = nil; ls1 = sp1->u.sketch->ls; while (!null(ls1)) { first = ls1->u.pair->car; length = first->u.line->length; if (length > lower && length < upper) { count++; ls = cons(first,ls); } ls1 = ls1->u.pair->cdr; } result->u.sketch->ls=ls; result->u.sketch->count=count; return result; } sexpr *filter_on_angle(sexpr *sp1, sexpr *sp2, sexpr *sp3) { sexpr *result; sexpr *ls1, *ls, *first; double lower, upper, theta; int count; if (!sketch(sp1) || !number(sp2) || !number(sp3)) { printf("filter-on-angle: Illegal arguments.\n"); longjmp(esc,1); } lower = sp2->u.x; upper = sp3->u.x; result = smalloc++; result->type = SKETCH; result->u.sketch = malloc(sizeof(sketch)); result->u.sketch->rows= sp1->u.sketch->rows; result->u.sketch->cols= sp1->u.sketch->cols; result->u.sketch->bixel_rows = 0; result->u.sketch->bixel_cols = 0; result->u.sketch->grid = NULL; count = 0; ls = nil; ls1 = sp1->u.sketch->ls; while (!null(ls1)) { first = ls1->u.pair->car; theta = (double) first->u.line->theta; if (angle_difference(theta, lower) <= PI && angle_difference(theta, upper) >= PI) { count++; ls = cons(first,ls); } ls1 = ls1->u.pair->cdr; } result->u.sketch->ls=ls; result->u.sketch->count=count; return result; } void generate_postscript_header(FILE *file, float ratio, char *title) { fprintf(file,"%%!PS-Adobe-2.0 EPSF-1.2\n"); fprintf(file,"%%%%Title: %s\n",title); fprintf(file,"%%%%Creator: University of New Mexico Scheme 2.0\n"); if (ratio < 1.0) fprintf(file,"%%%%BoundingBox: 0 %g 432 432\n",432*(1.0-ratio)); else fprintf(file,"%%%%BoundingBox: 0 0 %g 432\n",432/ratio); fprintf(file,"%%%%Pages: 1\n"); fprintf(file,"%%%%EndComments\n"); fprintf(file,"save\n"); fprintf(file,"/inch {72 mul} def\n"); fprintf(file,"/displine { /y2 exch def /x2 exch def /y1 exch def /x1 exch def\n"); fprintf(file,"x1 y1 moveto x2 y2 lineto stroke } def\n"); fprintf(file,"%%%%EndProlog\n"); fprintf(file,"%%%%%%Page: 1 1\n"); fprintf(file,"0 0 translate\n"); fprintf(file,"432 432 scale\n"); fprintf(file,"0.00115741 setlinewidth\n"); fprintf(file,"0.5 dup translate -90 rotate -0.5 dup translate\n"); if (ratio < 1.0) fprintf(file,"0 0 moveto 0 1 lineto %f 1 lineto %f 0 lineto closepath\n",ratio,ratio); else fprintf(file,"0 0 moveto 0 %f lineto 1 %f lineto 1 0 lineto closepath\n",1.0/ratio,1.0/ratio); fprintf(file,"stroke\n"); } sexpr *sketch2postscript(sexpr *sp1, sexpr *sp2) { sexpr *ls, *ln, *port; int rows, cols, count; char text[STRLEN]; float x0, y0, x1, y1, ratio, scale; FILE *file; if (!sketch(sp1)) { printf("sketch->postscript: First argument must be sketch.\n"); longjmp(esc,1); } if ((!string(sp2))) { printf("sketch->postscript: Illegal filename.\n"); longjmp(esc,1); } ls = sp1->u.sketch->ls; rows = sp1->u.sketch->rows; cols = sp1->u.sketch->cols; count = sp1->u.sketch->count; port = open_output_file(sp2); file = port->u.file; ratio = max((float) rows/cols,(float) cols/rows); sprintf(text,"#",rows,cols,count); generate_postscript_header(file, (float) rows/cols,text); fprintf(file,"0.00115741 setlinewidth\n"); scale = min(rows*ratio,cols*ratio); while(!null(ls)) { ln = car(ls); x0 = ln->u.line->x[0]/scale; y0 = ln->u.line->y[0]/scale; x1 = ln->u.line->x[1]/scale; y1 = ln->u.line->y[1]/scale; fprintf(file,"%f %f %f %f displine\n",y0,x0,y1,x1); ls = cdr(ls); } fprintf(file,"showpage\n"); close_output_port(port); return undefined; } char *convert_escape_sequences(char *text) { int i = 1, j = 0; char *converted_text = malloc(STRLEN*sizeof(char)); while (text[i] != '\0') { if (text[i-1] == '\\') { switch (text[i]) { case 'n': converted_text[j] = '\n'; i++; break; case 't': converted_text[j] = '\t'; i++; break; case 'v': converted_text[j] = '\v'; i++; break; case 'b': converted_text[j] = '\b'; i++; break; case 'r': converted_text[j] = '\r'; i++; break; case 'f': converted_text[j] = '\f'; i++; break; case 'a': converted_text[j] = '\a'; i++; break; case '\\': converted_text[j] = '\\'; i++; break; case '\?': converted_text[j] = '\?'; i++; break; case '\'': converted_text[j] = '\''; i++; break; default: converted_text[j] = text[i-1]; break; } } else converted_text[j] = text[i-1]; i++; j++; } converted_text[j++] = text[i-1]; converted_text[j] = '\0'; return converted_text; } sexpr *scheme_fprintf(sexpr *sp1, sexpr *sp2, sexpr *sp3) { char *format; char prefix[STRLEN]; sexpr *arg; int i, j, k; if (!output_port(sp1) || (!string(sp2)) || !vector(sp3)) { printf("fprintf: Illegal argument.\n"); longjmp(esc,1); } format = convert_escape_sequences(sp2->u.text); i = k = 0; while (format[i] != '\0') { j = 0; while (format[i] != '%' && format[i] != '\0') prefix[j++] = format[i++]; while (format[i] != 'e' && format[i] != 'f' && format[i] != 'g' && format[i] != 'c' && format[i] != 's' && format[i] != '\0') prefix[j++] = format[i++]; if (format[i] != '\0') prefix[j++] = format[i++]; prefix[j] = '\0'; if (k < sp3->u.vector->length) { arg = sp3->u.vector->vector[k++]; switch (arg->type) { case NUMBER: if ('e' <= prefix[j-1] && prefix[j-1] <= 'g') { fprintf(sp1->u.file,prefix,arg->u.x); } else { printf("fprintf: %c%c format used with number.\n",'%',prefix[j-1]); longjmp(esc,1); } break; case CHARACTER: if (prefix[j-1] == 'c') { fprintf(sp1->u.file,prefix,arg->u.c); } else { printf("fprintf: %c%c format used with character.\n",'%',prefix[j-1]); longjmp(esc,1); } break; case STRING: if (prefix[j-1] == 's') { fprintf(sp1->u.file,prefix,arg->u.text); } else { printf("fprintf: %c%c format used with string.\n",'%',prefix[j-1]); longjmp(esc,1); } break; default: printf("fprintf: No format for expression of this type.\n"); longjmp(esc,1); } } else fprintf(sp1->u.file,prefix); } if (k < sp3->u.vector->length) { printf("fprintf: Too many arguments for format string.\n"); longjmp(esc,1); } free(format); return undefined; } sexpr *scheme_fscanf(sexpr *sp1, sexpr *sp2) { float x; char *format, c; char prefix[STRLEN], text[STRLEN]; sexpr *result = nil; int i, j, k, status; if (!input_port(sp1) || (!string(sp2))) { printf("fscanf: Illegal argument.\n"); longjmp(esc,1); } format = convert_escape_sequences(sp2->u.text); i = k = 0; while (format[i] != '\0') { j = 0; while (format[i] != '%' && format[i] != '\0') prefix[j++] = format[i++]; while (format[i] != 'e' && format[i] != 'f' && format[i] != 'g' && format[i] != 'c' && format[i] != 's' && format[i] != '\0') prefix[j++] = format[i++]; if (format[i] == '\0') { free(format); return result; } prefix[j++] = format[i++]; prefix[j] = '\0'; switch (prefix[j-1]) { case 'e': case 'f': case 'g': status = fscanf(sp1->u.file,prefix,&x); if (status) { if (status != EOF) result = append(result,cons(num2exp((double) x),nil)); else { free(format); return eof_object; } } break; case 'c': status = fscanf(sp1->u.file,prefix,&c); if (status) { if (status != EOF) result = append(result,cons(char2exp(c),nil)); else { free(format); return eof_object; } } break; case 's': status = fscanf(sp1->u.file,prefix,text); if (status) { if (status != EOF) result = append(result,cons(str2exp(text),nil)); else { free(format); return eof_object; } } break; default: printf("fscanf: Unrecognized format: %c",prefix[j-1]); longjmp(esc,1); } } getc(sp1->u.file); free(format); return result; } sexpr *image2complex(sexpr *sp0) { int i, rows, cols, n; sexpr *sp1; if (!image(sp0)) { print_value_escape("image->complex_image: Argument is not an image: ",sp0); } rows = sp0->u.image->rows; cols = sp0->u.image->cols; n = rows*cols; sp1 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp1->u.complex_image->data[i] = Complex(sp0->u.image->data[i],0.0); } return sp1; } sexpr *register_images(sexpr *im1, sexpr *im2, sexpr *dy, sexpr *dx, sexpr *y0, sexpr *x0) { int i1, j1, i2, j2, index1, index2, count; int rows1, cols1, rows2, cols2, r0, c0, dr, dc; int u, v, max_u, max_v; int p1[256], p2[256], p12[256][256]; float p, h1, h2, h12, i12, max_i12; float *data1, *data2; float overlap; if (!image(im1) || !image(im2)) { printf("register-images: First two arguments must be images.\n"); longjmp(esc,1); } if (!number(y0) || !number(x0) || !number(dy) || !number(dx)) { printf("register-images: Search widths and offsets must be numbers.\n"); longjmp(esc,1); } im1 = image_normalize(im1); im2 = image_normalize(im2); data1 = im1->u.image->data; data2 = im2->u.image->data; rows1 = im1->u.image->rows; cols1 = im1->u.image->cols; rows2 = im2->u.image->rows; cols2 = im2->u.image->cols; for (i1 = 0; i1 < rows1; i1++) { index1 = i1*cols1; for (j1 = 0; j1 < cols1; j1++) { data1[index1 + j1] *= 255; } } for (i2 = 0; i2 < rows2; i2++) { index2 = i2*cols2; for (j2 = 0; j2 < cols2; j2++) { data2[index2 + j2] *= 255; } } r0 = (int) y0->u.x; c0 = (int) x0->u.x; dr = (int) dy->u.x; dc = (int) dx->u.x; if ((dr < 0) || (dc < 0)) { printf("register-images: Search widths must be positive.\n"); longjmp(esc,1); } max_u = max_v = 0; max_i12 = 0.0; for (u = r0 - dr; u <= r0 + dr; u++) { for (v = c0 - dc; v <= c0 + dc; v++) { for (i1 = 0; i1 < 256; i1++) { p1[i1] = p2[i1] = 0; for (j1 = 0; j1 < 256; j1++) { p12[i1][j1] = 0; } } count = 0; for (i1 = 0; i1 < rows1; i1++) { for (j1 = 0; j1 < cols1; j1++) { i2 = i1 + u; j2 = j1 + v; if ((i2 > 0) && (i2 < rows2) && (j2 > 0) && (j2 < cols2)) { index1 = i1*cols1 + j1; index2 = i2*cols2 + j2; p1[(int) data1[index1]]++; p2[(int) data2[index2]]++; p12[(int) data1[index1]][(int) data2[index2]]++; count++; } } } overlap = (float) count; h1 = h2 = h12 = 0.0; for (i1 = 0; i1 < 256; i1++) { h1 += (((p = p1[i1]/overlap) > 0) ? -p*log(p) : 0.0); h2 += (((p = p2[i1]/overlap) > 0) ? -p*log(p) : 0.0); for (i2 = 0; i2 < 256; i2++) { h12 += (((p = p12[i1][i2]/overlap) > 0) ? -p*log(p) : 0.0); } } i12 = h1 + h2 - h12; if (i12 > max_i12) { max_i12 = i12; max_u = u; max_v = v; } } } return cons(num2exp((float) max_u),cons(num2exp((float) max_v),nil)); } int power_of_two(int i) { int j; for (j = 1; j < 2*i; j *= 2) if (j == i) return 1; return 0; } sexpr *fft(sexpr *sp0, sexpr *direction) { int i; int N; float *data; float rootN; sexpr *sp1, *component; if (!vector(sp0) || !number(direction) || fabs(direction->u.x) != 1) { printf("fft: Illegal argument.\n"); longjmp(esc,1); } N=sp0->u.vector->length; if (!power_of_two(N)) { printf("fft: Length of vector must be integral power of two.\n"); longjmp(esc,1); } rootN = sqrt(N); data = malloc((2*N+1)*sizeof(float)); sp1 = make_vector(N); for (i = 0; i < N; i++) { component = sp0->u.vector->vector[i]; if (komplex(component)) { data[2*i+1] = component->u.z.r; data[2*i+2] = component->u.z.i; } else if (number(component)) { data[2*i+1] = component->u.x; data[2*i+2] = 0.0; } else print_value_escape("fft: Vector component not a number: ",component); } four1(data,N,(int) direction->u.x); for (i = 0; i < N; i++) { if (direction->u.x == 1) sp1->u.vector->vector[i] = complex2exp(Complex(data[2*i+1]/rootN,data[2*i+2]/rootN)); else sp1->u.vector->vector[i] = complex2exp(Complex(data[2*i+1]/rootN,data[2*i+2]/rootN)); } free(data); return sp1; } sexpr *image_fft(sexpr *sp0, sexpr* direction) { int i, j; int rows, cols, index; float *data_row, *data_col; float rootNM; sexpr *sp1; if ((!image(sp0) && !complex_image(sp0)) || !number(direction) || fabs(direction->u.x) != 1) { printf("fft-image: Illegal argument.\n"); longjmp(esc,1); } if (image(sp0)) sp0 = image2complex(sp0); cols = sp0->u.complex_image->cols; rows = sp0->u.complex_image->rows; if (!power_of_two(cols) || !power_of_two(rows)) { printf("fft-image: Image dimensions must be integral power of two.\n"); longjmp(esc,1); } rootNM = sqrt(rows*cols); data_row = malloc((2*cols+1)*sizeof(float)); data_col = malloc((2*rows+1)*sizeof(float)); sp1 = make_complex_image(rows,cols); for (i = 0; i < rows; i++) { for (j = 0; j < cols; j++) { index = cols*i + j; data_row[2*j+1] = sp0->u.complex_image->data[index].r; data_row[2*j+2] = sp0->u.complex_image->data[index].i; } four1(data_row,cols,(int) floor((float) direction->u.x)); for (j = 0; j < cols; j++) { index = cols*i + j; sp1->u.complex_image->data[index].r = data_row[2*j+1]; sp1->u.complex_image->data[index].i = data_row[2*j+2]; } } for (j = 0; j < cols; j++) { for (i = 0; i < rows; i++) { index = cols*i + j; data_col[2*i+1] = sp1->u.complex_image->data[index].r; data_col[2*i+2] = sp1->u.complex_image->data[index].i; } four1(data_col,rows,(int) floor((float) direction->u.x)); for (i = 0; i < rows; i++) { index = cols*i + j; sp1->u.complex_image->data[index].r = data_col[2*i+1]/rootNM; sp1->u.complex_image->data[index].i = data_col[2*i+2]/rootNM; } } free(data_row); free(data_col); return sp1; } sexpr *image_sum(sexpr *sp0) { int i, n; switch (sp0->type) { case IMAGE: { float *data, sum = 0; data = sp0->u.image->data; n = sp0->u.image->rows*sp0->u.image->cols; for (i = 0; i < n; i++) { sum += data[i]; } return num2exp(sum); } case COMPLEX_IMAGE: { fcomplex *cdata, csum = zero; cdata = sp0->u.complex_image->data; n = sp0->u.complex_image->rows*sp0->u.complex_image->cols; for (i = 0; i < n; i++) { csum = Cadd(csum,cdata[i]); } return complex2exp(csum); } default: print_value_escape("image-sum: Argument is not an image: ",sp0); } } sexpr *image_min(sexpr *sp0) { int i, n; float x, *data, minimum = FLT_MAX; if (!image(sp0)) { print_value_escape("image-min: Argument is not an image: ",sp0); } data = sp0->u.image->data; n = sp0->u.image->rows*sp0->u.image->cols; for (i = 0; i < n; i++) { x = data[i]; if (x < minimum) minimum = x; } return num2exp(minimum); } sexpr *image_max(sexpr *sp0) { int i, n; float x, *data, maximum = -FLT_MAX; if (!image(sp0)) { print_value_escape("image-max: Argument is not an image: ",sp0); } data = sp0->u.image->data; n = sp0->u.image->rows*sp0->u.image->cols; for (i = 0; i < n; i++) { x = data[i]; if (x > maximum) maximum = x; } return num2exp(maximum); } sexpr *image_normalize(sexpr *sp1) { int i, rows, cols, n; float value, scale; sexpr *sp2; float *data1, *data2; fcomplex *cdata1, *cdata2; float maximum = -FLT_MAX; float minimum = FLT_MAX; switch (sp1->type) { case IMAGE: rows = sp1->u.image->rows; cols = sp1->u.image->cols; sp2 = make_image(rows,cols); data1 = sp1->u.image->data; data2 = sp2->u.image->data; n = rows*cols; for (i=0; i < n; i++) { value = data1[i]; if (value > maximum) maximum = value; if (value < minimum) minimum = value; } scale = maximum-minimum; if (scale < 1.0) scale = 1.0; for (i=0; i < n; i++) data2[i] = (data1[i] - minimum)/scale; return sp2; case COMPLEX_IMAGE: rows = sp1->u.complex_image->rows; cols = sp1->u.complex_image->cols; sp2 = make_complex_image(rows,cols); cdata1 = sp1->u.complex_image->data; cdata2 = sp2->u.complex_image->data; n = rows*cols; for (i=0; i < n; i++) { value = Cabs(cdata1[i]); if (value > maximum) maximum = value; } if (maximum < 1.0) maximum = 1.0; scale = 1.0/maximum; for (i=0; i < n; i++) cdata2[i]=RCmul(scale,cdata1[i]); return sp2; default: printf("image-normalize: Argument must be image.\n"); longjmp(esc,1); } } sexpr *color_image_normalize(sexpr *sp1) { int i, rows, cols, n; float data, maximum, minimum, scale; sexpr *sp2; if (!color_image(sp1)) { printf("color-image-normalize: Argument must be color-image.\n"); longjmp(esc,1); } rows = sp1->u.image->rows; cols = sp1->u.image->cols; sp2 = make_color_image(rows,cols); maximum = -FLT_MAX; minimum = FLT_MAX; n = rows*cols; for (i=0; i < n*3; i++) { data = sp1->u.image->data[i]; if (data > maximum) maximum = data; if (data < minimum) minimum = data; } scale = maximum-minimum; if (scale < 1.0) scale = 1.0; for (i=0; i < n*3; i++) sp2->u.image->data[i]=(sp1->u.image->data[i] - minimum)/scale; return sp2; } sexpr *shrink(sexpr *sp1, sexpr *sp2) { int i, rows, cols, n; sexpr *sp3; float magnitude, threshold; fcomplex z; if (!number(sp2)) { printf("shrink: Second argument must be a number.\n"); longjmp(esc,1); } threshold = (float) sp2->u.x; switch (sp1->type) { case IMAGE: rows = sp1->u.image->rows; cols = sp1->u.image->cols; sp3 = make_image(rows,cols); n = rows*cols; for (i = 0; i < n; i++) { magnitude = sp1->u.image->data[i]; if (fabs(magnitude) < threshold) sp3->u.image->data[i] = 0; else sp3->u.image->data[i] = (magnitude > 0) ? magnitude - threshold : magnitude + threshold; } return sp3; case COMPLEX_IMAGE: rows = sp1->u.complex_image->rows; cols = sp1->u.complex_image->cols; sp3 = make_complex_image(rows,cols); n = rows*cols; for (i = 0; i < n; i++) { z = sp1->u.complex_image->data[i]; magnitude = Cmag(z); if (magnitude > 0.0) { if (magnitude > threshold) sp3->u.complex_image->data[i] = RCmul((magnitude - threshold)/magnitude,z); else sp3->u.complex_image->data[i] = zero; } } return sp3; default: printf("shrink: First argument must be an image.\n"); longjmp(esc,1); } } sexpr *image_rows(sexpr *sp) { switch (sp->type) { case IMAGE: case COLOR_IMAGE: return num2exp(sp->u.image->rows); case COMPLEX_IMAGE: return num2exp(sp->u.complex_image->rows); default: print_value_escape("image-rows: Argument is not an image: ",sp); } } sexpr *image_cols(sexpr *sp) { switch (sp->type) { case IMAGE: case COLOR_IMAGE: return num2exp(sp->u.image->cols); case COMPLEX_IMAGE: return num2exp(sp->u.complex_image->cols); default: print_value_escape("image-cols: Argument is not an image: ",sp); } } sexpr *real_image_complex_helper(fcomplex (func)(float,float), sexpr *sp1, sexpr *sp2) { int rows = sp2->u.image->rows; int cols = sp2->u.image->cols; int i, n = rows*cols; float x = (float) sp1->u.x; float *data2 = sp2->u.image->data; sexpr *sp3 = make_complex_image(rows,cols); fcomplex *cdata3 = sp3->u.complex_image->data; for (i = 0; i < n; i++) cdata3[i] = func((float) x, data2[i]); return sp3; } sexpr *image_real_complex_helper(fcomplex (func)(float,float), sexpr *sp1, sexpr *sp2) { int rows = sp1->u.image->rows; int cols = sp1->u.image->cols; int i, n = rows*cols; float *data1 = sp1->u.image->data; float x = (float) sp2->u.x; sexpr *sp3 = make_complex_image(rows,cols); fcomplex *cdata3 = sp3->u.complex_image->data; for (i = 0; i < n; i++) cdata3[i] = func(data1[i],x); return sp3; } sexpr *image_image_complex_helper(fcomplex (func)(float,float), sexpr *sp1, sexpr *sp2) { int rows = sp1->u.image->rows; int cols = sp1->u.image->cols; int i, n = rows*cols; float *data1 = sp1->u.image->data; float *data2 = sp2->u.image->data; sexpr *sp3 = make_complex_image(rows,cols); fcomplex *cdata3 = sp3->u.complex_image->data; check_image_sizes("complex-helper",rows,cols,sp2->u.image); for (i = 0; i < n; i++) cdata3[i] = func(data1[i],data2[i]); return sp3; } sexpr *complex_helper(fcomplex (func)(float,float), sexpr *sp1, sexpr *sp2) { switch (sp1->type) { case NUMBER: switch (sp2->type) { case NUMBER: return complex2exp(func((float) sp1->u.x, (float) sp2->u.x)); case IMAGE: return real_image_complex_helper(func,sp1,sp2); default: printf("complex-helper: Argument is not a number or image.\n"); longjmp(esc,1); } case IMAGE: switch (sp2->type) { case NUMBER: return real_image_complex_helper(func,sp1,sp2); case IMAGE: return image_image_complex_helper(func,sp1,sp2); default: printf("complex-helper: Incompatible argument types.\n"); longjmp(esc,1); } default: printf("complex-helper: Argument is not a number or image.\n"); longjmp(esc,1); } } fcomplex Polar(float rho, float theta) { return Complex(rho*cos(theta),rho*sin(theta)); } sexpr *polar2complex(sexpr *sp1, sexpr *sp2) { return complex_helper(Polar,sp1,sp2); } sexpr *make_complex(sexpr *sp1, sexpr *sp2) { return complex_helper(Complex,sp1,sp2); } sexpr *complex_image_real(sexpr *sp1) { int i, n, rows, cols; sexpr *sp2; if (!complex_image(sp1)) { printf("complex-image-real-part: Argument is not a complex image.\n"); longjmp(esc,1); } rows = sp1->u.complex_image->rows; cols = sp1->u.complex_image->cols; sp2 = make_image(rows,cols); n = rows*cols; for (i=0; i < rows*cols; i++) sp2->u.image->data[i] = (float) sp1->u.complex_image->data[i].r; return sp2; } sexpr *complex_image_imag(sexpr *sp1) { int i, n, rows, cols; sexpr *sp2; if (!complex_image(sp1)) { printf("complex-image-imag-part: Argument is not a complex image.\n"); longjmp(esc,1); } rows = sp1->u.complex_image->rows; cols = sp1->u.complex_image->cols; sp2 = make_image(rows,cols); n = rows*cols; for (i=0; i < n; i++) sp2->u.image->data[i] = (float) sp1->u.complex_image->data[i].i; return sp2; } sexpr *complex_image_angle(sexpr *sp1) { int i, rows, cols; sexpr *sp2; if (!complex_image(sp1)) { printf("complex-image-angle: Argument is not a complex image.\n"); longjmp(esc,1); } rows = sp1->u.complex_image->rows; cols = sp1->u.complex_image->cols; sp2 = make_image(rows,cols); for (i=0; i < rows*cols; i++) sp2->u.image->data[i] = atan2(sp1->u.complex_image->data[i].i,sp1->u.complex_image->data[i].r); return sp2; } sexpr *complex_image_magnitude(sexpr *sp1) { int i, rows, cols; sexpr *sp2; if (!complex_image(sp1)) { printf("complex-image-angle: Argument is not a complex image.\n"); longjmp(esc,1); } rows = sp1->u.complex_image->rows; cols = sp1->u.complex_image->cols; sp2 = make_image(rows,cols); for (i=0; i < rows*cols; i++) sp2->u.image->data[i] = (float) Cmag(sp1->u.complex_image->data[i]); return sp2; } sexpr *image_cat_helper_helper(sexpr *(*func)(), sexpr *sp1, sexpr *sp2) { return make_complex(func(complex_image_real(sp1),complex_image_real(sp2)), func(complex_image_imag(sp1),complex_image_imag(sp2))); } sexpr *image_cat_helper(sexpr *(*func)(), sexpr *sp1, sexpr *sp2) { switch (sp1->type) { case IMAGE: switch (sp2->type) { case IMAGE: return func(sp1,sp2); case COMPLEX_IMAGE: return image_cat_helper_helper(func,image2complex(sp1),sp2); default: return undefined; } case COMPLEX_IMAGE: switch (sp2->type) { case IMAGE: return image_cat_helper_helper(func,sp1,image2complex(sp2)); case COMPLEX_IMAGE: return image_cat_helper_helper(func,sp1,sp2); default: return undefined; } default: return undefined; } } sexpr *top2bottom(sexpr *sp1, sexpr *sp2) { int i, j, rows1, rows2, rows3, cols, size1, index; sexpr *sp3; cols = sp1->u.image->cols; check_image_sizes("top-to-bottom", -1, cols, sp2->u.image); rows1 = sp1->u.image->rows; rows2 = sp2->u.image->rows; rows3 = rows1 + rows2; sp3 = make_image(rows3,cols); size1 = rows1*cols; for (j = 0; j < cols; j++) { for (i = 0; i < rows1; i++) { index = cols*i + j; sp3->u.image->data[index] = sp1->u.image->data[index]; } for (i = 0; i < rows2; i++) { index = cols*i + j; sp3->u.image->data[index+size1] = sp2->u.image->data[index]; } } return sp3; } sexpr *left2right(sexpr *sp1, sexpr *sp2) { int i, j, cols1, cols2, cols3, rows, index1, index2, index3; sexpr *sp3; rows = sp1->u.image->rows; check_image_sizes("left-to-right", rows, -1, sp2->u.image); cols1 = sp1->u.image->cols; cols2 = sp2->u.image->cols; cols3 = cols1 + cols2; sp3 = make_image(rows,cols3); for (i = 0; i < rows; i++) { for (j = 0; j < cols1; j++) { index1 = cols1*i + j; index3 = cols3*i + j; sp3->u.image->data[index3] = sp1->u.image->data[index1]; } for (j = 0; j < cols2; j++) { index2 = cols2*i + j; index3 = cols3*i + cols1 + j; sp3->u.image->data[index3] = sp2->u.image->data[index2]; } } return sp3; } sexpr *scheme_top2bottom(sexpr *sp1, sexpr *sp2) { sexpr *sp3 = image_cat_helper(top2bottom,sp1,sp2); if (undefined(sp3)) { printf("top-to-bottom: Argument is not an image.\n"); longjmp(esc,1); } return sp3; } sexpr *scheme_left2right(sexpr *sp1, sexpr *sp2) { sexpr *sp3 = image_cat_helper(left2right,sp1,sp2); if (undefined(sp3)) { printf("left-to-right: Argument is not an image.\n"); longjmp(esc,1); } return sp3; } sexpr *matrix_product(sexpr *sp1, sexpr *sp2) { int i1, j1, j2, rows1, cols1, cols2; float sum, *data1, *data2, *data3; fcomplex csum, *cdata1, *cdata2, *cdata3; sexpr *sp3; switch (sp1->type) { case IMAGE: rows1 = sp1->u.image->rows; cols1 = sp1->u.image->cols; data1 = sp1->u.image->data; switch (sp2->type) { case IMAGE: check_image_sizes("matrix-product", cols1, -1, sp2->u.image); cols2 = sp2->u.image->cols; data2 = sp2->u.image->data; sp3 = make_image(rows1,cols2); data3 = sp3->u.image->data; for (i1 = 0; i1 < rows1; i1++) { for (j2 = 0; j2 < cols2; j2++) { sum = 0.0; for (j1 = 0; j1 < cols1; j1++) { sum += data1[i1*cols1 + j1]*data2[j1*cols2 + j2]; } data3[i1*cols2 + j2] = sum; } } return sp3; case COMPLEX_IMAGE: check_image_sizes("matrix-product", cols1, -1, (image *) sp2->u.complex_image); cols2 = sp2->u.complex_image->cols; cdata2 = sp2->u.complex_image->data; sp3 = make_complex_image(rows1,cols2); cdata3 = sp3->u.complex_image->data; for (i1 = 0; i1 < rows1; i1++) { for (j2 = 0; j2 < cols2; j2++) { csum = zero; for (j1 = 0; j1 < cols1; j1++) { csum = Cadd(csum,RCmul(data1[i1*cols1 + j1],cdata2[j1*cols2 + j2])); } cdata3[i1*cols2 + j2] = csum; } } return sp3; default: print_value_escape("matrix-product: Argument is not an image: ",sp2); } case COMPLEX_IMAGE: rows1 = sp1->u.complex_image->rows; cols1 = sp1->u.complex_image->cols; cdata1 = sp1->u.complex_image->data; switch (sp2->type) { case IMAGE: check_image_sizes("matrix-product", cols1, -1, sp2->u.image); cols2 = sp2->u.image->cols; data2 = sp2->u.image->data; sp3 = make_complex_image(rows1,cols2); cdata3 = sp3->u.complex_image->data; for (i1 = 0; i1 < rows1; i1++) { for (j2 = 0; j2 < cols2; j2++) { csum = zero; for (j1 = 0; j1 < cols1; j1++) { csum = Cadd(csum,RCmul(data2[j1*cols2 + j2],cdata1[i1*cols1 + j1])); } cdata3[i1*cols2 + j2] = csum; } } return sp3; case COMPLEX_IMAGE: check_image_sizes("matrix-product", cols1, -1, (image *) sp2->u.complex_image); cols2 = sp2->u.complex_image->cols; cdata2 = sp2->u.complex_image->data; sp3 = make_complex_image(rows1,cols2); cdata3 = sp3->u.complex_image->data; for (i1 = 0; i1 < rows1; i1++) { for (j2 = 0; j2 < cols2; j2++) { csum = zero; for (j1 = 0; j1 < cols1; j1++) { csum = Cadd(csum,Cmul(cdata1[i1*cols1 + j1],cdata2[j1*cols2 + j2])); } cdata3[i1*cols2 + j2] = csum; } } return sp3; default: print_value_escape("matrix-product: Argument is not an image: ",sp2); } default: print_value_escape("matrix-product: Argument is not an image: ",sp1); } } sexpr *times(sexpr *sp1, sexpr *sp2) { int i, n, rows, cols; float x; fcomplex z; sexpr *sp3; switch (sp1->type) { case NUMBER: x = sp1->u.x; switch (sp2->type) { case NUMBER: return num2exp(x*sp2->u.x); case COMPLEX: return complex2exp(RCmul(x,sp2->u.z)); case IMAGE: rows = sp2->u.image->rows; cols = sp2->u.image->cols; n = rows*cols; sp3 = make_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.image->data[i] = sp2->u.image->data[i]*x; } return sp3; case COMPLEX_IMAGE: rows = sp2->u.complex_image->rows; cols = sp2->u.complex_image->cols; n = rows*cols; sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = RCmul(x,sp2->u.complex_image->data[i]); } return sp3; default: print_value_escape("*: Argument is not a number or image: ",sp2); break; } case COMPLEX: z = sp1->u.z; switch (sp2->type) { case NUMBER: return complex2exp(RCmul(sp2->u.x,z)); case COMPLEX: return complex2exp(Cmul(z,sp2->u.z)); case IMAGE: rows = sp2->u.image->rows; cols = sp2->u.image->cols; n = rows*cols; sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = RCmul(sp2->u.image->data[i],z); } return sp3; case COMPLEX_IMAGE: rows = sp2->u.complex_image->rows; cols = sp2->u.complex_image->cols; n = rows*cols; sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Cmul(z,sp2->u.complex_image->data[i]); } return sp3; default: print_value_escape("*: Argument is not a number or image: ",sp2); break; } case IMAGE: rows = sp1->u.image->rows; cols = sp1->u.image->cols; n = rows*cols; switch (sp2->type) { case NUMBER: sp3 = make_image(rows,cols); x = sp2->u.x; for (i = 0; i < n; i++) { sp3->u.image->data[i] = sp1->u.image->data[i]*x; } return sp3; case COMPLEX: sp3 = make_complex_image(rows,cols); z = sp2->u.z; for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = RCmul(sp1->u.image->data[i],z); } return sp3; case IMAGE: check_image_sizes("*", rows, cols, sp2->u.image); sp3 = make_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.image->data[i] = sp1->u.image->data[i]*sp2->u.image->data[i]; } return sp3; case COMPLEX_IMAGE: check_image_sizes("*", rows, cols, (image *) sp2->u.complex_image); sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = RCmul(sp1->u.image->data[i],sp2->u.complex_image->data[i]); } return sp3; default: print_value_escape("*: Argument is not a number or image: ",sp2); break; } case COMPLEX_IMAGE: rows = sp1->u.complex_image->rows; cols = sp1->u.complex_image->cols; n = rows*cols; switch (sp2->type) { case NUMBER: sp3 = make_complex_image(rows,cols); x = sp2->u.x; for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = RCmul(x,sp1->u.complex_image->data[i]); } return sp3; case COMPLEX: sp3 = make_complex_image(rows,cols); z = sp2->u.z; for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Cmul(sp1->u.complex_image->data[i],z); } return sp3; case IMAGE: check_image_sizes("*", rows, cols, sp2->u.image); sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = RCmul(sp2->u.image->data[i],sp1->u.complex_image->data[i]); } return sp3; case COMPLEX_IMAGE: check_image_sizes("*", rows, cols, (image *) sp2->u.complex_image); sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Cmul(sp1->u.complex_image->data[i],sp2->u.complex_image->data[i]); } return sp3; default: print_value_escape("*: Argument is not a number or image: ",sp2); break; } default: print_value_escape("*: Argument is not a number or image: ",sp1); break; } } sexpr *plus(sexpr *sp1, sexpr *sp2) { int i, n, rows, cols; float x; fcomplex z; sexpr *sp3; switch (sp1->type) { case NUMBER: x = sp1->u.x; switch (sp2->type) { case NUMBER: return num2exp(x+sp2->u.x); case COMPLEX: return complex2exp(Cadd(Complex((float) x, 0.0),sp2->u.z)); case IMAGE: rows = sp2->u.image->rows; cols = sp2->u.image->cols; n = rows*cols; sp3 = make_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.image->data[i] = x+sp2->u.image->data[i]; } return sp3; case COMPLEX_IMAGE: rows = sp2->u.complex_image->rows; cols = sp2->u.complex_image->cols; n = rows*cols; sp3 = make_complex_image(rows,cols); z = Complex((float) x, 0.0); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Cadd(z,sp2->u.complex_image->data[i]); } return sp3; default: print_value_escape("+: Argument is not a number or image: ",sp2); break; } case COMPLEX: z = sp1->u.z; switch (sp2->type) { case NUMBER: return complex2exp(Cadd(z,Complex((float) sp2->u.x, 0.0))); case COMPLEX: return complex2exp(Cadd(z,sp2->u.z)); case IMAGE: rows = sp2->u.image->rows; cols = sp2->u.image->cols; n = rows*cols; sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Cadd(z,Complex(sp2->u.image->data[i], 0.0)); } return sp3; case COMPLEX_IMAGE: rows = sp2->u.complex_image->rows; cols = sp2->u.complex_image->cols; n = rows*cols; sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Cadd(z,sp2->u.complex_image->data[i]); } return sp3; default: print_value_escape("+: Argument is not a number or image: ",sp2); break; } case IMAGE: rows = sp1->u.image->rows; cols = sp1->u.image->cols; n = rows*cols; switch (sp2->type) { case NUMBER: sp3 = make_image(rows,cols); x = sp2->u.x; for (i = 0; i < n; i++) { sp3->u.image->data[i] = sp1->u.image->data[i]+x; } return sp3; case COMPLEX: sp3 = make_complex_image(rows,cols); z = sp2->u.z; for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Cadd(Complex(sp1->u.image->data[i], 0.0),z); } return sp3; case IMAGE: check_image_sizes("+", rows, cols, sp2->u.image); sp3 = make_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.image->data[i] = sp1->u.image->data[i]+sp2->u.image->data[i]; } return sp3; case COMPLEX_IMAGE: check_image_sizes("+", rows, cols, (image *) sp2->u.complex_image); sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Cadd(Complex(sp1->u.image->data[i],0.0), sp2->u.complex_image->data[i]); } return sp3; default: print_value_escape("+: Argument is not a number or image: ",sp2); break; } case COMPLEX_IMAGE: rows = sp1->u.complex_image->rows; cols = sp1->u.complex_image->cols; n = rows*cols; switch (sp2->type) { case NUMBER: sp3 = make_complex_image(rows,cols); z = Complex((float) sp2->u.x, 0.0); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Cadd(sp1->u.complex_image->data[i],z); } return sp3; case COMPLEX: sp3 = make_complex_image(rows,cols); z = sp2->u.z; for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Cadd(sp1->u.complex_image->data[i],z); } return sp3; case IMAGE: check_image_sizes("+", rows, cols, sp2->u.image); sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Cadd(sp1->u.complex_image->data[i],Complex(sp2->u.image->data[i], 0.0)); } return sp3; case COMPLEX_IMAGE: check_image_sizes("+", rows, cols, (image *) sp2->u.complex_image); sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Cadd(sp1->u.complex_image->data[i],sp2->u.complex_image->data[i]); } return sp3; default: print_value_escape("+: Argument is not a number or image: ",sp2); break; } default: print_value_escape("+: Argument is not a number or image: ",sp1); break; } } sexpr *minus(sexpr *sp1, sexpr *sp2) { int i, n, rows, cols; float x; fcomplex z; sexpr *sp3; switch (sp1->type) { case NUMBER: x = sp1->u.x; switch (sp2->type) { case NUMBER: return num2exp(x-sp2->u.x); case COMPLEX: return complex2exp(Csub(Complex((float) x, 0.0),sp2->u.z)); case IMAGE: rows = sp2->u.image->rows; cols = sp2->u.image->cols; n = rows*cols; sp3 = make_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.image->data[i] = x-sp2->u.image->data[i]; } return sp3; case COMPLEX_IMAGE: rows = sp2->u.complex_image->rows; cols = sp2->u.complex_image->cols; n = rows*cols; sp3 = make_complex_image(rows,cols); z = Complex((float) x, 0.0); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Csub(z,sp2->u.complex_image->data[i]); } return sp3; default: print_value_escape("-: Argument is not a number or image: ",sp2); break; } case COMPLEX: z = sp1->u.z; switch (sp2->type) { case NUMBER: return complex2exp(Csub(z,Complex((float) sp2->u.x, 0.0))); case COMPLEX: return complex2exp(Csub(z,sp2->u.z)); case IMAGE: rows = sp2->u.image->rows; cols = sp2->u.image->cols; n = rows*cols; sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Csub(z,Complex(sp2->u.image->data[i], 0.0)); } return sp3; case COMPLEX_IMAGE: rows = sp2->u.complex_image->rows; cols = sp2->u.complex_image->cols; n = rows*cols; sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Csub(z,sp2->u.complex_image->data[i]); } return sp3; default: print_value_escape("-: Argument is not a number or image: ",sp2); break; } case IMAGE: rows = sp1->u.image->rows; cols = sp1->u.image->cols; n = rows*cols; switch (sp2->type) { case NUMBER: sp3 = make_image(rows,cols); x = sp2->u.x; for (i = 0; i < n; i++) { sp3->u.image->data[i] = sp1->u.image->data[i]-x; } return sp3; case COMPLEX: sp3 = make_complex_image(rows,cols); z = sp2->u.z; for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Csub(Complex(sp1->u.image->data[i], 0.0),z); } return sp3; case IMAGE: check_image_sizes("-", rows, cols, sp2->u.image); sp3 = make_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.image->data[i] = sp1->u.image->data[i]-sp2->u.image->data[i]; } return sp3; case COMPLEX_IMAGE: check_image_sizes("-", rows, cols, (image *) sp2->u.complex_image); sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Csub(Complex(sp1->u.image->data[i],0.0), sp2->u.complex_image->data[i]); } return sp3; default: print_value_escape("-: Argument is not a number or image: ",sp2); break; } case COMPLEX_IMAGE: rows = sp1->u.complex_image->rows; cols = sp1->u.complex_image->cols; n = rows*cols; switch (sp2->type) { case NUMBER: sp3 = make_complex_image(rows,cols); z = Complex((float) sp2->u.x, 0.0); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Csub(sp1->u.complex_image->data[i],z); } return sp3; case COMPLEX: sp3 = make_complex_image(rows,cols); z = sp2->u.z; for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Csub(sp1->u.complex_image->data[i],z); } return sp3; case IMAGE: check_image_sizes("-", rows, cols, sp2->u.image); sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Csub(sp1->u.complex_image->data[i],Complex(sp2->u.image->data[i], 0.0)); } return sp3; case COMPLEX_IMAGE: check_image_sizes("-", rows, cols, (image *) sp2->u.complex_image); sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Csub(sp1->u.complex_image->data[i],sp2->u.complex_image->data[i]); } return sp3; default: print_value_escape("-: Argument is not a number or image: ",sp2); break; } default: print_value_escape("-: Argument is not a number or image: ",sp1); break; } } sexpr *divide(sexpr *sp1, sexpr *sp2) { int i, n, rows, cols; float x1, x2; fcomplex z1, z2; sexpr *sp3; switch (sp1->type) { case NUMBER: x1 = sp1->u.x; switch (sp2->type) { case NUMBER: x2 = sp2->u.x; if (x2 != 0.0) return num2exp(x1/x2); printf("/: Attempt to divide by zero.\n"); longjmp(esc,1); case COMPLEX: z1 = sp2->u.z; if (Cmag(z1) != 0.0) return complex2exp(Cdiv(Complex((float) x1, 0.0), z1)); printf("/: Attempt to divide by zero.\n"); longjmp(esc,1); case IMAGE: rows = sp2->u.image->rows; cols = sp2->u.image->cols; n = rows*cols; sp3 = make_image(rows,cols); for (i = 0; i < n; i++) { x2 = sp2->u.image->data[i]; if (x2 != 0) { sp3->u.image->data[i] = x1/x2; } else { printf("/: Attempt to divide by zero.\n"); longjmp(esc,1); } } return sp3; case COMPLEX_IMAGE: rows = sp2->u.complex_image->rows; cols = sp2->u.complex_image->cols; n = rows*cols; sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { z1 = sp2->u.complex_image->data[i]; if (Cmag(z1) != 0.0) { sp3->u.complex_image->data[i] = Cdiv(Complex((float) x1, 0.0), z1); } else { printf("/: Attempt to divide by zero.\n"); longjmp(esc,1); } } return sp3; default: print_value_escape("/: Argument is not a number or image: ",sp2); break; } case COMPLEX: z1 = sp1->u.z; switch (sp2->type) { case NUMBER: x2 = sp2->u.x; if (x2 != 0.0) return complex2exp(RCmul(1.0/x2,z1)); printf("/: Attempt to divide by zero.\n"); longjmp(esc,1); case COMPLEX: z2 = sp2->u.z; if (Cmag(z2) != 0.0) return complex2exp(Cdiv(z1,z2)); printf("/: Attempt to divide by zero.\n"); longjmp(esc,1); case IMAGE: rows = sp2->u.image->rows; cols = sp2->u.image->cols; n = rows*cols; sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { x2 = sp2->u.image->data[i]; if (x2 != 0.0) { sp3->u.complex_image->data[i] = RCmul(1.0/x2,z1); } else { printf("/: Attempt to divide by zero.\n"); longjmp(esc,1); } } return sp3; case COMPLEX_IMAGE: rows = sp2->u.complex_image->rows; cols = sp2->u.complex_image->cols; n = rows*cols; sp3 = make_complex_image(rows,cols); for (i = 0; i < n; i++) { z2 = sp2->u.complex_image->data[i]; if (Cmag(z2) != 0.0) { sp3->u.complex_image->data[i] = Cdiv(z1,z2); } else { printf("/: Attempt to divide by zero.\n"); longjmp(esc,1); } } return sp3; default: print_value_escape("/: Argument is not a number or image: ",sp2); break; } case IMAGE: rows = sp1->u.image->rows; cols = sp1->u.image->cols; n = rows*cols; switch (sp2->type) { case NUMBER: sp3 = make_image(rows,cols); x2 = sp2->u.x; if (x2 == 0.0) { printf("/: Attempt to divide by zero.\n"); longjmp(esc,1); } for (i = 0; i < n; i++) { sp3->u.image->data[i] = sp1->u.image->data[i]/x2; } return sp3; case COMPLEX: sp3 = make_complex_image(rows,cols); z2 = sp2->u.z; if (Cmag(z2) == 0.0) { printf("/: Attempt to divide by zero.\n"); longjmp(esc,1); } for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = Cdiv(Complex((float) sp1->u.image->data[i], 0.0), z2); } return sp3; case I