 /*
  * Khoros: $Id: gisutils.c,v 1.2 1992/03/20 23:35:12 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: gisutils.c,v 1.2 1992/03/20 23:35:12 dkhoros Exp $";
#endif

 /*
  * $Log: gisutils.c,v $
 * Revision 1.2  1992/03/20  23:35:12  dkhoros
 * VirtualPatch5
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1990, 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 too 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.
 *----------------------------------------------------------------------
 */

#include "unmcopyright.h"	 /* Copyright 1990 by UNM */

#include "vinclude.h"

/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
   >>>>                                                       <<<<
   >>>>     file name:  gisutils.c                            <<<<
   >>>>                                                       <<<<
   >>>>   description: utilities for GIS convertions          <<<<
   >>>>                                                       <<<<
   >>>>      routines: utm2arc()                              <<<<
   >>>>                arc2pixel()                            <<<<
   >>>>                pixel2arc()                            <<<<
   >>>>                readfixed()                            <<<<
   >>>>                nextline()                             <<<<
   >>>>                                                       <<<<
   >>>>                                                       <<<<
   >>>>                                                       <<<<
   >>>>                                                       <<<<
   >>>>                                                       <<<<
   >>>>                                                       <<<<
   >>>>                                                       <<<<
   >>>>                                                       <<<<
   >>>>                                                       <<<<
   >>>> modifications: added utm2arc, arc2pixel, pixel2arc    <<<<
   >>>>                10/7/90                                <<<<
   >>>>                added readfixed, nextline              <<<<
   >>>>                2/19/91                                <<<<
   >>>>                added assign_locations                 <<<<
   >>>>                2/23/91                                <<<<
   >>>>                                                       <<<<
   >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */

int utm2arc(); 
int arc2pixel();
int pixel2arc();
int ll_to_utm();
int dlg_line();
int dlg_point();
int readfixed();
void nextline();
int assign_locations();
int get_bit_pixel();
int put_bit_pixel();

/****************************************************************
*
* Routine: utm2arc
*
* Purpose: Given the coordinates of the four corners of a square in
*          both UTM and arc-seconds, this routines converts the UTM
*          coordinates of a point within the square to arc coordinates.
*
*   Input: The inputs are the coordinates of the four corners of
*          the square in both utm and arc.  Utm_west_edge is the
*          x utm coordinate of the west side of the square, utm_south_edge
*          is the y utm coordinate of the south side of the square, and
*          so forth.  Also input are x_utm and y_utm, which are the
*          utm coordinates of a point within the square.
*
*	   utm_west_edge - the utm x coordinate of the west edge
*	   utm_south_edge - the utm y coordinate of the south edge
*	   utm_east_edge - the utm x coordinate of the east edge
*	   utm_north_edge - the utm y coordinate of the north edge
*
*	   arc_west_edge - the longitude in seconds of the west edge 
*	   arc_south_edge - the latitude in seconds of the south edge 
*	   arc_east_edge - the longitude in seconds of the east edge 
*	   arc_north_edge - the latitude in seconds of the north edge 
*
*	   x_utm - the UTM x coordinate to be converted to an arc x coordinate
*	   y_utm - the UTM y coordinate to be converted to an arc y coordinate
*
*	   x_arc - the converted arc x coordinate
*	   y_arc - the converted arc y coordinate
*
*  Output: The outputs are x_arc and y_arc, which are the geographical
*          coordinates in arc which correspond to x_utm and y_utm.
*          This function returns 1 on success and 0 on failure.
*
* Comment: Note that the inputs and outputs are double.  All geographical
*          values are in arc-seconds.
*
*  Author:  Per Lysne  7/9/90
*
****************************************************************/
int
utm2arc (utm_west_edge, utm_south_edge, utm_east_edge, utm_north_edge,
             arc_west_edge, arc_south_edge, arc_east_edge, arc_north_edge,
             x_utm, y_utm, x_arc, y_arc)
double utm_west_edge, utm_south_edge, utm_east_edge, utm_north_edge;
double arc_west_edge, arc_south_edge, arc_east_edge, arc_north_edge;
double x_utm, y_utm;
double *x_arc, *y_arc;
{
    double utm_ns_run, utm_ew_run, arc_ns_run, arc_ew_run;
    double utm_x_to_west_edge, utm_y_to_south_edge;
    double arc_x_to_west_edge, arc_y_to_south_edge;

    /*
     * These are the delta x and y of the utm coordinates across the area.
     */
    utm_ns_run     = utm_north_edge - utm_south_edge;
    utm_ew_run     = utm_east_edge - utm_west_edge;

    /*
     * These are delta x and y of number of arc-seconds across the area.
     */
    arc_ns_run     = arc_north_edge - arc_south_edge;
    arc_ew_run     = arc_east_edge - arc_west_edge;
    if ((arc_ns_run == 0) || (arc_ew_run == 0)) {
        return (0);
    }

    /*
     * These are the utm offsets into the area from the south and west edges.
     */
    utm_x_to_west_edge = x_utm - utm_west_edge;
    utm_y_to_south_edge  = y_utm - utm_south_edge;

    /*
     * Set up a ratio and find the arc offsets from south and west edges;
     */
    arc_x_to_west_edge = (utm_x_to_west_edge / utm_ew_run) * arc_ew_run;
    arc_y_to_south_edge  = (utm_y_to_south_edge  / utm_ns_run) * arc_ns_run;

    /*
     * Add the offsets to the arc values of the south and west edges.
     */
    *x_arc = (float) (arc_west_edge + arc_x_to_west_edge);
    *y_arc = (float) (arc_south_edge  + arc_y_to_south_edge);

    return (1);
}



/****************************************************************
*
*    File: arc2pixel.c
*
* Routine: arc2pixel
*
* Purpose: Given the geographical coordinates of a square region, and
*          the size of an image to cover the region, this function
*          tells you which pixel in the image to plot corresponding
*          to a location within the geographical region defined.
*
*   Input: The first four inputs are define the geographical region.
*          Arc_west_edge is the arc-longitude of the west edge of the
*          region, arc_south_edge is the arc-latitude of the south
*          edge of the region, and so on.  Rows and cols give the
*          size of the image covering the region.  X_arc and y_arc
*          give the location of a point within the region which we
*          would like to plot in the image.
*
*	   arc_west_edge - the longitude in seconds of the west edge 
*	   arc_south_edge - the latitude in seconds of the south edge 
*	   arc_east_edge - the longitude in seconds of the east edge 
*	   arc_north_edge - the latitude in seconds of the north edge 
*
*	   rows - number of rows in the image covering the desired geog. area
*	   cols - number of cols in the image covering the desired geog. area
*
*	   x_arc - the converted arc x coordinate
*	   y_arc - the converted arc y coordinate
*
*	   x_pixel - the x coord of the pixel corresponding to the input
*		     geographical location returned
*	   y_pixel - the y coord of the pixel corresponding to the input
*		     geographical location returned
*
*  Output: The outputs are x_pixel and y_pixel which give the location
*          of the pixel in the image that we should plot corresponding
*          to the goegraphical location given by x_arc and y_arc.
*          This function returns 1 on success and 0 on failure.
*
* Comment: All geographical values are in arc-seconds and are double.
*
*  Author: Per Lysne  8/26/90
*
****************************************************************/
/* If we get an arc outside of this range it is wrong. */
#define MAX_X_ARC 60.0 * 60.0 * 180.0
#define MAX_Y_ARC 60.0 * 60.0 * 90.0

int
arc2pixel (arc_west_edge, arc_south_edge, arc_east_edge, arc_north_edge,
           rows, cols,
           x_arc, y_arc,
           x_pixel, y_pixel)
double arc_west_edge, arc_south_edge, arc_east_edge, arc_north_edge;
int rows, cols;
double x_arc, y_arc;
int *x_pixel, *y_pixel;
{
    double arc_ew_run, arc_ns_run, arc_x_offset, arc_y_offset;
    double pixel_x_offset, pixel_y_offset;

    /*
     * Error check the inputs.
     */
    if (arc_west_edge > MAX_X_ARC || arc_west_edge < -MAX_X_ARC) {
        fprintf (stderr, "error in arc2pixel - ");
        fprintf (stderr, "arc_west_edge is outside of allowable range\n");
        return (0);
    } 
    if (arc_east_edge > MAX_X_ARC || arc_east_edge < -MAX_X_ARC) {
        fprintf (stderr, "error in arc2pixel - ");
        fprintf (stderr, "arc_east_edge is outside of allowable range\n");
        return (0);
    }
    if (arc_south_edge > MAX_Y_ARC || arc_south_edge < -MAX_Y_ARC) {
        fprintf (stderr, "error in arc2pixel - ");
        fprintf (stderr, "arc_south_edge is outside of allowable range\n");
        return (0);
    }
    if (arc_north_edge > MAX_Y_ARC || arc_north_edge < -MAX_Y_ARC) {
        fprintf (stderr, "error in arc2pixel - ");
        fprintf (stderr, "arc_north_edge is outside of allowable range\n");
        return (0);
    }

    if (rows <= 0 || cols <= 0) {
        fprintf (stderr, "error in arc2pixel - ");
        fprintf (stderr, "either rows or cols inputs was < 0\n");
        return (0);
    }

    if (x_arc > arc_east_edge || x_arc < arc_west_edge) {
        printf ("warning in arc2pixel - ");
        printf ("x_arc is outside of the region defined\n");
    }
    if (y_arc > arc_north_edge || y_arc < arc_south_edge) {
        printf ("warning in arc2pixel - ");
        printf ("y_arc is outside of the region defined\n");
    }

    /*
     * These are delta x and y of number of arc-seconds across the area.
     */
    arc_ew_run = arc_east_edge - arc_west_edge;
    if (arc_ew_run <= 0.0) {
        fprintf (stderr, "error in arc2pixel - ");
        fprintf (stderr, "width of region given is <= 0\n");
        return (0);
    }
    arc_ns_run = arc_north_edge - arc_south_edge;
    if (arc_ns_run <= 0.0) {
        fprintf (stderr, "error in arc2pixel - ");
        fprintf (stderr, "height of region given is <= 0\n");
        return (0);
    }

    /*
     * Find the offsets in arc from the west and south edges of the area.
     */
    arc_x_offset = x_arc - arc_west_edge;
    arc_y_offset = y_arc - arc_south_edge;

    /*
     * Set up a ratio and find the pixel offsets from south and west edges;
     */
    pixel_x_offset = (arc_x_offset / arc_ew_run) * (double)cols;
    pixel_y_offset = (arc_y_offset / arc_ns_run) * (double)rows;

    /*
     * Add the offsets to the arc values of the south and west edges.  
     * Round off the float offsets to integer pixels, and remember that
     * in an image the 0th row is at the top so y_pixel is inverted.
     */
    if ((int)(pixel_x_offset-.5) == (int)pixel_x_offset) {
        *x_pixel = (int)pixel_x_offset +1;
    }
    else {
        *x_pixel = (int)pixel_x_offset;
    }
    if ((int)(pixel_y_offset-.5) == (int)pixel_y_offset) {
        *y_pixel = (int)pixel_y_offset +1;
    }
    else {
        *y_pixel = (int)pixel_y_offset;
    }
    *y_pixel = rows-(*y_pixel);

    /*
     * Make sure that the pixel x and y are within the image entered.
     */
    if (*x_pixel < 0) {
        *x_pixel=0;
    }
    if (*x_pixel > cols-1) {
        *x_pixel=cols-1;
    }
    if (*y_pixel < 0) {
        *y_pixel=0;
    }
    if (*y_pixel > rows-1) {
        *y_pixel=rows-1;
    }

    return (1);
}



/****************************************************************
*
*    File: pixel2arc.c
*
* Routine: pixel2arc
*
* Purpose: Given the geographical coordinates of a square region and
*          the size of an image that is to cover that regoin, this
*          function returns the geographical coordinates corresponding
*          to a pixel inside that image.
*
*   Input: The first four inputs are the geographical coordinates of
*          the edges of the square region in arc-secs.  Arc_west_edge
*          is the arc - longitude of the west edge of the area,
*          arc_south_edge is the arc - latitude of the south edge of
*          the area, and so on.  Rows and cols give the size of the
*          image which is covering the area.  X_pixel and Y_pixel are
*          the x and y which specify a pixel within the image.
*
*	   arc_west_edge - the longitude in seconds of the west edge 
*	   arc_south_edge - the latitude in seconds of the south edge 
*	   arc_east_edge - the longitude in seconds of the east edge 
*	   arc_north_edge - the latitude in seconds of the north edge 
*
*	   rows - number of rows in the image covering the desired geog. area
*	   cols - number of cols in the image covering the desired geog. area
*
*	   x_pixel - the x coord of the pixel to be converted to geog. coords.
*	   y_pixel - the y coord of the pixel to be converted to geog. coords.
*
*	   x_arc - the converted arc x coordinate returned
*	   y_arc - the converted arc y coordinate returned
*
*
*  Output: The outputs are x_arc and y_arc which give the arc lon and
*          lat of the pixel given by x_pixel and y_pixel.  This
*          function returns 1 on success and 0 on failure.
*
* Comment: All geographical values are in arc-seconds, and are double.
*
*  Author: Per Lysne  8/26/90
*
****************************************************************/
/* If we get an arc outside of this range it is wrong. */
#define MAX_X_ARC 60.0 * 60.0 * 180.0
#define MAX_Y_ARC 60.0 * 60.0 * 90.0

int
pixel2arc (arc_west_edge, arc_south_edge, arc_east_edge, arc_north_edge,
           rows, cols,
           x_pixel, y_pixel,
           x_arc, y_arc)
double arc_west_edge, arc_south_edge, arc_east_edge, arc_north_edge;
int rows, cols;
int x_pixel, y_pixel;
double *x_arc, *y_arc;
{
    double arc_ew_run, arc_ns_run, arc_x_offset, arc_y_offset;

    /*
     * Error check the inputs.
     */
    if (arc_west_edge > MAX_X_ARC || arc_west_edge < -MAX_X_ARC) {
        fprintf (stderr, "error in pixel2arc - ");
        fprintf (stderr, "arc_west_edge is outside of allowable range\n");
        return (0);
    }
    if (arc_east_edge > MAX_X_ARC || arc_east_edge < -MAX_X_ARC) {
        fprintf (stderr, "error in pixel2arc - ");
        fprintf (stderr, "arc_east_edge is outside of allowable range\n");
        return (0);
    }
    if (arc_south_edge > MAX_Y_ARC || arc_south_edge < -MAX_Y_ARC) {
        fprintf (stderr, "error in pixel2arc - ");
        fprintf (stderr, "arc_south_edge is outside of allowable range\n");
        return (0);
    }
    if (arc_north_edge > MAX_Y_ARC || arc_north_edge < -MAX_Y_ARC) {
        fprintf (stderr, "error in pixel2arc - ");
        fprintf (stderr, "arc_north_edge is outside of allowable range\n");
        return (0);
    }

    if (rows <= 0 || cols <= 0) {
        fprintf (stderr, "error in pixel2arc - ");
        fprintf (stderr, "either rows or cols inputs was < 0\n");
        return (0);
    }

    if (x_pixel < 0 || x_pixel >= cols) {
        fprintf (stderr, "error in pixel2arc - ");
        fprintf (stderr, "x_pixel was < 0 or greater than no. cols\n");
        return (0);
    }
    if (y_pixel < 0 || y_pixel >= rows) {
        fprintf (stderr, "error in pixel2arc - ");
        fprintf (stderr, "y_pixel was < 0 or greater than no. cols\n");
        return (0);
    }

    /*
     * These are delta x and y of number of arc-seconds across the area.
     */
    arc_ew_run = arc_east_edge - arc_west_edge;
    if (arc_ew_run <= 0.0) {
        fprintf (stderr, "error in pixel2arc - ");
        fprintf (stderr, "arc_east_edge <= arc_west_edge\n");
        return (0);
    }
    arc_ns_run = arc_north_edge - arc_south_edge;
    if (arc_ns_run <= 0.0) {
        fprintf (stderr, "error in pixel2arc - ");
        fprintf (stderr, "arc_north_edge <= arc_south_edge\n");
        return (0);
    }

    /*
     * Find the offsets in arc from the west and south edges of the area.
     */
    arc_x_offset = ((double)x_pixel/(double)cols) * arc_ew_run;
    arc_y_offset = ((double)(rows-y_pixel)/(double)rows) * arc_ns_run;

    /*
     * Add the offsets to the arc values of the south and west edges.  
     */
    *x_arc = arc_west_edge + arc_x_offset;
    *y_arc = arc_south_edge + arc_y_offset;

    return (1);
}



/****************************************************************
*
* Routine Name: ll_to_utm
*
* Purpose: This routine performs a trigonometric conversion between
*        geographical coordinates in latitude and longitude to utm
*        coordinates.
*
* Input: lon, lat - the longitude and latitude of a point to convert to utm
*        sec_flag - set to 1 to specify that the lon and lat are in seconds
*        zone_flag - to to 1 if the user is inputing the utm zone, otherwise
*                   it will be computed from the longitude
*
* Output: easting, northing - the converted utm coordinates
*         zone - the utm zone that the easting and northing belong to 
*
* Written By:  Per Lysne  3/1/91
*
****************************************************************/
#define degrees_to_radians 0.017453

int ll_to_utm (lon, lat, sec_flag, zone_flag, easting, northing, zone)
double lon, lat;
int sec_flag,zone_flag;
double *easting, *northing;
int *zone;
{
    double offset, central_meridian;
    double dlamr, phir;
    double phic, sinp, cosp, tanp, esqr, cosp3, dlam2, dlam3, dlam4;
    double false_easting, false_northing;
    double rp;
    double fac1, fac2, fac3, dm;
    double eta2, facx, facy;
    double x, y;

    /*
     * Use the Clark-66 spheroid.
     */
    double a = 6378206.0;
    double e = 0.082269;

    /*
     * If the seconds flag is set then convert the lon and lat from
     * seconds to degrees.
     */
    if (sec_flag == 1) {
        lon = lon / 3600.0;
        lat = lat / 3600.0;
    }

    /*
     * Make sure the input longitude is in the correct range.
     */
    if ((lon <= -180.0) || (lon > 180.0)) {
        fprintf (stderr, "error in ll_to_utm: lon must be in the range ");
        fprintf (stderr, "-180 < lon <= 180 or its equivalent in seconds\n");
        return (0);
    }

    /*
     * Make sure the input latitude is in the correct range.
     */
    if ((lat < -90.0) || (lat > 90.0)) {
        fprintf (stderr, "error in ll_to_utm: lat must be in the range ");
        fprintf (stderr, "-90 < lon < 90 or its equivalent in seconds\n");
        return (0);
    }

    /*
     * If the calling program is leaving it up to us to compute the utm
     * zone, zone_flg will be 0.
     */
    if (zone_flag == 0) {

        /*
         * Compute the zone that this longitude value requires, and the central
         * meridian for that zone.
         */
        if (lon <= 0) {
            *zone = 30 + (int)(lon/6.0);
        }
        else {
            *zone = 31 + (int)(lon/6.0);
        }
    }
    central_meridian = -183.0 + (double)(6*(*zone));

    /*
     * Find the offset of the lon from the central meridian of the zone.  Set
     * the lon to this offset because the rest of the alogrithim needs it
     * that way.
     */
    offset = lon - central_meridian;
    lon = offset;

    /*
     * Convert the lon and lat to radians.
     */
    dlamr = lon * degrees_to_radians;
    phir = lat * degrees_to_radians;

    /*
     * The rest of this function is taken from page 342 of Map Projections -
     * Theory and Applications, by Frederick Pearson II.
     */
    phic = fabs (phir);
    sinp = sin (phic);
    cosp = cos (phic);
    tanp = tan (phic);
    esqr = e*e;
    cosp3 = cosp*cosp*cosp;
    dlam2 = dlamr * dlamr;
    dlam3 = dlam2 * dlamr;
    dlam4 = dlam2 * dlam2;

    /*
     * Set the false easting and northing.
     */
    false_easting = 500000.0;
    false_northing = 0.0;
    if (phir < 0.0)
        false_northing = 1.0e+7;

    /*
     * Find the radius of curvature of the prime vertical.
     */
    rp = a /sqrt(1.0 - (esqr*sinp*sinp));

    /*
     * Find the distance along the meridian from the equator to the lat phir.
     */
    fac1 = (1.0 - (0.25*esqr) - (0.046875*esqr*esqr)) * phic;
    fac2 = ((0.375*esqr) + (0.093758*esqr*esqr)) * sin(2.0*phic);
    fac3 = (0.058594*esqr*esqr) * sin(4.0*phic);
    dm   = a * (fac1 - fac2 + fac3);

    /*
     * Values defined to be used in the next equation.
     */
    eta2 = (esqr * cosp * cosp)/(1.0 - esqr);
    facx = 1.0 - (tanp*tanp) + eta2;
    facy = 5.0 - (tanp*tanp) + (9.0*eta2);

    /*
     * Find the easting and northing.
     */
    x = rp * ((dlamr*cosp) + (dlam3*cosp3*facx/6.0));
    y = dm + rp * ((dlam2*sinp*cosp/2.0) + (dlam4*sinp*cosp3*facy/24.0));

    /*
     * Reflect about the equator - remember that we have been using the
     * absolute value of the latitude.
     */
    if (phir < 0.0)
       y = -y;

    /*
     * Adjust for the false easting and northing.
     */
    *easting = 0.9996 * x + false_easting;
    *northing = 0.9996 * y + false_northing;

    /*
     * Wasn't that easy!
     */
    return (1);
}



/****************************************************************
*
* Routine Name: dlg_line
*
* Purpose: To plot a line in an image band.
*
* Input: A pointer to the image band to plot in, the number of rows
*        and columns in the image band, and the x and y coordinates
*        of the beginning and end points of the line.
*
*	 band - pointer to the image band to plot line in.
*	 rows - number of rows in the image band.
*	 cols - number of cols in the image band.
*	 x1, y1 - starting coords of line in image
*	 x2, y2 - ending coords of line in image
*
* Output: The line is plotted in the image band.  Returns 1 on success
*         and 0 on failure.
*
* Written By:  Per Lysne  10/7/90
*
****************************************************************/
int
dlg_line (band, rows, cols, x1, y1, x2, y2)
char *band;
int rows, cols;
int x1, y1, x2, y2;
{
    int xinc, yinc, row, col;
    float slope, intercept;

    if (band == NULL) {
        fprintf (stderr, "error in dlg_line - bad image band\n");
        return (0);
    }

    if (rows <= 0) {
        fprintf (stderr, "error in dlg_line - bad number of rows\n");
        return (0);
    }
    if (cols <= 0) {
        fprintf (stderr, "error in dlg_line - bad number of columns\n");
        return (0);
    }

    if ((x1 < 0) || (x1 > cols-1)) {
        fprintf (stderr, "error in dlg_line - bad x1 entered\n");
        return (0);
    }

    if ((x2 < 0) || (x2 > cols-1)) {
        fprintf (stderr, "error in dlg_line - bad x2 entered\n");
        return (0);
    }

    if ((y1 < 0) || (y1 > rows-1)) {
        fprintf (stderr, "error in dlg_line - bad y1 entered\n");
        return (0);
    }

    if ((y2 < 0) || (y2 > rows-1)) {
        fprintf (stderr, "error in dlg_line - bad y2 entered\n");
        return (0);
    }

    xinc = 1;
    if (x2 < x1) {
        xinc = -1;
    }
    yinc = 1;
    if (y2 < y1) {
        yinc = -1;
    }

    if (x1 == x2) {
        for (row=y1; row!=y2; row=row+yinc) {
            dlg_point (band, rows, cols, x1, row);
        }
    }
    else {
        slope = (float)(y2-y1)/(float)(x2-x1);
        intercept = ((float)x1 * (float)y2 - (float)x2 * (float)y1) /
                    (float)(x1-x2);

        if ((slope < 1.0) && (slope > -1.0)) {
            for (col=x1; col!=x2; col=col+xinc) {
                row =  (int) ((float)col * slope +intercept);
                dlg_point (band, rows, cols, col, row);
            }
        }
        else {
            for (row=y1; row!=y2; row=row+yinc) {
                col = (int) (((float)row - intercept) / slope); 
                dlg_point (band, rows, cols, col, row);
            }
        }
    }
    dlg_point (band, rows, cols, x2, y2);

    return (1);
}



/****************************************************************
*
* Routine Name: dlg_point
*
* Purpose: To plot a point in an image band.
*
* Input: A pointer to the image band to plot in, the number of rows
*        and columns in the image band, and the x and y coordinates
*        of the point.
*
*	 band - pointer to the image band to plot point in.
*	 rows - number of rows in the image band.
*	 cols - number of cols in the image band.
*	 x, y - coords of point in image
*
* Output: The point is plotted in the image band.  Returns 1 on success
*         and 0 on failure.
*
* Written By:  Per Lysne  10/7/90
*
****************************************************************/
int
dlg_point (band, rows, cols, x, y) 
char *band;
int rows, cols;
int x, y;
{
    /*
     * Must convert the data to unsigned to avoid complaints about loosing
     * data from Saber.
     */
    unsigned char *imageband;

    if (band == NULL) {
        fprintf (stderr, "error in dlg_point - bad image band\n");
        return (0);
    }

    if (rows <= 0) {
        fprintf (stderr, "error in dlg_point - bad no. of rows\n");
        return (0);
    }

    if (cols <= 0) {
        fprintf (stderr, "error in dlg_point - bad no. of cols\n");
        return (0);
    }

    if ((x < 0) || (x > cols-1)) {
        fprintf (stderr, "error in dlg_point - bad x to plot\n");
        return (0);
    }

    if ((y < 0) || (y > rows-1)) {
        fprintf (stderr, "error in dlg_point - bad y to plot\n");
        return (0);
    }

    imageband = (unsigned char *) band;

    imageband[(y*cols+x)/8] = imageband[(y*cols+x)/8] | (1 << ((y*cols+x)%8));

    return (1);
}



/****************************************************************
*
* Routine: readfixed
*
* Purpose: This function reads fixed format data from a file.
*          For instance, if you need to read a six character integer
*          from a file and no delimiting characters are used, you
*          need this routines.  All this routine does is read the
*          specified number of characters from your file, and then
*          do a sscanf from the string to convert to the proper data
*          type.
*
*   Input: The inputs are a file pointer, the number, or length of
*          the fixed format to read, and a C format string to tell
*          the function what data type is to be read.  The pointer,
*          variable, tells the function where to put the number it
*          just read.
*
*	   file - open stream to an ascii data file
*	   number - length of field to read
*	   format - C format string for reading the data ie, "%f"
*	   variable - a pointer to the variable of the appropriate
*		      type to read the data into.
*    
*  Output: The output is the data read and is pointed to by 'variable'.
*          This function returns 1 upon success and 0 on failure.
*    
*  Author: Per Lysne  7/30/90
*    
****************************************************************/

#define STRINGLENGTH 256

int
readfixed (file, number, format, variable)
FILE *file;
int number;
char *format;
char *variable;
{
    int i, done;
    char string[STRINGLENGTH];

    if (file == NULL) {
        fprintf (stderr, "error in readfixed - bad file input\n");
        return (0);
    }

    if ((number>=STRINGLENGTH-2) || (number<=0)) {
       fprintf (stderr,"error in readfixed - bad no. of characters to read\n");
        return (0);
    }

    if (format == NULL) {
        fprintf (stderr, "error in readfixed - bad format string input\n");
        return (0);
    }
 
    if (variable == NULL) {
        fprintf (stderr, "error in readfixed - bad variable input\n");
        return (0);
    }

    for (i=0; i<=number; i++) string[i] = '\0';

    /*
     * Read number characters from the file into string.  If this
     * number of characters cannot be read, return with an error.
     */
    i=0;
    done=0;
    while (!done) {
        string[i] = fgetc(file);
        if (string[i] == '\n') {
            string[i]='\0';
            ungetc('\n', file);
            done=1;
        }
        if (i>=number-1) {
            done=1;
        }
        i++;
    }
    string[i] = '\0';

    /*
     * In some data, a 'D' is used to denote an exponent, but C requires
     * an 'e' or 'E' to denote an exponent.  This section changes all 'd'
     * characters to 'e' if the data is being read as a double.
     */
    if (strcmp (format, "%lf") == 0) {
        for (i=0; i<number; i++) {
            if (string[i] == 'D' || string[i] == 'd') {
                string[i] = 'e';
            }
        }
    }

    /*
     * Finally, read the data out of the string.
     */
    if (strcmp(format, "%s") == 0) {
        if (strncpy(variable, string, number) == NULL) {
            return (0);
        }
    }
    else {
        if (sscanf (string, format, variable) != 1) {
            return (0);
        }
    }

    /*
     * Return successfully.
     */
    return (1);
}



/****************************************************************
*
* Routine: nextline
*
* Purpose: This routine moves the file pointer to the beginning of
*          the next line of input.
*
*   Input: A file pointer.
*
*  Output: File pointer is moved.
*
*  Author: Per Lysne  12/24/90
*
****************************************************************/
void nextline(file)
FILE *file;
{
    char c;

    c=fgetc(file);
    while ((c!='\n') && (c!=EOF))
        c=fgetc(file);
}



/****************************************************************
*
* Routine: assign_locations
*
* Purpose: To assign explicit location data to an image.
*
*   Input: img - pointer to the viff image
*	   x_sw - the x coord of the left side of the image in the desired 
*		  reference system.
*	   y_sw - the y coord of the bottom of the image in the desired 
*		  reference system.
*	   x_ne - the x coord of the right side of the image in the desired 
*		  reference system.
*	   y_ne - the y coord of the top of the image in the desired 
*		  reference system.
*
*  Output: the input image has explicit locations are filled in 
*
*  Author: Per Lysne  2/23/91
*
****************************************************************/
int assign_locations (img, x_sw, y_sw, x_ne, y_ne)
struct xvimage *img;
double x_sw, y_sw, x_ne, y_ne;
{
    int i, j, rows, cols, index;
    float *x_loc_ptr, *y_loc_ptr;
    float x_base, y_base, delta_x, delta_y, x_inc, y_inc;
    float x_loc, y_loc;

    /*
     * Error Check the inputs.
     */
    if (img == NULL) {
        fprintf (stderr, "error in assign_locations: bad image input\n");
        return (0);
    }
    if (img->location == NULL) {
        fprintf (stderr, "error in assign_locations: bad locations input\n");
        return (0);
    }
    if (x_sw == x_ne) {
        fprintf (stderr, "error in assign_locations: x range must not be 0\n");
        return (0);
    }
    if (y_sw == y_ne) {
        fprintf (stderr, "error in assign_locations: y range must not be 0\n");
        return (0);
    }

    /*
     * Useful pointers.
     */
    rows = img->col_size;
    cols = img->row_size;
    x_loc_ptr = img->location;
    y_loc_ptr = &(img->location[rows*cols]);

    /*
     * Set up a ratio between the total x and y length of the image and
     * the dimension of a single pixel.
     */
    x_base = (float) x_sw;
    delta_x = (float) (x_ne - x_sw);
    x_inc = delta_x / (float) cols;
    y_base = (float) y_sw;
    delta_y = (float) (y_ne - y_sw);
    y_inc = delta_y / (float) rows;

    /*
     * Assign the locations.
     */
    index=0;
    for (i=rows-1; i>=0; i--) {
        y_loc = (float) i * y_inc + y_base;
        for (j=0; j<cols; j++) {
            x_loc = (float) j * x_inc + x_base;
            x_loc_ptr[index] = x_loc;
            y_loc_ptr[index] = y_loc;

            index++;
        }
    }

    return (1);
}



/****************************************************************
*
* Routine: cut_locations
*
* Purpose: This function is designed to be used to create overlay images.
*          If you have an image that has been created from a vector file,
*          such as dlg format, the image will contain some lines and a
*          lot of blank space between the lines.  The blank space means
*          that there is no data for the blank areas, and so there is no
*          point in storing them.  The input to this fuction is an image
*          in a structure xvimage.  The image must contain a full set of
*          explicit locations (i.e. it can't have gone through this
*          function before) and it must be evenly sampled on a rectangular
*          grid.  The output will be the same image except that the blank
*          areas, both image data and locations, will have been removed,
*          and only those areas containing data are saved.  Note that
*          the output can no longer be displayed as an image, but requires
*          a plotting package such as xprism3.
*
*   Input: A pointer to a xvimage structure - the image to be thinned.
*
*  Output: The same pointer - the image has been thinned.
*
*  Author: Per Lysne  3/4/91
*
****************************************************************/
int
cut_locations (img)
struct xvimage *img;
{
    unsigned char *imagedata, *new_imagedata;
    int index;
    int band, rows, cols, num_data_bands;
    int data_count, i, j, data, data_found;
    int image_count, location_count;
    float *location, *x_loc_ptr, *y_loc_ptr;
    float *new_location, *new_x_loc_ptr, *new_y_loc_ptr;

    /*
     * Error check the inputs.
     */
    if (img == NULL) {
        fprintf (stderr, "error in cut_locations: bad image input\n");
        return (0);
    }
    if (img->data_storage_type != VFF_TYP_BIT) {
        fprintf (stderr, "error in cut_locations: image must be bit\n");
        return (0);
    }
    if (img->location_type != VFF_LOC_EXPLICIT) {
        fprintf (stderr, "error in cut_locations: image must have locations\n");
        return (0);
    }
    if (img->location == NULL) {
        fprintf (stderr, "error in cut_locations: bad locations input\n");
        return (0);
    }

    /*
     * Some useful values.
     */
    rows=img->col_size;
    cols=img->row_size;
    num_data_bands=img->num_data_bands;
    imagedata= (unsigned char *) img->imagedata;
    location=img->location;
    x_loc_ptr=location;
    y_loc_ptr= &(location[rows*cols]);

    /*
     * Count the pixels containing data in the image.  If the same pixel
     * contains data in several different bands, it is only counted once.
     */
    data_count=0;
    for (i=0; i<rows; i++) {
        for (j=0; j<cols; j++) {

            data_found=0;
            for (band=0; band<num_data_bands; band++) {
                if (!get_bit_pixel(imagedata, band, rows, cols, j, i, &data)) {
                    fprintf (stderr,
                             "error in cut_locations: get_bit_pixel failed\n");
                    return (0);
                }

                if (data != 0)
                    data_found=1;
            }

            if (data_found)
                data_count++;
 
        }
    }
    
    /*
     * These are the amounts of space required for the image and location
     * data.
     */
    image_count = (int) ((data_count+7)/8);
    location_count = image_count*8;

    /*
     * This is the space that will be used for the final image data and
     * locations.
     */
    new_imagedata = (unsigned char *)
                    calloc (image_count, num_data_bands*sizeof(unsigned char));
    if (new_imagedata == NULL) {
        fprintf (stderr,
               "error in cut_locations: calloc failed for new imagedata\n");
        return (0);
    }

    new_location = (float *) calloc (location_count, 2*sizeof(float));
    if (new_location == NULL) {
        fprintf (stderr,
               "error in cut_locations: calloc failed for new location\n");
        return (0);
    }
    new_x_loc_ptr = new_location;
    new_y_loc_ptr = &(new_location[data_count]);

    /*
     * Transfer the data.
     */
    index=0;
    for (i=0; i<rows; i++) {
        for (j=0; j<cols; j++) {

            /*
             * First check to see if data is present for this pixel in
             * any one of the bands.
             */
            data_found=0;
            for (band=0; band<num_data_bands; band++) {
                if (!get_bit_pixel(imagedata, band, rows, cols, j, i, &data)) {
                    fprintf (stderr,
                             "error in cut_locations: get_bit_pixel failed\n");
                    return (0);
                }

                if (data != 0)
                    data_found=1;
            }

            /*
             * If data was found at this pixel then transfer the data
             * and the locations to the new image.
             */
            if (data_found) {
                for (band=0; band<num_data_bands; band++) {
                    if (!get_bit_pixel(imagedata,band,rows,cols,j,i,&data)) {
                        fprintf (stderr,
                              "error in cut_locations: get_bit_pixel failed\n");
                        return (0);
                    }

                    if (!put_bit_pixel(new_imagedata,band,1,location_count,
                    index,0,&data)){
                        fprintf (stderr,
                              "error in cut_locations: put_bit_pixel failed\n");
                        return (0);
                    }

                    new_x_loc_ptr[index] = x_loc_ptr[i*cols+j];
                    new_y_loc_ptr[index] = y_loc_ptr[i*cols+j];

                }
                index++;
            }
        }
    }

    free (img->imagedata);
    img->imagedata = (char *) new_imagedata;

    free (img->location);
    img->location = new_location;

    img->row_size=data_count;
    img->col_size=1;

    return (1);
}



/****************************************************************
*
* Routine: get_bit_pixel
*
* Purpose: This routine is used to address a pixel in a bit image.
*          The inputs are the data portion of an xvimage structure,
*          the band we want to access in the image, and the number
*          of rows and columns in the image.  Also input are the
*          x and y coordinates of the bit we want to retrieve
*          from the image.  The ouput is an integer, which will
*          be 0 or 1 depending on the value of the pixel we 
*          found.
*
*   Input: imagedata - a pointer to the data in the image
*          band - the band of the image we want to access
*          rows,cols - the number of rows and columns in the image
*          x,y - the coordinates of the pixel we want
*
*  Output: pixel the integer value of the bit at the given location
*
*  Author: Per Lysne  3/6/91
*
****************************************************************/
static int get_bit_pixel (imagedata, band, rows, cols, x, y, pixel)
unsigned char *imagedata;
int band, rows, cols;
int x, y;
int *pixel;
{
    unsigned char mask, data;
    int byte, offset;

    /*
     * Error check the inputs.
     */
    if (imagedata == NULL) {
        fprintf (stderr, "error in get_bit_pixel: bad imagedata input\n");
        return (0);
    }
    if (band < 0) {
        fprintf (stderr, "error in get_bit_pixel: bad band input: %d\n", band);
        return (0);
    }
    if (rows < 0) {
        fprintf (stderr, "error in get_bit_pixel: bad rows input: %d\n", rows);
        return (0);
    }
    if (cols < 0) {
        fprintf (stderr, "error in get_bit_pixel: bad cols input: %d\n", cols);
        return (0);
    }
    if ((x<0) || (x>cols)) {
        fprintf (stderr, "error in get_bit_pixel: bad x input: %d\n", x);
        return (0);
    }
    if ((y<0) || (y>rows)) {
        fprintf (stderr, "error in get_bit_pixel: bad y input: %d\n", y);
        return (0);
    }
    if (pixel == NULL) {
        fprintf (stderr, "error in get_bit_pixel: bad pixel input\n");
        return (0);
    }

    /*
     * Locate the byte which contains the pixel we want in the image.
     */
    byte = (int) (band*rows*cols/8) + (int) (y*cols/8) + (int) (x/8);
    offset = x%8;
    mask = 1<<offset;

    /*
     * Get the byte we want and AND it with the mask to isolate one pixel.
     */
    data = imagedata[byte];
    data = data & mask;

    /*
     * Set the return value depending on what was found in the pixel.
     */
    if (data == 0) {
        *pixel = 0;
    }
    else {
        *pixel = 1;
    }

    return (1);
}



/****************************************************************
*
* Routine: put_bit_pixel
*
* Purpose: This routine is used to address a pixel in a bit image.
*          The inputs are the data portion of an xvimage structure,
*          the band we want to access in the image, and the number
*          of rows and columns in the image.  Also input are the
*          x and y coordinates of the bit we want to place
*          from the image.  The pixel input is an integer that should
*          be 0 or 1 depending on the value of the pixel we 
*          want to set.
*
*   Input: imagedata - a pointer to the data in the image
*          band - the band of the image we want to access
*          rows,cols - the number of rows and columns in the image
*          x,y - the coordinates of the pixel we want
*          pixel -  the integer value to set the bit at the given location
*
*  Output: none
*
*  Author: Per Lysne  3/6/91
*
****************************************************************/
static int put_bit_pixel (imagedata, band, rows, cols, x, y, pixel)
unsigned char *imagedata;
int band, rows, cols;
int x, y;
int *pixel;
{
    unsigned char mask, data;
    int byte, offset;

    /*
     * Error check the inputs.
     */
    if (imagedata == NULL) {
        fprintf (stderr, "error in put_bit_pixel: bad imagedata input\n");
        return (0);
    }
    if (band < 0) {
        fprintf (stderr, "error in put_bit_pixel: bad band input: %d\n", band);
        return (0);
    }
    if (rows < 0) {
        fprintf (stderr, "error in put_bit_pixel: bad rows input: %d\n", rows);
        return (0);
    }
    if (cols < 0) {
        fprintf (stderr, "error in put_bit_pixel: bad cols input: %d\n", cols);
        return (0);
    }
    if ((x<0) || (x>cols)) {
        fprintf (stderr, "error in put_bit_pixel: bad x input: %d\n", x);
        return (0);
    }
    if ((y<0) || (y>rows)) {
        fprintf (stderr, "error in put_bit_pixel: bad y input: %d\n", y);
        return (0);
    }
    if (pixel == NULL) {
        fprintf (stderr, "error in put_bit_pixel: bad pixel input\n");
        return (0);
    }

    /*
     * Locate the byte which contains the pixel we want in the image.
     */
    byte = (int) (band*rows*cols/8) + (int) (y*cols/8) + (int) (x/8);
    offset = x%8;
    mask = 1<<offset;

    /* Get the byte containing the pixel. */
    data = imagedata[byte];

    /*
     * If we want to clear the bit, then invert the mask and AND it with
     * the data.  To set the bit, OR the byte with the mask.
     */
    if (*pixel == 0) {
        mask = mask ^ 255;
        data = data & mask;
    }
    else {
        data = data | mask;
    }

    /* Return the byte containing the pixel. */
    imagedata[byte] = data;

    return (1);
}
