coverage_vis.c

卫星仿真软件 卫星仿真软件 卫星仿真软件

void
update_display(int projection, grid * g, int coverage_flag, int tracks_flag,
	       int project_tracks_flag)
{
  Tk_PhotoHandle handle;

  /*
   * Re-write the entire image, unless we are just displaying ground
   * tracks, and there's not a new projection, in which case we can
   * simply write the foreground (which includes ground tracks).
   */
  int fgOnly = !new_foreground(projection, g, FALSE) && !coverage_flag &&
	       !tracks_flag && !project_tracks_flag;
  grid_and_foreground_to_image(projection, g, fgOnly);

  handle = Tk_FindPhoto(interp, imagename);
  if (handle) {
    /* double width of image to tile horizontally */
#if ( (TK_MAJOR_VERSION > 7) && (TK_MINOR_VERSION > 3) )
   /* changed for version 8.4, September 2002. */
    Tk_PhotoPutBlock(handle, get_image(), 0, 0, image->width*2, image->height,
        TK_PHOTO_COMPOSITE_OVERLAY);
#else
    Tk_PhotoPutBlock(handle, get_image(), 0, 0, image->width*2, image->height);
#endif
  } else {
    error("coverage_display: bad image name.");
  }

  /* display coverage percentages */
  set_coverage(projection, g);
}

/*
 * Clear the intensity array.
 */
void
clear_intensity(grid * g)
{
  unsigned int block;

  block = sizeof(int) * g->height * g->width;
  memset(g->data, 0, (size_t) block);
  memset(g->covered, 0, (size_t) sizeof(int) * g->height);
  g->count = 0;
}

/*
 * fill interval-decay with high value meaning many intervals since coverage
 */
void
clear_interval(grid * g)
{
  unsigned int block;

  block = sizeof(int) * g->height * g->width;
  memset(g->noaccess, num_noaccess_colors, (size_t) block);
}

/*
 * Increment everything on the interval grid towards decay.
 */
void
decay_interval(grid * g)
{
  unsigned int i;
  unsigned int image_size = g->height * g->width;

  for (i = 0; i < image_size; i++) {
  	(g->noaccess[i])++;
  }
}

/*
 * set_coverage - Display total percent coverage.
 */
static void
set_coverage(int projection, grid * g)
{
  char cmd[sizeof(format1) + 5];

  sprintf(cmd, format1, total_coverage(projection, TRUE, g));
  tcl_script(cmd);
  sprintf(cmd, format2, mean_coverage(projection, g));
  tcl_script(cmd);
}

/*
 * overlay_bitmap - Takes a bitmap and lies it on top of our image array.
 * h,w = height, width of the image array.
 */
static void
overlay_bitmap(unsigned int h, unsigned int w, unsigned char *b,
	       const char filename[])
{
  FILE *f;
  unsigned int pbmtype, image_size, width, height, i, j, c, mask;

  if (NULL == (f = fopen(filename, "r")))
    return;

  forward_over_comments(f);
  fscanf(f, "P%d", &pbmtype);
  forward_over_comments(f);
  fscanf(f, "%d", &width);
  forward_over_comments(f);
  fscanf(f, "%d", &height);
  /* skip one char */
  c = getc(f);

  if ((pbmtype == 4) && (width == w) && (height == h)) {
    image_size = h*w;
    for (i = 0; i < image_size; i += 8) {
      c = getc(f);
      mask = 0x80;
      for (j = 0; j < 8; j++) {
	if (c & mask) {
	  *(b++) = 1;
	} else {
	  b++;
	}
	mask = mask / 2;
      }
    }
  } else {
    fprintf(stderr, "Error in bitmap format: P%d %d %d\n",
		pbmtype, width, height);
  }

  fclose(f);
}

static void
draw_sin_map_edges(grid * g)
{
  SphericalCoordinates point;
  int left[2], right[2];
  unsigned int h = g->height;
  unsigned int w = g->width;
  unsigned int i;

  point.r = 1;
  point.theta = 0;
  point.phi = 0;

  for (i = 0; i < h; i++) {
    intensity_edges(left, right, &point, SINUSOIDAL, g);
    sin_foreground[(w * i + left[0])] = 1;
    sin_foreground[(w * i + right[0])] = 1;
    point.phi += PI / (h - 1);	/* Increment phi by one pixel */
  }
}

/*
 * new_foreground - Initialize foreground intensity image, including earth outline.
 * Returns TRUE if foreground must be re-drawn, FALSE otherwise.
 */
static int
new_foreground(int projection, grid * g, unsigned int force)
{
  unsigned int image_size;

  if (!force && (projection == current_proj))
    return (FALSE);

  if (current_proj != projection) {
    current_proj = projection;
    clear_intensity(g);
    clear_interval(g);
  }

  image_size = g->height * g->width;

  switch (projection) {
  case UNPROJECTED:
    Longitude_Center_Line = LONGITUDE_CENTER_LINE;
    memcpy(foreground, unp_foreground, image_size);
    break;
  case SINUSOIDAL:
    Longitude_Center_Line = LONGITUDE_CENTER_LINE_SINU;
    memcpy(foreground, sin_foreground, image_size);
    break;
  case SPHERICAL:
    Longitude_Center_Line = LONGITUDE_CENTER_LINE;
    memcpy(foreground, sph_foreground, image_size);
    break;
  case CYLINDRICAL:
  default:
    Longitude_Center_Line = LONGITUDE_CENTER_LINE;
    memcpy(foreground, cyl_foreground, image_size);
  }

  return (TRUE);
}

/*
 * Copy foreground (land masses and ground tracks, if any) only onto image.
 * The routine re-draws an image pixel only for corresponding foreground array
 * elements which are non-zero.
 */
static void
grid_and_foreground_to_image(int projection, grid * g, int fgOnly)
{
  unsigned int i, j, offset, ci;
  unsigned int w = g->width;
  unsigned int h = g->height;
  unsigned int num_colors_dec, num_noaccess_colors_dec;
  unsigned int data_value, map_threshold, send_coverage;
  int drawWidth = w;
  int middleWidth = w/2;
  unsigned int sinu_surround = FALSE;
  unsigned int sph_surround = FALSE;
  unsigned char *c;
  unsigned char *f = foreground;
  unsigned char fg, cross_fg;
  unsigned char *p = image->pixelPtr;
  unsigned char *q, *r;
  float lat, width, fw;

  /*
   * Do we need to draw a blank surround outside
   * the map?
   */
  if (!fgOnly) {
    if (projection == SINUSOIDAL) {
      sinu_surround = TRUE;
    }
    if (projection == SPHERICAL) {
      sph_surround = TRUE;
    }
  }

  fw = (float) w;
  num_colors_dec = num_colors - 1;
  num_noaccess_colors_dec = num_noaccess_colors - 1;

  send_coverage = geomview_flag && geomview_dynamic_texture_flag &&
                  texture_flag && (projection == CYLINDRICAL) &&
                  earth_on_flag && use_fancy_earth && earth_geom_exists;

   /*
    * If threshold is 0, we don't draw white background (0), but we
    * do draw map outline (1) above that. If 1, we will not draw
    * map outline. For 0/1, could have used || instead of MAX,
    * but MAX is clearer in intent.
    */
  map_threshold = MAX((send_coverage && !geomview_texture_with_map), !map_flag);

  for (i = 0; i < h; i++) {
    lat = 90.0 - 180.0 * i / h;
    if (sinu_surround) {
      width = (fw * cos(DEG_TO_RAD * lat));
      drawWidth = (int) (0.5 + (fw - width) / 2);
    }

    if (sph_surround) {
      width = fw/4 - sqrt(fw*fw/16-fw*fw*lat*lat/16/90/90);
      drawWidth = (int) width;
      middleWidth = w/2 - drawWidth;
    }

    offset = i * w;
    q = p;

    for (j = 0; j < w; j++, p += PIXELSIZE) {
      fg = *(f++);

      if (sinu_surround || sph_surround) {
	q = p; /* updating for cross wrap-around */
	if ((j <= drawWidth) || j >= (w - drawWidth)) {
	  p[0] = outline_colors[0];
	  p[1] = outline_colors[1];
	  p[2] = outline_colors[2];
	  continue;
	}
	if (sph_surround && (j > middleWidth) && (j < w - middleWidth)) {
	  p[0] = outline_colors[0];
	  p[1] = outline_colors[1];
	  p[2] = outline_colors[2];
	  continue;
	}
      }

      if (fg > map_threshold) {
	/* skipping Earth outline if map is to be sent */
	c = outline_colors + PIXELSIZE * fg;

	p[0] = c[0];
	p[1] = c[1];
	p[2] = c[2];

        /*
	 * Indicate sub-satellite point by drawing four arms of cross.
	 */
	if ((fg == GROUND_TRACKS_CI) || (fg == SPECIAL_GROUND_TRACKS_CI)) {
	  cross_fg = fg + 1;

	  if (j > 0) {
	    /* draw left arm of cross */
	    *(f-2) = cross_fg; /* set retroactively for future updates. */
	    r = p - PIXELSIZE;
	    r[0] = c[0];
	    r[1] = c[1];
	    r[2] = c[2];
	  } else {
	    /* left arm wrap around - set pixel at j=w */
	    /* not well-matched to sinusoidal */
	    *(f+w-2) = cross_fg; /* ready to draw at end of this line */
	  }
	  if (i > 0) {
	    /* draw top arm of cross */
	    *(f-w-1) = cross_fg; /* set retroactively for future updates. */
	    r = p - w * PIXELSIZE;
	    r[0] = c[0];
	    r[1] = c[1];
	    r[2] = c[2];
	  }
	  if (j < w-1) {
	    /* ready to draw right arm of cross, next j iteration */
	    *f = cross_fg;
	  } else {
	    /* right arm wrap around - set pixel at j=0 i same value */
	    /* not well-matched to sinusoidal */
	    *(f-j) = cross_fg;
	    q[0] = c[0];
	    q[1] = c[1];
	    q[2] = c[2];
	  }
	  if (i < h-1) {
	    /* ready to draw bottom of cross, next i iteration/row pass */
	    *(f-1 + w) = cross_fg;
	  }
	}

	continue;
      }

      if (!fgOnly) {
	data_value = g->data[offset + j];
	if (data_value || !no_access_flag) {
	  ci = MIN(data_value, num_colors_dec);
	  c = colors + PIXELSIZE * ci;
	} else {
	  ci = MIN(g->noaccess[offset + j], num_noaccess_colors_dec);
	  c = noaccess_colors + PIXELSIZE * ci;
	}
	p[0] = c[0];
	p[1] = c[1];
	p[2] = c[2];
      }
    }
  }

  if (send_coverage) {
#ifndef NO_ZLIB
    write_image_to_file_gz(temp_scratchfile);
#else
    write_image_to_file(temp_scratchfile);
#endif /* NO_ZLIB */
    coverage_send_scratchfile();
  }

  if (!map_threshold || !map_flag) {
    /* map outline is already drawn; no need to fill in for coverage window. */
    return;
  }

  f = foreground;
  p = image->pixelPtr;
  drawWidth = h*w;
  c = outline_colors + PIXELSIZE;
  for (i = 0; i <drawWidth; i++, p+=PIXELSIZE) {
    fg = *(f++);
    if (fg == 1) {
      p[0] = c[0];
      p[1] = c[1];
      p[2] = c[2];
    }
  }
}

void
coverage_send_static_file(void)
{
  fprintf(gv_out, "(read geometry { define earth_h {\n %s \n} } )", static_description);
}

void
coverage_send_scratchfile(void)
{
  fprintf(gv_out, "(read geometry { define earth_h {\n %s \n} } )", dynamic_description);
}

/*
 * tracks_to_foreground - place dot on the foreground map which corresponds
 * to the subsatellite point (nadir).  Do this for every satellite
 * in the list.
 */
void
tracks_to_foreground(const Satellite_list SL, int projection, grid * g,
		     unsigned int color_index, const CentralBody * pcb)
{
  unsigned int color;
  Satellite_list sl = SL;

  while (sl) {
    if (sl->s->tag == 1) {
      /* inelegant choice of special color */
      color = SPECIAL_GROUND_TRACKS_CI;
    } else {
      color = color_index;
    }
    track_to_foreground(sl->s, projection, g, color, pcb);
    sl = sl->next;
  }
}

/*
 * track_to_foreground - routine to increment the foreground map at the point
 * underneath the satellite.  For a single satellite.
 *
 * grid_index = color index to insert into the foreground array
 */
void
track_to_foreground(Satellite s, int projection, grid * g,
		    unsigned int color_index, const CentralBody * pcb)
{
  LatLon l;
  int coords[2];

  spherical_to_lat_lon(&l, &(s->x_S), s->t, pcb);
  latlon_to_grid_index(coords, &l, projection, g);
  foreground[g->width * coords[1] + coords[0]] = color_index;

}

/*
 * Given lat/lon coordinates of point on central body, return corresponding
 * index into grid.
 *        grid_index[0]=column index, grid_index[1]=row index.
 */
static void
latlon_to_grid_index(int grid_index[2], LatLon * pl, int projection, grid * g)
{
  double proj[2];
  LatLon l;

  /*
   * The vertical center-line of the grid corresponds to longitude =
   * Longitude_Center_Line, whereas both 2-d projections place longitude = 0
   * along this line.  So shift the input longitude as required.
   */
  l.lat = pl->lat;
  l.lon = pl->lon - Longitude_Center_Line;

  if (l.lon < -180.0) {
    l.lon += 360.0;
  } else if (l.lon > 180.0) {
    l.lon -= 360.0;
  }

  switch (projection) {
  case UNPROJECTED:
    project_latlon_unprojected(proj, &l);
    break;
  case SINUSOIDAL:
    project_latlon_sinusoidal(proj, &l);
    break;
  case SPHERICAL:
    project_latlon_spherical(proj, &l);
    break;
  case CYLINDRICAL:
  default:
    project_latlon_cylindrical(proj, &l);
  }

  /*
   * It's assumed that proj[0] ranges from -180 to 180, and
   *                   proj[1] ranges from -90 to 90,  no matter what.
   */
  grid_index[0] = ((int) (0.5 + (proj[0] + 180.0) * (g->width) / 360.0));
  grid_index[1] = ((int) (0.5 + proj[1] * (g->height - 1) / 180.0));
}