/*
################################################################
# PDPVISUAL.C
#
# This file uses the IJG code to read and write JPEG image files.
# The JPEG code is based on an example given in the IJG 
# distribution.
#
# PDPVISUAL reads a current full disk JPEG file 'current.jpg'
#(could also be user defined and defined as an argument to main),
# which has to be P-angle corrected. It then calculates
# the current ephemeris, center the image, remap the image to
# a standard solar size (360 pixel diameter), and adds a 
# Stonyhurst grid. 
#
#    Created: 1997-11-27 by Anders Johannesson
#
################################################################ 
*/

#include <stdio.h>
#include "jpeglib.h"
#include <setjmp.h>
#include "soleph.h"

unsigned char *image;	/* Points to large array of grayscale data */
unsigned char *outim;   /* Points to new array of grayscale data */
unsigned char timetext[100][26]; /* GRAPHICAL TIME */
unsigned char datetext[96][26]; /* GRAPHICAL DATE */
int image_height=0;	/* Number of rows in image */
int image_width=0;	/* Number of columns in image */
long outim_height=385;	/* Number of rows in outim */
long outim_width=384;	/* Number of columns in outim */
int disk_diameter=360;  /* Number of pixels in remapped diameter */
double dc_level=170.0;  /* Images normalized to this disk center level */

/*
###########################################################
                    WRITE A JPEG FILE 
########################################################### 
*/
GLOBAL(void)
write_JPEG_file (char * filename, int quality)
{

  struct jpeg_compress_struct cinfo;
  struct jpeg_error_mgr jerr;
  FILE * outfile;		/* target file */
  JSAMPROW row_pointer[1];	/* pointer to JSAMPLE row[s] */
  int row_stride;		/* physical row width in image buffer */

  /* STEP 1: ALLOCATE AND INITIALIZE JPEG COMPRESSION OBJECT */
  cinfo.err = jpeg_std_error(&jerr);
  jpeg_create_compress(&cinfo);

  /* STEP 2: SPECIFY DATA DESTINATION (EG, A FILE) */
  if ((outfile = fopen(filename, "wb")) == NULL) {
    fprintf(stderr, "can't open %s\n", filename);
    exit(1);
  }
  jpeg_stdio_dest(&cinfo, outfile);

  /* STEP 3: SET PARAMETERS FOR COMPRESSION */
  cinfo.image_width = outim_width; 	/* image width and height, in pixels */
  cinfo.image_height = outim_height;
  cinfo.input_components = 1;		/* # of color components per pixel */
  cinfo.in_color_space = JCS_GRAYSCALE; /* colorspace of input image */
  jpeg_set_defaults(&cinfo);
  jpeg_set_quality(&cinfo, quality, TRUE /* limit to baseline-JPEG values */);

  /* STEP 4: START COMPRESSOR */
  jpeg_start_compress(&cinfo, TRUE);

  /* STEP 5: WHILE (SCAN LINES REMAIN TO BE WRITTEN) */
  row_stride = outim_width*cinfo.input_components;  /* JSAMPLEs per row in image_buffer */
  while (cinfo.next_scanline < cinfo.image_height) {
    row_pointer[0] = &outim[cinfo.next_scanline * row_stride]; 
    (void) jpeg_write_scanlines(&cinfo, row_pointer, 1);
  }

  /* STEP 6: FINISH COMPRESSION */
  jpeg_finish_compress(&cinfo);
  fclose(outfile);

  /* STEP 7: RELEASE JPEG COMPRESSION OBJECT */
  jpeg_destroy_compress(&cinfo);

}

/*
 * ERROR HANDLING:
 *
 * The JPEG library's standard error handler (jerror.c) is divided into
 * several "methods" which you can override individually.  This lets you
 * adjust the behavior without duplicating a lot of code, which you might
 * have to update with each future release.
 *
 * Our example here shows how to override the "error_exit" method so that
 * control is returned to the library's caller when a fatal error occurs,
 * rather than calling exit() as the standard error_exit method does.
 *
 * We use C's setjmp/longjmp facility to return control.  This means that the
 * routine which calls the JPEG library must first execute a setjmp() call to
 * establish the return point.  We want the replacement error_exit to do a
 * longjmp().  But we need to make the setjmp buffer accessible to the
 * error_exit routine.  To do this, we make a private extension of the
 * standard JPEG error handler object.  (If we were using C++, we'd say we
 * were making a subclass of the regular error handler.)
 *
 * Here's the extended error handler struct:
 */

struct my_error_mgr {
  struct jpeg_error_mgr pub;	/* "public" fields */

  jmp_buf setjmp_buffer;	/* for return to caller */
};

typedef struct my_error_mgr * my_error_ptr;

/*
 * Here's the routine that will replace the standard error_exit method:
 */

METHODDEF(void)
my_error_exit (j_common_ptr cinfo)
{
  /* cinfo->err really points to a my_error_mgr struct, so coerce pointer */
  my_error_ptr myerr = (my_error_ptr) cinfo->err;

  /* Always display the message. */
  /* We could postpone this until after returning, if we chose. */
  (*cinfo->err->output_message) (cinfo);

  /* Return control to the setjmp point */
  longjmp(myerr->setjmp_buffer, 1);
}


/*
###########################################################
 READ A JPEG FILE 
########################################################### 
*/
GLOBAL(int)
read_JPEG_file (char * filename)
{

  struct jpeg_decompress_struct cinfo;
  struct my_error_mgr jerr;
  FILE * infile;		/* source file */
  JSAMPARRAY buffer;		/* Output row buffer */
  int row_stride, i;		/* physical row width in output buffer */
  long pixel_number;

  /* OPEN FILE */
  if ((infile = fopen(filename, "rb")) == NULL) {
    fprintf(stderr, "can't open %s\n", filename);
    return(-1);
  }

  /* STEP 1: ALLOCATE AND INITIALIZE JPEG DECOMPRESSION OBJECT */
  cinfo.err = jpeg_std_error(&jerr.pub);
  jerr.pub.error_exit = my_error_exit;
  if (setjmp(jerr.setjmp_buffer)) {
    jpeg_destroy_decompress(&cinfo);
    fclose(infile);
    return 0;
  }
  jpeg_create_decompress(&cinfo);

  /* STEP 2: SPECIFY DATA SOURCE (EG, A FILE) */
  jpeg_stdio_src(&cinfo, infile);

  /* STEP 3: READ FILE PARAMETERS WITH JPEG_READ_HEADER() */
  (void) jpeg_read_header(&cinfo, TRUE);
  
  /* STEP 4: SET PARAMETERS FOR DECOMPRESSION */

  /* In this example, we don't need to change any of the defaults set by
   * jpeg_read_header(), so we do nothing here.
   */

  /* STEP 5: START DECOMPRESSOR */
  (void) jpeg_start_decompress(&cinfo);
  row_stride = cinfo.output_width * cinfo.output_components;
  buffer = (*cinfo.mem->alloc_sarray)
		((j_common_ptr) &cinfo, JPOOL_IMAGE, row_stride, 1);
		
  /* ALLOCATE IMAGE ARRAY */
  pixel_number=(long)cinfo.output_width * (long)cinfo.output_height;
  image_width=cinfo.output_width;
  image_height=cinfo.output_height; 
  image = (unsigned char *) malloc(pixel_number);

  /* STEP 6: WHILE (SCAN LINES REMAIN TO BE READ) */
  while (cinfo.output_scanline < cinfo.output_height) {
    (void) jpeg_read_scanlines(&cinfo, buffer, 1);
    for (i=0; i<cinfo.output_width; i++)
    {
        image[(cinfo.output_scanline-1)*row_stride+i] = (char)*(buffer[0]+i); 
    }
  }

  /* STEP 7: FINISH DECOMPRESSION */
  (void) jpeg_finish_decompress(&cinfo);
  
  /* STEP 8: RELEASE JPEG DECOMPRESSION OBJECT */
  jpeg_destroy_decompress(&cinfo);
  fclose(infile);

  return(0);
}


/*
###########################################################
 GET TIME TEXT (ASSUMES KODAK 4.2 IMAGES)
########################################################### 
*/

void get_time_text()
{
    
   short row, col;
   
   /* EXTRACT */
   for (row=0; row<26; row++)
   {
      for (col=180; col<280; col++)
      {
         timetext[col-180][row] = image[row*image_width+col];
      }
   }
   
   /* EXTRACT */
   for (row=0; row<26; row++)
   {
      for (col=85; col<180; col++)
      {
         datetext[col-85][row] = image[row*image_width+col];
      }
   }
   
}

/*
###########################################################
 PUT TIME TEXT 
########################################################### 
*/

void put_time_text()
{

   short row, col;
   
   /* INSERT */
   for (row=0; row<26; row++)
   {
      for (col=0; col<100; col++)
      {
          outim[row*outim_width+col] = timetext[col][row];
      }
   }
   
   /* INSERT */
   for (row=0; row<26; row++)
   {
      for (col=0; col<96; col++)
      {
          outim[row*outim_width+col+(outim_width-96-2)] = datetext[col][row];
      }
   }
   
}

/*
###########################################################
 NORMALIZE IMAGE 
########################################################### 
*/

int normalize()
{
   
   double maxlevel=0.0, darklevel=0.0, pixel;
   long counter, row, col;
      
   /* FIND DISK CENTER BRIGHTNESS */
   counter = 0;
   for (row=(image_height/2-image_height/10); row<(image_height/2+image_height/10); row++)
   {
      for (col=(image_width/2-image_width/10); col<(image_width/2+image_width/10); col++)
      {
         maxlevel+= (double) image[row*image_width+col];
         counter++;       
      }
   }
   maxlevel = maxlevel / (double) counter;
   
   /* FIND DARK BACKGROUND */
   counter=0;
   for (row=0; row<image_height; row++)
   {
      for (col=image_width-8; col<image_width-5; col++)
      {
         darklevel+= (double) image[row*image_width+col];
         counter++;
      }
   }
   darklevel = darklevel / (double) counter;
   
   /* NORMALIZE */
   for (row=0; row<image_height; row++)
   {
      for (col=0; col<image_width; col++)
      {
         pixel = (double) image[row*image_width+col];
         if ((maxlevel-darklevel) > 0.0)
            pixel = (pixel - darklevel)/(maxlevel-darklevel)*dc_level;
         else
            pixel=0.0;
         if (pixel < 0.0)
            pixel=0.0;
         if (pixel > 255.0)
            pixel=255.0;
         image[row*image_width+col] = (unsigned char) (pixel+0.5);    
         counter++;       
      }
   }
   
   return(0); 
}

/*
###########################################################
 FIND CENTER OF IMAGE 
########################################################### 
*/

int center_im(long *x1,long *x2,long *y1,long* y2)
{
    long detector,i,j,number_of_done,offset=1,check_sum=0; /* 5 */
    int show_graphics=0;
    long threshold;
    
    /* CALCULATE THRESHOLD */
    threshold=(long) (0.50*dc_level);
    
    /*# FIND CENTERING (HORIZONTAL) */
    number_of_done=0;
    for (i=0; i<image_width; i++)
    {
        detector=0;
        for (j=image_height/5; j<(image_height-image_height/5); j++)
        {
            if( image[j*image_width+i] > threshold )
                detector=detector+1;
        }
        if ((detector > 10) && (number_of_done == 0))
        {
           *x1=i-offset;
           for (j=0; j<image_height; j++)
           {
               if (show_graphics == 1)
               {
                  image[j*image_width+i-offset] = 255;
                  image[j*image_width+i-offset-1] = 255;
               }
           }
           number_of_done=1;
        } 
        if ((detector < 11) && (number_of_done == 1)) 
        {
           *x2=i+offset;
           for (j=0; j<image_height; j++)
           {
               if (show_graphics == 1)
               {
                  image[j*image_width+i+offset] = 255;
                  image[j*image_width+i+offset+1] = 255;
               }
           }
           number_of_done=2;
        }
    }
    check_sum=number_of_done;
     
    /*# FIND CENTERING (VERTICAL) */
    number_of_done=0;
    for (j=0; j<image_height; j++)
    {
        detector=0;
        for (i=image_width/5; i<(image_width-image_width/5); i++)
        {
            if( image[j*image_width+i] > threshold )
                detector=detector+1;
        }
        if ((detector > 10) && (number_of_done == 0))
        {
           *y1=j-offset;
           for (i=0; i<image_width; i++)
           {
               if (show_graphics == 1)
               {
                  image[(j-offset)*image_width+i] = 255;
                  image[(j-offset-1)*image_width+i] = 255;
               }
           }
           number_of_done=1;
        } 
        if ((detector < 11) && (number_of_done == 1)) 
        {
           *y2=j+offset;
           for (i=0; i<image_width; i++)
           {
               if (show_graphics == 1)
               {
                  image[(j+offset)*image_width+i] = 255;
                  image[(j+offset+1)*image_width+i] = 255;
               }
           }
           number_of_done=2;
        }
     }       
     
     check_sum=check_sum+number_of_done;
     
     if ( (check_sum == 4)  && ((*x2-*x1) > 100) && ((*y2-*y1) > 100) )
         return(0);
     else
         return(-1);

}

/*
###########################################################
 REMAP IMAGE 
########################################################### 
*/

int remap_im(long x1,long x2,long y1,long y2)
{

   unsigned char new_pixel;
   int row,col,startx,starty, old_row, old_col;
   long pixel_number;
   float scalex, scaley;
    
   /*# ALLOCATE SPACE FOR OUTPUT IMAGE */
   pixel_number= (long)outim_height * (long)outim_width;
   outim = (unsigned char *) malloc(pixel_number);
   
   /*# FILL IN PIXELS IN OUTPUT IMAGE */
   startx=(outim_width-disk_diameter)/2;
   starty=(outim_height-disk_diameter)/2;
   scaley=(float)(y2-y1+1)/(float)disk_diameter;
   scalex=(float)(x2-x1+1)/(float)disk_diameter;
   for (row=0; row<outim_height; row++)
   {
      for (col=0; col<outim_width; col++)
      {
         old_row=(int)(scaley*row+y1-starty*scaley);
         old_col=(int)(scalex*col+x1-startx*scalex);
         new_pixel=image[old_row*image_width+old_col];
         outim[row*outim_width+col]=new_pixel;       
      }
   } 
   
   return(0);  
   
}

/*
###########################################################
 SHARPEN IMAGE 
########################################################### 
*/

int sharpen_im(short ksize)
{

   unsigned char *smooth;
   double new_pixel;
   long pixel_number;
   short xkernel, ykernel, newrow, newcol, row, col;
   double blacklevel=100.0;
    
   /*# ALLOCATE SPACE FOR OUTPUT IMAGE */
   pixel_number= (long)outim_height * (long)outim_width;
   smooth = (unsigned char *) malloc(pixel_number);
    
   /*# MAKE A BOXCAR SMOOTHED VERSION OF THE IMAGE */
   for (row=0; row<outim_height; row++)
   {
      for (col=0; col<outim_width; col++)
      {
         new_pixel=0.0;
         for (xkernel=-ksize; xkernel<ksize; xkernel++)
         {
            for (ykernel=-ksize; ykernel<ksize; ykernel++)
            {
               newrow=row+ykernel;
               if (newrow < 0) newrow=0;
               if (newrow >= outim_height) newrow=outim_height-1;
               newcol=col+xkernel;
               if (newcol < 0) newcol=0;
               if (newcol >= outim_width) newcol=outim_width-1;
               new_pixel+= (double) outim[newrow*outim_width+newcol];
            }
         }
         smooth[row*outim_width+col]=(unsigned char)(new_pixel/((2*ksize+1)*(2*ksize+1)));              
      }
   }
   
   /*# SHARPEN IMAGE (DIFF MASK) */
   for (row=0; row<outim_height; row++)
   {
      for (col=0; col<outim_width; col++)
      {
          new_pixel=3.0*outim[row*outim_width+col]-2.0*smooth[row*outim_width+col];        
          new_pixel=(new_pixel-blacklevel)/(255.0-blacklevel)*255.0;
          if (new_pixel < 0) new_pixel=0;
          if (new_pixel >255) new_pixel=255;
          outim[row*outim_width+col] = (unsigned char) new_pixel;
      }
   }
   
   free(smooth);
   
   return(0);
   
}

/*
###########################################################
 MAKE EMPTY IMAGE
########################################################### 
*/

void make_empty()
{

   long pixel_number;
   short row, col;
    
   /*# ALLOCATE SPACE FOR OUTPUT IMAGE */
   pixel_number= (long)outim_height * (long)outim_width;
   outim = (unsigned char *) malloc(pixel_number);
   
   for (row=0; row<outim_height; row++)
   {
      for (col=0; col<outim_width; col++)
      {
         outim[row*outim_width+col]=128;       
      }
   } 
   
}
   
   
/*
###########################################################
 OVERLAY STONYHURST GRID
########################################################### 
*/

int overlay_stonyhurst(double B0)
/* P-ANGLE ASSUMED TO BE CORRECTED */

{

   short  row, col, error_code;
   double S,P=0.0,L0=0.0,x,y,L,B,BSTEP,LSTEP;
   long cenx, ceny;
   
   /*# GRID STEPS IN B AND L */
   BSTEP=10.0;
   LSTEP=10.0;
   
   /*# DISK OFFSET IN OUTPUT IMAGE */
   cenx=outim_width/2;
   ceny=outim_height/2;
   
   /*# SEMI DIAMETER OF THE SUN IN PIXELS */
   S= ((double) disk_diameter)/2.0;

   /*# FOR STEPS IN LATITUDE DRAW GRID */
   for (B=-90.0; B<=90.0; B+=BSTEP) 
   {
      for (L=-180.0; L<=180.0; L+=1.0)
      {
         error_code=helio_to_cart(S,B0,P,L0,B,L,&x,&y); 
         if (error_code != FARSUN)
         {
            row=(short)(y+0.5)+ceny;
            row=outim_height-row-1; /* COORDINATE SYSTEM REVERSED */
            col=(short)(x+0.5)+cenx; 
            if ((B > -0.5) && (B < 1.0)) 
               outim[row*outim_width+col]=0; 
            else
               outim[row*outim_width+col]=255;      
         }
      }
   }
   
   /*# FOR STEPS IN LONGITUDE DRAW GRID */
   for (B=-90.0; B<=90.0; B+=1.0) 
   {
      for (L=-180.0; L<=180.0; L+=LSTEP)
      {
         error_code=helio_to_cart(S,B0,P,L0,B,L,&x,&y); 
         if (error_code != FARSUN)
         {
            row=(short)(y+0.5)+ceny;
            row=outim_height-row-1; /* COORDINATE SYSTEM REVERSED */
            col=(short)(x+0.5)+cenx; 
            if ((L > -0.5) && (L < 0.5)) 
               outim[row*outim_width+col]=0; 
            else
               outim[row*outim_width+col]=255;      
         }
      }
   }
      
   return(0); 
}


/*
###########################################################
 MAIN
########################################################### 
 Usage: pdpvisual image.jpg  (default is: current.jpg)
###########################################################
*/
void main(int argc, char *argv[])
{
 
 /* TAG INPUT PARAMETER. OM SAKNAS ANVAEND CURRENT  OM SAKNAS ANVAEND SVART */
 
   long x1, x2, y1, y2;
   int error;
   struct tm *Current_Time_UT;
   time_t Current_Time;
   TimeType Etime;
   EphType *Ephem;
   double B0;
   char *filename;
   
   /* GET IMAGE NAME */
   if (argc > 1) 
      sprintf(filename,"%s",argv[1]);
   else
      filename="current.jpg";
      
   /* GET CURRENT TIME */
   time(&Current_Time);
   Current_Time_UT=gmtime(&Current_Time);
   Etime.iy=Current_Time_UT->tm_year+1900;
   Etime.im=Current_Time_UT->tm_mon+1;
   Etime.id=Current_Time_UT->tm_mday;
   Etime.ih=Current_Time_UT->tm_hour;
   Etime.imi=Current_Time_UT->tm_min;
   Etime.is=Current_Time_UT->tm_sec; 
   
   /* CALCULATE EPHEMERID (USE ONLY B0-ANGLE) */
   soleph(Etime,Ephem);
   B0=Ephem->b;

   /* READ INPUT JPEG FILE */
   error=read_JPEG_file (filename);
   
   if (error == -1)
   {
      make_empty();
   }
   else
   {
   
      /* GET TIME TEXT */
      get_time_text();
   
      /* NORMALIZE IMAGE */
      error=normalize();
   
      /* CENTER DISK */
      error=center_im(&x1,&x2,&y1,&y2);
      
      /* CENTERING OK? */
      if (error == -1)
      {
         make_empty();
      }
      else
      {
   
         /* REMAP TO STANDARD SIZE */
         error=remap_im(x1,x2,y1,y2);
   
         /* SHARPEN THE IMAGE */
         error=sharpen_im(7);
         
      }
      
   }
   
   /* OVERLAY STONYHURST GRID */
   error=overlay_stonyhurst(B0);
   
   /* PUT TIME TEXT */
   put_time_text();
   
   /* WRITE OUTPUT JPEG */
   write_JPEG_file ("newsun.jpg", 75);
   
}