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

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

/*
    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);
	    }
    }
}