tomwalters@0: /* tomwalters@0: Copyright (c) Applied Psychology Unit, Medical Research Council. 1988, 1989 tomwalters@0: =========================================================================== tomwalters@0: tomwalters@0: Permission to use, copy, modify, and distribute this software without fee tomwalters@0: is hereby granted for research purposes, provided that this copyright tomwalters@0: notice appears in all copies and in all supporting documentation, and that tomwalters@0: the software is not redistributed for any fee (except for a nominal shipping tomwalters@0: charge). Anyone wanting to incorporate all or part of this software in a tomwalters@0: commercial product must obtain a license from the Medical Research Council. tomwalters@0: tomwalters@0: The MRC makes no representations about the suitability of this tomwalters@0: software for any purpose. It is provided "as is" without express or implied tomwalters@0: warranty. tomwalters@0: tomwalters@0: THE MRC DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING tomwalters@0: ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL THE tomwalters@0: A.P.U. BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY tomwalters@0: DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN tomwalters@0: AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF tomwalters@0: OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. tomwalters@0: */ tomwalters@0: tomwalters@0: /* tomwalters@0: spiral.c tomwalters@0: ======== tomwalters@0: tomwalters@0: Tap SAI frames and draw each one as a spiral in the given window. tomwalters@0: */ tomwalters@0: tomwalters@0: /* tomwalters@0: Modifications made to the routines draw_spiral, map_ribbon, map_spiral tomwalters@0: so that the nwidth section of the SAI image is taken into account while tomwalters@0: doing the calculations, (in image.c, doColumn()) and while plotting, tomwalters@0: only the positive section of the time axis is drawn. tomwalters@0: Added a new field "nwid" in the _draw_state object. tomwalters@0: 17 March, 1995. A Jay Datta tomwalters@0: */ tomwalters@0: tomwalters@0: #include tomwalters@0: #include tomwalters@0: tomwalters@0: #include "windows.h" tomwalters@0: #include "stitch.h" tomwalters@0: #include "source.h" tomwalters@0: #include "draw.h" tomwalters@0: #include "ops.h" tomwalters@0: #include "spiral.h" tomwalters@0: tomwalters@0: #ifndef lint tomwalters@0: static char *sccs_id = "@(#)spiral.c 1.9 Mike Allerhand (MRC-APU) 6/1/91" ; tomwalters@0: #endif tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * Misc definitions tomwalters@0: ****************************************************************************/ tomwalters@0: #define TWOPI 6.28318530717 tomwalters@0: #define log2(x) (log((double)x)/log(2.0)) tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * Model parameters. tomwalters@0: * Defaults: tomwalters@0: * form_spl = 'arch' tomwalters@0: * dotsize_spl = 2 tomwalters@0: * axis_spl = off tomwalters@0: * zero_spl = 4.072 tomwalters@0: * dotthresh_spl = 50 tomwalters@0: ****************************************************************************/ tomwalters@0: int frameheight, framewidth, nwid; /* Dimensions of each sai */ tomwalters@0: tomwalters@0: char form_spl = 'A' ; /* Form of spiral, 'A' or 'L' */ tomwalters@0: int dotsize_spl = 2 ; /* Dotsize in pixels */ tomwalters@0: int axis_spl = 0 ; /* Flag for spiral axis */ tomwalters@0: double zero_spl = 4.072 ; /* For rotation and circuit removal */ tomwalters@0: int dotthresh_spl = 25 ; /* Min peak height */ tomwalters@0: tomwalters@0: /*=================== COORDINATE SYSTEM OF THE DATA =======================*/ tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * Size and location, (program arguments), of calibrated box, in tomwalters@0: * coord system of the data. tomwalters@0: * range = range or size of both axes, (ie. "magnification" of plot). tomwalters@0: * xorg,yorg = real plot coordinates at centre of screen. tomwalters@0: ****************************************************************************/ tomwalters@0: static float range; tomwalters@0: static float xorg; tomwalters@0: static float yorg; tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * Limits of calibrated box in coord system of the data. tomwalters@0: ****************************************************************************/ tomwalters@0: static float Xd0, Yd0; /* Bottom-left corner */ tomwalters@0: static float Xd1, Yd1; /* Top-right corner */ tomwalters@0: tomwalters@0: tomwalters@0: /*===================== PIXEL COORDINATE SYSTEM ==========================*/ tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * Dimensions of window (box plus border) in pixel coord system. tomwalters@0: ****************************************************************************/ tomwalters@0: static int width,height; tomwalters@0: static int border; /* border = min(width,height) / 6; */ tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * Limits of on-screen calibrated box in pixel coord system. tomwalters@0: ****************************************************************************/ tomwalters@0: static int Xp0, Yp0; /* Bottom-left corner (border, height-border) */ tomwalters@0: static int Xp1, Yp1; /* Top-right corner (width-border, border) */ tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * Dimensions of calibrated box in pixel coord system tomwalters@0: ****************************************************************************/ tomwalters@0: static int Wb; /* Box width Wb = Xp1-Xp0 */ tomwalters@0: static int Hb; /* Box height Hb = Yp0-Yp1 */ tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * The parameters for the calibrated box in the pixel coord system are tomwalters@0: * initialized for a given window size (width, height, in pixels), tomwalters@0: * as follows: tomwalters@0: * tomwalters@0: * border = min(width,height) / 6; tomwalters@0: * Xp0 = border; tomwalters@0: * Xp1 = width - border; tomwalters@0: * Yp0 = border; tomwalters@0: * Yp1 = height - border; tomwalters@0: * Wb = Xp1-Xp0; tomwalters@0: * Hb = Yp1-Yp0; tomwalters@0: ****************************************************************************/ tomwalters@0: tomwalters@0: /*================= COORDINATE SYSTEM TRANSFORMATIONS =====================*/ tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * Transformations from coord system of the data to pixel coord system. tomwalters@0: * For example, if the origin in the data coord system is (0,0), then tomwalters@0: * the origin in the pixel coord system is (px(0),py(0)). tomwalters@0: ****************************************************************************/ tomwalters@0: #define px(x) ((int)(Xp0+Wb/2+(x-xorg)*(Wb/range))) tomwalters@0: #define py(y) ((int)(Yp0+Hb/2+(y-yorg)*(Hb/range))) tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * Screen Position (in pixel coord system): test position is within box. tomwalters@0: ****************************************************************************/ tomwalters@0: #define Xinbox(x) ((Xp0<=x && x<=Xp1) ? 1 : 0) tomwalters@0: #define Yinbox(y) ((Yp0<=y && y<=Yp1) ? 1 : 0) tomwalters@0: #define inbox(x,y) ((Xinbox(x) && Yinbox(y)) ? 1 : 0) tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * This routine is called once per SAI frame. The SAI data is in array "frame". tomwalters@0: * The SAI frame is formatted row-wise, starting with the lowest-frequency row. tomwalters@0: * The state structure is: tomwalters@0: * struct _draw_state { tomwalters@0: * Source source; tomwalters@0: * WindowObject window; tomwalters@0: * int min, max; tomwalters@0: * int framewidth, frameheight; dimensions of the SAI tomwalters@0: * int nwid; tomwalters@0: * long frames; tomwalters@0: * long framenumber; tomwalters@0: * void (*interceptor)(); function which stores window for replay tomwalters@0: * void (*drawer)(); drawing function tomwalters@0: * }; tomwalters@0: * tomwalters@0: * This struct is defined in stitch/draw.h tomwalters@0: * It is initialized by routine "SourceDraw" in stitch/draw.c, which is called tomwalters@0: * from gen.c. Note that this takes place after the SaiEntry routine has been tomwalters@0: * called, so the sai dimensions set there will apply. tomwalters@0: * tomwalters@0: ****************************************************************************/ tomwalters@0: void draw_spiral( state, image ) tomwalters@0: struct _draw_state *state ; tomwalters@0: short *image ; tomwalters@0: { tomwalters@0: static int first=1; /* flag for first-call of this routine */ tomwalters@0: tomwalters@0: /* Coordinate arrays for axis and spiral points */ tomwalters@0: static short *Xaxis; tomwalters@0: static short *Yaxis; tomwalters@0: static int naxis = 0; /* number of points plotted along spiral axis */ tomwalters@0: /* (naxis = (100 * nsamples) +1) */ tomwalters@0: static short *Xspiral; tomwalters@0: static short *Yspiral; tomwalters@0: int npoints = 0; /* number of samples in pulse ribbon and spiral */ tomwalters@0: tomwalters@0: static short *Xdot; tomwalters@0: static short *Ydot; tomwalters@0: tomwalters@0: static short *frame ; tomwalters@0: tomwalters@0: frameheight = state->frameheight; tomwalters@0: framewidth= (state->framewidth)+(state->nwid); /* -nwid */ tomwalters@0: /* fprintf(stderr, "THE NWID VALUE is %d\n", state->nwid); */ tomwalters@0: nwid=-state->nwid; tomwalters@0: /* fprintf(stderr, "framewidth=%d, frameheight=%d\n", frameheight, framewidth); */ tomwalters@0: /* Initialize parameters of pixel coord system (see spiral.h) */ tomwalters@0: width = Width( state->window ) ; tomwalters@0: height = Height( state->window ) ; tomwalters@0: tomwalters@0: border = 0; tomwalters@0: tomwalters@0: Xp0 = border; tomwalters@0: Xp1 = width - border; tomwalters@0: Yp0 = border; tomwalters@0: Yp1 = height - border; tomwalters@0: Wb = Xp1-Xp0; tomwalters@0: Hb = Yp1-Yp0; tomwalters@0: tomwalters@0: /* Initialize parameters of data coord system */ tomwalters@0: if (form_spl == 'A') range = 2 * (TWOPI * log2((double)(framewidth+nwid)) - TWOPI * zero_spl); tomwalters@0: if (form_spl == 'L') range = 2 * (framewidth+nwid); tomwalters@0: xorg = 0; tomwalters@0: yorg = 0; tomwalters@0: tomwalters@0: /* Allocate space (first-time call only). */ tomwalters@0: if (first) { tomwalters@0: first = 0; tomwalters@0: frame = NewArray(short, frameheight*(framewidth+nwid), "frame spiral.c" ); tomwalters@0: if (axis_spl) { tomwalters@0: /* Allocate space and compute (X,Y) coords for spiral axis. */ tomwalters@0: /* Resolution of points on axis is 100 times greater, so that */ tomwalters@0: /* inner circuits look smooth. Add 1 as time starts from zero.*/ tomwalters@0: naxis = (100*framewidth+nwid)+1; tomwalters@0: Xaxis = NewArray(short, naxis, "Xaxis, spiral.c" ); tomwalters@0: Yaxis = NewArray(short, naxis, "Yaxis, spiral.c" ); tomwalters@0: naxis = gen_axis(Xaxis,Yaxis,naxis,form_spl,zero_spl); tomwalters@0: } tomwalters@0: if (dotsize_spl) { tomwalters@0: /* Allocate space for coords of spiral points.*/ tomwalters@0: Xspiral = NewArray(short, frameheight * framewidth+nwid, "Xspiral spiral.c" ); tomwalters@0: Yspiral = NewArray(short, frameheight * framewidth+nwid, "Yspiral spiral.c" ); tomwalters@0: Xdot = NewArray(short, frameheight * framewidth+nwid, "Xdot spiral.c" ); tomwalters@0: Ydot = NewArray(short, frameheight * framewidth+nwid, "Ydot spiral.c" ); tomwalters@0: } tomwalters@0: } tomwalters@0: tomwalters@0: CopyArray( image, frame, (framewidth+nwid)*frameheight ) ; tomwalters@0: tomwalters@0: /* Map one SAI frame onto a pulse ribbon, and return the number of 1's */ tomwalters@0: npoints = map_ribbon(frame,dotthresh_spl); tomwalters@0: tomwalters@0: /* Map 1's in a pulse ribbon onto points in spiral coords. */ tomwalters@0: /* Return the number of points in the plot window. */ tomwalters@0: npoints = map_spiral(frame,Xspiral,Yspiral,form_spl,zero_spl); tomwalters@0: tomwalters@0: /* Plot spiral points and axis in window. */ tomwalters@0: if (dotsize_spl) tomwalters@0: plotdots(state->window,Xspiral,Yspiral,Xdot,Ydot,npoints,dotsize_spl-1); tomwalters@0: if (axis_spl) tomwalters@0: Draw(state->window,Xaxis,Yaxis,naxis); tomwalters@0: tomwalters@0: return ; tomwalters@0: } tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: Map an SAI frame onto a "pulse ribbon". tomwalters@0: tomwalters@0: Convert a "frame" into a binary array where a `1' corresponds to a peak in tomwalters@0: the SAI frame, and everwhere else is `0'. This binary array is called a tomwalters@0: pulse ribbon. tomwalters@0: tomwalters@0: Each SAI is a "frame". A frame is a matrix of "frameheight" rows by tomwalters@0: "framewidth" columns. Each row corresponds to a frequency channel, and tomwalters@0: frameheight corresponds to "chans" in image.c. The columns correspond to tomwalters@0: the number of time samples in the "image_duration", and framewidth tomwalters@0: corresponds to framewidth (or imagewidth/chans) in image.c. tomwalters@0: tomwalters@0: The format of data in the input file, (short *frame), is "unmultiplexed". tomwalters@0: This is an SAI frame in row-wise format, with the row corresponding tomwalters@0: to the lowest-frequency channel comming first. tomwalters@0: tomwalters@0: Every point in the sai frame is a coordinate pair: tomwalters@0: (f,t); f=0,1,...,frameheight-1, t=0,1,...,framewidth-1. tomwalters@0: (The lowest frequency channel has the lowest channel number f). tomwalters@0: (The time-origin is the rightmost point, framewidth-1). tomwalters@0: tomwalters@0: Peaks in the SAI frame are maxima along the time axis. Peaks are independent tomwalters@0: of adjacent frequency channels. tomwalters@0: A peak is defined as the highest point between zeroes, which also exceeds a tomwalters@0: threshold. (A zero-crossing between peaks is guaranteed, since the cochleagram tomwalters@0: signal has been half-wave rectified). tomwalters@0: tomwalters@0: The number of peaks found is returned. tomwalters@0: ****************************************************************************/ tomwalters@0: map_ribbon(frame,threshold) tomwalters@0: short *frame; tomwalters@0: int threshold; tomwalters@0: { tomwalters@0: short p, pmax; tomwalters@0: int i, j, k, jk, npoints=0; tomwalters@0: tomwalters@0: /* for each SAI row (each frequency channel) in turn */ tomwalters@0: for (i=0, k=0 ; i=0 ; j--) { /* SSS */ tomwalters@0: /* while zero, advance to next peak */ tomwalters@0: while (j>=(0-1) && frame[j+k]<=0) { /* SSS */ tomwalters@0: frame[j+k]=0; tomwalters@0: j--; tomwalters@0: } tomwalters@0: /* while greater than zero, search for max peak */ tomwalters@0: for (pmax=0 ; j>=(0-1) && (p=frame[j+k])>0 ; j--) { /* SSS */ tomwalters@0: frame[j+k]=0; tomwalters@0: if (p>pmax) { tomwalters@0: pmax=p; jk=j+k; tomwalters@0: } tomwalters@0: } tomwalters@0: /* record max peak if it exceeds the threshold */ tomwalters@0: if (pmax > threshold) { tomwalters@0: frame[jk]=1; tomwalters@0: npoints++; tomwalters@0: } tomwalters@0: } tomwalters@0: return npoints; tomwalters@0: } tomwalters@0: tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * Map 1's in pulse ribbon onto (x,y) spiral coords, and store in the tomwalters@0: * given arrays if the points will appear in the plot window. tomwalters@0: * Return the number of mapped points within the plot window. tomwalters@0: ****************************************************************************/ tomwalters@0: map_spiral(frame,Xarray,Yarray,form_spl,zero_spl) tomwalters@0: short *frame; tomwalters@0: short *Xarray, *Yarray; tomwalters@0: char form_spl; tomwalters@0: float zero_spl; tomwalters@0: { tomwalters@0: float t, f, x, y; tomwalters@0: float r0, r, theta; tomwalters@0: int i, j, k, X, Y; tomwalters@0: int points=0; tomwalters@0: tomwalters@0: /* for each SAI row (each frequency channel) in turn */ tomwalters@0: for (i=0, k=0 ; i=0 ; j--) /* SSS */ tomwalters@0: if (frame[j+k] == 1) { /* for each `1' in pulse ribbon */ tomwalters@0: f = i; tomwalters@0: t = (framewidth-1) - j ; /* time origin is on the right */ tomwalters@0: tomwalters@0: /* calculate polar coords of point, (r,theta): */ tomwalters@0: /* radius r at outer edge of spiral cycle, (ie at min f) */ tomwalters@0: /* radius r0, at same theta, at min f on previous cycle */ tomwalters@0: if (t > 1) theta = TWOPI * (log2(t) - zero_spl); tomwalters@0: else theta = 0; /* special case of 1st cycle */ tomwalters@0: if (form_spl == 'A') { /* spiral == ARCHEMEDIAN */ tomwalters@0: r = theta; tomwalters@0: if (r < 0) r = 0; tomwalters@0: if (r > 0) r0 = r - TWOPI; tomwalters@0: else r0 = 0; /* special case of 1st cycle */ tomwalters@0: } tomwalters@0: if (form_spl == 'L') { /* spiral == LOGARITHMIC */ tomwalters@0: r = t; tomwalters@0: if (t > 1) r0 = r/2; tomwalters@0: else r0 = 0; /* special case of 1st cycle */ tomwalters@0: } tomwalters@0: /* scale radius r for current frequency f (>= min f) */ tomwalters@0: /* scale such that spiral axis occupies its own channel */ tomwalters@0: r = r - ((f+1)/(frameheight+1)) * (r-r0); tomwalters@0: /* convert polar to cartesian coords */ tomwalters@0: x = r*cos(theta); tomwalters@0: y = r*sin(theta); tomwalters@0: /* Transform into pixel coord system and store point */ tomwalters@0: X = px(x); tomwalters@0: Y = py(y); tomwalters@0: if (inbox(X,Y)) { tomwalters@0: Xarray[points] = X; tomwalters@0: Yarray[points++] = Y; tomwalters@0: } tomwalters@0: } tomwalters@0: tomwalters@0: return points; tomwalters@0: } tomwalters@0: tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * Generate points (X,Y coords) for spiral axis and store in given arrays. tomwalters@0: * Return the number of points stored. tomwalters@0: ****************************************************************************/ tomwalters@0: int gen_axis(Xarray,Yarray,naxis,form_spl,zero_spl) tomwalters@0: short *Xarray, *Yarray; tomwalters@0: int naxis; tomwalters@0: char form_spl; tomwalters@0: float zero_spl; tomwalters@0: { tomwalters@0: int t; tomwalters@0: float r, theta; tomwalters@0: int X, Y; tomwalters@0: float x, y; tomwalters@0: int points=0; tomwalters@0: tomwalters@0: /* Calculate axis points in data coord system */ tomwalters@0: for (t=100 ; t < naxis ; t++) { tomwalters@0: theta = TWOPI * (log2((float)t/100) - zero_spl); tomwalters@0: if (theta < 0) theta = 0; tomwalters@0: if (form_spl == 'A') /* spiral == Achemedian */ tomwalters@0: r = theta; tomwalters@0: else /* spiral == Logarithmic */ tomwalters@0: r = (float)t/100; tomwalters@0: x = r*cos(theta); tomwalters@0: y = r*sin(theta); tomwalters@0: /* Transform into pixel coord system and store */ tomwalters@0: X = px(x); tomwalters@0: Y = py(y); tomwalters@0: if (inbox(X,Y)) { tomwalters@0: Xarray[points] = X; tomwalters@0: Yarray[points++] = Y; tomwalters@0: } tomwalters@0: tomwalters@0: } tomwalters@0: return points; tomwalters@0: } tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: /**************************************************************************** tomwalters@0: * Plot dots of size d. tomwalters@0: * tomwalters@0: * Plotting in a given window uses the glib routines Draw and Plot. tomwalters@0: * These routines, and associated structures, are declared in windows.h. tomwalters@0: * The WindowObject structure has an "entries" pointer to a WindowEntries tomwalters@0: * structure, and this holds pointers to a set of functions which operate on tomwalters@0: * the window. The window functions, and the routine newDisplayWindow which tomwalters@0: * allocates a new WindowObject, are defined in X.c. tomwalters@0: * tomwalters@0: * void Draw(win,X,Y,points) tomwalters@0: * Takes two arrays of shorts, one with x-coords the other with y-coords. tomwalters@0: * The number of x's = the number of y's = "points". tomwalters@0: * The user must ensure that the data coords are scaled into the pixel tomwalters@0: * coord system. This is as follows: (1,1) is bottom-left corner, and tomwalters@0: * (width,height) is top-right corner. tomwalters@0: * void Plot(win,X,Y,points) tomwalters@0: * Same as Draw, but plots points. tomwalters@0: ****************************************************************************/ tomwalters@0: plotdots(win,Xspiral,Yspiral,X,Y,npoints,d) tomwalters@0: WindowObject win; tomwalters@0: short *Xspiral, *Yspiral; tomwalters@0: short *X, *Y; tomwalters@0: int npoints, d; tomwalters@0: { tomwalters@0: int x, y, i; tomwalters@0: tomwalters@0: if (d==0) /* dotsize = 1 */ tomwalters@0: Plot(win,Xspiral,Yspiral,npoints); tomwalters@0: else { /* dotsize > 1 */ tomwalters@0: for (x=(-d) ; x<=d ; x++) tomwalters@0: for (y=(-d) ; y<=d ; y++) { tomwalters@0: for (i=0 ; i