/* 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 IMAGE: check_image_sizes("/", rows, cols, sp2->u.image); sp3 = make_image(rows,cols); for (i = 0; i < n; i++) { x2 = sp2->u.image->data[i]; if (x2 != 0.0) { sp3->u.image->data[i] = sp1->u.image->data[i]/x2; } else { printf("/: Attempt to divide by zero.\n"); longjmp(esc,1); } } 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++) { z2 = sp2->u.complex_image->data[i]; if (Cmag(z2) != 0.0) { sp3->u.complex_image->data[i] = Cdiv(Complex((float) sp1->u.image->data[i], 0.0), 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 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); x2 = sp2->u.x; if (x2 == 0.0) { printf("/: Attempt to divide by zero.\n"); longjmp(esc,1); } x2 = 1.0/x2; for (i = 0; i < n; i++) { sp3->u.complex_image->data[i] = RCmul(x2,sp1->u.complex_image->data[i]); } 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(sp1->u.complex_image->data[i],z2); } return sp3; case IMAGE: check_image_sizes("/", rows, cols, sp2->u.image); 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,sp1->u.complex_image->data[i]); } else { printf("/: Attempt to divide by zero.\n"); longjmp(esc,1); } } 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++) { z2 = sp2->u.complex_image->data[i]; if (Cmag(z2) != 0.0) { sp3->u.complex_image->data[i] = Cdiv(sp1->u.complex_image->data[i],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; } default: print_value_escape("/: Argument is not a number or image: ",sp1); break; } } sexpr *setcar(sexpr *s1, sexpr *s2) { if (!pair(s1)) { print_value_escape("Attempt to set-car! of non-pair: ",s1); } s1->u.pair->car = s2; return undefined; } sexpr *setcdr(sexpr *s1, sexpr *s2) { if (!pair(s1)) { print_value_escape("Attempt to set-cdr! of non-pair: ",s2); } s1->u.pair->cdr = s2; return undefined; } sexpr *make_number_list(int n) { if (n == 0) return nil; return cons(num2exp((double) n),make_number_list(n-1)); } sexpr *make_complex_number_list(int n) { if (n == 0) return nil; return cons(complex2exp(Complex(0.0,(float) n)),make_complex_number_list(n-1)); } sexpr *compile_app(sexpr *sp) { return cps2vector(car(compile(sp,cons(bytecode$halt,nil),global_env_vars))); } sexpr *image2list(sexpr *sp0) { int i, n; float *data; fcomplex *cdata; sexpr *sp1 = nil; switch (sp0->type) { case IMAGE: n = sp0->u.image->rows*sp0->u.image->cols; data = sp0->u.image->data; for (i = 0; i < n; i++) sp1 = cons(num2exp(data[i]),sp1); return sp1; case COMPLEX_IMAGE: n = sp0->u.complex_image->rows*sp0->u.complex_image->cols; cdata = sp0->u.complex_image->data; for (i = 0; i < n; i++) sp1 = cons(complex2exp(cdata[i]),sp1); return sp1; default: printf("image->list: Illegal argument.\n"); longjmp(esc,1); } } sexpr *scheme_image_map(sexpr *sp0, sexpr *sp1) { int i, j, k, n; int rows, cols, index, cflag; vector *vp; float **images; fcomplex **complex_images; float *data2; fcomplex *cdata3; sexpr *code, *result, *args, *sp2, *sp3, *sp4; if (!vector(sp1)) { printf("image-map: Illegal argument.\n"); longjmp(esc,1); } cflag = 0; sp4 = sp1->u.vector->vector[0]; n = sp1->u.vector->length; switch (sp0->type) { case PRIMITIVE: if (sp0->u.primitive->arity != n) { printf("image-map: First argument must be a procedure of %d arguments.\n",n); longjmp(esc,1); } break; case VIRGIN: case CLOSURE: if (sp0->u.closure->n != n && abs(sp0->u.closure->n+1) > n) { printf("image-map: First argument must be a procedure of %d arguments.\n",n); longjmp(esc,1); } break; default: printf("image-map: First argument must be a procedure.\n"); longjmp(esc,1); } images = (float **) malloc(n*sizeof(float *)); complex_images = (fcomplex **) malloc(n*sizeof(fcomplex *)); args = make_complex_number_list(n); vp = list2vector(args)->u.vector; switch (sp4->type) { case IMAGE: rows = sp4->u.image->rows; cols = sp4->u.image->cols; images[0] = sp4->u.image->data; vp->vector[0]->type = NUMBER; break; case COMPLEX_IMAGE: rows = sp4->u.complex_image->rows; cols = sp4->u.complex_image->cols; complex_images[0] = sp4->u.complex_image->data; break; default: printf("image-map: Illegal argument.\n"); free(images); longjmp(esc,1); } for (i = 1; i < n; i++) { sp4 = sp1->u.vector->vector[i]; switch (sp4->type) { case IMAGE: check_image_sizes("image-map", rows, cols, sp4->u.image); images[i] = sp4->u.image->data; vp->vector[i]->type = NUMBER; break; case COMPLEX_IMAGE: check_image_sizes("image-map", rows, cols, (image *) sp4->u.complex_image); complex_images[i] = sp4->u.complex_image->data; break; default: printf("image-map: Illegal argument.\n"); free(images); longjmp(esc,1); } } sp2 = make_image(rows,cols); data2 = sp2->u.image->data; sp3 = make_complex_image(rows,cols); cdata3 = sp3->u.complex_image->data; code = compile_app(cons(sp0,args)); index = 0; for (i = 0; i < rows; i++) { for (j = 0; j < cols; j++) { for (k = 0; k < n; k++) { if (number(vp->vector[k])) { vp->vector[k]->u.x = images[k][index]; } else { vp->vector[k]->u.z = complex_images[k][index]; } } result = virtual_machine(nil,code,global_env_vals); switch (result->type) { case NUMBER: data2[index] = result->u.x; cdata3[index++] = Complex((float) result->u.x, 0.0); break; case COMPLEX: cflag = 1; cdata3[index++] = result->u.z; break; default: printf("image-map: Error at location (%d,%d).\n",i,j); longjmp(esc,1); } } } free(images); return cflag ? sp3 : sp2; } sexpr *image_crop(sexpr *image, sexpr *row0, sexpr *col0, sexpr *rows, sexpr *cols) { int input_cols, i, j, output_rows, output_cols, r0, c0, n, index; sexpr *sp; if (!image(image)) { printf("image-crop: First argument must be an image.\n"); longjmp(esc,1); } if (!number(row0) || !number(col0) || !number(rows) || !number(cols)) { printf("image-crop: Second through fifth arguments must be numbers.\n"); longjmp(esc,1); } input_cols = image->u.image->cols; n = input_cols*image->u.image->rows; r0 = row0->u.x; c0 = col0->u.x; output_rows = (int) floor((float) rows->u.x); output_cols = (int) floor((float) cols->u.x); if ((r0 < 0) || (c0 < 0) || (output_rows <= 0) || (output_cols <= 0)) { printf("image-crop: Number arguments must be positive.\n"); longjmp(esc,1); } sp = make_image(output_rows,output_cols); for (i=0; i < output_rows; i++) { for (j=0; j < output_cols; j++) { index = (i + r0)*input_cols + j + c0; if ((index >= 0) && (index < n)) { sp->u.image->data[i*output_cols + j] = image->u.image->data[index]; } else { printf("image-crop: Incompatible output dimensions.\n"); longjmp(esc,1); } } } return sp; } sexpr *image_pad(sexpr *image, sexpr *rows, sexpr *cols) { int input_rows, input_cols, i, j, output_rows, output_cols, n; sexpr *sp; if (!image(image)) { printf("image-pad: First argument must be an image.\n"); longjmp(esc,1); } if (!number(rows) || !number(cols)) { printf("image-pad: Row and column size must be numbers.\n"); longjmp(esc,1); } input_rows = image->u.image->rows; input_cols = image->u.image->cols; output_rows = (int) floor((float) rows->u.x); output_cols = (int) floor((float) cols->u.x); n = output_rows*output_cols; if ((input_rows < 0) || (input_cols < 0) || (output_rows <= 0) || (output_cols <= 0)) { printf("image-pad: Number arguments must be positive.\n"); longjmp(esc,1); } if (output_rows < input_rows || output_cols < input_cols) { printf("image-pad: Output sizes must exceed input sizes.\n"); longjmp(esc,1); } sp = make_image(output_rows,output_cols); for (i=0; i < n; i++) sp->u.image->data[i] = 0.0; for (i=0; i < input_rows; i++) { for (j=0; j < input_cols; j++) { sp->u.image->data[i*output_cols + j] = image->u.image->data[i*input_cols + j]; } } return sp; } sexpr *scheme_make_image(sexpr *sp1, sexpr *sp2, sexpr *sp3, sexpr *sp4) { int i, j, cflag, not_filter; int rows, cols, fold_rows, fold_cols; int n, index; sexpr *sp5, *sp6, *code, *result; sexpr *value1 = num2exp(0); sexpr *value2 = num2exp(0); cflag = 0; if (!number(sp1) || !number(sp2)) { printf("make-image: First two arguments must be numbers.\n"); longjmp(esc,1); } if (!primitive(sp3) && !virgin(sp3) && !closure(sp3)) { print_value_escape("make-image: Third argument must be a procedure: ",sp3); } if (!boolean(sp4)) { print_value_escape("make-image: Fourth argument must be boolean: ",sp4); } not_filter = sp4->u.i; rows = (int) floor((float) sp1->u.x); cols = (int) floor((float) sp2->u.x); sp5 = make_image(rows,cols); sp6 = make_complex_image(rows,cols); n = rows*cols; fold_rows = rows/2; fold_cols = cols/2; if (primitive(sp3)) { sexpr *(*func)() = sp3->u.primitive->func; if (sp3->u.primitive->arity != 2) { printf("make-image: Third argument must be a procedure of two numbers.\n"); longjmp(esc,1); } index = 0; for (i = 0; i < rows; i++) { value1->u.x = (not_filter || i <= fold_rows) ? i : i - rows; for (j = 0; j < cols; j++) { value2->u.x = (not_filter || j <= fold_cols) ? j : j - cols; result = func(value1,value2); switch (result->type) { case NUMBER: sp5->u.image->data[index] = result->u.x; sp6->u.complex_image->data[index++] = Complex(result->u.x,0.0); break; case COMPLEX: cflag = 1; sp6->u.complex_image->data[index++] = result->u.z; break; default: printf("make-image: Error at location (%d,%d).\n",i,j); longjmp(esc,1); } } } } else { if (sp3->u.closure->n != 2 && sp3->u.closure->n >= 0) { printf("make-image: Third argument must be a procedure of two arguments.\n"); longjmp(esc,1); } code = compile_app(cons(sp3,list2(value1,value2))); index = 0; for (i = 0; i < rows; i++) { value1->u.x = (not_filter || i <= fold_rows) ? i : i - rows; for (j = 0; j < cols; j++) { value2->u.x = (not_filter || j <= fold_cols) ? j : j - cols; result = virtual_machine(nil,code,global_env_vals); switch (result->type) { case NUMBER: sp5->u.image->data[index] = result->u.x; sp6->u.complex_image->data[index++] = Complex(result->u.x,0.0); break; case COMPLEX: cflag = 1; sp6->u.complex_image->data[index++] = result->u.z; break; default: printf("make-image: Error at location (%d,%d).\n",i,j); longjmp(esc,1); } } } } return cflag ? sp6 : sp5; } sexpr *rgb2color(sexpr *r, sexpr* g, sexpr *b) { int i, rows, cols, n; sexpr *sp; if (!image(r) || !image(g) || !image(b)) { printf("rgb->color-image: Argument is not an image.\n"); longjmp(esc,1); } rows = r->u.image->rows; cols = r->u.image->cols; check_image_sizes("rgb->color-image", rows, cols, g->u.image); check_image_sizes("rgb->color-image", rows, cols, b->u.image); sp = make_color_image(rows,cols); n = rows*cols; for (i=0; i < n; i++) { sp->u.image->data[i*3+0] = r->u.image->data[i]; sp->u.image->data[i*3+1] = g->u.image->data[i]; sp->u.image->data[i*3+2] = b->u.image->data[i]; } return sp; } sexpr *make_hot_image(sexpr *sp1) { int i, rows, cols, n; float value, maximum, minimum; sexpr *sp2; if (!image(sp1)) { printf("make-hot-image: Argument is not an 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; i++) { value = sp1->u.image->data[i]; if (value > maximum) maximum = value; if (value < minimum) minimum = value; } for (i=0; i < n; i++) { value = sp1->u.image->data[i]; sp2->u.image->data[i*3+0]=rhot(minimum,maximum,value); sp2->u.image->data[i*3+1]=ghot(minimum,maximum,value); sp2->u.image->data[i*3+2]=bhot(minimum,maximum,value); } return sp2; } /* Richardson, John L., Visualizing quantum scattering on the CM-2 supercomputer, Computer Physics Communications 63 (1991), pp. 84-94. http://www.hallym.ac.kr/~physics/education/TIPTOP/VLAB/QmSct/complex.html */ sexpr *complex2color(sexpr *sp1) { int i, rows, cols, n; float x, y, radius, d, a, b, R, G, B, scale, maximum, minimum; sexpr *sp2; if (!complex_image(sp1)) { printf("complex-image->color-image: Argument is not a complex 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; i++) { x = sp1->u.complex_image->data[i].r; y = sp1->u.complex_image->data[i].i; if (x > maximum) maximum = x; if (x < minimum) minimum = x; if (y > maximum) maximum = y; if (y < minimum) minimum = y; } scale = 2.0/(maximum-minimum); for (i=0; i < n; i++) { x = scale*(sp1->u.complex_image->data[i].r); y = scale*(sp1->u.complex_image->data[i].i); radius = sqrt(x*x+y*y); a = 0.40824829046386301636 * x; b = 0.70710678118654752440 * y; d = 1.0/(1.0 + radius*radius); R = 0.5 + 0.81649658092772603273 * x * d; G = 0.5 - d * ( a - b ); B = 0.5 - d * ( a + b ); d = 0.5 - radius*d; if (radius < 1) d = -d; sp2->u.image->data[i*3+0] = R + d; sp2->u.image->data[i*3+1] = G + d; sp2->u.image->data[i*3+2] = B + d; } return sp2; } sexpr *color_image_red(sexpr *sp1) { int i, rows, cols; sexpr *sp2; if (!color_image(sp1)) { printf("color-image-red: Argument is not a color image.\n"); longjmp(esc,1); } rows = sp1->u.image->rows; cols = sp1->u.image->cols; sp2 = make_image(rows,cols); for (i=0; i < rows*cols; i++) sp2->u.image->data[i] = sp1->u.image->data[i*3]; return sp2; } sexpr *color_image_green(sexpr *sp1) { int i, rows, cols; sexpr *sp2; if (!color_image(sp1)) { printf("color-image-green: Argument is not a color image.\n"); longjmp(esc,1); } rows = sp1->u.image->rows; cols = sp1->u.image->cols; sp2 = make_image(rows,cols); for (i=0; i < rows*cols; i++) sp2->u.image->data[i] = sp1->u.image->data[i*3+1]; return sp2; } sexpr *color_image_blue(sexpr *sp1) { int i, rows, cols; sexpr *sp2; if (!color_image(sp1)) { printf("color-image-blue: Argument is not a color image.\n"); longjmp(esc,1); } rows = sp1->u.image->rows; cols = sp1->u.image->cols; sp2 = make_image(rows,cols); for (i=0; i < rows*cols; i++) sp2->u.image->data[i] = sp1->u.image->data[i*3+2]; return sp2; } sexpr *rgb2hsi(sexpr *spr, sexpr *spg, sexpr *spb) { sexpr *sph, *sps, *spi; float v1, v2, r, g, b; int rows, cols, j, n; if (!image(spr) || !image(spg) || !image(spb)) { printf("rgb->hsi: Argument is not an image.\n"); longjmp(esc,1); } rows = spr->u.image->rows; cols = spr->u.image->cols; check_image_sizes("rgb->hsi",rows,cols,spg->u.image); check_image_sizes("rgb->hsi",rows,cols,spb->u.image); sph = make_image(rows,cols); sps = make_image(rows,cols); spi = make_image(rows,cols); n = rows*cols; for (j=0; j < n; j++) { r = spr->u.image->data[j]; g = spg->u.image->data[j]; b = spb->u.image->data[j]; /* v1 = (2.0*b-r-g)/2.449489742783178; v2 = (r-g)/2.449489742783178; */ v1 = (2.0*r-g-b)/2.449489742783178; v2 = (g-b)/2.449489742783178; if (v1 != 0) sph->u.image->data[j] = atan2(v2,v1); else sph->u.image->data[j] = 0; sps->u.image->data[j] = sqrt(v1*v1 + v2*v2); spi->u.image->data[j] = (r+g+b)/3; } return cons(sph,list2(sps,spi)); } sexpr *hsi2rgb(sexpr *sph, sexpr *sps, sexpr *spi) { sexpr *spr, *spg, *spb; float v1, v2, h, s, i; int rows, cols, j, n; if (!image(sph) || !image(sps) || !image(spi)) { printf("hsi->rgb: Argument is not an image.\n"); longjmp(esc,1); } rows = sph->u.image->rows; cols = sph->u.image->cols; check_image_sizes("hsi->rgb",rows,cols,sps->u.image); check_image_sizes("hsi->rgb",rows,cols,spi->u.image); spr = make_image(rows,cols); spg = make_image(rows,cols); spb = make_image(rows,cols); n = rows*cols; for (j=0; j < n; j++) { h = sph->u.image->data[j]; s = sps->u.image->data[j]; i = spi->u.image->data[j]; v1 = 2.449489742783178*s*cos(h)/3; v2 = 2.449489742783178*s*sin(h)/2; /* spr->u.image->data[j] = i - v1/2 + v2; spg->u.image->data[j] = i - v1/2 - v2; spb->u.image->data[j] = i + v1; */ spg->u.image->data[j] = i - v1/2 + v2; spb->u.image->data[j] = i - v1/2 - v2; spr->u.image->data[j] = i + v1; } return cons(spr,list2(spg,spb)); } sexpr *upsample_cols(sexpr *sp1) { int i, j, cols; int rows1, rows2, index1, index2; sexpr *sp2; switch (sp1->type) { case IMAGE: cols = sp1->u.image->cols; rows1 = sp1->u.image->rows; rows2 = rows1*2; sp2 = make_image(rows2,cols); for (i = 0; i < rows2; i++) { for (j = 0; j < cols; j++) { index2 = cols*i + j; if (even(i)) { index1 = cols*i/2 + j; sp2->u.image->data[index2] = sp1->u.image->data[index1]; } else sp2->u.image->data[index2] = 0; } } return sp2; case COMPLEX_IMAGE: cols = sp1->u.complex_image->cols; rows1 = sp1->u.complex_image->rows; rows2 = rows1*2; sp2 = make_complex_image(rows2,cols); for (i = 0; i < rows2; i++) { for (j = 0; j < cols; j++) { index2 = cols*i + j; if (even(i)) { index1 = cols*i/2 + j; sp2->u.complex_image->data[index2] = sp1->u.complex_image->data[index1]; } else sp2->u.complex_image->data[index2] = zero; } } return sp2; default: printf("upsample-cols: Argument must be an image.\n"); longjmp(esc,1); } } sexpr *upsample_rows(sexpr *sp1) { int i, j, rows; int cols1, cols2, index1, index2; sexpr *sp2; switch (sp1->type) { case IMAGE: rows = sp1->u.image->rows; cols1 = sp1->u.image->cols; cols2 = cols1*2; sp2 = make_image(rows,cols2); for (j = 0; j < cols2; j++) { for (i = 0; i < rows; i++) { index2 = cols2*i + j; if (even(j)) { index1 = cols1*i + j/2; sp2->u.image->data[index2] = sp1->u.image->data[index1]; } else sp2->u.image->data[index2] = 0; } } return sp2; case COMPLEX_IMAGE: rows = sp1->u.complex_image->rows; cols1 = sp1->u.complex_image->cols; cols2 = cols1*2; sp2 = make_complex_image(rows,cols2); for (j = 0; j < cols2; j++) { for (i = 0; i < rows; i++) { index2 = cols2*i + j; if (even(j)) { index1 = cols1*i + j/2; sp2->u.complex_image->data[index2] = sp1->u.complex_image->data[index1]; } else sp2->u.complex_image->data[index2] = zero; } } return sp2; default: printf("upsample-rows: Argument must be an image.\n"); longjmp(esc,1); } } sexpr *downsample_cols(sexpr *sp1) { int i1, i2, j, cols; int rows1, rows2, index1, index2; sexpr *sp2; switch (sp1->type) { case IMAGE: cols = sp1->u.image->cols; rows1 = sp1->u.image->rows; if (rows1 <= 1) { printf("downsample-cols: Too few rows.\n"); longjmp(esc,1); } rows2 = even(rows1) ? rows1/2 : rows1/2+1; sp2 = make_image(rows2,cols); for (i1 = 0, i2 = 0; i1 < rows1; i1 += 2, i2++) { for (j = 0; j < cols; j++) { index1 = cols*i1 + j; index2 = cols*i2 + j; sp2->u.image->data[index2] = sp1->u.image->data[index1]; } } return sp2; case COMPLEX_IMAGE: cols = sp1->u.complex_image->cols; rows1 = sp1->u.complex_image->rows; if (rows1 <= 1) { printf("downsample-cols: Too few rows.\n"); longjmp(esc,1); } rows2 = even(rows1) ? rows1/2 : rows1/2+1; sp2 = make_complex_image(rows2,cols); for (i1 = 0, i2 = 0; i1 < rows1; i1 += 2, i2++) { for (j = 0; j < cols; j++) { index1 = cols*i1 + j; index2 = cols*i2 + j; sp2->u.complex_image->data[index2] = sp1->u.complex_image->data[index1]; } } return sp2; default: printf("downsample-cols: Argument is not an image.\n"); longjmp(esc,1); } } sexpr *downsample_rows(sexpr *sp1) { int i, j1, j2, rows; int cols1, cols2, index1, index2; sexpr *sp2; switch (sp1->type) { case IMAGE: cols1 = sp1->u.image->cols; rows = sp1->u.image->rows; if (cols1 <= 1) { printf("downsample-rows: Too few columns.\n"); longjmp(esc,1); } cols2 = even(cols1) ? cols1/2 : cols1/2+1; sp2 = make_image(rows,cols2); for (j1 = 0, j2 = 0; j1 < cols1; j1 += 2, j2++) { for (i = 0; i < rows; i++) { index1 = cols1*i + j1; index2 = cols2*i + j2; sp2->u.image->data[index2] = sp1->u.image->data[index1]; } } return sp2; case COMPLEX_IMAGE: cols1 = sp1->u.complex_image->cols; rows = sp1->u.complex_image->rows; if (cols1 <= 1) { printf("downsample-rows: Too few columns.\n"); longjmp(esc,1); } cols2 = even(cols1) ? cols1/2 : cols1/2+1; sp2 = make_complex_image(rows,cols2); for (j1 = 0, j2 = 0; j1 < cols1; j1 += 2, j2++) { for (i = 0; i < rows; i++) { index1 = cols1*i + j1; index2 = cols2*i + j2; sp2->u.complex_image->data[index2] = sp1->u.complex_image->data[index1]; } } return sp2; default: printf("downsample-rows: Argument is not an image.\n"); longjmp(esc,1); } } sexpr *image_transpose(sexpr *sp1) { int i, j, rows, cols; int index1, index2; sexpr *sp2; switch (sp1->type) { case IMAGE: rows = sp1->u.image->rows; cols = sp1->u.image->cols; sp2 = make_image(cols,rows); for (i = 0; i < rows; i++) { for (j = 0; j < cols; j++) { index1 = cols*i + j; index2 = rows*j + i; sp2->u.image->data[index2] = sp1->u.image->data[index1]; } } return sp2; case COMPLEX_IMAGE: rows = sp1->u.complex_image->rows; cols = sp1->u.complex_image->cols; sp2 = make_complex_image(cols,rows); for (i = 0; i < rows; i++) { for (j = 0; j < cols; j++) { index1 = cols*i + j; index2 = rows*j + i; sp2->u.complex_image->data[index2] = sp1->u.complex_image->data[index1]; } } return sp2; default: printf("image-transpose: Argument is not an image.\n"); longjmp(esc,1); } } sexpr *array2complex(sexpr *sp1) { int i, j, rows, cols, index; sexpr *row, *sp2; rows = sp1->u.vector->length; cols = sp1->u.vector->vector[0]->u.vector->length; sp2 = make_complex_image(rows,cols); for (i = 0; i < rows; i++) { row = sp1->u.vector->vector[i]; /* Need to check row length */ for (j = 0; j < cols; j++) { index = cols*i + j; switch (row->u.vector->vector[j]->type) { case NUMBER: sp2->u.complex_image->data[index] = Complex(row->u.vector->vector[j]->u.x,0.0); break; case COMPLEX: sp2->u.complex_image->data[index] = row->u.vector->vector[j]->u.z; break; default: printf("array->complex: Non-number at location (%d, %d).\n",i,j); longjmp(esc,1); } } } return sp2; } sexpr *array2image(sexpr *sp1) { int i, j, rows, cols, index; sexpr *row, *sp2; if (!vector(sp1) || !vector(sp1->u.vector->vector[0])) { print_value_escape("array->image: Argument is not an array: ",sp1); } rows = sp1->u.vector->length; cols = sp1->u.vector->vector[0]->u.vector->length; for (i = 1; i < rows; i++) { if (!vector(sp1->u.vector->vector[i]) || sp1->u.vector->vector[i]->u.vector->length != cols) { print_value_escape("array->image: Argument is not an array: ",sp1); } } sp2 = make_image(rows,cols); for (i = 0; i < rows; i++) { row = sp1->u.vector->vector[i]; for (j = 0; j < cols; j++) { index = cols*i + j; switch (row->u.vector->vector[j]->type) { case NUMBER: sp2->u.image->data[index] = (float) row->u.vector->vector[j]->u.x; break; case COMPLEX: return array2complex(sp1); default: printf("array->image: Non-number at location (%d, %d).\n",i,j); longjmp(esc,1); } } } return sp2; } sexpr *list2vector(sexpr *ls) { int i, len; sexpr *sp; if (!pair(ls) && !null(ls)) { printf("list->vector: Argument is not a list.\n"); longjmp(esc,1); } len = length(ls); sp = smalloc++; sp->type = VECTOR; sp->u.vector = malloc(sizeof(vector)); sp->u.vector->length = len; sp->u.vector->vector = malloc(len*sizeof(sexpr)); for (i = 0; i < len; i++) { sp->u.vector->vector[i] = car(ls); ls = cdr(ls); } return sp; } sexpr *vector2list(sexpr *sp1) { int i, len; sexpr *sp2; if (!vector(sp1)) { printf("vector->list: Argument is not a vector.\n"); longjmp(esc,1); } len = sp1->u.vector->length; sp2 = smalloc++; *sp2 = *nil; for (i = len-1; i >= 0; i--) sp2 = cons(sp1->u.vector->vector[i],sp2); return sp2; } sexpr *convolve(sexpr *sp1, sexpr *sp2) { int i1, j1, i2, j2, mn; int rows, cols, m, n, m2, n2, index1, index2; float sum, *data1, *data2, *data3; fcomplex csum, *cdata1, *cdata2, *cdata3; sexpr *sp3; switch (sp1->type) { case IMAGE: rows = sp1->u.image->rows; cols = sp1->u.image->cols; data1 = sp1->u.image->data; switch (sp2->type) { case VECTOR: return convolve(sp1,array2image(sp2)); case IMAGE: m = sp2->u.image->rows; n = sp2->u.image->cols; data2 = sp2->u.image->data; m2 = even(m) ? m/2-1 : (int) floor(m/2); n2 = even(n) ? n/2-1 : (int) floor(n/2); sp3 = make_image(rows,cols); data3 = sp3->u.image->data; for (i1 = 0; i1 < rows; i1++) { for (j1 = 0; j1 < cols; j1++) { sum = 0; for (i2 = 0; i2 < m; i2++) { for (j2 = 0; j2 < n; j2++) { index1 = cols*mod(i1+i2-m2,rows) + mod(j1+j2-n2,cols); index2 = n*(m-i2-1) + (n-j2-1); sum += data1[index1]*data2[index2]; } } data3[cols*i1 + j1] = sum; } } return sp3; case COMPLEX_IMAGE: m = sp2->u.complex_image->rows; n = sp2->u.complex_image->cols; mn = m*n; cdata2 = sp2->u.complex_image->data; m2 = even(m) ? m/2-1 : (int) floor(m/2); n2 = even(n) ? n/2-1 : (int) floor(n/2); sp3 = make_complex_image(rows,cols); cdata3 = sp3->u.complex_image->data; for (i1 = 0; i1 < rows; i1++) { for (j1 = 0; j1 < cols; j1++) { csum = zero; for (i2 = 0; i2 < m; i2++) { for (j2 = 0; j2 < n; j2++) { index1 = cols*mod(i1+i2-m2,rows) + mod(j1+j2-n2,cols); index2 = n*(m-i2-1) + (n-j2-1); csum = Cadd(csum,RCmul(data1[index1],Conjg(cdata2[index2]))); } } cdata3[cols*i1 + j1] = csum; } } return sp3; default: printf("convolve: Incompatible kernel.\n"); longjmp(esc,1); } case COMPLEX_IMAGE: rows = sp1->u.image->rows; cols = sp1->u.image->cols; cdata1 = sp1->u.complex_image->data; switch (sp2->type) { case VECTOR: return convolve(sp1,array2image(sp2)); case IMAGE: m = sp2->u.complex_image->rows; n = sp2->u.complex_image->cols; data2 = sp2->u.image->data; m2 = even(m) ? m/2-1 : (int) floor(m/2); n2 = even(n) ? n/2-1 : (int) floor(n/2); sp3 = make_complex_image(rows,cols); cdata3 = sp3->u.complex_image->data; for (i1 = 0; i1 < rows; i1++) { for (j1 = 0; j1 < cols; j1++) { csum = zero; for (i2 = 0; i2 < m; i2++) { for (j2 = 0; j2< n; j2++) { index1 = cols*mod(i1+i2-m2,rows) + mod(j1+j2-n2,cols); index2 = n*(m-i2-1) + (n-j2-1); csum = Cadd(csum,RCmul(data2[index2],cdata1[index1])); } } cdata3[cols*i1 + j1] = csum; } } return sp3; case COMPLEX_IMAGE: m = sp2->u.complex_image->rows; n = sp2->u.complex_image->cols; mn = m*n; cdata2 = sp2->u.complex_image->data; m2 = even(m) ? m/2-1 : (int) floor(m/2); n2 = even(n) ? n/2-1 : (int) floor(n/2); sp3 = make_complex_image(rows,cols); cdata3 = sp3->u.complex_image->data; for (i1 = 0; i1 < rows; i1++) { for (j1 = 0; j1 < cols; j1++) { csum = zero; for (i2 = 0; i2 < m; i2++) { for (j2 = 0; j2 < n; j2++) { index1 = cols*mod(i1+i2-m2,rows) + mod(j1+j2-n2,cols); index2 = n*(m-i2-1) + (n-j2-1); csum = Cadd(csum,Cmul(cdata1[index1],Conjg(cdata2[index2]))); } } cdata3[cols*i1 + j1] = csum; } } return sp3; default: printf("convolve: Incompatible kernel.\n"); longjmp(esc,1); } default: printf("convolve: First argument must be an image.\n"); longjmp(esc,1); } } sexpr *outline(sexpr *sp1, sexpr *sp2, sexpr *sp3) { int i, j, index0, index1; int rows, cols, edge, non_edge; float *data1, *data4; sexpr *sp4; if (!image(sp1)) { print_value_escape("outline: First argument must be image: ",sp1); } if (!number(sp2) || !number(sp3)) { printf("outline: Second and third arguments must be numbers.\n"); } rows = sp1->u.image->rows; cols = sp1->u.image->cols; data1 = sp1->u.image->data; edge = sp2->u.x; non_edge = sp3->u.x; sp4 = make_image(rows,cols); data4 = sp4->u.image->data; for (i = 0; i < rows; i++) { index0 = i*cols; for (j = 0; j < cols; j++) { data4[index0+j] = non_edge; } } for (i = 0; i < rows-1; i++) { index0 = i*cols; index1 = (i+1)*cols; for (j = 0; j < cols; j++) { if (data1[index0+j] + data1[index1 + j] == 1) { data4[index0+j] = edge; } } } index0 = (rows-1)*cols; for (j = 0; j < cols; j++) { if (data1[index0+j] + data1[j] == 1) { data4[index0+j] = edge; } } for (i = 0; i < rows; i++) { index0 = i*cols; for (j = 1; j < cols-1; j++) { if (data1[index0+j] + data1[index0 + j + 1] == 1) { data4[index0+j] = edge; } } } for (i = 0; i < rows; i++) { index0 = i*cols; if (data1[index0+cols-1] + data1[index0] == 1) { data4[index0+cols-1] = edge; } } return sp4; } sexpr *perimeters(sexpr *sp1) { int i, j, index0, index1, index2; int rows, cols, maximum, n; float *data1, *data2; sexpr *sp2; if (!image(sp1)) { print_value_escape("perimeters: Argument must be image: ",sp1); } rows = sp1->u.image->rows; cols = sp1->u.image->cols; data1 = sp1->u.image->data; n = rows*cols; maximum = 0; for (i = 0; i < n; i++) if (data1[i] > maximum) maximum = data1[i]; sp2 = make_image(1,maximum+1); data2 = sp2->u.image->data; for (i = 0; i <= maximum; i++) data2[i] = 0.0; for (i = 0; i < rows-1; i++) { index0 = i*cols; index1 = (i+1)*cols; for (j = 0; j < cols; j++) { index2 = (int) data1[index0 + j]; if (index2 > 0) { if (data1[index1 + j] == 0) { data2[index2] = data2[index2] + 1; } } index2 = (int) data1[index1 + j]; if (index2 > 0) { if (data1[index0 + j] == 0) { data2[index2] = data2[index2] + 1; } } } } index0 = (rows-1)*cols; for (j = 0; j < cols; j++) { index2 = (int) data1[index0 + j]; if (index2 > 0) { if (data1[j] == 0) { data2[index2] = data2[index2] + 1; } } index2 = (int) data1[j]; if (index2 > 0) { if (data1[index0 + j] == 0) { data2[index2] = data2[index2] + 1; } } } for (i = 0; i < rows; i++) { index0 = i*cols; for (j = 1; j < cols-1; j++) { index2 = data1[index0 + j]; if (index2 > 0) { if (data1[index0 + j + 1] == 0) { data2[index2] = data2[index2] + 1; } } index2 = data1[index0 + j + 1]; if (index2 > 0) { if (data1[index0 + j] == 0) { data2[index2] = data2[index2] + 1; } } } } for (i = 0; i < rows; i++) { index0 = i*cols; index2 = (int) data1[index0 + cols - 1]; if (index2 > 0) { if (data1[index0 + cols] == 0) { data2[index2] = data2[index2] + 1; } } index2 = (int) data1[index0 + cols]; if (index2 > 0) { if (data1[index0 + cols - 1] == 0) { data2[index2] = data2[index2] + 1; } } } return image2array(sp2)->u.vector->vector[0]; } sexpr *areas(sexpr *sp1) { int i, index; int rows, cols, maximum, n; float *data1, *data2; sexpr *sp2; if (!image(sp1)) { print_value_escape("areas: Argument must be image: ",sp1); } rows = sp1->u.image->rows; cols = sp1->u.image->cols; data1 = sp1->u.image->data; n = rows*cols; maximum = 0; for (i = 0; i < n; i++) if (data1[i] > maximum) maximum = data1[i]; sp2 = make_image(1,maximum+1); data2 = sp2->u.image->data; for (i = 0; i <= maximum; i++) data2[i] = 0.0; for (i = 0; i < n; i++) { index = (int) data1[i]; if (index > 0) { data2[index] = data2[index] + 1; } } return image2array(sp2)->u.vector->vector[0]; } sexpr *centers_of_mass(sexpr *sp1) { int i, j, index; int rows, cols, maximum, n; int *r, *c, *a; float *data; sexpr *sp2; if (!image(sp1)) { print_value_escape("centers-of-mass: Argument must be image: ",sp1); } rows = sp1->u.image->rows; cols = sp1->u.image->cols; data = sp1->u.image->data; n = rows*cols; maximum = 0; for (i = 0; i < n; i++) if (data[i] > maximum) maximum = data[i]; maximum++; r = (int *) malloc(maximum*sizeof(int)); c = (int *) malloc(maximum*sizeof(int)); a = (int *) malloc(maximum*sizeof(int)); for (i = 0; i < maximum; i++) { r[i] = 0; c[i] = 0; a[i] = 0; } sp2 = make_vector(maximum); for (i = 0; i < rows; i++) { for (j = 0; j < cols; j++) { index = (int) data[i*cols + j]; r[index] += i; c[index] += j; a[index]++; } } for (i = 0; i < maximum; i++) { sp2->u.vector->vector[i] = cons(num2exp((float) r[i]/(float) a[i]),cons(num2exp((float) c[i]/(float) a[i]),nil)); } return sp2; } sexpr *bounding_boxes(sexpr *sp1) { int i, j, index; int rows, cols, maximum, n; int *rmin, *rmax, *cmin, *cmax; float *data; sexpr *sp2; if (!image(sp1)) { print_value_escape("bounding-boxes: Argument must be image: ",sp1); } rows = sp1->u.image->rows; cols = sp1->u.image->cols; data = sp1->u.image->data; n = rows*cols; maximum = 0; for (i = 0; i < n; i++) if (data[i] > maximum) maximum = data[i]; maximum++; rmin = (int *) malloc(maximum*sizeof(int)); rmax = (int *) malloc(maximum*sizeof(int)); cmin = (int *) malloc(maximum*sizeof(int)); cmax = (int *) malloc(maximum*sizeof(int)); for (i = 0; i < maximum; i++) { rmin[i] = 0; rmax[i] = 0; cmin[i] = 0; cmax[i] = 0; } sp2 = make_vector(maximum); for (i = 0; i < rows; i++) { for (j = 0; j < cols; j++) { index = (int) data[i*cols + j]; rmax[index] = i; cmax[index] = j; } } for (i = rows-1; i >= 0; i--) { for (j = cols-1; j >= 0; j--) { index = (int) data[i*cols + j]; rmin[index] = i; cmin[index] = j; } } for (i = 0; i < maximum; i++) { sp2->u.vector->vector[i] = cons(num2exp((float) rmin[i]),cons(num2exp((float) cmin[i]),cons(num2exp((float) rmax[i]),cons(num2exp((float) cmax[i]), nil)))); } return sp2; } sexpr *label(sexpr *sp1) { int i, j, index0, index1, n; int rows, cols; int a, b, c; float *data1, *data2; int *alias; sexpr *sp2; if (!image(sp1)) { print_value_escape("label: First argument must be image: ",sp1); } rows = sp1->u.image->rows; cols = sp1->u.image->cols; data1 = sp1->u.image->data; c = 1; sp2 = make_image(rows,cols); data2 = sp2->u.image->data; n = rows*cols; for (i = 0; i < n; i++) data2[i] = 0.0; alias = malloc(n*sizeof(int)); for (i = 0; i < n; i++) alias[i] = i; for (i = 0; i < rows; i++) { for (j = 0; j < cols; j++) { if (data1[cols*i + j] > 0) { a = (int) data2[cols*mod(i-1,rows) + j]; b = (int) data2[cols*i + mod(j-1,cols)]; if (!a && !b) { data2[cols*i + j] = (float) c++; } else if (a && b) { if (alias[a] == alias[b]) { data2[cols*i + j] = (float) alias[a]; } else { /* use lower and set higher alias to lower */ if (alias[a] < alias[b]) { data2[cols*i + j] = (float) alias[a]; alias[b] = alias[a]; } else { data2[cols*i + j] = (float) alias[b]; alias[a] = alias[b]; } } } else if (a) { data2[cols*i + j] = (float) alias[a]; } else { data2[cols*i + j] = (float) alias[b]; } } } } /* find lowest label in each equivalence class */ for (i = 0; i < n; i++) { index0 = (int) data2[i]; if (index0 > 0) { index1 = alias[index0]; while (index0 != index1) { index0 = index1; index1 = alias[index0]; } data2[i] = (float) index1; } } /* renumber labels */ for (i = 1, j = 1; i < (int) c; i++) { if (alias[i] == i) alias[i] = j++; } /* re-label image */ for (i = 0; i < n; i++) { data2[i] = (float) alias[(int) data2[i]]; } free(alias); return sp2; } sexpr *distance_transform(sexpr *sp1) { int i, j, k, l, index0, index1; int rows, cols, m, n; float sum, x, *data1, *data2, *kdata1, *kdata2; sexpr *sp2; sexpr *kernel1 = make_image(5,5); sexpr *kernel2 = make_image(5,5); if (!image(sp1)) { print_value_escape("distance-transform: Argument must be an image: ",sp1); } rows = sp1->u.image->rows; cols = sp1->u.image->cols; n = (int) sqrt(rows*rows+cols*cols)/2.0; data1 = sp1->u.image->data; kdata1 = kernel1->u.image->data; kdata2 = kernel2->u.image->data; kdata1[0] = kdata2[24] = 2.8284; kdata1[1] = kdata2[23] = 2.2; kdata1[2] = kdata2[22] = 2.0; kdata1[3] = kdata2[21] = 2.2; kdata1[4] = kdata2[20] = 2.8284; kdata1[5] = kdata2[19] = 2.2; kdata1[6] = kdata2[18] = 1.4; kdata1[7] = kdata2[17] = 1.0; kdata1[8] = kdata2[16] = 1.4; kdata1[9] = kdata2[15] = 2.2; kdata1[10] = kdata2[14] = 2.0; kdata1[11] = kdata2[13] = 1.0; kdata1[12] = kdata2[12] = 0.0; sp2 = outline(sp1,num2exp(0.0),num2exp(10e7)); data2 = sp2->u.image->data; for (m = 0; m < n; m++) { for (i = 0; i < rows; i++) { for (j = 0; j < cols; j++) { sum = 10e7; for (k = 2; k < 5; k++) { index0 = cols*mod(i+k-2,rows); index1 = 5*k; for (l = 0; l < 5; l++) { if (index1 + l >= 12) { x = data2[index0+mod(j+l-2,cols)]+kdata2[index1 + l]; sum = min(sum,x); } } } data2[cols*i + j] = sum; } } for (i = rows-1; i >= 0; i--) { for (j = cols-1; j >= 0; j--) { sum = 10e7; for (k = 0; k < 3; k++) { index0 = cols*mod(i+k-2,rows); index1 = 5*k; for (l = 0; l < 5; l++) { if (index1 + l <= 12) { x = data2[index0+mod(j+l-2,cols)]+kdata1[index1 + l]; sum = min(sum,x); } } } data2[cols*i + j] = sum; } } } n = rows*cols; for (i = 0; i < n; i++) data2[i] *= data1[i]; return sp2; } sexpr *set_image_display_scale(sexpr *scale) { if (!number(scale)) { printf("set-image-scale!: Illegal argument.\n"); longjmp(esc,1); } image_display_scale = scale->u.x; return undefined; } sexpr *set_graphic_default_color(sexpr *color) { if (!pair(color)) { printf("set-graphic-color!: Illegal argument.\n"); longjmp(esc,1); } graphic$default_color = color; return undefined; } sexpr *draw_transparent(sexpr *p1, sexpr *sp) { sexpr *len = cadr(sp->u.sexpr); sexpr *p2; double x, y, heading; heading = p1->u.plumber->heading; x = p1->u.plumber->x + (len->u.x/300.0)*cos(heading); y = p1->u.plumber->y + (len->u.x/300.0)*sin(heading); p2 = make_plumber(x, y, heading); return p2; } sexpr *draw_bend(sexpr *p1, sexpr *sp) { sexpr *second = cadr(sp->u.sexpr); sexpr *p2; double x, y, heading; x = p1->u.plumber->x; y = p1->u.plumber->y; heading = p1->u.plumber->heading; p2 = make_plumber(x,y,angle_difference(heading,deg2rad(second->u.x))); return p2; } sexpr *graphic2postscript(sexpr *sp1, sexpr *sp2) { sexpr *port, *p; char text[STRLEN]; if (!graphic(sp1)) { printf("graphic->postscript: First argument must be graphic object.\n"); longjmp(esc,1); } if ((!string(sp2))) { printf("graphic->postscript: Illegal filename.\n"); longjmp(esc,1); } port = open_output_file(sp2); sprintf(text,"#",gtype(sp1)->u.symbol->name); generate_postscript_header(port->u.file,1.0,text); p = make_plumber(1.0,-0.6667,HALFPI); draw_plumber_postscript(p,graphic$default_color,port->u.file); p = draw_graphic_postscript(p,sp1,port->u.file); draw_plumber_postscript(p,graphic$default_color,port->u.file); fprintf(port->u.file,"showpage\n"); close_output_port(port); return undefined; } sexpr *draw_graphic_postscript(sexpr *p, sexpr *sp, FILE *file) { sexpr *g, *first; g = sp->u.sexpr; if (null(g)) return p; first = car(g); if (eq(first,symbol$straight)) return draw_straight_postscript(p,sp,file); else if (eq(first,symbol$spot)) return draw_spot_postscript(p,sp,file); else if (eq(first,symbol$transparent)) return draw_transparent(p,sp); else if (eq(first,symbol$adjoin)) return draw_adjoin_postscript(p,sp,file); else if (eq(first,symbol$adorn)) { draw_adorn_postscript(p,sp,file); return p; } else if (eq(first,symbol$bend)) return draw_bend(p,sp); else if (eq(first,symbol$text)) return draw_text_postscript(p,sp,file); else { printf("Unrecognized graphic type.\n"); return p; } } sexpr *draw_text_postscript(sexpr *p, sexpr *sp, FILE *file) { sexpr *text = cadr(sp->u.sexpr); sexpr *color = cadr(cdr(sp->u.sexpr)); sexpr *size = cadr(cddr(sp->u.sexpr)); sexpr *r = car(color); sexpr *g = cadr(color); sexpr *b = car(cddr(color)); if (!string(text)) { printf("graphic->postscript: Illegal text argument.\n"); longjmp(esc,1); } if (!number(r) || !number(g) || !number(b)) { printf("graphic->postscript: Illegal color argument.\n"); longjmp(esc,1); } if (!number(size)) { printf("graphic->postscript: Illegal size argument.\n"); longjmp(esc,1); } fprintf(file,"%f %f moveto\n",(1.0-p->u.plumber->y)/2.0,p->u.plumber->x/2.0); fprintf(file,"90 rotate\n"); fprintf(file,"/Helvetica findfont\n"); fprintf(file,"%f scalefont\n",size->u.x*0.02); fprintf(file,"setfont\n"); if (r->u.x + g->u.x + b->u.x > 764) fprintf(file,"%f %f %f setrgbcolor\n",0.0,0.0,0.0); else fprintf(file,"%f %f %f setrgbcolor\n",r->u.x/128,g->u.x/128,b->u.x/128); fprintf(file,"(%s) show\n",text->u.text); fprintf(file,"-90 rotate\n"); return p; } sexpr *draw_plumber_postscript(sexpr *p, sexpr *color, FILE *file) { sexpr *r, *g, *b; double c, s; float x0, y0, x1, y1, x2, y2; c = cos(p->u.plumber->heading-PI/2.0)/60.0; s = sin(p->u.plumber->heading-PI/2.0)/60.0; r = car(color); g = cadr(color); b = car(cddr(color)); if (!number(r) || !number(g) || !number(b)) { printf("graphic->postscript: Illegal color argument.\n"); longjmp(esc,1); } x0 = p->u.plumber->x-C30*c+S30*s; y0 = 1.0-p->u.plumber->y+C30*s+S30*c; x1 = p->u.plumber->x-s; y1 = 1.0-p->u.plumber->y-c; x2 = p->u.plumber->x+C30*c+S30*s; y2 = 1.0-p->u.plumber->y-C30*s+S30*c; fprintf(file,"%f setlinewidth\n",1/864.0); if (r->u.x + g->u.x + b->u.x > 764) fprintf(file,"%f %f %f setrgbcolor\n",0.0,0.0,0.0); else fprintf(file,"%f %f %f setrgbcolor\n",r->u.x/128,g->u.x/128,b->u.x/128); fprintf(file,"%f %f moveto\n",y0/2.0,x0/2.0); fprintf(file,"%f %f lineto\n",y1/2.0,x1/2.0); fprintf(file,"%f %f lineto\n",y2/2.0,x2/2.0); fprintf(file,"%f %f lineto\n",y0/2.0,x0/2.0); fprintf(file,"stroke\n"); return p; } sexpr *draw_straight_postscript(sexpr *p1, sexpr *sp, FILE *file) { sexpr *len = cadr(sp->u.sexpr); sexpr *color = car(cddr(sp->u.sexpr)); sexpr *width = cadr(cddr(sp->u.sexpr)); sexpr *p2; sexpr *r, *g, *b; double x, y, heading; heading = p1->u.plumber->heading; x = p1->u.plumber->x + (len->u.x/300.0)*cos(heading); y = p1->u.plumber->y + (len->u.x/300.0)*sin(heading); p2 = make_plumber(x, y, heading); if (!number(len)) { printf("graphic->postscript: Illegal length argument.\n"); longjmp(esc,1); } r = car(color); g = cadr(color); b = car(cddr(color)); if (!number(r) || !number(g) || !number(b)) { printf("graphic->postscript: Illegal color argument.\n"); longjmp(esc,1); } if (!number(width)) { printf("graphic->postscript: Illegal width argument.\n"); longjmp(esc,1); } fprintf(file,"%f setlinewidth\n",width->u.x/864.0); if (r->u.x + g->u.x + b->u.x > 764) fprintf(file,"%f %f %f setrgbcolor\n",0.0,0.0,0.0); else fprintf(file,"%f %f %f setrgbcolor\n",r->u.x/128,g->u.x/128,b->u.x/128); fprintf(file,"%f %f %f %f displine\n",(1.0-p1->u.plumber->y)/2.0,(p1->u.plumber->x)/2.0,(1.0-y)/2.0,x/2.0); return p2; } sexpr *draw_spot_postscript(sexpr *p, sexpr *sp, FILE *file) { sexpr *width = cadr(sp->u.sexpr); sexpr *color = car(cddr(sp->u.sexpr)); sexpr *r, *g, *b; r = car(color); g = cadr(color); b = car(cddr(color)); if (!number(r) || !number(g) || !number(b)) { printf("graphic->postscript: Illegal color argument.\n"); longjmp(esc,1); } if (!number(width)) { printf("graphic->postscript: Illegal width argument.\n"); longjmp(esc,1); } fprintf(file,"newpath\n"); if (r->u.x + g->u.x + b->u.x > 764) fprintf(file,"%f %f %f setrgbcolor\n",0.0,0.0,0.0); else fprintf(file,"%f %f %f setrgbcolor\n",r->u.x/128,g->u.x/128,b->u.x/128); fprintf(file,"%f %f %f 0.0 360.0 arc\n",(1.0-p->u.plumber->y)/2.0,(p->u.plumber->x)/2.0,width->u.x/432.0); fprintf(file,"fill\n"); return p; } sexpr *draw_adjoin_postscript(sexpr *p, sexpr *sp, FILE *file) { sexpr *g=sp->u.sexpr; if (null(g)) return p; return draw_graphic_postscript(draw_graphic_postscript(p,gcar(sp),file),gcdr(sp),file); } sexpr *draw_adorn_postscript(sexpr *p, sexpr *sp, FILE *file) { sexpr *g=sp->u.sexpr; if (null(g)) return p; draw_graphic_postscript(p,gcar(sp),file); draw_graphic_postscript(p,gcdr(sp),file); return p; } #ifdef GL sexpr *draw_graphic(sexpr *p, sexpr *sp) { sexpr *g, *first; g = sp->u.sexpr; if (null(g)) return p; first = car(g); if (eq(first,symbol$straight)) return draw_straight(p,sp); else if (eq(first,symbol$spot)) return draw_spot(p,sp); else if (eq(first,symbol$transparent)) return draw_transparent(p,sp); else if (eq(first,symbol$adjoin)) return draw_adjoin(p,sp); else if (eq(first,symbol$adorn)) { draw_adorn(p,sp); return p; } else if (eq(first,symbol$bend)) return draw_bend(p,sp); else if (eq(first,symbol$text)) return draw_text(p,sp); else { printf("Unrecognized graphic type.\n"); return p; } } sexpr *draw_text(sexpr *p, sexpr *sp) { char *c; sexpr *text = cadr(sp->u.sexpr); sexpr *color = cadr(cdr(sp->u.sexpr)); sexpr *size = cadr(cddr(sp->u.sexpr)); sexpr *r, *g, *b; double s; if (!string(text)) { printf("draw-text: Illegal text argument.\n"); longjmp(esc,1); } r = car(color); g = cadr(color); b = car(cddr(color)); if (!number(r) || !number(g) || !number(b)) { printf("draw-text: Illegal color argument.\n"); longjmp(esc,1); } if (!number(size)) { printf("draw-text: Illegal size argument.\n"); longjmp(esc,1); } if (size->u.x < 2.0) { glColor3f(r->u.x/128,g->u.x/128,b->u.x/128); glRasterPos2f(p->u.plumber->x, p->u.plumber->y); for (c = text->u.text; *c; c++) glutBitmapCharacter(GLUT_BITMAP_8_BY_13, *c); } else { s=size->u.x*0.0003; glPushMatrix(); glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); glEnable(GL_BLEND); glEnable(GL_LINE_SMOOTH); glLineWidth(1.0); glTranslatef(p->u.plumber->x, p->u.plumber->y, 0.0); glScalef(s, s, s); glColor3f(r->u.x/128,g->u.x/128,b->u.x/128); for (c = text->u.text; *c; c++) glutStrokeCharacter(GLUT_STROKE_ROMAN, *c); glPopMatrix(); } return p; } sexpr *draw_plumber(sexpr *p, sexpr *color) { sexpr *r, *g, *b; double x, y, c, s; c = cos(p->u.plumber->heading-PI/2.0); s = sin(p->u.plumber->heading-PI/2.0); x = p->u.plumber->x; y = p->u.plumber->y; r = car(color); g = cadr(color); b = car(cddr(color)); if (!number(r) || !number(g) || !number(b)) { printf("draw-plumber: Illegal color argument.\n"); longjmp(esc,1); } glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); glEnable(GL_BLEND); glEnable(GL_LINE_SMOOTH); glLineWidth(1.0); glBegin(GL_LINE_LOOP); glColor3f(r->u.x/128,g->u.x/128,b->u.x/128); glVertex2f(x-(C30*c - S30*s)/60.0,y+(-C30*s - S30*c)/60.0); glVertex2f(x-s/60.0,y+c/60.0); glVertex2f(x-(-C30*c - S30*s)/60.0,y+(C30*s - S30*c)/60.0); glEnd(); return p; } sexpr *draw_straight(sexpr *p1, sexpr *sp) { sexpr *len = cadr(sp->u.sexpr); sexpr *color = car(cddr(sp->u.sexpr)); sexpr *width = cadr(cddr(sp->u.sexpr)); sexpr *p2; sexpr *r, *g, *b; double x, y, heading; if (!number(len)) { printf("draw-straight: Illegal length argument.\n"); longjmp(esc,1); } heading = p1->u.plumber->heading; x = p1->u.plumber->x + (len->u.x/300.0)*cos(heading); y = p1->u.plumber->y + (len->u.x/300.0)*sin(heading); p2 = make_plumber(x, y, heading); r = car(color); g = cadr(color); b = car(cddr(color)); if (!number(r) || !number(g) || !number(b)) { printf("draw-straight: Illegal color argument.\n"); longjmp(esc,1); } if (!number(width)) { printf("draw-straight: Illegal width argument.\n"); longjmp(esc,1); } glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); glEnable(GL_BLEND); glEnable(GL_LINE_SMOOTH); glLineWidth(width->u.x); glBegin(GL_LINES); glColor3f(r->u.x/128,g->u.x/128,b->u.x/128); glVertex2f(p1->u.plumber->x,p1->u.plumber->y); glVertex2f(p2->u.plumber->x,p2->u.plumber->y); glEnd(); return p2; } sexpr *draw_spot(sexpr *p, sexpr *sp) { sexpr *width = cadr(sp->u.sexpr); sexpr *color = car(cddr(sp->u.sexpr)); sexpr *r, *g, *b; if (!number(width)) { printf("draw-spot: Illegal width argument.\n"); longjmp(esc,1); } r = car(color); g = cadr(color); b = car(cddr(color)); if (!number(r) || !number(g) || !number(b)) { printf("draw-spot: Illegal color argument.\n"); longjmp(esc,1); } glPointSize(width->u.x); glBegin(GL_POINTS); glColor3f(r->u.x/128,g->u.x/128,b->u.x/128); glVertex2f(p->u.plumber->x,p->u.plumber->y); glEnd(); return p; } sexpr *draw_adjoin(sexpr *p, sexpr *sp) { sexpr *g=sp->u.sexpr; if (null(g)) return p; return draw_graphic(draw_graphic(p,gcar(sp)),gcdr(sp)); } sexpr *draw_adorn(sexpr *p, sexpr *sp) { sexpr *g=sp->u.sexpr; if (null(g)) return p; draw_graphic(p,gcar(sp)); draw_graphic(p,gcdr(sp)); return p; } #endif sexpr *scheme_read(sexpr *port) { char c; sexpr *input; if (!input_port(port)) { print_value_escape("read: Argument is not an input-port: ",port); } if ((c = getc(port->u.file)) != EOF) { ungetc(c,port->u.file); input = parse(port->u.file,' '); if (pair(input)) return car(input); } return undefined; } sexpr *scheme_write(sexpr *sp, sexpr *port) { if (!output_port(port)) { print_value_escape("write: Argument is not an output-port: ",port); } display_or_print(sp,PRINT,TOPLEVEL,port->u.file); return undefined; } int eq(sexpr *s1, sexpr *s2) { if (s1->type == s2->type) { switch(s1->type) { case NIL: return 1; case NUMBER: return (s1->u.x == s2->u.x); case TOKEN: case CHARACTER: return (s1->u.c == s2->u.c); case STRING: return !strcmp(s1->u.text,s2->u.text); case BOOLEAN: return (s1->u.i == s2->u.i); case VIRGIN: case CLOSURE: case PRIMITIVE: case VECTOR: case INPUT_PORT: case OUTPUT_PORT: case PAIR: case SYMBOL: case GRAPHIC: case IMAGE: case LINE: case SKETCH: return (s1 == s2); case COMPLEX: return (s1->u.z.r == s2->u.z.r && s1->u.z.i == s2->u.z.i); default: printf("eq: Illegal type.\n"); return 0; } } else return 0; } sexpr *image_map(double (*func)(double), sexpr *sp1) { int rows = sp1->u.image->rows; int cols = sp1->u.image->cols; int i, n = rows*cols; float *data1 = sp1->u.image->data; sexpr *sp2 = make_image(rows,cols); float *data2 = sp2->u.image->data; for (i = 0; i < n; i++) data2[i] = (float) func((double) data1[i]); return sp2; } sexpr *complex_image_map(fcomplex (*func)(fcomplex), sexpr *sp1) { int rows = sp1->u.complex_image->rows; int cols = sp1->u.complex_image->cols; int i, n = rows*cols; fcomplex *cdata1 = sp1->u.complex_image->data; sexpr *sp2 = make_complex_image(rows,cols); fcomplex *cdata2 = sp2->u.complex_image->data; for (i = 0; i < n; i++) cdata2[i] = func(cdata1[i]); return sp2; } sexpr *complex_image_map_real_result(float (*func)(fcomplex), sexpr *sp1) { int rows = sp1->u.complex_image->rows; int cols = sp1->u.complex_image->cols; int i, n = rows*cols; fcomplex *cdata1 = sp1->u.complex_image->data; sexpr *sp2 = make_image(rows,cols); float *data2 = sp2->u.image->data; for (i = 0; i < n; i++) data2[i] = func(cdata1[i]); return sp2; } sexpr *sine(sexpr *s) { switch (s->type) { case NUMBER: return num2exp(sin(s->u.x)); case COMPLEX: return complex2exp(Csin(s->u.z)); case IMAGE: return image_map(sin,s); case COMPLEX_IMAGE: return complex_image_map(Csin,s); default: printf("sin: Argument is not a number or image.\n"); longjmp(esc,1); } } sexpr *cosine(sexpr *s) { switch (s->type) { case NUMBER: return num2exp(cos(s->u.x)); case COMPLEX: return complex2exp(Ccos(s->u.z)); case IMAGE: return image_map(cos,s); case COMPLEX_IMAGE: return complex_image_map(Ccos,s); default: printf("cos: Argument is not a number or image.\n"); longjmp(esc,1); } } sexpr *tangent(sexpr *s) { switch (s->type) { case NUMBER: return num2exp(tan(s->u.x)); case IMAGE: return image_map(tan,s); case COMPLEX: case COMPLEX_IMAGE: return divide(sine(s),cosine(s)); default: printf("tan: Argument is not a number or image.\n"); longjmp(esc,1); } } sexpr *arcsine(sexpr *s) { switch (s->type) { case NUMBER: return num2exp(asin(s->u.x)); case COMPLEX: return complex2exp(Casin(s->u.z)); case IMAGE: return image_map(asin,s); case COMPLEX_IMAGE: return complex_image_map(Casin,s); default: printf("asin: Argument is not a number or image.\n"); longjmp(esc,1); } } sexpr *arccosine(sexpr *s) { switch (s->type) { case NUMBER: return num2exp(acos(s->u.x)); case COMPLEX: return complex2exp(Cacos(s->u.z)); case IMAGE: return image_map(acos,s); case COMPLEX_IMAGE: return complex_image_map(Cacos,s); default: printf("acos: Argument is not a number or image.\n"); longjmp(esc,1); } } sexpr *logarithm(sexpr *s) { switch(s->type) { case NUMBER: return num2exp(log(s->u.x)); case COMPLEX: return complex2exp(Clog(s->u.z)); case IMAGE: return image_map(log,s); case COMPLEX_IMAGE: return complex_image_map(Clog,s); default: printf("log: Argument is not a number or image.\n"); longjmp(esc,1); } } sexpr *exponential(sexpr *s) { switch(s->type) { case NUMBER: return num2exp(exp(s->u.x)); case COMPLEX: return complex2exp(Cexp(s->u.z)); case IMAGE: return image_map(exp,s); case COMPLEX_IMAGE: return complex_image_map(Cexp,s); default: printf("exp: Argument is not a number or image.\n"); longjmp(esc,1); } } sexpr *image_map_complex_result(fcomplex (*func)(fcomplex), sexpr *sp1) { int rows = sp1->u.image->rows; int cols = sp1->u.image->cols; int i, n = rows*cols; float *data1 = sp1->u.image->data; sexpr *sp2 = make_complex_image(rows,cols); fcomplex *cdata2 = sp2->u.complex_image->data; for (i = 0; i < n; i++) cdata2[i] = func(Complex(data1[i],0.0)); return sp2; } sexpr *square_root(sexpr *s) { switch(s->type) { case NUMBER: if (s->u.x >= 0) { return num2exp(sqrt(s->u.x)); } else { return complex2exp(Complex(0.0,sqrt(fabs(s->u.x)))); } case COMPLEX: return complex2exp(Csqrt(s->u.z)); case IMAGE: if (image_min(s)->u.x >= 0) { return image_map(sqrt,s); } else { return image_map_complex_result(Csqrt, s); } case COMPLEX_IMAGE: return complex_image_map(Csqrt,s); default: printf("sqrt: Argument is not a number or image.\n"); longjmp(esc,1); } } sexpr *scheme_abs(sexpr *s) { switch(s->type) { case NUMBER: return num2exp(fabs(s->u.x)); case COMPLEX: return num2exp(Cabs(s->u.z)); case IMAGE: return image_map(fabs,s); case COMPLEX_IMAGE: return complex_image_map_real_result(Cabs,s); default: printf("abs: Argument is not a number or image.\n"); longjmp(esc,1); } } sexpr *scheme_floor(sexpr *s) { switch(s->type) { case NUMBER: return num2exp(floor(s->u.x)); case IMAGE: return image_map(floor,s); default: printf("floor: Argument is not a real number or image.\n"); longjmp(esc,1); } } sexpr *scheme_ceil(sexpr *s) { switch(s->type) { case NUMBER: return num2exp(ceil(s->u.x)); case IMAGE: return image_map(ceil,s); default: printf("ceil: Argument is not a real number or image.\n"); longjmp(esc,1); } } double random_func(double x) { return floor(x*random()/RAND_MAX); } sexpr *scheme_random(sexpr *s) { switch(s->type) { case NUMBER: return num2exp((int) random_func(s->u.x)); case IMAGE: return image_map(random_func,s); default: printf("random: Argument is not a real number or image.\n"); longjmp(esc,1); } } double angle_func(double x) { return (x < 0) ? PI : 0.0; } float complex_angle_func(fcomplex z) { return atan2(z.i,z.r); } sexpr *real_part(sexpr *s) { switch (s->type) { case NUMBER: return s; case COMPLEX: return num2exp(s->u.z.r); case IMAGE: return s; case COMPLEX_IMAGE: return complex_image_real(s); default: printf("real-part: Argument is not a number or image.\n"); longjmp(esc,1); } } double zero_func(double x) { return 0.0; } sexpr *imag_part(sexpr *s) { switch (s->type) { case NUMBER: return num2exp(0); case COMPLEX: return num2exp(s->u.z.i); case IMAGE: return image_map(zero_func,s); case COMPLEX_IMAGE: return complex_image_imag(s); default: printf("imag-part: Argument is not a number or image.\n"); longjmp(esc,1); } } sexpr *angle(sexpr *s) { switch (s->type) { case NUMBER: return num2exp(angle_func(s->u.x)); case COMPLEX: return num2exp(complex_angle_func(s->u.z)); case IMAGE: return image_map(angle_func,s); case COMPLEX_IMAGE: return complex_image_map_real_result(complex_angle_func,s); default: printf("angle: Argument is not a number or image.\n"); longjmp(esc,1); } } sexpr *conjugate(sexpr *s) { switch (s->type) { case NUMBER: return s; case COMPLEX: return complex2exp(Conjg(s->u.z)); case IMAGE: return s; case COMPLEX_IMAGE: return complex_image_map(Conjg,s); default: printf("conjugate: Argument is not a number or image.\n"); longjmp(esc,1); } } sexpr *expt(sexpr *s1, sexpr *s2) { switch(s1->type) { case NUMBER: switch(s2->type) { case NUMBER: return num2exp(pow(s1->u.x,s2->u.x)); default: return exponential(times(s2,logarithm(s1))); } default: return exponential(times(s2,logarithm(s1))); } } #define torf(x) ((x) ? true : false); sexpr *eqp(sexpr *s1, sexpr *s2) { return torf(eq(s1,s2)); } sexpr *real_equals_image(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_image(rows,cols); float *data3 = sp3->u.image->data; for (i = 0; i < n; i++) data3[i] = (data2[i] == x) ? 1 : 0; return sp3; } sexpr *image_equals_image(sexpr *sp1, sexpr *sp2) { int rows = sp2->u.image->rows; int cols = sp2->u.image->cols; int i, n = rows*cols; float *data1 = sp1->u.image->data; float *data2 = sp2->u.image->data; sexpr *sp3 = make_image(rows,cols); float *data3 = sp3->u.image->data; check_image_sizes("=", rows, cols, sp2->u.image); for (i = 0; i < n; i++) data3[i] = (data1[i] == data2[i]) ? 1 : 0; return sp3; } sexpr *zero_imag_part(sexpr *sp) { return equals(imag_part(sp),num2exp(0.0)); } sexpr *equals(sexpr *sp1, sexpr *sp2) { switch (sp1->type) { case NUMBER: switch (sp2->type) { case NUMBER: return torf(sp1->u.x == sp2->u.x); case COMPLEX: return torf((sp1->u.x == sp2->u.z.r) && (sp2->u.z.i == 0.0)); case IMAGE: return real_equals_image(sp1,sp2); case COMPLEX_IMAGE: return times(real_equals_image(sp1,real_part(sp2)),zero_imag_part(sp2)); default: printf("=: Incompatible arguments.\n"); longjmp(esc,1); } case COMPLEX: switch (sp2->type) { case NUMBER: return torf((sp1->u.z.r == sp2->u.x) && (sp1->u.z.i == 0.0)); case COMPLEX: return torf((sp1->u.z.r == sp2->u.z.r) && (sp1->u.z.i == sp2->u.z.i)); case IMAGE: if (sp1->u.z.i == 0.0) { return equals(real_part(sp1),sp2); } else { return times(sp2,num2exp(0.0)); } case COMPLEX_IMAGE: return times(equals(real_part(sp1),real_part(sp2)),equals(imag_part(sp1),imag_part(sp2))); default: printf("=: Incompatible arguments.\n"); longjmp(esc,1); } case IMAGE: switch (sp2->type) { case NUMBER: case COMPLEX: return equals(sp2,sp1); case IMAGE: return image_equals_image(sp1,sp2); case COMPLEX_IMAGE: return times(image_equals_image(sp1,real_part(sp2)),zero_imag_part(sp2)); default: printf("=: Incompatible arguments.\n"); longjmp(esc,1); } case COMPLEX_IMAGE: switch (sp2->type) { case NUMBER: case COMPLEX: case IMAGE: return equals(sp2,sp1); case COMPLEX_IMAGE: return times(equals(real_part(sp1),real_part(sp2)),equals(imag_part(sp1),imag_part(sp2))); default: printf("=: Incompatible arguments.\n"); longjmp(esc,1); } default: print_value_escape("=: Argument is not a number or image: ",sp1); } } sexpr *real_image_helper(double (func)(double,double), sexpr *sp1, sexpr *sp2) { int rows = sp2->u.image->rows; int cols = sp2->u.image->cols; int i, n = rows*cols; double x = sp1->u.x; float *data2 = sp2->u.image->data; sexpr *sp3 = make_image(rows,cols); float *data3 = sp3->u.image->data; for (i = 0; i < n; i++) data3[i] = (float) func(x, (double) data2[i]); return sp3; } sexpr *image_real_helper(double (func)(double,double), 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; double x = sp2->u.x; sexpr *sp3 = make_image(rows,cols); float *data3 = sp3->u.image->data; for (i = 0; i < n; i++) data3[i] = (float) func((double) data1[i], x); return sp3; } sexpr *image_image_helper(double (func)(double,double), 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_image(rows,cols); float *data3 = sp3->u.image->data; check_image_sizes("real-helper",rows,cols,sp2->u.image); for (i = 0; i < n; i++) data3[i] = (float) func((double) data1[i],(double) data2[i]); return sp3; } sexpr *real_helper(double (func)(double,double), sexpr *sp1, sexpr *sp2) { switch (sp1->type) { case NUMBER: switch (sp2->type) { case NUMBER: return torf(func(sp1->u.x,sp2->u.x)); case IMAGE: return real_image_helper(func,sp1,sp2); default: printf("real-helper: Incompatible argument types.\n"); longjmp(esc,1); } case IMAGE: switch (sp2->type) { case NUMBER: return image_real_helper(func,sp1,sp2); case IMAGE: return image_image_helper(func,sp1,sp2); default: printf("real-helper: Incompatible argument types.\n"); longjmp(esc,1); } default: printf("real-helper: Argument is not a real number or image.\n"); longjmp(esc,1); } } double lt_func(double x, double y) { return (x < y) ? 1 : 0; } double gt_func(double x, double y) { return (x > y) ? 1 : 0; } double leq_func(double x, double y) { return (x <= y) ? 1 : 0; } double geq_func(double x, double y) { return (x >= y) ? 1 : 0; } double min_func(double x, double y) { return (x <= y) ? x : y; } double max_func(double x, double y) { return (x >= y) ? x : y; } sexpr *lt(sexpr *sp1, sexpr *sp2) { return real_helper(lt_func,sp1,sp2); } sexpr *gt(sexpr *sp1, sexpr *sp2) { return real_helper(gt_func,sp1,sp2); } sexpr *leq(sexpr *sp1, sexpr *sp2) { return real_helper(leq_func,sp1,sp2); } sexpr *geq(sexpr *sp1, sexpr *sp2) { return real_helper(geq_func,sp1,sp2); } sexpr *scheme_min(sexpr *sp1, sexpr *sp2) { if (!number(sp1) || !number(sp2)) return real_helper(min_func,sp1,sp2); return num2exp((sp1->u.x >= sp2->u.x) ? sp2->u.x : sp1->u.x); } sexpr *scheme_max(sexpr *sp1, sexpr *sp2) { if (!number(sp1) || !number(sp2)) return real_helper(max_func,sp1,sp2); return num2exp((sp1->u.x >= sp2->u.x) ? sp1->u.x : sp2->u.x); } sexpr *positivep(sexpr *sp) { return gt(sp,num2exp(0.0)); } sexpr *negativep(sexpr *sp) { return lt(sp,num2exp(0.0)); } sexpr *zerop(sexpr *sp) { return equals(sp,num2exp(0.0)); } /* http://www.astro.princeton.edu/~esirko/idl_html_help/A17.html atan(Zy, Zx) = -i alog((Zx + iZy)/sqrt(Zx^2 + Zy^2)) */ fcomplex complex_complex_atan2(fcomplex z1, fcomplex z2) { return Cmul(minus_root_minus_one,Cexp(Cdiv(Cadd(z2,Cmul(root_minus_one,z1)),Csqrt(Cadd(Cmul(z1,z1),Cmul(z2,z2)))))); } fcomplex complex_real_atan2(fcomplex z, float x) { double u = z.r/x; double v = z.i/x; return Cdiv(Csub(Clog(Complex(1-v,u)),Clog(Complex(1+v,-u))),two_root_minus_one); } fcomplex real_complex_atan2(float x, fcomplex z) { double m = z.r*z.r + z.i*z.i; double u = x*z.r/m; double v = x*z.i/m; return Cdiv(Csub(Clog(Complex(1+v,u)),Clog(Complex(1-v,-u))),two_root_minus_one); } sexpr *not_likely_to_be_called(sexpr *sp1, sexpr *sp2) { int rows = sp1->u.complex_image->rows; int cols = sp1->u.complex_image->cols; int i, n = rows*cols; fcomplex *cdata1 = sp1->u.complex_image->data; fcomplex *cdata2 = sp2->u.complex_image->data; sexpr *sp3 = make_complex_image(rows,cols); fcomplex *cdata3 = sp3->u.complex_image->data; check_image_sizes("atan",rows,cols,(image *) sp2->u.complex_image); for (i = 0; i < n; i++) cdata3[i] = complex_complex_atan2(cdata1[i],cdata2[i]); return sp3; } sexpr *scalar2image(sexpr *sp1, sexpr *sp2) { sexpr *z = complex2exp(zero); return plus(times(sp2,z),plus(sp1,z)); } sexpr *real2complex(sexpr *sp) { return make_complex(sp,num2exp(0.0)); } sexpr *arctangent(sexpr *sp1, sexpr *sp2) { switch (sp1->type) { case NUMBER: switch (sp2->type) { case NUMBER: return num2exp(atan2(sp1->u.x,sp2->u.x)); case COMPLEX: if (Cmag(sp2->u.z) == 0.0) { printf("atan: Second argument cannot have zero magnitude.\n"); longjmp(esc,1); } return complex2exp(real_complex_atan2(sp1->u.x,sp2->u.z)); case IMAGE: return real_image_helper(atan2,sp1,sp2); case COMPLEX_IMAGE: return not_likely_to_be_called(scalar2image(sp1,sp2),sp2); default: printf("atan: Illegal or unsupported argument.\n"); longjmp(esc,1); } case COMPLEX: switch (sp2->type) { case NUMBER: if (sp2->u.x == 0.0) { printf("atan: Second argument cannot be zero.\n"); longjmp(esc,1); } return complex2exp(complex_real_atan2(sp1->u.z,sp2->u.x)); case COMPLEX: return complex2exp(complex_complex_atan2(sp1->u.z,sp2->u.z)); case IMAGE: return not_likely_to_be_called(scalar2image(sp1,sp2),real2complex(sp2)); case COMPLEX_IMAGE: return not_likely_to_be_called(scalar2image(sp1,sp2),sp2); default: printf("atan: Illegal or unsupported argument type.\n"); longjmp(esc,1); } case IMAGE: switch (sp2->type) { case NUMBER: return image_real_helper(atan2,sp1,sp2); case COMPLEX: return not_likely_to_be_called(real2complex(sp1),scalar2image(sp2,sp1)); case IMAGE: return image_image_helper(atan2,sp1,sp2); case COMPLEX_IMAGE: return not_likely_to_be_called(real2complex(sp1),sp2); default: printf("atan: Illegal or unsupported argument type.\n"); longjmp(esc,1); } case COMPLEX_IMAGE: switch (sp2->type) { case NUMBER: return not_likely_to_be_called(sp1,scalar2image(sp2,sp1)); case COMPLEX: return not_likely_to_be_called(sp1,scalar2image(sp2,sp1)); case IMAGE: return not_likely_to_be_called(sp1,real2complex(sp2)); case COMPLEX_IMAGE: return not_likely_to_be_called(sp1,sp2); default: printf("atan: Illegal or unsupported argument type.\n"); longjmp(esc,1); } default: printf("atan: Illegal or unsupported argument type.\n"); longjmp(esc,1); } } sexpr *gilgalad() { printf("Gil-galad was an Elven king.\n"); printf("Of him the harpers sadly sing:\n"); printf("the last whose realm was fair and free\n"); printf("between the Mountains and the Sea.\n\n"); printf("His sword was long, his lance was keen,\n"); printf("his shining helm afar was seen;\n"); printf("the countless stars of heaven's field\n"); printf("were mirrored in his silver shield.\n\n"); printf("But long ago he rode away,\n"); printf("and where he dwelleth none can say;\n"); printf("for into darkness fell his star\n"); printf("in Mordor where the shadows are.\n"); return undefined; } sexpr *char2token(char c) { sexpr *sp = malloc(sizeof(sexpr)); sp->type = TOKEN; sp->u.c = c; return sp; } sexpr *char2exp(char c) { sexpr *sp = smalloc++; sp->type = CHARACTER; sp->u.c = c; return sp; } sexpr *num2exp(double x) { sexpr *sp = smalloc++; sp->type = NUMBER; sp->u.x = x; return sp; } sexpr *complex2exp(fcomplex z) { sexpr *sp = smalloc++; sp->type = COMPLEX; sp->u.z = z; return sp; } sexpr *num2str(sexpr *sp) { char *c; if (!number(sp)) { print_value_escape("number->string: Argument is not a real number: ",sp); } c = (char *) malloc(STRLEN); sprintf(c,"%g",sp->u.x); return str2exp(c); } sexpr *str2num(sexpr *sp) { float x; if (!string(sp)) { print_value_escape("string->number: Argument is not a string: ",sp); } sscanf(sp->u.text,"%g",&x); return num2exp(x); } sexpr *symbol2exp(sexpr *sp) { if (!symbol(sp)) { print_value_escape("symbol->string: Argument is not a symbol: ",sp); } return str2exp(sp->u.symbol->name); } sexpr *str2exp(char *text) { sexpr *sp = smalloc++; sp->type = STRING; sp->u.text = malloc((strlen(text) + 1)*sizeof(char)); strcpy(sp->u.text,text); return sp; } sexpr *exp2symbol(sexpr *sp) { int m, n = strlen(sp->u.text); if (!string(sp)) { print_value_escape("string->symbol: Argument not a string: ",sp); } m = 0; if (islegal1st(sp->u.text[0])) while (islegal2nd(sp->u.text[m++])); if (m < n) { char *name = malloc(STRLEN*sizeof(char)); sprintf(name,"|%s|",sp->u.text); return str2collectsymbol(name); } return str2collectsymbol(sp->u.text); } int length(sexpr *ls) { int x = 0; while (pair(ls)) { x++; ls = ls->u.pair->cdr; } if (null(ls)) return x; if (symbol(ls)) return -(x+1); return -1000; } sexpr *make_vector(int len) { sexpr *sp; sp = smalloc++; sp->type = VECTOR; sp->u.vector = malloc(sizeof(vector)); sp->u.vector->length = len; sp->u.vector->vector = malloc(len*sizeof(sexpr)); return sp; } sexpr *scheme_make_vector(sexpr *size, sexpr *init) { int i, len; sexpr *sp; if (!number(size)) { printf("make-vector: Illegal argument.\n"); longjmp(esc,1); } len = (int) size->u.x; sp = smalloc++; sp->type = VECTOR; sp->u.vector = malloc(sizeof(vector)); sp->u.vector->length = len; sp->u.vector->vector = malloc(len*sizeof(sexpr)); for (i=0; i < len; i++) sp->u.vector->vector[i]=init; return sp; } sexpr *vector_set(sexpr *sp, sexpr *i, sexpr *value) { if (!vector(sp) || !number(i)) { printf("vector-set!: Illegal argument.\n"); longjmp(esc,1); } if (i->u.x < 0 || i->u.x >= sp->u.vector->length) { printf("vector-set!: Index out of range: %g\n",i->u.x); longjmp(esc,1); } sp->u.vector->vector[(int) floor((float) i->u.x)]=value; return undefined; } sexpr *vector_ref(sexpr *sp, sexpr *i) { if (!vector(sp) || !number(i)) { printf("vector-ref: Illegal argument.\n"); longjmp(esc,1); } if (i->u.x < 0 || i->u.x >= sp->u.vector->length) { printf("vector-ref: Index out of range: %g\n",i->u.x); longjmp(esc,1); } return sp->u.vector->vector[(int) floor((float) i->u.x)]; } sexpr *vector_length(sexpr *sp) { if (!vector(sp)) { print_value_escape("vector-length: Argument is not a vector: ",sp); } return num2exp(sp->u.vector->length); } sexpr *string_eq(sexpr *sp1, sexpr *sp2) { if (!string(sp1)) { print_value_escape("string=?: Argument is not a string: ",sp1); } if (!string(sp2)) { print_value_escape("string=?: Argument is not a string: ",sp2); } return torf(!strcmp(sp1->u.text,sp2->u.text)); } sexpr *string_lt(sexpr *sp1, sexpr *sp2) { if (!string(sp1)) { print_value_escape("stringu.text,sp2->u.text) < 0); } sexpr *string_gt(sexpr *sp1, sexpr *sp2) { if (!string(sp1)) { print_value_escape("string>?: Argument is not a string: ",sp1); } if (!string(sp2)) { print_value_escape("string>?: Argument is not a string: ",sp2); } return torf(strcmp(sp1->u.text,sp2->u.text) > 0); } sexpr *string_leq(sexpr *sp1, sexpr *sp2) { if (!string(sp1)) { print_value_escape("string<=?: Argument is not a string: ",sp1); } if (!string(sp2)) { print_value_escape("string<=?: Argument is not a string: ",sp2); } return torf(strcmp(sp1->u.text,sp2->u.text) <= 0); } sexpr *string_geq(sexpr *sp1, sexpr *sp2) { if (!string(sp1)) { print_value_escape("string>=?: Argument is not a string: ",sp1); } if (!string(sp2)) { print_value_escape("string>=?: Argument is not a string: ",sp2); } return torf(strcmp(sp1->u.text,sp2->u.text) >= 0); } sexpr *char_eq(sexpr *sp1, sexpr *sp2) { if (!character(sp1)) { print_value_escape("char=?: Argument is not a character: ",sp1); } if (!character(sp2)) { print_value_escape("char=?: Argument is not a character: ",sp2); } return torf(sp1->u.c == sp2->u.c); } sexpr *char_lt(sexpr *sp1, sexpr *sp2) { if (!character(sp1)) { print_value_escape("charu.c < sp2->u.c); } sexpr *char_gt(sexpr *sp1, sexpr *sp2) { if (!character(sp1)) { print_value_escape("char>?: Argument is not a character: ",sp1); } if (!character(sp2)) { print_value_escape("char>?: Argument is not a character: ",sp2); } return torf(sp1->u.c > sp2->u.c); } sexpr *char_leq(sexpr *sp1, sexpr *sp2) { if (!character(sp1)) { print_value_escape("char<=?: Argument is not a character: ",sp1); } if (!character(sp2)) { print_value_escape("char<=?: Argument is not a character: ",sp2); } return torf(sp1->u.c <= sp2->u.c); } sexpr *char_geq(sexpr *sp1, sexpr *sp2) { if (!character(sp1)) { print_value_escape("char>=?: Argument is not a character: ",sp1); } if (!character(sp2)) { print_value_escape("char>=?: Argument is not a character: ",sp2); } return torf(sp1->u.c >= sp2->u.c); } sexpr *string_append(sexpr *sp1, sexpr *sp2) { char *c; if (!string(sp1)) { print_value_escape("string-append: Argument is not a string: ",sp1); } if (!string(sp2)) { print_value_escape("string-append: Argument is not a string: ",sp2); } c = (char *) malloc(strlen(sp1->u.text)+strlen(sp2->u.text)+1); strcpy(c,sp1->u.text); return str2exp(strcat(c,sp2->u.text)); } sexpr *string_ref(sexpr *text, sexpr *i) { if (!string(text)) { print_value_escape("string-ref: Argument is not a string: ",text); } if (!number(i)) { print_value_escape("string-ref: Argument is not a number: ",i); } if ((i->u.x < 0) || (i->u.x > strlen(text->u.text))) { printf("string-ref: Index out of range: %g",i->u.x); } return char2exp(text->u.text[(int) floor((float) i->u.x)]); } sexpr *substring(sexpr *text1, sexpr *i, sexpr *j) { int len, start, end, pos; char *text2; if (!string(text1)) { print_value_escape("substring: Argument is not a string: ",text1); } if (!number(i) || !number(j)) { printf("substring: Index is not a number.\n"); longjmp(esc,1); } len = strlen(text1->u.text); start = (int) floor((float) i->u.x); end = (int) floor((float) j->u.x); if ((start < 0) || (start > end) || (start > len)) { printf("substring: Start index out of range: %d", start); longjmp(esc,1); } if ((end < 0) || (end > len)) { printf("substring: End index out of range: %d", end); longjmp(esc,1); } text2 = (char *) malloc(sizeof(char)*(end-start)+1); for (pos=0; pos < (end-start); pos++) text2[pos]=text1->u.text[start+pos]; text2[end-start]= '\0'; return str2exp(text2); } sexpr *string_length(sexpr *text) { if (!string(text)) { print_value_escape("string-length: Argument is not a string: ",text); } return num2exp(strlen(text->u.text)); } sexpr *char2str(sexpr *sp) { char *text; if (!character(sp)) { print_value_escape("character->string: Argument is not a character: ",sp); } text = (char *) malloc(2); text[0] = sp->u.c; text[1] = '\0'; return str2exp(text); } int notquote(char c) { return (c != '\"'); } int notstile(char c) { return (c != '|'); } sexpr *lex(FILE *file, int stop) { char c, d, *name; sexpr *sp; double x, y; c = getc(file); if (c == stop || c == EOF) return nil; else if (isspace(c)) { if (stop == ' ') return nil; else return lex(file,stop); } else if (c=='(') return cons(token$left,lex(file,stop)); else if (c==')') return cons(token$right,lex(file,stop)); else if (c=='\'') return cons(token$quote,lex(file,stop)); else if (c=='`') return cons(token$quasiquote,lex(file,stop)); else if (c==',') { d = getc(file); if (d=='@') return cons(token$unquote_splicing,lex(file,stop)); else { ungetc(d,file); return cons(token$unquote,lex(file,stop)); } } else if (c=='#') { d = getc(file); if (d=='!') { eatline(file, c); return lex(file,stop); } else if (d=='t') return cons(true,lex(file,stop)); else if (d=='f') return cons(false,lex(file,stop)); else if (d=='\\') return backslash(file,stop); else if (d=='(') { ungetc(d,file); return cons(token$pound_sign,lex(file,stop)); } else { printf("Illegal character, '%c', following '#'.\n",d); eatline(file, d); longjmp(esc,1); } } else if (c==';') { eatline(file, c); return lex(file,stop); } else if (isdigit(c)) { ungetc(c,file); x = eat_number(file,stop); d = getc(file); switch(d) { case 'i': return cons(complex2exp(Complex(0.0,x)),lex(file,stop)); case '@': y = eat_number(file,stop); return cons(polar2complex(num2exp(x),num2exp(y)),lex(file,stop)); case '+': y = eat_number(file,stop); if (getc(file) == 'i') return cons(complex2exp(Complex(x,y)),lex(file,stop)); else { printf("Malformed number.\n"); longjmp(esc,1); } case '-': y = eat_number(file,stop); if (getc(file) == 'i') return cons(complex2exp(Complex(x,-y)),lex(file,stop)); else { printf("Malformed number.\n"); longjmp(esc,1); } default: ungetc(d,file); return cons(num2exp(x),lex(file,stop)); } } else if (c=='\"') { name = malloc(STRLEN*sizeof(char)); eatit(notquote,name,file,'\"'); d = getc(file); sp = str2exp(name); free(name); return cons(sp,lex(file,stop)); } else if (c=='|') { char *between=malloc(STRLEN*sizeof(char)); name = malloc(STRLEN*sizeof(char)); eatit(notstile,between,file,'|'); d = getc(file); sprintf(name,"|%s|",between); sp = str2collectsymbol(name); free(between); free(name); return cons(sp,lex(file,stop)); } else if (c=='.') { d = getc(file); if (!isdigit(d)) { ungetc(d,file); return cons(token$dot,lex(file,stop)); } else { ungetc(d,file); ungetc(c,file); x = eat_number(file,stop); d = getc(file); switch(d) { case 'i': return cons(complex2exp(Complex(0.0,x)),lex(file,stop)); case '+': y = eat_number(file,stop); if (getc(file) == 'i') return cons(complex2exp(Complex(x,y)),lex(file,stop)); else { printf("Malformed number.\n"); longjmp(esc,1); } case '-': y = eat_number(file,stop); if (getc(file) == 'i') return cons(complex2exp(Complex(x,-y)),lex(file,stop)); else { printf("Malformed number.\n"); longjmp(esc,1); } default: ungetc(d,file); return cons(num2exp(x),lex(file,stop)); } } } else if (c=='-') { d=getc(file); if (isdigit_or_pt(d)) { ungetc(d,file); x = eat_number(file,stop); d = getc(file); switch(d) { case 'i': return cons(complex2exp(Complex(0.0,-x)),lex(file,stop)); case '+': y = eat_number(file,stop); if (getc(file) == 'i') return cons(complex2exp(Complex(-x,y)),lex(file,stop)); else { printf("Malformed number.\n"); longjmp(esc,1); } case '-': y = eat_number(file,stop); if (getc(file) == 'i') return cons(complex2exp(Complex(-x,-y)),lex(file,stop)); else { printf("Malformed number.\n"); longjmp(esc,1); } default: ungetc(d,file); return cons(num2exp(-x),lex(file,stop)); } } else if (islegal1st(d)) { ungetc(d,file); ungetc(c,file); name = malloc(STRLEN*sizeof(char)); eatit(islegal2nd,name,file,stop); sp = str2collectsymbol(name); free(name); return cons(sp,lex(file,stop)); } else { ungetc(d,file); return cons(symbol$hyphen,lex(file,stop)); } } else if (islegal1st(c)) { ungetc(c,file); name = malloc(STRLEN*sizeof(char)); eatit(islegal2nd,name,file,stop); sp = str2collectsymbol(name); free(name); return cons(sp,lex(file,stop)); } else longjmp(esc,1); } sexpr *read_newline(FILE *file, int stop) { char c2, c3, c4, c5, c6, c7; c2 = getc(file); if (c2 != 'e') { ungetc(c2,file); return(cons(char2exp('n'),lex(file,stop))); } else { c3 = getc(file); if (c3 != 'w') return lex_error(file, c3); else { c4 = getc(file); if (c4 != 'l') return lex_error(file, c4); else { c5 = getc(file); if (c5 != 'i') return lex_error(file, c5); else { c6 = getc(file); if (c6 != 'n') return lex_error(file, c6); else { c7 = getc(file); if (c7 == 'e') return(cons(char2exp('\n'),lex(file,stop))); else return lex_error(file, c7); } } } } } } sexpr *read_space(FILE *file, int stop) { char c2, c3, c4, c5; c2 = getc(file); if (c2 != 'p') { ungetc(c2,file); return(cons(char2exp('s'),lex(file,stop))); } else { c3 = getc(file); if (c3 != 'a') return lex_error(file, c3); else { c4 = getc(file); if (c4 != 'c') return lex_error(file, c4); else { c5 = getc(file); if (c5 == 'e') return(cons(char2exp(' '),lex(file,stop))); else return lex_error(file, c5); } } } } sexpr *read_tab(FILE *file, int stop) { char c2, c3; c2 = getc(file); if (c2 != 'a') { ungetc(c2,file); return(cons(char2exp('t'),lex(file,stop))); } else { c3 = getc(file); if (c3 == 'b') return(cons(char2exp('\t'),lex(file,stop))); else return lex_error(file, c3); } } sexpr *backslash(FILE *file, int stop) { char c1; c1 = getc(file); switch (c1) { case 'n': return read_newline(file,stop); break; case 's': return read_space(file,stop); break; case 't': return read_tab(file,stop); break; default: if (c1 != stop) return(cons(char2exp(c1),lex(file,stop))); else return lex_error(file, c1); } } sexpr *list2(sexpr *s1, sexpr *s2) { return cons(s1,cons(s2,nil)); } sexpr *ploop0(sexpr *ls, FILE *file, int stop) { sexpr *foo; if (null(ls)) return ls; foo = ploop1(ls,nil,0,file,stop); return cons(car(foo),ploop0(cdr(foo),file,stop)); } sexpr *ploop1(sexpr *ls, sexpr *acc, int d, FILE *file, int stop) { sexpr *sp; if (d==0 && !null(acc)) return cons(acc,ls); if (null(ls)) return ploop1(lex(file,stop),acc,d,file,stop); sp = car(ls); if (eq(sp,token$dot) || symbol(sp) || number(sp) || character(sp) || string(sp) || boolean(sp) || komplex(sp)) { return ploop1(cdr(ls),cons(sp,acc),d,file,stop); } if (eq(sp,token$left)) return ploop1(cdr(ls),cons(sp,acc),d+1,file,stop); if (eq(sp,token$right)) return ploop1(cdr(ls),ploop2(acc,nil,file,stop),d-1,file,stop); if (eq(sp,token$quote)) { sexpr *foo = ploop1(cdr(ls),nil,0,file,stop); return ploop1(cdr(foo),cons(list2(symbol$quote,caar(foo)),acc),d,file,stop); } if (eq(sp,token$quasiquote)) { sexpr *foo = ploop1(cdr(ls),nil,0,file,stop); return ploop1(cdr(foo),cons(list2(symbol$quasiquote,caar(foo)),acc),d,file,stop); } if (eq(sp,token$unquote)) { sexpr *foo = ploop1(cdr(ls),nil,0,file,stop); return ploop1(cdr(foo),cons(list2(symbol$unquote,caar(foo)),acc),d,file,stop); } if (eq(sp,token$unquote_splicing)) { sexpr *foo = ploop1(cdr(ls),nil,0,file,stop); return ploop1(cdr(foo),cons(list2(symbol$unquote_splicing,caar(foo)),acc),d,file,stop); } if (eq(sp,token$pound_sign)) { sexpr *foo = ploop1(cdr(ls),nil,0,file,stop); return ploop1(cdr(foo),cons(list2vector(caar(foo)),acc),d,file,stop); } printf("Parsing error.\n"); longjmp(esc,1); } sexpr *ploop2(sexpr *ls, sexpr *acc, FILE *file, int stop) { sexpr *sp; if (null(ls)) { printf("Parsing error.\n"); longjmp(esc,1); } sp = car(ls); if (eq(sp,token$left)) return cons(acc,cdr(ls)); if (eq(sp,token$dot)) { if (null(cdr(acc))) return ploop2(cdr(ls),car(acc),file,stop); printf("Parsing error.\n"); longjmp(esc,1); } return ploop2(cdr(ls),cons(sp,acc),file,stop); } sexpr *parse(FILE *file, int stop) { sexpr *s1, *s2; s1 = lex(file,stop); s2 = ploop0(s1,file,stop); return map(car,s2); } sexpr *scheme_exit() { exit(0); } sexpr *gensym() { char *name = malloc(sizeof(char)*10); sprintf(name,"g%d",gensyms); gensyms++; return(str2collectsymbol(name)); } sexpr *nullp(sexpr *sp) { return torf(null(sp)); } sexpr *eof_objectp(sexpr *sp) { return torf(eof_object(sp)); } sexpr *pairp(sexpr *sp) { return torf(pair(sp)); } sexpr *numberp(sexpr *sp) { return torf(number(sp)); } sexpr *characterp(sexpr *sp) { return torf(character(sp)); } sexpr *stringp(sexpr *sp) { return torf(string(sp)); } sexpr *symbolp(sexpr *sp) { return torf(symbol(sp)); } sexpr *procedurep(sexpr *sp) { return torf(virgin(sp) || closure(sp) || primitive(sp)); } sexpr *primitivep(sexpr *sp) { return torf(primitive(sp)); } sexpr *booleanp(sexpr *sp) { return torf(boolean(sp)); } sexpr *vectorp(sexpr *sp) { return torf(vector(sp)); } sexpr *input_portp(sexpr *sp) { return torf(input_port(sp)); } sexpr *output_portp(sexpr *sp) { return torf(output_port(sp)); } sexpr *voidp(sexpr *sp) { return torf(undefined(sp)); } sexpr *graphicp(sexpr *sp) { return torf(graphic(sp)); } sexpr *imagep(sexpr *sp) { return torf(image(sp)); } sexpr *complex_imagep(sexpr *sp) { return torf(complex_image(sp)); } sexpr *color_imagep(sexpr *sp) { return torf(color_image(sp)); } sexpr *linep(sexpr *sp) { return torf(line(sp)); } sexpr *sketchp(sexpr *sp) { return torf(sketch(sp)); } sexpr *undefinedp(sexpr *sp) { return torf(undefined(sp)); } sexpr *complexp(sexpr *sp) { return torf(komplex(sp)); } sexpr *realp(sexpr *sp) { return torf(number(sp) || (komplex(sp) && (sp->u.z.i == 0))); } int list(sexpr *sp) { if (null(sp)) return 1; else if (pair(sp)) return list(cdr(sp)); else return 0; } sexpr *display(sexpr *sp, sexpr *port) { if (undefined(sp)) return sp; if (!output_port(port)) { print_value_escape("display: Argument is not an output-port:",port); } display_or_print(sp,DISPLAY,TOPLEVEL,port->u.file); return undefined; } sexpr *print(sexpr *sp) { if (undefined(sp)) return sp; print2stdout(sp); printf("\n"); return undefined; } void pair_print(sexpr *sp, int print_flag, FILE *file) { if (null(sp)) fprintf(file,"()"); else if (!pair(sp)) display_or_print(sp,print_flag,EMBEDDED,file); else { fprintf(file,"("); pair_print(car(sp),print_flag,file); pair_cdr_print(cdr(sp),print_flag,file); fprintf(file,")"); } } void pair_cdr_print(sexpr *sp, int print_flag, FILE *file) { if (!null(sp)) { if (pair(sp)) { sexpr *rest = cdr(sp); fprintf(file," "); pair_print(car(sp),print_flag,file); if (!null(rest)) { if (pair(rest)) pair_cdr_print(rest,print_flag,file); else { fprintf(file," . "); pair_print(rest,print_flag,file); } } } else { fprintf(file," . "); display_or_print(sp,print_flag,EMBEDDED,file); } } } void vector_print(sexpr *sp, int print_flag, FILE *file) { int i, l; if (sp == global_vals) { fprintf(file,"#"); } else { l = sp->u.vector->length; fprintf(file,"#("); for (i = 0; i < l; i++) { display_or_print(sp->u.vector->vector[i],print_flag,EMBEDDED,file); if (i < l-1) fprintf(file," "); } fprintf(file,")"); } } sexpr *newline(sexpr *port) { if (!output_port(port)) { print_value_escape("newline: Argument is not an output-port: ",port); } fprintf(port->u.file,"\n"); return undefined; } sexpr *scheme_void() { return undefined; } sexpr *scheme_error() { longjmp(esc,1); return nil; } sexpr *scheme_system(sexpr *sp0) { if (!string(sp0)) { print_value_escape("system: Argument is not a string: ",sp0); } return(num2exp(system(sp0->u.text))); } sexpr *assoc(sexpr *symbol, sexpr *frame) { int i, n = frame->u.vector->length; sexpr **data = frame->u.vector->vector; for (i = 0; i < n; i++) if (symbol == data[i]) return num2exp(i); return undefined; } sexpr *add_new_global_def(sexpr *symbol, sexpr *value) { int n = global_vars->u.vector->length; global_vars->u.vector->vector[n] = symbol; global_vals->u.vector->vector[n] = value; global_vars->u.vector->length = n+1; global_vals->u.vector->length = n+1; return num2exp(n); } sexpr *add_def_compile(sexpr *symbol, sexpr *env) { sexpr *frame = env->u.pair->car; int i, n = frame->u.vector->length; sexpr **data = frame->u.vector->vector; sexpr **extended_data; for (i = 0; i < n; i++) if (symbol == data[i]) return num2exp(i); extended_data = malloc((n+1)*sizeof(sexpr)); frame->u.vector->length = n+1; for (i = 0; i < n; i++) extended_data[i]=data[i]; extended_data[n] = symbol; frame->u.vector->vector = extended_data; free(data); return num2exp(n); } void add_def(int i, sexpr *value, sexpr *env) { sexpr *frame = env->u.pair->car; int n = frame->u.vector->length; sexpr **data = frame->u.vector->vector; sexpr **extended_data; if (i < n) { data[i] = value; return; } else { extended_data = malloc((n+1)*sizeof(sexpr)); frame->u.vector->length = n+1; for (i = 0; i < n; i++) extended_data[i]=data[i]; extended_data[n] = value; frame->u.vector->vector = extended_data; free(data); } } void add_def_macro(sexpr *symbol, sexpr *value, sexpr *env) { sexpr *frame = env->u.pair->car; int i, n = frame->u.vector->length; sexpr **data = frame->u.vector->vector; sexpr **extended_data; frame->u.vector->length = n+2; for (i = 0; i < n; i += 2) { if (symbol==data[i]) { data[i+1] = value; return; } } extended_data = malloc((n+2)*sizeof(sexpr)); frame->u.vector->length = n+2; for (i = 0; i < n; i++) extended_data[i]=data[i]; extended_data[n] = symbol; extended_data[n+1] = value; frame->u.vector->vector = extended_data; free(data); } void add_macro(sexpr *symbol, sexpr *expander) { sexpr *code = expander->u.closure->compiled; int i, n = code->u.vector->length; sexpr **data = code->u.vector->vector; sexpr **extended_data = malloc((n+2)*sizeof(sexpr)); code->u.vector->length = n+2; extended_data[0] = bytecode$frame; extended_data[1] = make_vector(1); extended_data[1]->u.vector->vector[0] = bytecode$halt; for (i=0; i < n; i++) extended_data[2+i] = data[i]; code->u.vector->vector = extended_data; free(data); symbol->u.symbol->macro = expander->u.closure; /* so garbage collector can see expanders */ add_def_macro(symbol,expander,macros); } sexpr *list2frame(int n, sexpr *vals) { int i; sexpr *frame, **data; /* (lambda args ... ) => -1 */ if (n == -1) { frame = make_vector(1); data = frame->u.vector->vector; data[0] = vals; return frame; } /* (lambda (x y . z) ... ) => -3 */ if (n < 0) { n = fabs(n); frame = make_vector(n); data = frame->u.vector->vector; for (i = 0; i < n-1; i++) { data[i] = car(vals); vals = vals->u.pair->cdr; } data[n-1] = vals; return frame; } frame = make_vector(n); data = frame->u.vector->vector; /* (lambda (x y) ... ) => 2 */ for (i = 0; i < n; i++) { data[i] = car(vals); vals = vals->u.pair->cdr; } return frame; } sexpr *frame4apply(int n, int r0, int r, sexpr **args) { int i, m; sexpr *frame, **data, *ls; m = r-r0; /* (lambda (x y . z) ... ) => -3 */ if (n < 0) { n = abs(n); frame = make_vector(n); data = frame->u.vector->vector; /* extra args become named args */ if (m <= n) { for (i = 0; i < m-1; i++) data[i] = args[r-i-1]; ls = args[r-m]; for (i = m-1; i < n-1; i++) { data[i] = car(ls); ls = cdr(ls); } data[n-1] = ls; } else { /* named args become extra args */ for (i = 0; i < n-1; i++) data[i] = args[r-i-1]; ls = args[r-m]; for (i = m-1; i > n-1; i--) ls = cons(args[r-i],ls); data[n-1] = ls; } return frame; } if (m > n) { printf("apply: Wrong number of arguments for user-defined function.\n"); longjmp(esc,1); } /* (lambda (x y) ... ) => 2 */ frame = make_vector(n); data = frame->u.vector->vector; /* extra args become named args */ for (i = 0; i < m-1; i++) data[i] = args[r-i-1]; ls = args[r-m]; for (i = m-1; i < n; i++) { data[i] = car(ls); ls = cdr(ls); } return frame; } sexpr *frame4funcall(int n, int r0, int r, sexpr **args) { int i; sexpr *frame, **data, *ls; /* (lambda args ... ) => -1 */ if (n == -1) { frame = make_vector(1); data = frame->u.vector->vector; ls = nil; for (i = r0; i < r; i++) ls = cons(args[i],ls); data[0] = ls; return frame; } /* (lambda (x y . z) ... ) => -3 */ if (n < 0) { n = abs(n); frame = make_vector(n); data = frame->u.vector->vector; if (r-r0 < n-1) { printf("funcall: Wrong number of arguments for user-defined function.\n"); longjmp(esc,1); } for (i = 0; i < n-1; i++) data[i] = args[r-i-1]; ls = nil; for (i = r0; i < r-n+1; i++) ls = cons(args[i],ls); data[n-1] = ls; return frame; } if (r-r0 != n) { printf("funcall: Wrong number of arguments for user-defined function.\n"); longjmp(esc,1); } /* (lambda (x y) ... ) => 2 */ frame = make_vector(n); data = frame->u.vector->vector; for (i = 0; i < n; i++) data[i] = args[r-i-1]; return frame; } sexpr *reduce_pred(sexpr *pred, sexpr *ls) { sexpr *rest; sexpr *(*func)() = pred->u.primitive->func; if (null(ls)) longjmp(esc,1); rest = cdr(ls); if (null(rest)) return true; if (!isfalse(func(car(ls),car(rest)))) return reduce_pred(pred,rest); return false; } sexpr *reduce(sexpr *proc, sexpr *id, sexpr *ls) { sexpr *rest; sexpr *(*func)() = proc->u.primitive->func; if (null(ls)) return id; rest = cdr(ls); if (null(rest)) return func(car(ls),id); return func(car(ls),reduce(proc,id,rest)); } sexpr *map(sexpr *(*proc)(sexpr *), sexpr *ls) { if (null(ls)) return ls; return cons(proc(car(ls)),map(proc,cdr(ls))); } sexpr *scheme_map(sexpr *proc, sexpr *ls) { return map(proc->u.primitive->func,ls); } sexpr *expand_let(sexpr *sp) { sexpr *vars, *vals, *body; sexpr *first, *second; first = cdr(sp); second = car(first); vars = map(car,second); vals = map(cadr,second); body = cdr(first); return cons(cons(symbol$lambda,cons(vars,body)),vals); } sexpr *append(sexpr *ls1, sexpr *ls2) { if (null(ls1)) return ls2; return cons(car(ls1),append(cdr(ls1),ls2)); } sexpr *reverse(sexpr *ls) { if (!pair(ls) && !null(ls)) { print_value_escape("reverse: Argument is not a list: ",ls); } if (null(ls)) return ls; return append(reverse(cdr(ls)),cons(car(ls),nil)); } sexpr *expand_let_star(sexpr *sp) { sexpr *vars, *vals, *body; sexpr *rest, *bindings, *first_binding, *rest_bindings; rest = cdr(sp); bindings = car(rest); first_binding = car(bindings); rest_bindings = cdr(bindings); if (null(rest_bindings)) { vars = cons(car(first_binding),nil); vals = cdr(first_binding); body = cdr(rest); return cons(cons(symbol$lambda,cons(vars,body)),vals); } else { vars = cons(car(first_binding),nil); vals = cons(cadr(first_binding),nil); body = cons(symbol$let_star,cons(rest_bindings,cdr(rest))); return cons(cons(symbol$lambda,cons(vars,cons(body,nil))),vals); } } sexpr *make_set_list(sexpr *vars, sexpr *vals) { if (null(vars)) return nil; return cons(cons(symbol$set_bang,cons(car(vars),cons(car(vals),nil))),make_set_list(cdr(vars),cdr(vals))); } sexpr *do_nothing(sexpr *sp) { return nil; } sexpr *expand_letrec(sexpr *sp) { sexpr *vars, *vals, *body, *setbangs, *dummy; sexpr *first, *second; first=cdr(sp); second=car(first); vars = map(car,second); vals = map(cadr,second); dummy = map(do_nothing,second); body = first->u.pair->cdr; setbangs = make_set_list(vars,vals); return cons(cons(symbol$lambda,cons(vars,append(setbangs,body))),dummy); } sexpr *not(sexpr *sp) { return torf(isfalse(sp)); } continuation *make_continuation(int s, int r) { continuation *kp = (continuation *) malloc(sizeof(continuation)); kp->args = malloc(r*sizeof(sexpr *)); kp->cv_stack = malloc(s*sizeof(sexpr *)); kp->e_stack = malloc(s*sizeof(sexpr *)); kp->lu_stack = malloc(s*sizeof(sexpr *)); kp->pc_stack = malloc(s*sizeof(sexpr **)); kp->r_stack = malloc(s*sizeof(int)); return kp; } sexpr *virtual_machine(sexpr *a, sexpr *cv, sexpr *e) { sexpr **global_data = global_vals->u.vector->vector; int r = 0, s = 0; sexpr *args[STACKSIZE]; sexpr *lu_stack[STACKSIZE]; sexpr *cv_stack[STACKSIZE]; sexpr **pc_stack[STACKSIZE]; sexpr *e_stack[STACKSIZE]; int r_stack[STACKSIZE]; int i, j, s0; sexpr *sp, *lu = nil; closure *cp; continuation *kp; sexpr **data; sexpr *(*func)(); sexpr **pc = cv->u.vector->vector; virtual_machine: switch((*pc)->type) { case HALT: return a; case REFER: ++pc; j = (int) (*pc)->u.x; ++pc; i = (int) (*pc)->u.x; a = lu->u.vector->vector[j-1]->u.vector->vector[i]; ++pc; goto virtual_machine; case LOCAL_REFER: ++pc; i = (int) (*pc)->u.x; a = e->u.pair->car->u.vector->vector[i]; ++pc; goto virtual_machine; case GLOBAL_REFER: ++pc; i = (int) (*pc)->u.x; a = global_data[i]; ++pc; goto virtual_machine; case CONSTANT: ++pc; a = *pc; ++pc; goto virtual_machine; case ARGUMENT: args[r++] = a; ++pc; goto virtual_machine; case TEST: if ((a->type == BOOLEAN) && (!a->u.i)) { ++pc; ++pc; } else { ++pc; cv = *pc; pc = cv->u.vector->vector; } goto virtual_machine; case APPLY: switch (a->type) { case PRIMITIVE: func = a->u.primitive->func; i = a->u.primitive->arity; sp = frame4apply(i,r_stack[s-1],r,args); data = sp->u.vector->vector; switch(i) { case 0: a = func(); break; case 1: a = func(data[0]); break; case 2: a = func(data[0],data[1]); break; case 3: a = func(data[0],data[1],data[2]); break; case 4: a = func(data[0],data[1],data[2],data[3]); break; case 5: a = func(data[0],data[1],data[2],data[3],data[4]); break; case 6: a = func(data[0],data[1],data[2],data[3],data[4],data[5]); break; case 7: a = func(data[0],data[1],data[2],data[3],data[4],data[5],data[6]); break; case 8: a = func(data[0],data[1],data[2],data[3],data[4],data[5],data[6],data[7]); break; case 9: a = func(data[0],data[1],data[2],data[3],data[4],data[5],data[6],data[7],data[8]); } s--; cv = cv_stack[s]; pc = pc_stack[s]; lu = lu_stack[s]; e = e_stack[s]; r = r_stack[s]; goto virtual_machine; case VIRGIN: cv = a->u.closure->compiled; pc = cv->u.vector->vector; lu = a->u.closure->lookup = list2vector(a->u.closure->env); sp = frame4apply(a->u.closure->n,r_stack[s-1],r,args); e = cons(sp,a->u.closure->env); r = r_stack[s-1]; goto virtual_machine; case CLOSURE: cv = a->u.closure->compiled; pc = cv->u.vector->vector; lu = a->u.closure->lookup; sp = frame4apply(a->u.closure->n,r_stack[s-1],r,args); e = cons(sp,a->u.closure->env); r = r_stack[s-1]; goto virtual_machine; default: print_value_escape("apply: Attempt to apply non-function: ",a); } case FUNCALL: switch (a->type) { case PRIMITIVE: func = a->u.primitive->func; i = a->u.primitive->arity; if (i != r-r_stack[s-1]) { print_value_escape("apply: Wrong number of arguments for primitive function: ",a); } switch(i) { case 0: a = func(); break; case 1: a = func(args[r-1]); break; case 2: a = func(args[r-1],args[r-2]); break; case 3: a = func(args[r-1],args[r-2],args[r-3]); break; case 4: a = func(args[r-1],args[r-2],args[r-3],args[r-4]); break; case 5: a = func(args[r-1],args[r-2],args[r-3],args[r-4],args[r-5]); break; case 6: a = func(args[r-1],args[r-2],args[r-3],args[r-4],args[r-5],args[r-6]); break; case 7: a = func(args[r-1],args[r-2],args[r-3],args[r-4],args[r-5],args[r-6],args[r-7]); break; case 8: a = func(args[r-1],args[r-2],args[r-3],args[r-4],args[r-5],args[r-6],args[r-7],args[r-8]); break; case 9: a = func(args[r-1],args[r-2],args[r-3],args[r-4],args[r-5],args[r-6],args[r-7],args[r-8],args[r-9]); } s--; cv = cv_stack[s]; pc = pc_stack[s]; lu = lu_stack[s]; e = e_stack[s]; r = r_stack[s]; goto virtual_machine; case VIRGIN: cv = a->u.closure->compiled; pc = cv->u.vector->vector; lu = a->u.closure->lookup = list2vector(a->u.closure->env); sp = frame4funcall(a->u.closure->n,r_stack[s-1],r,args); e = cons(sp,a->u.closure->env); r = r_stack[s-1]; goto virtual_machine; case CLOSURE: cv = a->u.closure->compiled; pc = cv->u.vector->vector; lu = a->u.closure->lookup; sp = frame4funcall(a->u.closure->n,r_stack[s-1],r,args); e = cons(sp,a->u.closure->env); r = r_stack[s-1]; goto virtual_machine; default: print_value_escape("apply: Attempt to apply non-function: ",a); } case FRAME: ++pc; cv = *pc; cv_stack[s] = cv; pc_stack[s] = cv->u.vector->vector; lu_stack[s] = lu; e_stack[s] = e; r_stack[s] = r; s++; ++pc; goto virtual_machine; case RETURN: s--; cv = cv_stack[s]; pc = pc_stack[s]; lu = lu_stack[s]; e = e_stack[s]; r = r_stack[s]; goto virtual_machine; case CLOSE: cp = (closure *) malloc(sizeof(closure)); sp = smalloc++; sp->type = VIRGIN; sp->u.closure = cp; a = sp; ++pc; cp->n = (int) (*pc)->u.x; ++pc; cp->lookup = nil; cp->compiled = *pc; cp->env = e; ++pc; goto virtual_machine; case ASSIGN: ++pc; j = (int) (*pc)->u.x; ++pc; i = (int) (*pc)->u.x; lu->u.vector->vector[j-1]->u.vector->vector[i] = a; a = undefined; ++pc; goto virtual_machine; case LOCAL_ASSIGN: ++pc; i = (int) (*pc)->u.x; e->u.pair->car->u.vector->vector[i] = a; a = undefined; ++pc; goto virtual_machine; case GLOBAL_ASSIGN: ++pc; i = (int) (*pc)->u.x; global_data[i] = a; a = undefined; ++pc; goto virtual_machine; case CONTI: cp = (closure *) malloc(sizeof(closure)); sp = smalloc++; sp->type = CLOSURE; sp->u.closure = cp; cp->n = 1; cp->compiled = make_vector(2); cp->compiled->u.vector->vector[0] = bytecode$nuate; a = sp; sp = smalloc++; sp->type = CONTINUATION; cp->compiled->u.vector->vector[1] = sp; cp->env = global_env_vals; cp->lookup = nil; kp = make_continuation(s,r); sp->u.continuation = kp; for (i = 0; i < r; i++) kp->args[i] = args[i]; for (i = 0; i < s; i++) { kp->cv_stack[i] = cv_stack[i]; kp->pc_stack[i] = pc_stack[i]; kp->lu_stack[i] = lu_stack[i]; kp->e_stack[i] = e_stack[i]; kp->r_stack[i] = r_stack[i]; } kp->s = s; ++pc; goto virtual_machine; case NUATE: ++pc; kp = (*pc)->u.continuation; s = kp->s; for (i = 0; i < s; i++) { cv_stack[i] = kp->cv_stack[i]; pc_stack[i] = kp->pc_stack[i]; lu_stack[i] = kp->lu_stack[i]; e_stack[i] = kp->e_stack[i]; r_stack[i] = kp->r_stack[i]; } a = e->u.pair->car->u.vector->vector[0]; s--; cv = cv_stack[s]; pc = pc_stack[s]; lu = lu_stack[s]; e = e_stack[s]; r = r_stack[s]; for (i = 0; i < r; i++) args[i] = kp->args[i]; goto virtual_machine; case THROW: cp = (closure *) malloc(sizeof(closure)); sp = smalloc++; sp->type = CLOSURE; sp->u.closure = cp; cp->n = 1; cp->compiled = make_vector(2); cp->compiled->u.vector->vector[0] = bytecode$catch; cp->compiled->u.vector->vector[1] = num2exp(s); cp->env = global_env_vals; cp->lookup = nil; a = sp; ++pc; goto virtual_machine; case CATCH: ++pc; s0 = (*pc)->u.i; if (s0 > s) { printf("Attempt to return to escape continuation which no longer exists.\n"); longjmp(esc,1); } else s = s0; a = e->u.pair->car->u.vector->vector[0]; cv = cv_stack[s]; pc = pc_stack[s]; e = e_stack[s]; r = r_stack[s]; goto virtual_machine; case DEFINE: ++pc; add_def((int) (*pc)->u.x,a,e); a = undefined; ++pc; goto virtual_machine; default: print_value_escape("Unrecognized virtual machine instruction: ",*pc); } } sexpr *eval(sexpr *x) { return virtual_machine(nil,cps2vector(car(compile(expand(x),cons(bytecode$halt,nil),global_env_vars))),global_env_vals); } sexpr *scheme_compile(sexpr *x) { return cps2vector(car(compile(expand(x),cons(bytecode$halt,nil),global_env_vars))); } sexpr *expand_quasiquote_list(sexpr *sp) { if (!pair(sp)) return cons(symbol$quote,cons(cons(sp,nil),nil)); else { sexpr *first = sp->u.pair->car; if (symbol(first)) { if (first==symbol$unquote) return cons(primitive$cons,cons(cadr(sp),cons(nil,nil))); else if (first==symbol$unquote_splicing) return cadr(sp); else if (first==symbol$quasiquote) return expand_quasiquote_list(expand_quasiquote(cadr(sp))); } return cons(primitive$cons,cons(cons(primitive$append,cons(expand_quasiquote_list(car(sp)),cons(expand_quasiquote(cdr(sp)),nil))),cons(nil,nil))); } } sexpr *expand_quasiquote(sexpr *sp) { if (!pair(sp)) return cons(symbol$quote,cons(sp,nil)); else { sexpr *first = sp->u.pair->car; if (symbol(first)) { if (first==symbol$unquote) return cadr(sp); else if (first==symbol$unquote_splicing) { printf("expand: unquote-splicing not inside list.\n"); longjmp(esc,1); } else if (first==symbol$quasiquote) return expand_quasiquote(expand_quasiquote(cadr(sp))); } return cons(primitive$append,cons(expand_quasiquote_list(car(sp)),cons(expand_quasiquote(cdr(sp)),nil))); } } sexpr *expand(sexpr *sp) { sexpr *first; if (self_evaluating(sp) || symbol(sp)) return sp; first = car(sp); if (symbol(first)) { if (first==symbol$quote) return sp; if (first==symbol$quasiquote) return expand(expand_quasiquote(cadr(sp))); if (first==symbol$let) return expand(expand_let(sp)); if (first==symbol$let_star) return expand(expand_let_star(sp)); if (first==symbol$letrec) return expand(expand_letrec(sp)); if (first==symbol$lambda) return cons(symbol$lambda,cons(cadr(sp),map(expand,cddr(sp)))); if (first->u.symbol->macro != NULL) { sexpr *cv = first->u.symbol->macro->compiled; sexpr *frame = list2frame(first->u.symbol->macro->n,cdr(sp)); sexpr *env = cons(frame,first->u.symbol->macro->env); sp = virtual_machine(nil,cv,env); return expand(sp); } } return cons(expand(car(sp)),map(expand,cdr(sp))); } sexpr *gtype(sexpr *g) { if (!graphic(g)) { print_value_escape("graphic-type: Argument is not a graphic object: ",g); } if (pair(g->u.sexpr)) return car(g->u.sexpr); else return symbol$nil; } sexpr *gargs(sexpr *g) { if (!graphic(g)) { print_value_escape("graphic-args: Argument is not a graphic object: ",g); } return cdr(g->u.sexpr); } sexpr *gnull(sexpr *g) { if (!graphic(g)) { print_value_escape("graphic-null?: Argument is not a graphic object: ",g); } return torf(null(g->u.sexpr)); } sexpr *gcar(sexpr *g1) { sexpr *g2; if (!graphic(g1)) { print_value_escape("gcar: Argument is not a graphic object: ",g1); } g2 = gtype(g1); if (eq(g2,symbol$adjoin) || eq(g2,symbol$adorn)) { g2 = smalloc++; g2->type = GRAPHIC; g2->u.sexpr = cadr(g1->u.sexpr); return g2; } print_value_escape("gcar: Argument is not an adjoin or adorn: ",g2); } sexpr *gcdr(sexpr *g1) { sexpr *g2; sexpr *rest; if (!graphic(g1)) { print_value_escape("gcdr: Argument is not a graphic object: ",g1); } g2 = gtype(g1); if (eq(g2,symbol$adjoin) || eq(g2,symbol$adorn)) { rest = cddr(g1->u.sexpr); g2 = smalloc++; g2->type = GRAPHIC; if (null(rest)) g2->u.sexpr = nil; else g2->u.sexpr = cons(car(g1->u.sexpr),rest); return g2; } print_value_escape("gcdr: Argument is not an adjoin or adorn: ",g2); } sexpr *adjoin(sexpr *g1, sexpr *g2) { sexpr *g3; if (!graphic(g1)) { print_value_escape("adjoin: Argument is not a graphic object: ",g1); } if (!graphic(g2)) { print_value_escape("adjoin: Argument is not a graphic object: ",g2); } if (null(g1->u.sexpr)) return g2; if (null(g2->u.sexpr)) return g1; g3 = smalloc++; g3->type = GRAPHIC; g3->u.sexpr = cons(symbol$adjoin,cons(g1->u.sexpr,cons(g2->u.sexpr,nil))); return g3; } sexpr *adorn(sexpr *g1, sexpr *g2) { sexpr *g3; if (!graphic(g1)) { print_value_escape("adorn: Argument is not a graphic object: ",g1); } if (!graphic(g2)) { print_value_escape("adorn: Argument is not a graphic object: ",g2); } if (null(g1->u.sexpr)) return g2; if (null(g2->u.sexpr)) return g1; g3 = smalloc++; g3->type = GRAPHIC; g3->u.sexpr = cons(symbol$adorn,cons(g1->u.sexpr,cons(g2->u.sexpr,nil))); return g3; } sexpr *make_plumber(double x, double y, double heading) { sexpr *p; p = smalloc++; p->type = PLUMBER; p->u.plumber=(plumber *) malloc(sizeof(plumber)); p->u.plumber->x = x; p->u.plumber->y = y; p->u.plumber->heading = heading; return p; } sexpr *straight(sexpr *len) { sexpr *g; if (!number(len)) { print_value_escape("straight: Length argument is not a number: ",len); } g = smalloc++; g->type = GRAPHIC; g->u.sexpr = cons(symbol$straight,cons(len,cons(graphic$default_color,cons(num2exp(1.0),nil)))); return g; } sexpr *fancy_straight(sexpr *len, sexpr *color, sexpr *width) { sexpr *g; if (!number(len)) { print_value_escape("straight: Length argument is not a number: ",len); } if (!pair(color)) { print_value_escape("straight: Color argument is not a pair: ",color); } if (!number(width)) { print_value_escape("straight: Width argument is not a number: ",width); } g = smalloc++; g->type = GRAPHIC; g->u.sexpr = cons(symbol$straight,cons(len,cons(color,cons(width,nil)))); return g; } sexpr *spot(sexpr *width) { sexpr *g; if (!number(width)) { print_value_escape("spot: Width argument is not a number: ",width); } g = smalloc++; g->type = GRAPHIC; g->u.sexpr = cons(symbol$spot,cons(width,cons(graphic$default_color,nil))); return g; } sexpr *fancy_spot(sexpr *width, sexpr *color) { sexpr *g; if (!number(width)) { print_value_escape("spot: Width argument is not a number: ",width); } if (!pair(color)) { print_value_escape("spot: Color argument is not a pair: ",color); } g = smalloc++; g->type = GRAPHIC; g->u.sexpr = cons(symbol$spot,cons(width,cons(color,nil))); return g; } sexpr *transparent(sexpr *len) { sexpr *g; if (!number(len)) { print_value_escape("transparent: Length argument is not a number: ",len); } g = smalloc++; g->type = GRAPHIC; g->u.sexpr = cons(symbol$transparent,cons(len,nil)); return g; } sexpr *bend(sexpr *angle) { sexpr *g; if (!number(angle)) { print_value_escape("bend: Argument is not a number: ",angle); } g = smalloc++; g->type = GRAPHIC; g->u.sexpr = cons(symbol$bend,cons(angle,nil)); return g; } sexpr *text(sexpr *text) { sexpr *g; if (!string(text)) { print_value_escape("text: Argument is not a string: ",text); } g = smalloc++; g->type = GRAPHIC; g->u.sexpr = cons(symbol$text,cons(text,cons(graphic$default_color,cons(num2exp(1.0),nil)))); return g; } sexpr *fancy_text(sexpr *text, sexpr *color, sexpr *size) { sexpr *g; if (!string(text)) { print_value_escape("text: Argument is not a string: ",text); } if (!pair(color)) { print_value_escape("text: Color argument is not a pair: ",color); } if (!number(size)) { print_value_escape("text: Size argument is not a number: ",size); } g = smalloc++; g->type = GRAPHIC; g->u.sexpr = cons(symbol$text,cons(text,cons(color,cons(size,nil)))); return g; } sexpr *gcolor(sexpr *g, sexpr *color) { sexpr *first; if (!graphic(g)) { print_value_escape("gcolor: Argument is not a graphic object: ",g); } if (!pair(color)) { print_value_escape("gcolor: Color argument is not a pair: ",color); } if (null(g->u.sexpr)) return g; first = car(g->u.sexpr); if (eq(first,symbol$adjoin)) return adjoin(gcolor(gcar(g),color),gcolor(gcdr(g),color)); else if (eq(first,symbol$adorn)) return adorn(gcolor(gcar(g),color),gcolor(gcdr(g),color)); else if (eq(first,symbol$straight)) return fancy_straight(cadr(g->u.sexpr),color,cadr(cddr(g->u.sexpr))); else if (eq(first,symbol$spot)) return fancy_spot(cadr(g->u.sexpr),color); else if (eq(first,symbol$text)) return fancy_text(cadr(g->u.sexpr),color,cadr(cddr(g->u.sexpr))); else return g; } sexpr *gscale(sexpr *g, sexpr *lscale, sexpr *wscale ) { sexpr *first; if (!graphic(g)) { print_value_escape("gscale: Argument is not a graphic object: ",g); } if (!number(lscale)) { print_value_escape("gscale: Length-scale argument is not a number: ",lscale); } if (!number(wscale)) { print_value_escape("gscale: Width-scale argument is not a number: ",wscale); } if (null(g->u.sexpr)) return g; first = car(g->u.sexpr); if (eq(first,symbol$adjoin)) return adjoin(gscale(gcar(g),lscale,wscale),gscale(gcdr(g),lscale,wscale)); else if (eq(first,symbol$adorn)) return adorn(gscale(gcar(g),lscale,wscale),gscale(gcdr(g),lscale,wscale)); else if (eq(first,symbol$straight)) return fancy_straight(times(cadr(g->u.sexpr),lscale),car(cddr(g->u.sexpr)),times(cadr(cddr(g->u.sexpr)),wscale)); else if (eq(first,symbol$transparent)) return transparent(times(cadr(g->u.sexpr),lscale)); else if (eq(first,symbol$text)) return fancy_text(cadr(g->u.sexpr),car(cddr(g->u.sexpr)),times(cadr(cddr(g->u.sexpr)),lscale)); else return g; } sexpr *gmirror(sexpr *g) { sexpr *first; if (!graphic(g)) { print_value_escape("gmirror: Argument is not a graphic object: ",g); } if (null(g->u.sexpr)) return g; first = car(g->u.sexpr); if (eq(first,symbol$adjoin)) return adjoin(gmirror(gcar(g)),gmirror(gcdr(g))); else if (eq(first,symbol$adorn)) return adorn(gmirror(gcar(g)),gmirror(gcdr(g))); else if (eq(first,symbol$bend)) return bend(minus(num2exp(0.0),cadr(g->u.sexpr))); else return g; } sexpr *greverse(sexpr *g) { sexpr *first; if (!graphic(g)) { print_value_escape("greverse: Argument is not a graphic object: ",g); } if (null(g->u.sexpr)) return g; first = car(g->u.sexpr); if (eq(first,symbol$adjoin)) return adjoin(greverse(gcdr(g)),greverse(gcar(g))); else if (eq(first,symbol$adorn)) return adorn(greverse(gcar(g)),greverse(gcdr(g))); else return g; } sexpr *make_primitive(sexpr *(*func)(), int arity, char *name) { sexpr *sp = malloc(sizeof(sexpr)); primitive *primitive = malloc(sizeof(primitive)); sp->type = PRIMITIVE; sp->u.primitive = primitive; primitive->arity= arity; primitive->func = func; primitive->name = name; return sp; } sexpr *cps2vector_helper(sexpr *code) { if (null(code)) return nil; else { sexpr *first = car(code); if (first==bytecode$constant || first==bytecode$close) return cons(first,cons(cadr(code),cps2vector_helper(cddr(code)))); else return cons(cps2vector(first),cps2vector_helper(cdr(code))); } } sexpr *cps2vector(sexpr *code) { if (pair(code)) return list2vector(cps2vector_helper(code)); else return code; } sexpr *compile_symbol(sexpr *symbol, sexpr *next, sexpr *env) { sexpr *frame, *rest_env, *result; int j = 0; rest_env = env; frame = rest_env->u.pair->car; while (frame != global_vars) { result = assoc(symbol,frame); if (!undefined(result)) { if (j == 0) { return cons(append(list2(bytecode$local_refer,result),next),env); } else { return cons(append(cons(bytecode$refer,list2(num2exp(j),result)),next),env); } } rest_env = rest_env->u.pair->cdr; frame = rest_env->u.pair->car; j++; } result = assoc(symbol,global_vars); if (undefined(result)) result = add_new_global_def(symbol,undefined); return cons(append(list2(bytecode$global_refer,result),next),env); } sexpr *compile_set_bang(sexpr *symbol, sexpr *value, sexpr *next, sexpr *env) { sexpr *frame, *rest_env, *result; int j = 0; rest_env = env; frame = rest_env->u.pair->car; while (frame != global_vars) { result = assoc(symbol,frame); if (!undefined(result)) { if (j == 0) { return compile(value,append(list2(bytecode$local_assign,result),next),env); } else { return compile(value,append(cons(bytecode$assign,list2(num2exp(j),result)),next),env); } } rest_env = rest_env->u.pair->cdr; frame = rest_env->u.pair->car; j++; } result = assoc(symbol,global_vars); if (undefined(result)) result = add_new_global_def(symbol,undefined); return compile(value,append(list2(bytecode$global_assign,result),next),env); } sexpr *compile_body(sexpr *body, sexpr *next, sexpr *env) { if (!null(body)) { sexpr *extended_env = cdr(compile(car(body),next,env)); return car(compile(car(body),compile_body(cdr(body),next,extended_env),env)); } else return next; } sexpr *compile(sexpr *x, sexpr *next, sexpr *env) { if (symbol(x)) return compile_symbol(x,next,env); if (self_evaluating(x)) return cons(append(list2(bytecode$constant,x),next),env); if (pair(x)) { sexpr *first = car(x); sexpr *args, *c; if (first==symbol$quote) return cons(append(list2(bytecode$constant,cadr(x)),next),env); if (first==symbol$lambda) { sexpr *vars = cadr(x); int len = length(vars); c = compile_body(cddr(x),cons(bytecode$return,nil),cons(list2frame(len,vars),env)); return cons(append(cons(bytecode$close,list2(num2exp(len),c)),next),env); } if (first==symbol$begin) { if (null(cdr(x))) return cons(compile(undefined,next,env),env); else return cons(compile_body(cdr(x),next,env),env); } if (first==symbol$if) { sexpr *test, *consequence, *alternative; test = cdr(x); consequence = cdr(test); alternative = cdr(consequence); if (!null(alternative)) { alternative = car(compile(cadr(consequence),next,env)); consequence = car(compile(car(consequence),next,env)); } else { consequence = car(compile(car(consequence),next,env)); alternative = car(compile(undefined,next,env)); } test = car(test); return compile(test,append(list2(bytecode$test,consequence),alternative),env); } if (first==symbol$define) { sexpr *symbol = cdr(x), *value = cadr(symbol); return compile(value,append(list2(bytecode$define,add_def_compile(car(symbol),env)),next),env); } if (first==symbol$set_bang) { sexpr *symbol = cdr(x), *value = cadr(symbol); return compile_set_bang(car(symbol),value,next,env); } if (first==symbol$callcc) { sexpr *esc = cadr(x); c = cons(bytecode$conti,append(cons(bytecode$argument,nil),car(compile(esc,cons(bytecode$funcall,nil),env)))); if (car(next) != bytecode$return) return cons(append(list2(bytecode$frame,next),c),env); else return cons(c,env); } if (first==symbol$callec) { sexpr *esc = cadr(x); c = cons(bytecode$throw,append(cons(bytecode$argument,nil),car(compile(esc,cons(bytecode$funcall,nil),env)))); if (car(next) != bytecode$return) return cons(append(list2(bytecode$frame,next),c),env); else return cons(c,env); } if (first==symbol$define_macro) { sexpr *var = cdr(x); add_macro(car(var),eval(cadr(var))); return compile(undefined,next,env); } if (first==symbol$apply) { args = cddr(x); c = car(compile(cadr(x),cons(bytecode$apply,nil),env)); while (!null(args)) { c = car(compile(car(args),append(cons(bytecode$argument,nil),c),env)); args = cdr(args); } if (car(next) != bytecode$return) return cons(append(list2(bytecode$frame,next),c),env); else return cons(c,env); } args = cdr(x); c = car(compile(car(x),cons(bytecode$funcall,nil),env)); while (!null(args)) { c = car(compile(car(args),append(cons(bytecode$argument,nil),c),env)); args = cdr(args); } if (car(next) != bytecode$return) return cons(append(list2(bytecode$frame,next),c),env); else return cons(c,env); return cons(append(list2(bytecode$frame,next),c),env); } printf("Compilation error.\n"); longjmp(esc,1); } sexpr *print_closure(sexpr *sp) { if (!virgin(sp) && !closure(sp)) { printf("print-closure: Argument must be closure.\n"); longjmp(esc,1); } printf("n = %d\n",sp->u.closure->n); printf("lookup = "); print(sp->u.closure->lookup); printf("compiled = "); print(sp->u.closure->compiled); printf("env = "); print(sp->u.closure->env); return undefined; } sexpr *print_macro(sexpr *sp) { if (!symbol(sp)) { printf("print-macro: Argument must be symbol.\n"); longjmp(esc,1); } printf("n = %d\n",sp->u.symbol->macro->n); printf("compiled = "); print(sp->u.symbol->macro->compiled); printf("env = "); print(sp->u.symbol->macro->env); return undefined; } #define install_constant(name,constant) add_new_global_def(str2symbol(name),constant) #define install_function(name,function,arity) add_new_global_def(str2symbol(name),make_primitive(function,arity,name)) sexpr *make_global_env() { /* Install Scheme primitives */ install_function("car",car,1); install_function("cdr",cdr,1); add_new_global_def(str2symbol("cons"),primitive$cons); install_function("eq?",eqp,2); install_function("positive?",positivep,1); install_function("negative?",negativep,1); install_function("zero?",zerop,1); install_function("scheme:=",equals,2); install_function("scheme:<",lt,2); install_function("scheme:>",gt,2); install_function("scheme:<=",leq,2); install_function("scheme:>=",geq,2); install_function("sin",sine,1); install_function("cos",cosine,1); install_function("tan",tangent,1); install_function("asin",arcsine,1); install_function("acos",arccosine,1); install_function("scheme:atan",arctangent,2); install_function("log",logarithm,1); install_function("exp",exponential,1); install_function("expt",expt,2); install_function("sqrt",square_root,1); install_function("abs",scheme_abs,1); install_function("floor",scheme_floor,1); install_function("ceil",scheme_ceil,1); install_function("real-part",real_part,1); install_function("imag-part",imag_part,1); install_function("angle",angle,1); install_function("magnitude",scheme_abs,1); install_function("conjugate",conjugate,1); install_function("polar->complex",polar2complex,2); install_function("complex",make_complex,2); install_function("conjugate",conjugate,1); install_function("scheme:min",scheme_min,2); install_function("scheme:max",scheme_max,2); install_function("random",scheme_random,1); install_function("scheme:reduce",reduce,3); install_function("scheme:reduce-pred",reduce_pred,2); install_function("scheme:map",scheme_map,2); add_new_global_def(str2symbol("scheme:append"),primitive$append); install_function("reverse",reverse,1); install_function("null?",nullp,1); install_function("symbol?",symbolp,1); install_function("char?",characterp,1); install_function("string?",stringp,1); install_function("number?",numberp,1); install_function("boolean?",booleanp,1); install_function("procedure?",procedurep,1); install_function("primitive?",primitivep,1); install_function("pair?",pairp,1); install_function("vector?",vectorp,1); install_function("complex?",complexp,1); install_function("real?",realp,1); install_function("not",not,1); install_function("delete!",delete,2); install_function("set-car!",setcar,2); install_function("set-cdr!",setcdr,2); install_function("load",load,1); install_function("exit",scheme_exit,0); install_function("scheme:read",scheme_read,1); install_function("scheme:display",display,2); install_function("scheme:write",scheme_write,2); install_function("scheme:newline",newline,1); install_function("scheme:error",scheme_error,0); install_function("void",scheme_void,0); install_function("void?",voidp,1); install_function("gensym",gensym,0); install_function("character->string",char2str,1); install_function("number->string",num2str,1); install_function("string->number",str2num,1); install_function("symbol->string",symbol2exp,1); install_function("string->symbol",exp2symbol,1); install_function("scheme:string-append",string_append,2); install_function("string-ref",string_ref,2); install_function("substring",substring,3); install_function("string-length",string_length,1); install_function("scheme:string=?",string_eq,2); install_function("scheme:string?",string_gt,2); install_function("scheme:string<=?",string_leq,2); install_function("scheme:string>=?",string_geq,2); install_function("scheme:char=?",char_eq,2); install_function("scheme:char?",char_lt,2); install_function("scheme:char<=?",char_leq,2); install_function("scheme:char>=?",char_geq,2); install_function("list->vector",list2vector,1); install_function("vector->list",vector2list,1); install_function("scheme:make-vector",scheme_make_vector,2); install_function("vector-set!",vector_set,3); install_function("vector-ref",vector_ref,2); install_function("vector-length",vector_length,1); install_function("open-input-file",open_input_file,1); install_function("open-output-file",open_output_file,1); install_function("close-input-port",close_input_port,1); install_function("close-output-port",close_output_port,1); install_function("input-port?",input_portp,1); install_function("output-port?",output_portp,1); install_function("eof-object?",eof_objectp,1); install_function("scheme:read-char",read_char,1); install_function("scheme:write-char",write_char,2); install_function("scheme:peek-char",peek_char,1); install_function("scheme:+",plus,2); install_function("scheme:-",minus,2); install_function("scheme:*",times,2); install_function("scheme:/",divide,2); install_function("fscanf",scheme_fscanf,2); install_function("scheme:fprintf",scheme_fprintf,3); install_function("undefined?",undefinedp,1); install_function("expand",expand,1); install_function("eval",eval,1); install_function("compile",scheme_compile,1); install_function("print-closure",print_closure,1); install_function("print-macro",print_macro,1); install_function("system",scheme_system,1); install_function("gil-galad",gilgalad,0); /* Install Plumbing graphics constants */ install_constant("scheme:graphic-nil",graphic$nil); install_constant("scheme:graphic-white",color(255.0,255.0,255.0)); install_constant("scheme:graphic-black",color(0.0,0.0,0.0)); install_constant("scheme:graphic-red",color(255.0,0.0,0.0)); install_constant("scheme:graphic-green",color(0.0,255.0,0.0)); install_constant("scheme:graphic-blue",color(0.0,0.0,255.0)); install_constant("scheme:graphic-orange",color(255.0,64.0,0.0)); /* Install Scheme constants */ install_constant("scheme:current-input-port",current_input_port); install_constant("scheme:current-output-port",current_output_port); /* Install Boldt line grouping primitives */ install_function("scheme:zero-crossings",zerocrossings,4); install_function("scheme:filter-on-contrast",filter_on_contrast,3); install_function("scheme:filter-on-length",filter_on_length,3); install_function("filter-on-angle",filter_on_angle,3); install_function("po-link!",po_link,2); install_function("po-replace!",po_replace,2); install_function("scheme:po-merge!",po_merge,2); install_function("svd",test_svd,1); install_function("make-sketch",make_sketch,3); install_function("sketch?",sketchp,1); install_function("sketch-lines",sketch_lines,1); install_function("sketch-rows",sketch_rows,1); install_function("sketch-cols",sketch_cols,1); install_function("sketch-count",sketch_count,1); install_function("sketch-union",sketch_union,2); install_function("scheme:make-line",make_line,5); install_function("line?",linep,1); install_function("line-x0",line_x0,1); install_function("line-y0",line_y0,1); install_function("line-x1",line_x1,1); install_function("line-y1",line_y1,1); install_function("line-contrast",line_contrast,1); install_function("line-theta",line_theta,1); install_function("line-length",line_length,1); install_function("line-links",line_links,1); install_function("set-line-contrast!",set_line_contrast,2); install_function("sketch->postscript",sketch2postscript,2); #ifdef GL install_function("show-link-graph",show_link_graph,1); #endif /* Install Plumbing graphics primitives */ install_function("scheme:fancy-straight",fancy_straight,3); install_function("scheme:spot",spot,1); install_function("scheme:fancy-spot",fancy_spot,2); install_function("transparent",transparent,1); install_function("graphic?",graphicp,1); install_function("graphic-type",gtype,1); install_function("graphic-args",gargs,1); install_function("graphic-null?",gnull,1); install_function("scheme:text",text,1); install_function("scheme:fancy-text",fancy_text,3); install_function("gcar",gcar,1); install_function("gcdr",gcdr,1); install_function("gmirror",gmirror,1); install_function("greverse",greverse,1); install_function("gcolor",gcolor,2); install_function("scheme:gscale",gscale,3); install_function("set-graphic-color!",set_graphic_default_color,1); install_function("scheme:adjoin",adjoin,2); install_function("scheme:adorn",adorn,2); install_function("bend",bend,1); install_function("scheme:straight",straight,1); install_function("graphic->postscript",graphic2postscript,2); /* Install image processing primitives */ install_function("set-image-display-scale!",set_image_display_scale,1); install_function("scheme:image-map",scheme_image_map,2); install_function("read-image",readpgm,1); install_function("scheme:write-image",writepgm,2); install_function("scheme:image-crop",image_crop,5); install_function("scheme:image-pad",image_pad,3); install_function("scheme:make-image",scheme_make_image,4); install_function("make-hot-image",make_hot_image,1); install_function("read-color-image",readppm,1); install_function("scheme:write-color-image",writeppm,2); install_function("rgb->color-image",rgb2color,3); install_function("color-image-red",color_image_red,1); install_function("color-image-green",color_image_green,1); install_function("color-image-blue",color_image_blue,1); install_function("rgb->hsi",rgb2hsi,3); install_function("hsi->rgb",hsi2rgb,3); install_function("color-image?",color_imagep,1); install_function("complex-image?",complex_imagep,1); install_function("image->complex-image",image2complex,1); install_function("complex-image->color-image",complex2color,1); install_function("image?",imagep,1); install_function("image-rows",image_rows,1); install_function("image-cols",image_cols,1); install_function("image-ref",image_ref,3); install_function("image-set!",image_set,4); install_function("image-transpose",image_transpose,1); install_function("scheme:top-to-bottom",scheme_top2bottom,2); install_function("scheme:left-to-right",scheme_left2right,2); install_function("upsample-rows",upsample_rows,1); install_function("upsample-cols",upsample_cols,1); install_function("downsample-rows",downsample_rows,1); install_function("downsample-cols",downsample_cols,1); install_function("convolve",convolve,2); install_function("array->image",array2image,1); install_function("image->array",image2array,1); install_function("shrink",shrink,2); install_function("image-normalize",image_normalize,1); install_function("color-image-normalize",color_image_normalize,1); install_function("scheme:fft",fft,2); install_function("scheme:image-fft",image_fft,2); install_function("scheme:make-covariance-matrix",make_covariance_matrix,1); install_function("matrix-product",matrix_product,2); install_function("image-sum",image_sum,1); install_function("image-max",image_max,1); install_function("image-min",image_min,1); install_function("label",label,1); install_function("distance-transform",distance_transform,1); install_function("scheme:outline",outline,3); install_function("perimeters",perimeters,1); install_function("areas",areas,1); install_function("centers-of-mass",centers_of_mass,1); install_function("bounding-boxes",bounding_boxes,1); install_function("image->list",image2list,1); install_function("scheme:register-images",register_images,6); return undefined; } unsigned hash(char *name) { unsigned hashval; for (hashval = 0; *name != '\0'; name++) hashval = *name + 31 * hashval; return hashval % HASHSIZE; } node *symbol_table_lookup(char *name) { node *np; for (np=hashtab[hash(name)]; np != NULL; np=np->next) if (!strcmp(name,np->sp->u.symbol->name)) return np; return NULL; } void symbol_table_delete(char *name) { node *n1, *n2; int hashval = hash(name); n2=hashtab[hashval]; if (!strcmp(name,n2->sp->u.symbol->name)) { hashtab[hashval] = n2->next; free(n2); } else { n1=n2; n2=n2->next; while (n2 != NULL) { if (!strcmp(name,n2->sp->u.symbol->name)) { n1->next = n2->next; free(n2); break; } else { n1=n2; n2=n2->next; } } } } sexpr *str2symbol(char *name) { node *np; unsigned hashval; np = symbol_table_lookup(name); if (np == NULL) { np = (node *) malloc(sizeof(node)); np->sp = malloc(sizeof(sexpr)); np->sp->type = SYMBOL; hashval = hash(name); np->sp->u.symbol = malloc(sizeof(symbol)); np->sp->u.symbol->name = malloc((strlen(name) + 1)*sizeof(char)); strcpy(np->sp->u.symbol->name,name); np->sp->u.symbol->macro = (closure *) NULL; np->next = hashtab[hashval]; hashtab[hashval] = np; } return np->sp; } sexpr *str2collectsymbol(char *name) { node *np; unsigned hashval; np = symbol_table_lookup(name); if (np == NULL) { np = (node *) malloc(sizeof(node)); np->sp = smalloc++; np->sp->type = SYMBOL; hashval = hash(name); np->sp->u.symbol = malloc(sizeof(symbol)); np->sp->u.symbol->name = malloc((strlen(name) + 1)*sizeof(char)); strcpy(np->sp->u.symbol->name,name); np->sp->u.symbol->macro = (closure *) NULL; np->next = hashtab[hashval]; hashtab[hashval] = np; } return np->sp; } sexpr *copy(sexpr *start, sexpr *end, sexpr *sp0) { int i, len; sexpr *sp1; if (sp0 < start || sp0 >= end) return sp0; if (copied(sp0)) return sp0->u.sexpr; sp1 = smalloc++; *sp1 = *sp0; if (symbol(sp0)) { node *np = symbol_table_lookup(sp0->u.symbol->name); np->sp = sp1; } sp0->type = COPIED; sp0->u.sexpr = sp1; switch (sp1->type) { case PAIR: sp1->u.pair->car = copy(start,end,sp1->u.pair->car); sp1->u.pair->cdr = copy(start,end,sp1->u.pair->cdr); break; case SYMBOL: if (sp1->u.symbol->macro != NULL) { sp1->u.symbol->macro->compiled = copy(start,end,sp1->u.symbol->macro->compiled); sp1->u.symbol->macro->env = copy(start,end,sp1->u.symbol->macro->env); } break; case VIRGIN: case CLOSURE: sp1->u.closure->lookup = copy(start,end,sp1->u.closure->lookup); sp1->u.closure->compiled = copy(start,end,sp1->u.closure->compiled); sp1->u.closure->env = copy(start,end,sp1->u.closure->env); break; case VECTOR: len = sp1->u.vector->length; for (i=0; i < len; i++) { sp1->u.vector->vector[i] = copy(start,end,sp1->u.vector->vector[i]); } break; case GRAPHIC: sp1->u.sexpr = copy(start,end,sp1->u.sexpr); break; case LINE: sp1->u.line->links[0] = copy(start,end,sp1->u.line->links[0]); sp1->u.line->links[1] = copy(start,end,sp1->u.line->links[1]); break; case SKETCH: len = sp1->u.sketch->bixel_rows*sp1->u.sketch->bixel_cols; for (i=0; i < len; i++) { sp1->u.sketch->grid[i] = copy(start,end,sp1->u.sketch->grid[i]); } sp1->u.sketch->ls = copy(start,end,sp1->u.sketch->ls); break; case CONTINUATION: len = sp1->u.continuation->s; for (i = 0; i < len; i++) { *(sp1->u.continuation->pc_stack[i]) = copy(start,end,*(sp1->u.continuation->pc_stack[i])); sp1->u.continuation->cv_stack[i] = copy(start,end,sp1->u.continuation->cv_stack[i]); sp1->u.continuation->e_stack[i] = copy(start,end,sp1->u.continuation->e_stack[i]); sp1->u.continuation->lu_stack[i] = copy(start,end,sp1->u.continuation->lu_stack[i]); } len = sp1->u.continuation->r_stack[len-1]; for (i = 0; i < len; i++) { sp1->u.continuation->args[i] = copy(start,end,sp1->u.continuation->args[i]); } } return sp1; } void free_garbage(sexpr *start, sexpr *end) { sexpr *sp; for (sp = start; sp < end; sp++) { if (!copied(sp)) { switch (sp->type) { case PAIR: free(sp->u.pair); break; case SYMBOL: symbol_table_delete(sp->u.symbol->name); free(sp->u.symbol->name); free(sp->u.symbol); break; case VIRGIN: case CLOSURE: free(sp->u.closure); break; case VECTOR: free(sp->u.vector->vector); free(sp->u.vector); break; case PLUMBER: free(sp->u.plumber); break; case IMAGE: case COLOR_IMAGE: free(sp->u.image->data); free(sp->u.image); break; case COMPLEX_IMAGE: free(sp->u.complex_image->data); free(sp->u.complex_image); break; case LINE: free(sp->u.line); break; case SKETCH: free(sp->u.sketch->grid); free(sp->u.sketch); break; case CONTINUATION: free(sp->u.continuation->args); free(sp->u.continuation->r_stack); free(sp->u.continuation->pc_stack); free(sp->u.continuation->cv_stack); free(sp->u.continuation->e_stack); free(sp->u.continuation->lu_stack); free(sp->u.continuation); break; } } } } void garbage_collect() { sexpr *start, *end; end = smalloc; /* switch copy spaces */ if (gc) { start = (sexpr *) smalloc1; smalloc = (sexpr *) smalloc0; gc = 0; } else { start = (sexpr *) smalloc0; smalloc = (sexpr *) smalloc1; gc = 1; } /* initiate copying from all roots */ global_vars = copy(start, end, global_vars); global_vals = copy(start, end, global_vals); global_env_vars = copy(start, end, global_env_vars); global_env_vals = copy(start, end, global_env_vals); primitive$cons = copy(start, end, primitive$cons); primitive$append = copy(start, end, primitive$append); macros = copy(start, end, macros); graphic$default_color = copy(start, end, graphic$default_color); /* free non-sexpr stuff pointed to by uncopied sexpr's */ free_garbage(start, end); } void print_table() { int i; node *np; printf("\n"); for (i = 0; i < HASHSIZE; i++) { np = hashtab[i]; while (np != NULL) { print(np->sp); np = np->next; } } }