diff model/spiral.c @ 0:5242703e91d3 tip

Initial checkin for AIM92 aimR8.2 (last updated May 1997).

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