sattools/src/jpgstack.c

#include <stdio.h>
#include <stdlib.h>
#include <jpeglib.h>
#include <math.h>
#include <cpgplot.h>

#define NMAX 1024

struct image
{
  int nx, ny, nz;
  float *z;
};
struct image read_jpg (char *filename);
void write_jpg (char *filename, struct image img);


int
main (int argc, char *argv[])
{
  int i, j, flag, n;
  struct image avg, max, raw;

  for (flag = 0, i = 1, n = 0; i < argc; i++, n++)
    {
      printf ("%d %s\n", i, argv[i]);
      // Read image
      raw = read_jpg (argv[i]);

      // If first image, initialize
      if (flag == 0)
	{
	  avg.nx = raw.nx;
	  avg.ny = raw.ny;
	  avg.nz = raw.nz;
	  avg.z =
	    (float *) malloc (sizeof (float) * avg.nx * avg.ny * avg.nz);
	  for (j = 0; j < avg.nx * avg.ny * avg.nz; j++)
	    avg.z[j] = 0.0;

	  max.nx = raw.nx;
	  max.ny = raw.ny;
	  max.nz = raw.nz;
	  max.z =
	    (float *) malloc (sizeof (float) * max.nx * max.ny * max.nz);
	  for (j = 0; j < max.nx * max.ny * max.nz; j++)
	    max.z[j] = 0.0;

	  flag = 1;
	}

      // Add values
      for (j = 0; j < avg.nx * avg.ny * avg.nz; j++)
	avg.z[j] += raw.z[j];

      // Add values
      for (j = 0; j < avg.nx * avg.ny * avg.nz; j++)
	if (raw.z[j] > max.z[j])
	  max.z[j] = raw.z[j];

      // Free
      free (raw.z);
    }

  // Average;
  for (j = 0; j < avg.nx * avg.ny * avg.nz; j++)
    avg.z[j] /= (float) n;

  // Write
  write_jpg ("average.jpg", avg);
  write_jpg ("maximum.jpg", max);

  // Free
  free (avg.z);
  free (max.z);

  return 0;
}

struct image
read_jpg (char *filename)
{
  int i = 0, j, k, l, m;
  unsigned long location = 0;
  struct image img;
  struct jpeg_decompress_struct cinfo;
  struct jpeg_error_mgr jerr;
  JSAMPROW row_pointer[1];
  unsigned char *raw_image = NULL;
  FILE *file;

  // Open file
  file = fopen (filename, "rb");
  if (!file)
    perror ("Error opening file");

  // Get header info, decompress
  cinfo.err = jpeg_std_error (&jerr);
  jpeg_create_decompress (&cinfo);
  jpeg_stdio_src (&cinfo, file);
  jpeg_read_header (&cinfo, TRUE);
  jpeg_start_decompress (&cinfo);

  // Allocate memory
  raw_image =
    (unsigned char *) malloc (cinfo.output_width * cinfo.output_height *
			      cinfo.num_components);

  // Read image, one scan at a time
  row_pointer[0] =
    (unsigned char *) malloc (cinfo.output_width * cinfo.num_components);
  while (cinfo.output_scanline < cinfo.image_height)
    {
      jpeg_read_scanlines (&cinfo, row_pointer, 1);
      for (i = 0; i < cinfo.image_width * cinfo.num_components; i++)
	raw_image[location++] = row_pointer[0][i];
    }
  // wrap up decompression, destroy objects, free pointers and close open files
  jpeg_finish_decompress (&cinfo);
  jpeg_destroy_decompress (&cinfo);

  // Copy image to image struct
  img.nx = cinfo.image_width;
  img.ny = cinfo.image_height;
  img.nz = cinfo.num_components;
  img.z = (float *) malloc (sizeof (float) * img.nx * img.ny * img.nz);

  // Fill image
  for (i = 0; i < img.nx; i++)
    {
      for (j = 0; j < img.ny; j++)
	{
	  for (k = 0; k < img.nz; k++)
	    {
	      l = img.nz * (i + img.nx * j) + k;
	      img.z[l] = (float) raw_image[l];
	    }
	}
    }

  // Free allocated memory
  free (row_pointer[0]);
  free (raw_image);

  // Close file
  fclose (file);

  return img;
}

// Write jpg
void
write_jpg (char *filename, struct image img)
{
  int i, j, k, l, m;
  struct jpeg_compress_struct cinfo;
  struct jpeg_error_mgr jerr;
  JSAMPROW row_pointer[1];
  FILE *outfile;
  unsigned char *raw_image = NULL;

  outfile = fopen (filename, "wb");
  cinfo.err = jpeg_std_error (&jerr);
  jpeg_create_compress (&cinfo);
  jpeg_stdio_dest (&cinfo, outfile);
  cinfo.image_width = img.nx;
  cinfo.image_height = img.ny;
  cinfo.input_components = 3;
  cinfo.in_color_space = JCS_RGB;
  jpeg_set_defaults (&cinfo);
  jpeg_start_compress (&cinfo, TRUE);

  // Allocate memory
  raw_image =
    (unsigned char *) malloc (cinfo.image_width * cinfo.image_height *
			      cinfo.input_components);

  // Fill image
  for (i = 0; i < img.nx; i++)
    {
      for (j = 0; j < img.ny; j++)
	{
	  for (k = 0; k < img.nz; k++)
	    {
	      l = img.nz * (i + img.nx * j) + k;
	      raw_image[l] = (unsigned char) img.z[l];
	    }
	}
    }

  while (cinfo.next_scanline < cinfo.image_height)
    {
      row_pointer[0] =
	&raw_image[cinfo.next_scanline * cinfo.image_width *
		   cinfo.input_components];
      jpeg_write_scanlines (&cinfo, row_pointer, 1);
    }
  jpeg_finish_compress (&cinfo);
  jpeg_destroy_compress (&cinfo);
  fclose (outfile);

  return;
}