//---------------------------------------------------------------------------

#include "global.h"
#include "stdio.h"
#include "imatrix.h"
#include "fits/fitsio.h"
#include "swephexp.h"

/***************************************************************************
                   ************* LoadImage **************
   this function is used only in the jpg wolf
   the function load the data from a BITMAP data structure to the matrix
   data structure. In fits wolf, the function "LoadFitsFile" is used.
*/
int ImageMatrix::LoadImage(BITMAP bitmap,short inverted)
{  int a,b,x,y;
   unsigned long color_sum=0, color_num=0;
   baseline_peaks=NULL;
   width=cd_bitmapwidth[bitmap];
   height=cd_bitmapheight[bitmap];
   /* allocate memory for the image's pixels */
   data=new COLOR_TYPE *[width];
   if (!data) return(0); /* memory allocation failed */
   for (a=0;a<width;a++)
   {  data[a]=new COLOR_TYPE[height];
      if (!data[a]) return (0);   /* memory allocation failed */
   }
   /* allocate memory for the image's pixels identified as orange */
   starspix=new COLOR_TYPE *[width];
   if (!starspix) return(0); /* memory allocation failed */
   for (a=0;a<width;a++)
   {  starspix[a]=new COLOR_TYPE[height];
      if (!starspix[a]) return (0);   /* memory allocation failed */
      for (b=0;b<height;b++) starspix[a][b]=0; /* initialize */
   }
   /* load the picture */
   for (y=0;y<height;y++)
     for (x=0;x<width;x++)
     {  data[x][y]=(COLOR_TYPE)(0.3333*cd_rawget(bitmap,x,y,0)+0.3333*cd_rawget(bitmap,x,y,1)+0.3333*cd_rawget(bitmap,x,y,2));
        if (inverted) data[x][y]=255-data[x][y];
        starspix[x][y]=data[x][y];
        /* calculate the average color of the image's background */
        if (x>width*0.25 && x<width*0.75 && y>height*0.25 && y<height*0.75 && data[x][y]<200)
        {  color_sum=color_sum+data[x][y];
           color_num=color_num+1;
        }
     }
   AvgColor=(color_sum/color_num);   /* set the average color */
   return(1);
}
/* **************************************************************************
                 *************** FreeImage *****************
 free the ImageMatrix object with all memory allocated                   
 this function should be called when the object is not needed any longer 
*/
void ImageMatrix::FreeImage()
{  int a;
   /* free "starpix" array */
   if (starspix)
   {  for (a=0;a<width;a++)
        delete starspix[a];
      delete starspix;
   }
   /* free "data" array */
   if (data)
   {  for (a=0;a<width;a++)
        delete data[a];
      delete data;
   }
   /* free the baseline_peaks */
   if (baseline_peaks) delete baseline_peaks;
   /* initialize values */
   data=NULL;
   starspix=NULL;
   baseline_peaks=NULL;
}
/**************************************************************************
            ******************* LoadFits ********************
   loads a fits file into the ImageMatrix 
   input: fits_file_name - char * - a valid name of a fits image file.
   returned value: integer - 1 if successful.
   The function opens the fits file specified by "fits_file_name" and loads
   its data to the matrix data structure.
   It also calculates the midean pixel value and stores it in "AvgColor" field.
*/

/* sort_function
   this function is used for the qsort function call that sorts the pixels' 
   values by their brightness so the midean value can be found 
   the function gets to integer values and returns 1 if a>b or -1 if b>a
   the function returns 0 if the two integers are equal.
*/

int sort_function(const void *a, const void *b)
{  return((*((unsigned short *)(a))>*((unsigned short *)(b)))*2-1);
}

int ImageMatrix::LoadFits(char *fits_file_name)
{  int a,b;
   COLOR_TYPE val;
   unsigned long color_sum=0, color_num=0;
   fitsfile *fptr;  /* pointer to the fits file. defined in fitsio.h */
   int status,ii,jj;
   int image_type,anynull,simple,bitpix,extend;
   long pcount,gcount;
   COLOR_TYPE *short_array;
   long nelements;
   int naxis = 2;
   long naxes[naxis];
   char value[FLEN_VALUE],comment[FLEN_COMMENT],card[FLEN_CARD];

   strcpy(name,fits_file_name);
   JulianDate=GetFileDate(fits_file_name);

   /* open the fits file */
   status=0;    /* very important! fitsio won't work when status is not initialized */
   fits_open_image(&fptr,fits_file_name,0,&status);
   if (status!=0)            /* failed in openning the fits file */
   {  printf("Can not open file '%s'\n",fits_file_name);
      return(0);
   }
   fits_read_imghdr(fptr,MAX_COMPRESS_DIM,&simple,&bitpix,&naxis,naxes,&pcount,&gcount,&extend,&status);
   if (status!=0)            /* failed to read the header */
   {  printf("Can not open image header. file '%s'\n",fits_file_name); 
      return(0); 
   }
   
   /* read keywords */
   ffgkey(fptr,"EXPTIME", value, comment, &status);
   exposure_time=atof(value);
   status=0;  /* in case the keyword wasn't found */

   /* allocate memory for the image */
   baseline_peaks=NULL;
   width=naxes[0];
   height=naxes[1];
   data=new COLOR_TYPE *[width];
   if (!data) return(0);        /* memory allocation failed */
   for (a=0;a<width;a++)
   {  data[a]=new COLOR_TYPE[height];
      if (!data[a]) return(0);  /* failed to allocate memory */
   }
   
   /* allocate memory for the pixels identified as stars */  
   starspix=new COLOR_TYPE *[width];
   if (!starspix) return(0);     /* failed to allocate memory */
   for (a=0;a<width;a++)
   {  starspix[a]=new COLOR_TYPE[height];
      if (!starspix[a])          /* failed to allocate memory */
      {  printf("Out of memory\n");
         return(0);
      } 
      for (b=0;b<height;b++) starspix[a][b]=0;  /* initialize */
   }
   
   /* read the image and load the pixels' values into "data" and "starpix" arrays */
   short_array=new COLOR_TYPE[width*height];
   fits_read_img(fptr,TUSHORT,1,width*height,NULL,short_array,NULL,&status);
   fits_close_file(fptr,&status);

   /* load the fits image into the data structure */
   for (jj = 0;jj < height; jj++)
     for (ii = 0;ii < width; ii++)
     {  val=short_array[(jj)*width + ii];
        if (INVERTED_IMAGE==1) val=65535-val;
        data[ii][height-jj-1]=val;
        starspix[ii][height-jj-1]=data[ii][height-jj-1];
        /* add the colors for getting the average color */
//        if ((ii>width*0.25 && ii<width*0.75 && jj>height*0.25 && jj<height*0.75) && (val<3000))
//        {  color_sum=color_sum+val;
//           color_num=color_num+1;
//        }
     }
//data[400][700]=1900;data[400][699]=1908;
//data[399][700]=2000;data[399][699]=2102;
//data[401][700]=2050;data[401][699]=2052;
//data[400][690]=1890;data[400][689]=1905;
//data[399][690]=2100;data[399][689]=1894;
//data[401][690]=2100;data[401][689]=1894;

   /* get the midean pixel value:
      Sorting the whole pixels is CPU consuming so only the middle 40000 pixels are sorted.
      The midean value of the 40000 pixels at the middle of the image proved to be very close to the "true" midean value 
   */
   qsort(&(short_array[(int)(width*height/2)-20000]),40000,sizeof(unsigned short),sort_function);
   AvgColor=short_array[(int)(width*height/2)];   /* the midean value is the value at the entry that is at the middel of the sorted range */
   delete short_array;
//   if (color_num==0) color_num=1;  /* prevent division by zero           */
//   AvgColor=(color_sum/color_num); /* set the average color of the image */
//   if (AvgColor==0) AvgColor=1;    /* prevent division by zero           */
   return(1);                        /* loading fits image aucceeded         */
}

/* ***********************************************************************
   ***************************** SaveFits ********************************
   purpose: save the data hold in "data" array to a fits file
   parameters: fits_file_name - the name of the new fits file 
*/ 
int ImageMatrix::SaveFits(char *fits_file_name)
{  fitsfile *fptr;       /* pointer to the FITS file; defined in fitsio.h */
   int status, ii, jj;
   int image_type,anynull,simple,bitpix,extend;
   long pcount,gcount;
   unsigned short *short_array;
   long  nelements,fpixel = 1,naxis = 2, exposure;
   long naxes[2];
   COLOR_TYPE val;

   naxes[0]=width;   /* set the image's width  */
   naxes[1]=height;  /* set the image's height */
   status=0;         /* improtant!             */
   remove(fits_file_name); /* make sure the file does not exist (or else it will fail creating the file) */
   fits_create_file(&fptr, fits_file_name, &status);  /* create a new file */
   if (status) return(0);
   fits_create_img(fptr, USHORT_IMG, naxis, naxes, &status); /* USHORT_IMG */
   if (status) return(0);
   exposure=1500;
   if (status) return(0);
   short_array=new unsigned short [height*width];
   /* save the image data into the fits file */
   for (jj = 0;jj < height;jj++)
     for (ii = 0;ii < width;ii++)
       {  val=data[ii][height-jj-1];
       if (INVERTED_IMAGE==1) val=65535-val;
       short_array[(jj)*width + ii]=val;
       }
   nelements = naxes[0] * naxes[1];
   fits_write_img(fptr, TUSHORT, 1, nelements, short_array, &status);  /* TUSHORT */
   fits_close_file(fptr, &status);
   delete short_array;
   return(1);  /* successfuly saved */
}

/************************************************************************
           ***************** GetPicture ****************
   purpose: Load the image into a bitmap image
   input:  int - a handle to a BITMAP structue
   output: int - 1 if succeeded
   The function copies the data stored in the matrix to a bitmap structure
   so it can later be saved to the disk as a .bmp file.
*/
int ImageMatrix::GetPicture(BITMAP h_bitmap)
{  int x,y,val;
   double AvgSig;
   AvgSig=sqrt(AvgColor);
   for (y=0;y<height;y++)
     for (x=0;x<width;x++)
     {  val=data[x][y];
//        if (val-AvgColor>20*AvgSig) val=255;
//        else
//	if (val<AvgColor) val=70-(int)((((double)(AvgColor-val))/AvgColor)*70);
//        else val=70+(int)((((double)(val-AvgColor))/(20*AvgSig))*185);
        if (val<=AvgColor*0.4) val=0;
        else
	if (val>=1.4*AvgColor) val=255;  /* 3 */
        else
	val=(unsigned char)(((val-AvgColor*0.4)*255)/(1.0*AvgColor)); /* 2.5 */
        if (val>255) val=255;
        cd_rawassign(h_bitmap,x,y,0,(unsigned char)val);
        cd_rawassign(h_bitmap,x,y,1,(unsigned char)val);
        cd_rawassign(h_bitmap,x,y,2,(unsigned char)val);
     }
   return(1);
}
/**************************************************************************
         ****************** ClosestStar *******************
   purpose: finds he closest star to the given x,y points.
   inputs: x,y - coodinates within the matrix.
   output: the distance of the closest star (peak) identified in the image.
   The function goes over the list of "star_points" and get the closest star 
   to the x,y coordinates in means of Euclidean distance.
*/
double ImageMatrix::ClosestStar(long x,long y)
{  int star_point_index,min_star,a;
   double dist,min_dist;
   star_point_index=0;
   min_dist=100000;       // min_dist <- infinite value
   while (star_points[star_point_index].x>=0)
   {
      {  dist=sqrt(pow(star_points[star_point_index].x-x,2)+pow(star_points[star_point_index].y-y,2));  // get the Eucldean distance
         if (dist<min_dist)
         {  min_dist=dist;
            min_star=star_point_index;
         }
      }
      star_point_index++;
   }
   return(min_dist);
}

/* *******************************************************************
   **************************** FindArea *****************************
   purpose: Find an area covered with values which are higher than the
   background.
   mark area */
pixel_item *ImageMatrix::AddMarkPix(pixel_item *p_tail,short x,short y,long *counter)
{  pixel_item *p_new_item;
   p_new_item=new pixel_item;
   p_new_item->x=x;
   p_new_item->y=y;
   p_new_item->p_next=NULL;
   p_tail->p_next=p_new_item;
   starspix[x][y]=0;
   *counter=*counter+1;
   return(p_new_item);
}

long ImageMatrix::FindArea(short x,short y,short value,short *min_x,short *max_x,short *min_y,short *max_y)
{ pixel_item *p_head,*p_tail,*p_temp;
  long counter=0;
  if (starspix[x][y]<value) return(0);
  p_head=new pixel_item;
  p_head->x=x;
  p_head->y=y;
  p_head->p_next=NULL;
  p_tail=p_head;
  while (p_head)
  {  x=p_head->x;
     y=p_head->y;
     if (x<*min_x) *min_x=x;
     if (x>*max_x) *max_x=x;
     if (y<*min_y) *min_y=y;
     if (y>*max_y) *max_y=y;
     if (x<width-1) if (starspix[x+1][y]>=value) p_tail=AddMarkPix(p_tail,x+1,y,&counter);
     if (x>0) if (starspix[x-1][y]>=value) p_tail=AddMarkPix(p_tail,x-1,y,&counter);
     if (x>0 && y>0) if (starspix[x-1][y-1]>=value) p_tail=AddMarkPix(p_tail,x-1,y-1,&counter);
     if (x<width-1 && y>0) if (starspix[x+1][y-1]>=value) p_tail=AddMarkPix(p_tail,x+1,y-1,&counter);
     if (y>0) if (starspix[x][y-1]>=value) p_tail=AddMarkPix(p_tail,x,y-1,&counter);
     if (x>0 && y<height-1) if (starspix[x-1][y+1]>=value) p_tail=AddMarkPix(p_tail,x-1,y+1,&counter);
     if (x<width-1 && y<height-1) if (starspix[x+1][y+1]>=value) p_tail=AddMarkPix(p_tail,x+1,y+1,&counter);
     if (y<height-1) if (starspix[x][y+1]>=value) p_tail=AddMarkPix(p_tail,x,y+1,&counter);
     starspix[p_head->x][p_head->y]=0; /* mark the pixel */
     p_temp=p_head;
     if (p_head->p_next) p_head=p_head->p_next;
     else p_head=NULL;
     delete p_temp;
   }
   return(counter);
}

//-----------------------------------------------------------------------
/* IsStar 
   purpose: return 1 if the x,y coordinates are star cooridinates 
   inputs: x,y coordinates of a pixel that should be checked.
   output: 1 if the x,y coordinates are a peak (a star).
   The function "checks" around the pixel to find out whether it is 
   substainsially brighter than all its surounding pixels.
*/

unsigned short ImageMatrix::IsStar(short x,short y)
{  int eqnum=0,lsnum=0,lssum=0;
   long stars_color,stars_background;
   double sqrt_AvgColor;
#ifdef FITS
   if (x<5 || x>width-5 || y<5 || y>height-5) return(0);
   /* check minimum color */
   // stars_color=GetStarsColor(x,y);
   stars_color=(data[x][y]+data[x-1][y]+data[x+1][y]+data[x][y-1]+data[x][y+1])/5;
   if (stars_color<AvgColor*1.05) return(0);
   /* make sure no memory access violation will occur */
   /* check difference from the background */
#define min(a,b) (((a)<(b)) ? (a):(b))
   sqrt_AvgColor=sqrt(AvgColor);
   if (data[x][y]-min(min(data[x-4][y],data[x-3][y]),data[x-2][y]) < 5*sqrt_AvgColor) lssum++;
   if (data[x][y]-min(min(data[x+4][y],data[x+3][y]),data[x+2][y]) < 5*sqrt_AvgColor) lssum++;
   if (data[x][y]-min(min(data[x][y-4],data[x][y-3]),data[x][y-2]) < 5*sqrt_AvgColor) lssum++;
   if (data[x][y]-min(min(data[x][y+4],data[x][y+3]),data[x][y+2]) < 5*sqrt_AvgColor) lssum++;
   if (lssum>0) return(0);
   /* check if the star's environment is also higher than the average background */
   //   if ((data[x-1][y]+data[x+1][y]+data[x][y-1]+data[x][y+1])/4<=1.12*AvgColor) return(0);
   /* check if higher than its surrounding */
   if (data[x-1][y]<=data[x][y])     eqnum++; /* left         */
   if (data[x+1][y]<=data[x][y])     eqnum++; /* right        */
   if (data[x][y-1]<=data[x][y])     eqnum++; /* top          */
   if (data[x][y+1]<=data[x][y])     eqnum++; /* bottom       */
   if (data[x-1][y-1]<=data[x][y])   eqnum++; /* upper left   */
   if (data[x-1][y+1]<=data[x][y])   eqnum++; /* lower left   */
   if (data[x+1][y-1]<=data[x][y])   eqnum++; /* right upper  */
   if (data[x+1][y+1]<=data[x][y])   eqnum++; /* lower right  */
   if (eqnum<8) return(0); /* true only if all the surrounding pixel are less bright */
   stars_background=GetStarsBackground(x,y);
   if (stars_background==0) return(0); /* probably a bad pixel outside the horizon */
   if (data[x][y]-stars_background<15*sqrt(stars_background)) return(0); /* ignore peaks less than 15 sigmas */
   if (stars_background>AvgColor*2) return(0); /* probably a cloud */
//   if (stars_background==0 || (float)stars_color/(float)stars_background<1.2
//       || (((float)(5*stars_color-data[x][y]))/4)/(float)stars_background<1.1) return(0); /* make sure it's not a bad pixel */
   return(1);
#else
   if (data[x][y]<MIN_COLOR) return(0);
   if (sqrt(pow(x-MID_X,2)+pow(y-MID_Y,2))<CENTER_COLOR_DIST && data[x][y]<MIN_CENTER_COLOR) return(0);
   /* check if the star is different from it's area */
   lssum=0; 
   if (data[x][y]-data[x-4][y]<MIN_DEF_COLOR) lssum++;
   if (data[x][y]-data[x+4][y]<MIN_DEF_COLOR) lssum++;
   if (data[x][y]-data[x][y+4]<MIN_DEF_COLOR) lssum++;
   if (data[x][y]-data[x][y-4]<MIN_DEF_COLOR) lssum++;
   if (lssum>0) return(0);
//   if (((float)(data[x][y]+data[x+1][y]+data[x][y+1]+data[x][y-1]))/5>MIN_AVERAGE_COLOR) return(0);
   /* check if it is one pixel only */
//   lssum=0;
//   lssum=lssum+data[x][y]-data[x-1][y];
//   lssum=lssum+data[x][y]-data[x+1][y];
//   lssum=lssum+data[x][y]-data[x][y-1];
//   lssum=lssum+data[x][y]-data[x][y+1];
//   if (lssum/4>12) return(0);
   lssum=0;
   /* check if it is a star */
   if (data[x-2][y]<=data[x][y])            /* left       */
   {  lsnum++;
      lssum=lssum+data[x][y]-data[x-2][y];
   }
   if (data[x-2][y]==data[x][y]) eqnum++;   /* left       */
   if (data[x+2][y]<=data[x][y])            /* right      */
   {  lsnum++;
      lssum=lssum+data[x][y]-data[x+2][y];
   }
   if (data[x+2][y]==data[x][y]) eqnum++;   /* right      */
   if (data[x][y-2]<=data[x][y])            /* up         */
   {  lsnum++;
      lssum=lssum+data[x][y]-data[x][y-2];
   }
   if (data[x][y-2]==data[x][y]) eqnum++;   /* up         */
   if (data[x][y+2]<=data[x][y])            /* down       */
   {  lsnum++;
      lssum=lssum+data[x][y]-data[x][y+2];
   }
   if (data[x][y+2]==data[x][y]) eqnum++;   /* down       */
   if (data[x-2][y-2]<=data[x][y])          /* up left    */
   {  lsnum++;
      lssum=lssum+data[x][y]-data[x-2][y-2];
   }
   if (data[x-2][y-2]==data[x][y]) eqnum++; /* up left    */
   if (data[x+2][y-2]<=data[x][y])          /* up right   */
   {  lsnum++;
      lssum=lssum+data[x][y]-data[x+2][y-2];
   }
   if (data[x+2][y-2]==data[x][y]) eqnum++; /* up right   */
   if (data[x-2][y+2]<=data[x][y])          /* down left  */
   {  lsnum++;
      lssum=lssum+data[x][y]-data[x-2][y+2];
   }
   if (data[x-2][y+2]==data[x][y]) eqnum++; /* down left  */
   if (data[x+2][y+2]<=data[x][y])          /* down right */
   {  lsnum++;
      lssum=lssum+data[x][y]-data[x+2][y+2];
   }
   if (data[x+2][y+2]==data[x][y]) eqnum++; /* down right */
   if (lsnum==8 && eqnum<=2 && lssum/8>30) return(1);
   if (lsnum==8 && eqnum<=1 && lssum/8>25) return(1);
   if (lsnum==8 && eqnum==0 && lssum/8>15) return(1);
   return(0);
#endif
}

/*************************************************************************************/
/* ********************************* GetStarBackground ******************************
   purpose: Gets the background count for a specific pixel (usually a star pixel)
   input: x,y: The x,y coordinates of the pixel
   comments: The background count is computed by finding the midean value of the area
             sourounding the star.
             this function also uses the "sort_function" function used by "LoadFits" for
             sorting 
*/
long ImageMatrix::GetStarsBackground(long x,long y)
{  long xc,yc,bg;
   unsigned short half_width_size=10;
   unsigned short pixel_counts[2*10*2*10];
   if (x<=half_width_size || y<=half_width_size || x>=width-half_width_size || y>=height-half_width_size) return(0);
   /* get the pixel counts around the star and put them in an array */
   for (yc=y-half_width_size;yc<y+half_width_size;yc++)
     for (xc=x-half_width_size;xc<x+half_width_size;xc++)
       pixel_counts[(half_width_size*2)*(yc-(y-half_width_size))+(xc-(x-half_width_size))]=data[xc][yc];  
   /* sort the array for finding the midean value */
   qsort(pixel_counts,2*half_width_size*2*half_width_size,sizeof(unsigned short),sort_function);  
   bg=pixel_counts[(long)((float)(half_width_size*2*half_width_size*2)/2)];
   if (bg==0) bg=1;
   return(bg);
}

/*************************************************************************************/
/* *********************************** GetStarSurface **********************************
   purpose: Gets size of the star surface (number of pixels) and the avegare and maximum
            edge.
   input: x,y: The x,y coordinates of the pixel
          top_pix_count: the number of top pixel to average from (possibly 5,9 or 16)
*/
void ImageMatrix::GetStarSurface(long x,long y,long base_x,long base_y,unsigned short background,unsigned short sqrt_background,unsigned short *pixel_count,unsigned short *edge_count,unsigned short *sum_edge,unsigned short *maximum_edge,COLOR_TYPE *pixel_counts)
{  long edge_pix_value=0;
#define MAX_DIST_FROM_PEAK 3
   if (*pixel_count>250) return;
   if (x<base_x-MAX_DIST_FROM_PEAK) return;  /* stop when X pixels away from the peak of the star */
   if (x>base_x+MAX_DIST_FROM_PEAK) return;  /* stop when X pixels away from the peak of the star */
   if (y<base_y-MAX_DIST_FROM_PEAK) return;  /* stop when X pixels away from the peak of the star */
   if (y>base_y+MAX_DIST_FROM_PEAK) return;  /* stop when X pixels away from the peak of the star */
   if (x<0 || x>=width || y<0 || y>=height) return;
   if (starspix[x][y]==1) return;  /* the pixel was already counted */
   if (x<5 || x>width-5) return;
   if (y<5 || y>height-5) return;
   starspix[x][y]=1;
   pixel_counts[(*pixel_count)++]=data[x][y];

   edge_pix_value=data[x][y]-background;
   /* check if this is an "edge" pixel */
   if ((data[x-1][y]-background<9*sqrt_background) && starspix[x-1][y]==0)
   {  if (edge_pix_value>*maximum_edge) *maximum_edge=edge_pix_value;
      (*sum_edge)+=edge_pix_value;
      (*edge_count)++;
      starspix[x-1][y]=1;
   }
   else GetStarSurface(x-1,y,base_x,base_y,background,sqrt_background,pixel_count,edge_count,sum_edge,maximum_edge,pixel_counts);

   /* check if this is an "edge" pixel */
   if (data[x+1][y]-background<9*sqrt_background && starspix[x+1][y]==0)
   {  if (edge_pix_value>*maximum_edge) *maximum_edge=edge_pix_value;
     (*sum_edge)+=edge_pix_value;
     (*edge_count)++;
     starspix[x+1][y]=1;
   }
   else GetStarSurface(x+1,y,base_x,base_y,background,sqrt_background,pixel_count,edge_count,sum_edge,maximum_edge,pixel_counts);

   /* check if this is an "edge" pixel */
   if (data[x][y-1]-background<9*sqrt_background && starspix[x][y-1]==0)
   {  if (edge_pix_value>*maximum_edge) *maximum_edge=edge_pix_value;
     (*sum_edge)+=edge_pix_value;
     (*edge_count)++;
     starspix[x][y-1]=1;
   }
   else GetStarSurface(x,y-1,base_x,base_y,background,sqrt_background,pixel_count,edge_count,sum_edge,maximum_edge,pixel_counts);

   /* check if this is an "edge" pixel */
   if (data[x][y+1]-background<9*sqrt_background && starspix[x][y+1]==0)
   {  if (edge_pix_value>*maximum_edge) *maximum_edge=edge_pix_value;
      (*sum_edge)+=edge_pix_value;
      (*edge_count)++;
      starspix[x][y+1]=1;
   }
   else GetStarSurface(x,y+1,base_x,base_y,background,sqrt_background,pixel_count,edge_count,sum_edge,maximum_edge,pixel_counts);
   return;
}


/*************************************************************************************/
/* *********************************** GetStarsColor **********************************
   purpose: Gets the foreground count for a specific pixel (usually a star pixel)
   input: x,y: The x,y coordinates of the pixel
          top_pix_count: the number of top pixel to average from (possibly 5,9 or 16)
*/
long ImageMatrix::GetStarsColor(long x,long y,short top_pix_count,unsigned short background)
{  long xc,yc,num,sum;
   COLOR_TYPE pixel_count[256];
   unsigned short surface_size=0,sum_edge=0,edge_count=0,edges_average=0,max_edge=0;
   double res;
   unsigned short sqrt_background=(unsigned short)(sqrt(background));
   
   if (x<8 || x>width-8 || y<8 || y>height-8) return(0);  /* too close to the edge */
   /* check if the peak is a cosmic ray hit */
   for (xc=0;xc<256;xc++) pixel_count[xc]=0;
   for (xc=x-6;xc<x+6;xc++)
     for (yc=y-6;yc<y+6;yc++) starspix[xc][yc]=0;
   GetStarSurface(x,y,x,y,background,sqrt_background,&surface_size,&edge_count,&sum_edge,&max_edge,pixel_count);
   if (edge_count>0) edges_average=sum_edge/edge_count;
   /* ignore cosmic ray hits */
   if (!strstr(name,"clear")) /* don't activate this algorithm for clear frames */
   {  res=GamaRaysModel.CalculateRules(surface_size,edges_average/sqrt_background,max_edge/sqrt_background);
      if (res>=0.7) return(0);   /* detected as a cosmic ray hit */
   }

   /* initialize the pixels' count array */
   if (surface_size>256) surface_size=256;
   qsort(pixel_count,surface_size,sizeof(COLOR_TYPE),sort_function);
   sum=num=0;
   while (num<top_pix_count && surface_size-num-1>=0 && pixel_count[surface_size-num-1]>0)
   { sum=sum+pixel_count[surface_size-num-1];
     num++;
   }
   if (num<top_pix_count) sum+=background*(top_pix_count-num); /* if the nearby pixels are lower than the background - add the background */
   num=top_pix_count;  /* devide always with the with the number of top pixels e.g. C16 - devide always by 16*/
   if (num==0) return(0);  /* prevent division by zero possible error */
   else return((long)(sum/num));
}


/* *************************************************************************** */
/* ******************************  GetStarPoints  **************************** 
   purpose: find a list of all peaks in the image. Each peak is expected to 
   be a star or some other astronomic object 
   inputs: star_points - an array that will store the peaks. This array should be
                         empty when calling the function, but when the function
                         ends, it contain all the stars (peaks) that were found.
           stripe_point - the same like "star_points" but for stripes.
           max_dist     - the maximum distance from the center of the image. This 
                          parameter is used so pixels outside the hemisphere will
                          be ignored.
    The fits version of the function ignores stripes.
    The function searches for peaks and for each peak checks if it is not too far
    from the center of the image.   
*/
long ImageMatrix::GetStarsPoints(star_point *star_points,stripe_point *stripes,double max_dist)
{  short x,y,min_x,max_x,min_y,max_y,value;
   int pix_num,a,add,star_point_index=0,stripe_index=0;
   unsigned short star_color;
   double dist,stripe_length;
   min_x=32000;
   max_x=-1;
   min_y=32000;
   max_y=-1;
   /* find the very bright stars and stripes */
   for (y=1;y<height-1;y++)
     for (x=1;x<width-1;x++)
     {  if (x>current_location.bad.x1 && x<current_location.bad.x2
	   && y>current_location.bad.y1 && y<current_location.bad.y2)
             continue; /* ignore the bad area */
        pix_num=FindArea(x,y,(long)(AvgColor+40*sqrt(AvgColor)),&min_x,&max_x,&min_y,&max_y);
        if (pix_num>4)
        {  dist=sqrt(pow(((max_x+min_x)/2)-MID_X,2)+pow(((max_y+min_y)/2)-MID_Y,2));
           stripe_length=sqrt(pow(max_x-min_x,2)+pow(max_y-min_y,2));
           if (pix_num<300 && pix_num>4 && (dist+0.5*(stripe_length-2)<max_dist*0.8))
	   {  /* check if it is a stripe */
              if (stripes && pix_num<0.2*((0.5*stripe_length)*(0.5*stripe_length)*3.14))
	      {  /* this is a stripe */
	         stripes[stripe_index].x=(max_x+min_x)/2;
	         stripes[stripe_index].y=(max_y+min_y)/2;
	         stripes[stripe_index].min_x=min_x;
	         stripes[stripe_index].max_x=max_x;
	         stripes[stripe_index].min_y=min_y;
	         stripes[stripe_index].max_y=max_y;
	         if (stripe_index<STRIPES_NUM-2) stripe_index++;
	      }
#ifndef FITS
           /* a "normal" star */
           /* check if the star was already added */
           add=1;
           a=star_point_index-1;
           while(a>=0)
	   {  dist=sqrt(pow(star_points[a].x-x,2)+pow(star_points[a].y-y,2));
	      if (dist<3 || (dist<6 && star_points[a].bright_star==1))
	      {  add=0;
	         break;
	      }
	      a--;
	   }
           /* add the star */
           if (add && star_point_index<MAX_STARS_POINTS)
	   {  star_points[star_point_index].x=(max_x+min_x)/2;
	      star_points[star_point_index].y=(max_y+min_y)/2;
	      star_points[star_point_index].bright_star=1;  /* a bright star */
	      star_points[star_point_index].color=data[x][y];
	      star_points[star_point_index].TimesAvg=(float)data[x][y]/(float)AvgColor;
	      star_points[star_point_index].cloud=0;
	      star_points[star_point_index].cloud_opacity=0;
	      star_points[star_point_index].cloud_glow=0;
              star_points[star_point_index].star=NULL;
	      star_point_index++;
	   }
#endif
	   }
        }
        /* initialize the minimum and maximum values */
        min_x=32000;
        max_x=-1;
        min_y=32000;
        max_y=-1;
     }
   /* search for peaks */
   {  for (y=1;y<height-1;y++)
        for (x=1;x<width-1;x++)
//        if (starspix[x][y]>0) /* this pixel wasn't calculated in the bright star search */
	  {  if (star_point_index>=MAX_STARS_POINTS-1) break;  /* don't exceed the allocated memory buffer */
	     if (x>current_location.bad.x1 && x<current_location.bad.x2
	         && y>current_location.bad.y1 && y<current_location.bad.y2)
	           continue; /* ignore the bad area */
           /* check if the star is not too far from the middle */
           dist=sqrt(pow(x-MID_X,2)+pow(y-MID_Y,2));
           if (dist<max_dist && IsStar(x,y))        /* if the pixel is a peak and not too far from the center of the image */
	   {  unsigned short background;
              /* check colors and cosmic rays */
	      add=1;
              background=GetStarsBackground(x,y);
	      star_color=GetStarsColor(x,y,5,background);
              if (star_color==0) add=0;  /* the star is detected as a cosmic ray hit */
              a=star_point_index-1;
//#ifndef FITS
              /* make sure that a star with two peaks won't be counted as two stars */
              while (a>=0)
              {  dist=sqrt(pow(star_points[a].x-x,2)+pow(star_points[a].y-y,2));
                 if (dist<3)
                 { add=0;
                   break;
                 }
                 a--;
              }
//#endif
              /* add the star */
              if (add)
	      { /* add the star's attributes */
		 /* the color attribute is the avarage of the 5 pixel counts around the pixel's location */
                 star_points[star_point_index].x=x;
                 star_points[star_point_index].y=y;
                 star_points[star_point_index].color=star_color;
		 star_points[star_point_index].C1=GetStarsColor(x,y,1,background);
		 star_points[star_point_index].C5=star_color;//GetStarsColor(x,y,5);
		 star_points[star_point_index].C9=GetStarsColor(x,y,9,background);
		 star_points[star_point_index].C16=GetStarsColor(x,y,16,background);
		 star_points[star_point_index].C25=GetStarsColor(x,y,25,background);
                 star_points[star_point_index].background=background;
                 star_points[star_point_index].bright_star=NOT_BRIGHT;
                 star_points[star_point_index].cloud=0;
		 star_points[star_point_index].star=NULL;
                 star_point_index++;
              }
              /* initialize the minimum and maximum values (for the next iteration) */
              min_x=32000;
              max_x=-1;
              min_y=32000;
              max_y=-1;
           }
        }
   }
   /* set a sign for the last star point */
   star_points[star_point_index].x=-1;
   star_points[star_point_index].y=-1;
   /* set a sign for the last stripe point */
   if (stripes)
   {  stripes[stripe_index].x=-1;
      stripes[stripe_index].y=-1;
   }
   return(star_point_index-1);  /* return the total number of star points found */
}

/* *************************************************************
   ********************** IsCloudArea ************************** 
   return 1 if the area is a cloud or 0 if the area is not a cloud 
   inputs: x,y coordinates.
   output: 1 if the x,y coordinates are inside a cloud.
   The function checks the area around the x,y coordinates and checks the average 
   of the pixels' values in that area. If the average is too high - the area is a cloud.
*/

int ImageMatrix::IsCloudArea(int x, int y)
{  long xc,yc,white_pix_num=0;
   if (x<CLOUD_AREA_SIZE) x=CLOUD_AREA_SIZE;
   if (y<CLOUD_AREA_SIZE) y=CLOUD_AREA_SIZE;
   if (x>width-CLOUD_AREA_SIZE)  x=width-CLOUD_AREA_SIZE-1;
   if (y>height-CLOUD_AREA_SIZE) y=height-CLOUD_AREA_SIZE-1;
   for (yc=y-CLOUD_AREA_SIZE;yc<y+CLOUD_AREA_SIZE;yc++)
     for(xc=x-CLOUD_AREA_SIZE;xc<x+CLOUD_AREA_SIZE;xc++)
       if (data[xc][yc]-AvgColor>sqrt(AvgColor)*MIN_CLOUD_VALUE) white_pix_num++;
   if (white_pix_num>(CLOUD_AREA_SIZE*2*CLOUD_AREA_SIZE*2)/2) return(1);   /* the area is a cloud */
   return(0);     /* the area is not a cloud */
}


/* ************************************************************************
                     ********** DrawLabels ************                  
purpose: Draws a list of labesl (using the "labesl" data structure on the image
         This is done after all labels are collected.
parameters: bitmap - BITMAP: hanlde of the bitmap to be drawn.
            resize - short: the resize scale of the output image.
            right - short: number of pixels to type the label at the right of the object
            below - short: number of pixel to write the label below the object
*/
void ImageMatrix::DrawLabels(BITMAP bitmap, short resize,short right, short below,FONT font)
{  long LabelIndex=0,x,y,star_type;
   unsigned char aqua_style[10],stars_style[10],const_style[10],planets_style[10],sky_objects_style[10],white_style[10],unexpected_style[10],*style;
   char display_name[50];

   /* set the styles to be used */
   cd_setstyle(aqua_style,0,0,0,0,0,0,0,0);                  /* no color                        */
   cd_setstyle(stars_style,0,0,255,255,255,0,1,0);           /* stars' names font color style   */
   cd_setstyle(const_style,0,155,0,255,255,255,1,0);         /* constellations names font color */
   cd_setstyle(planets_style,255,0,0,0,255,255,1,0);         /* planet names font color         */
   cd_setstyle(sky_objects_style,255,0,255,255,255,255,1,0); /* sky object (from external file) */
   cd_setstyle(white_style,255,255,255,255,255,255,1,0);     /* white font color - for image name, time etc' */
   cd_setstyle(unexpected_style,255,255,0,255,255,255,1,0);  /* unexpected object */

   while (labels[LabelIndex].x>0)
   {  labels[LabelIndex].text[49]='\0';   /* prevent memory glitches. make sure copy block is limitted */
      strcpy(display_name,labels[LabelIndex].text);
      x=labels[LabelIndex].x;
      y=labels[LabelIndex].y;
      if (resize!=100)
      {  x=(long)(x*resize)/100;
         y=(long)(y*resize)/100;
      }
      star_type=labels[LabelIndex].type;
      /* select the font according to the type of the object */
      if (star_type==STAR_TYPE_CONSTELLATION) style=const_style;
      if (star_type==STAR_TYPE_BRIGHTSTAR) style=stars_style;
      if (star_type==STAR_TYPE_PLANET) style=planets_style;
      if (star_type==STAR_TYPE_SKY_OBJECT) 
      {  cd_rect(bitmap,x-1,y-1,x+1,y+1,1,sky_objects_style,aqua_style);
         style=sky_objects_style;  /* make sure the label will be at the same color */
      }
      if (star_type==STAR_TYPE_UNEXPECTED)
      {  cd_rect(bitmap,x-15,y-15,x+15,y+15,1,unexpected_style,aqua_style);
         goto end_loop; /* don't draw the label for unexpected objects */
      }
      /* draw the label */
      // cd_plot(bitmap,x,y,style);
      cd_drawtext(bitmap,font,x+right,y+below,style,display_name,NO,NO,NO);
end_loop:
      LabelIndex++;
   }
} 

/**************************************************************************
         ***************** LoadBaselinePeaks *******************
This function loads the peaks from an external file that was pre-prepared
based on some known "good" image.
The data is stored in "baseline_peaks" attribute and being used to compare
the peaks of the current image (stars_points) with this list of peaks in 
order to detect "unexpected" stars or magnitude changes.
inputs - filename the filename of the baseline peaks.
output - 1 if succeeded.
*/
int ImageMatrix::LoadBaselinePeaks(char *filename)
{  int peak_counter=0;
   FILE *peaks_file;
   char *p_one_line,one_line[256];      /* this string should load one line of the text file */
   if (!(peaks_file=fopen(filename,"r"))) return(0);    /* open the peaks text file */  
   baseline_peaks=new star_point[MAX_BASELINE_PEAKS];

   /* read the peaks from the file */ 
   p_one_line=fgets(one_line,255,peaks_file);
   while (p_one_line)
   {  if (strncmp(p_one_line,"//",2)!=0)  /* check if this line is a comment line */
      {  p_one_line=&(one_line[14]); /* skip the name of the star */  
         p_one_line=strtok(p_one_line," \t\n");
	 baseline_peaks[peak_counter].x=atoi(p_one_line);  /* read the x coordinate */
	 p_one_line=strtok(NULL," \t\n");
	 baseline_peaks[peak_counter].y=atoi(p_one_line);  /* read the y coordinate */
	 p_one_line=strtok(NULL," \t\n");                  /* step on the distance from the center */
	 p_one_line=strtok(NULL," \t\n");                  /* step on the distance from the zenith */ 
	 p_one_line=strtok(NULL," \t\n"); 
         baseline_peaks[peak_counter].color=atoi(p_one_line);
         p_one_line=strtok(NULL," \t\n");
	 baseline_peaks[peak_counter].background=atoi(p_one_line);
         baseline_peaks[peak_counter].cloud=0;             /* initlially - not a cloud */
	 baseline_peaks[peak_counter].star=NULL;           /* initially - no stars assigned to points */
         peak_counter++;
      }
      p_one_line=fgets(one_line,255,peaks_file);
   }

   /* set the last peak that indicates on end of the list */
   baseline_peaks[peak_counter].x=-10000;
   baseline_peaks[peak_counter].y=-10000;

   /* exit the function */
   fclose(peaks_file);
   return(1);
}

/* ************************************************************************
         ******************* FindClouds3 *********************
   purpose: detect areas covered with light or heavy clouds and mark them in the image.
   input - star_points - a list of all stars appears in the image.
           baseline_image_name - the canonical image file name.
           bitmap - the output image.
   The last star_point should have negative values in x and y fields.
*/

/* Do a linear interpolation using the closest baseline peaks that were found in the area of the pixel. */

short interpolate(long type,long x,long y,star_point *baseline_peaks,double *dists,long *closests)
{ long close_index,opoq=0;
  long left_x=-1,right_x=-1,top_y=-1,bottom_y=-1,min_dist;
  double dist,min_dist_left_x=100000,min_dist_right_x=100000,min_dist_top_y=100000,min_dist_bottom_y=100000;
  long directions=0;
  if (type==1)
  {
     for (close_index=0;close_index<10;close_index++)
       if (closests[close_index]>-1 && baseline_peaks[closests[close_index]].cloud!=-1)
       { dist=pow(pow(baseline_peaks[closests[close_index]].x-x,2)+pow(baseline_peaks[closests[close_index]].y-y,2),0.5);
//	 if (dist>60) continue;
         /* check if there are stars in all directions from that pixel*/
         if (baseline_peaks[closests[close_index]].x<=x) directions=directions | 1; /* left found */
         if (baseline_peaks[closests[close_index]].x>=x) directions=directions | 2; /* right found */
         if (baseline_peaks[closests[close_index]].y<=y) directions=directions | 4; /* top found */
         if (baseline_peaks[closests[close_index]].y>=y) directions=directions | 8; /* bottom found */

	 /* find the closest left point */
	 if (baseline_peaks[closests[close_index]].x<x && dist<min_dist_left_x &&
	 abs(baseline_peaks[closests[close_index]].y-y)<x-baseline_peaks[closests[close_index]].x)
	 {  min_dist_left_x=dist;
	    left_x=close_index;
	 }
	 /* find the closest right point */
	 if (baseline_peaks[closests[close_index]].x>=x && dist<min_dist_right_x &&
	    abs(baseline_peaks[closests[close_index]].y-y)<baseline_peaks[closests[close_index]].x-x)
	 {  min_dist_right_x=dist;
	    right_x=close_index;
	 }
	 /* find the closest top point */
	 if (baseline_peaks[closests[close_index]].y<y && dist<min_dist_top_y &&
	    abs(baseline_peaks[closests[close_index]].x-x)<y-baseline_peaks[closests[close_index]].y)
	 {  min_dist_top_y=dist;
	    top_y=close_index;
	 }
	 /* find the closest bottom point */
	 if (baseline_peaks[closests[close_index]].y>=y && dist<min_dist_bottom_y &&
	    abs(baseline_peaks[closests[close_index]].x-x)<baseline_peaks[closests[close_index]].y-y)
	 {  min_dist_bottom_y=dist;
	    bottom_y=close_index;
	 }
       }
     if (directions!=15) return(255); /* one direction has no stars at all */
     if (right_x==-1) right_x=left_x;
     if (left_x==-1) left_x=right_x;
     if (top_y==-1) top_y=bottom_y;
     if (bottom_y==-1) bottom_y=top_y;
     if (left_x==-1 || top_y==-1) return(255);
     return((unsigned short)((dists[left_x]*baseline_peaks[closests[left_x]].cloud
                       + dists[right_x]*baseline_peaks[closests[right_x]].cloud
                       + dists[top_y]*baseline_peaks[closests[top_y]].cloud
			     + dists[bottom_y]*baseline_peaks[closests[bottom_y]].cloud)
			    / (dists[left_x]+dists[right_x]+dists[top_y]+dists[bottom_y])));
   }

   /* N nearest neighbors without interpolation */
   if (type==2)
   {  double close_sum=0.0,close_num=0.0;
      /* count the non detected stars */
      for (close_index=0;close_index<5;close_index++)
        if (closests[close_index]>-1
        && baseline_peaks[closests[close_index]].cloud>-1 /* ignore dim stars that didn't appear */
        && dists[close_index]<40)                         /* ignore points that are too far      */
        {  close_sum+=baseline_peaks[closests[close_index]].cloud*(dists[close_index]/10);
           close_num+=dists[close_index]/10;
        }
      if (close_num==0) return(255); /* no stars where detected at this area at all */
      else return((unsigned short)(close_sum/close_num));
   }
}

/* a sort function that is given as a parameter to the qsort function */
int sort_by_x(const void *star_point1, const void *star_point2)
{  const star_point *p1=(const star_point *)star_point1;
   const star_point *p2=(const star_point *)star_point2;
   if ((*p1).x>(*p2).x) return(1);
   else return(-1);
}

int ImageMatrix::FindClouds3(site_info *location, star_point* StarPoints,char *baseline_image_name,BITMAP bitmap,short find_changes)
{  ImageMatrix *baseline_image;
   short star_point_index,baseline_point_index,matched_point_index,star_points_num=0;
   long x,y;
   short *sorted_by_x,last_x_pos;
   double dist,scale_item,min_dist,max_dist;
   FONT font;
   FILE *font_file;
   unsigned char font_buffer[10],cloud_style[10];
   /* count the stars peaks and initialize all peaks to transients 
      (until they are matched with a canonical peak) */
   while(StarPoints[star_points_num].x>0) StarPoints[star_points_num++].bright_star=TRANSIENT;
   /* load the frame at the closest sidereal time */
   baseline_image=new ImageMatrix;
   if (!(baseline_image->LoadFits(baseline_image_name)))
   {  delete baseline_image;
      return(0);  /* failed to open the canonical frame */
   }
   baseline_peaks=new star_point[MAX_BASELINE_PEAKS];
   sorted_by_x=new short[1024+50];
   for (baseline_point_index=0;baseline_point_index<MAX_BASELINE_PEAKS;baseline_point_index++)
     baseline_peaks[baseline_point_index].x=32000;
   altaz2xy(MAX_DIST_FROM_CENTER,90,&x,&y);
   baseline_image->GetStarsPoints(baseline_peaks,NULL,sqrt(pow(x-MID_X,2)+pow(y-MID_Y,2)));
   qsort(baseline_peaks,MAX_BASELINE_PEAKS,sizeof(star_point),sort_by_x);
   /* arrange the points in an array, which its index is the x coordinate value */
   last_x_pos=-1;
   baseline_point_index=1; /* the first is the -1 point */
   while (baseline_point_index<MAX_BASELINE_PEAKS)
   {  if (baseline_peaks[baseline_point_index].x!=last_x_pos)
      {  int a;
         if (baseline_peaks[baseline_point_index].x>=1024) break;/* > 1024 means garbadge information(probably 32000) */
         for (a=last_x_pos+1;a<=baseline_peaks[baseline_point_index].x;a++)
           sorted_by_x[a]=baseline_point_index;
         last_x_pos=baseline_peaks[baseline_point_index].x;
      }
      baseline_point_index++;
   }
   for (int a=last_x_pos+1;a<1024+50;a++) sorted_by_x[a]=baseline_point_index-1;
   /* get the opacity of each star point */
   baseline_point_index=1; /* the first is -1 */
   while (baseline_peaks[baseline_point_index].x<32000 && baseline_point_index<MAX_BASELINE_PEAKS)
   {  min_dist=10000.0;
      star_point_index=0;
      int cosmic_ray_hit=0;
      cosmic_ray_hit=(GetStarsColor(baseline_peaks[baseline_point_index].x,baseline_peaks[baseline_point_index].y,5,baseline_peaks[baseline_point_index].background)==0); 
      while (star_points[star_point_index].x>-1)
      {  /* check if this is the same star point */
	     dist=sqrt(pow(star_points[star_point_index].x-baseline_peaks[baseline_point_index].x,2)+pow(star_points[star_point_index].y-baseline_peaks[baseline_point_index].y,2));
         if (dist<min_dist)
         {  min_dist=dist;
	    matched_point_index=star_point_index;
         }
         star_point_index++;
      }
      if (cosmic_ray_hit)  /* compare CRs in the baseline image with transients in the given image */
      {  if (min_dist<=4) star_points[matched_point_index].bright_star=NOT_BRIGHT; /* if there is a match - don't alert. This is due to possible noise in the CR detection */
         baseline_peaks[baseline_point_index].cloud=-1;   /* don't detect clouds using CR hits */
         goto end_loop1;
      }
 
      if (min_dist>4)
      { /* the star is not visible at all */
        if (baseline_peaks[baseline_point_index].background>0 && baseline_peaks[baseline_point_index].color-baseline_peaks[baseline_point_index].background>baseline_image->AvgColor*3)
	    baseline_peaks[baseline_point_index].cloud=255;
        else baseline_peaks[baseline_point_index].cloud=-1;
      }
      else
      {  double star_sigma,canonical_sigma;
	 star_sigma=sqrt(star_points[matched_point_index].background);
         canonical_sigma=sqrt(baseline_peaks[baseline_point_index].background);      
	 /* check if the star is significantly brighter than in the canonical frame */
         if ((star_points[matched_point_index].color-star_points[matched_point_index].background)/star_sigma>20+2.4*(baseline_peaks[baseline_point_index].color-baseline_peaks[baseline_point_index].background)/canonical_sigma) 
	   star_points[matched_point_index].bright_star=BRIGHT_STAR;
         else 
          star_points[matched_point_index].bright_star=NOT_BRIGHT; /* this star is matched with a canonical peak and is not too bright */

        /* initialize */
	baseline_peaks[baseline_point_index].cloud=0;
        star_points[matched_point_index].cloud=0;
	/*  skyglow_current  = skyglow_canonical  x cloud opacity + cloud glow
            starglow_current = starglow_canonical x cloud opacity + cloud glow
            equations results:
            cloud_opacity = (skyglow_current - starglow_current) / (skyglow_canonical - starglow_canonical)
            cloud_glow = starglow_current - starglow_canonical x cloud_opacity
	*/
        if (baseline_peaks[baseline_point_index].background-baseline_peaks[baseline_point_index].color==0)
          baseline_peaks[baseline_point_index].cloud=-1;  /* prevent division by zero error */
        else 
        {  baseline_peaks[baseline_point_index].star=star_points[matched_point_index].star;
           baseline_peaks[baseline_point_index].cloud_opacity=(double)(star_points[matched_point_index].background-star_points[matched_point_index].color)/(double)(baseline_peaks[baseline_point_index].background-baseline_peaks[baseline_point_index].color);
           baseline_peaks[baseline_point_index].cloud_glow=star_points[matched_point_index].color - (unsigned short)(baseline_peaks[baseline_point_index].color * baseline_peaks[baseline_point_index].cloud_opacity);
	   star_points[matched_point_index].cloud_opacity=baseline_peaks[baseline_point_index].cloud_opacity;
	   star_points[matched_point_index].cloud_glow=baseline_peaks[baseline_point_index].cloud_glow;
           if (star_points[matched_point_index].cloud_opacity<0.99 && star_points[matched_point_index].cloud_glow<AvgColor*9.75 && star_points[matched_point_index].background>0.75*AvgColor)
	   {  /* the parameter of the pow function determines the effect of the visibility of the clouds */
	      star_points[matched_point_index].cloud=(unsigned char)(255*pow((0.99-star_points[matched_point_index].cloud_opacity)/0.99,1.00));        /* the star is covered by a cloud */
              baseline_peaks[baseline_point_index].cloud=star_points[matched_point_index].cloud;
           }
        }
      }
end_loop1:
      baseline_point_index++;
   }
   /* initialize the stars pix array */
   /* this array holds the cloud value (opacity) for each pixel */
   for (y=0;y<height;y++)
     for (x=0;x<width;x++)
       starspix[x][y]=0;

   /* find the opacity of a pixel by the opacities of the closets stars */
   /* draw the scale for each pixel  */
   max_dist=AltModel90.CalcRules(15);
   for (y=0;y<height;y+=5)
     for (x=0;x<width;x+=5)
     {
       if (sqrt(pow(x-MID_X,2)+pow(y-MID_Y,2))<max_dist)
       {  /* find the closest baseline peaks to that pixel */
          long closests[10],close_index;
          double dists[10];
          long closest,second_closest,third_closest;
          double min_dist1=10000,min_dist2=10000,min_dist3=10000;
          long opoq;
          for (close_index=0;close_index<10;close_index++)    /* initialize the closests star's array */
          {  dists[close_index]=100000.0;
             closests[close_index]=-1;
          }
          if (x-60>0) baseline_point_index=sorted_by_x[x-60]; /* find the first star point to start with */
          else baseline_point_index=1;  /* start from the first point (after the -1 point) */
	     while(baseline_point_index<=sorted_by_x[x+60])
          {  /* find the closest visible star */
             dist=sqrt(pow(x-baseline_peaks[baseline_point_index].x,2)+pow(y-baseline_peaks[baseline_point_index].y,2));
             if (dist<=60.0)  /* save time by not searching in every distance */
               if (baseline_peaks[baseline_point_index].cloud!=-1)
		 if (pow(baseline_peaks[baseline_point_index].x-MID_X,2)+pow(baseline_peaks[baseline_point_index].y-MID_Y,2)<pow(max_dist+25,2))
               for (close_index=0;close_index<10;close_index++)
               if (dist<dists[close_index])
               {  memmove(&(dists[close_index+1]),&(dists[close_index]),sizeof(double)*(9-close_index));
                  memmove(&(closests[close_index+1]),&(closests[close_index]),sizeof(long)*(9-close_index));
                  dists[close_index]=dist;
                  closests[close_index]=baseline_point_index;
                  break;
               }
             baseline_point_index++;
         }

         /* interpolation (first parameter can be 1 or 2 to switch neighbor and linear interpolation) */
         opoq=interpolate(1,x,y,baseline_peaks,dists,closests);
         if (opoq>255) opoq=255;     /* make sure the value does not excceed 255 */
         if (opoq<20 && opoq>3) opoq=opoq+20;
         /* set the opacity for all the pixels within the 5X5 cluster */
         for (long b=y;b<y+5;b++)
         for (long a=x;a<x+5;a++) starspix[a][b]=opoq;
       }
    }
    /* draw the pixels */
    for (y=(long)(MID_Y-max_dist);y<MID_Y+max_dist;y++)
      for(x=(long)(MID_X-max_dist);x<MID_X+max_dist;x++)
    if (sqrt(pow(x-MID_X,2)+pow(y-MID_Y,2))<max_dist)
    {  unsigned char opoq;
       unsigned char cloud_style[10];
       long xc,yc,sum=0,num=0;
       for (yc=y-20;yc<y+20;yc++)
         for(xc=x-20;xc<x+20;xc++)
         {  sum+=starspix[xc][yc];
            num++;
         }
         if (num>0) opoq=(unsigned char)(sum/num);
         else opoq=0;
         if (opoq>0)
         {  cd_setstyle(cloud_style,0,0,255,opoq,opoq,opoq,1,0);
       	       cd_plot(bitmap,x,y,cloud_style);
         }
    }
   
    /* ************************************************* */   
    /* printout points brighter than the canonical frame */
    /* ************************************************* */

    if (find_changes)
    {  char image_name[255];
       long star_color,star_background;
       max_dist=AltModel90.CalcRules(25); /* due to ground-based objects around CONCAM, search for transients only above 25 degrees */
       star_point_index=0;
       /* get the name of the image */
       strcpy(image_name,name);
       if (strrchr(name,'/')) strcpy(image_name,&((strrchr(name,'/'))[1]));  /* get only the file name */
       while (star_points[star_point_index].x>0)
       {  int planet_index;
          /* first check if theis star is too close to a planet */
          if (star_points[star_point_index].bright_star!=NOT_BRIGHT)
          {  /* search the bright planets to check if the point is near one of them */
	     swe_set_topo(location->latitude,location->longitude,location->height);
             for (planet_index=1;planet_index<=PLANETS_NUM-4;planet_index++)  /* calculate all planets but the sun and the moon */
             {  double xx[6],ra,dec,zd,az,rar,decr,dist,min_dist,max_dist;
                int res;
                char serr[200];
                res=swe_calc_ut(JulianDate,planet_index, SEFLG_TOPOCTR | SEFLG_SWIEPH | SEFLG_SPEED | SEFLG_EQUATORIAL ,xx, serr);
                equ2hor(JulianDate,DELTA_T,0,0,location,xx[0]/15,xx[1],0,&zd,&az,&rar,&decr);
                res=altaz2xy(90-zd,az,&x,&y);  /* calculate the x,y coordinates */
                if (res==1)  /* only if the star is not beyond the horizon */  	       
                {  dist=sqrt(pow(star_points[star_point_index].x-x,2)+pow(star_points[star_point_index].y-y,2));
                   if (dist<15) goto end_loop;  /* the star is too close to a bright planet so it is probably not a "real" star */
                }
             }
             /* write the information to the std output stream */
             x=star_points[star_point_index].x;
             y=star_points[star_point_index].y;
             star_color=star_points[star_point_index].color;
             star_background=star_points[star_point_index].background;
             if ((star_color-star_background)>30*sqrt((double)(star_background)) /* alert only if the transient is brighter than 30 sigmas */
	     && (data[x][y]-star_background>30*sqrt(star_background))
	     && (pow(x-MID_X,2)+pow(y-MID_Y,2)<max_dist*max_dist)
/* && (pow(x-MID_X,2)+pow(y-MID_Y,2)>150*150) */
	     && (star_color-baseline_image->data[x][y]>30*sqrt(baseline_image->data[x][y])
	     && (data[x][y]-star_background)>1.4*(baseline_image->data[x][y]-star_background))
             && (data[x][y]-baseline_image->data[x][y]>30*sqrt(baseline_image->data[x][y])))
	     {  unsigned short bs_color,bs_background;
	        char dir[128],jpg_name[64];
	        bs_background=baseline_image->GetStarsBackground(x,y);
	        bs_color=baseline_image->GetStarsColor(x,y,5,0.5*bs_background);
                if ((star_color-star_background>30*sqrt(bs_background)+(bs_color-bs_background))
                && (data[x][y]-star_background>30*sqrt(bs_background)+(bs_color-bs_background)))   /* also should be 30 sigmas over the canonical color-background */
                {
                  double alt,az,ra,dec;
                  long C1_B,C5_B,C9_B,C16_B,C25_B;
                  xy2altaz(x,y,&alt,&az);
 	          altaz2radec(&current_location,alt,az,&ra,&dec);
                  C1_B=star_points[star_point_index].C1-star_points[star_point_index].background;
                  C5_B=star_points[star_point_index].C5-star_points[star_point_index].background;
                  C9_B=star_points[star_point_index].C9-star_points[star_point_index].background;
                  C16_B=star_points[star_point_index].C16-star_points[star_point_index].background;
                  C25_B=star_points[star_point_index].C25-star_points[star_point_index].background;
                  if (star_points[star_point_index].bright_star==TRANSIENT) printf("transient:<transient>");
                  else printf("luminance change:<luminance_change>");
                  strcpy(dir,"http://concam.net/");
                  strncat(dir,image_name,2);
                  strcat(dir,"/");
                  strncat(dir,image_name,8);
                  strcpy(jpg_name,image_name);
                  *strchr(jpg_name,'.')='\0';
                  strcat(jpg_name,".jpg");
//                  printf(" <a href=\"%s\">%s</a> jd=<jd>%.6f</jd> ra=~<ra>%2.2f</ra> dec=~<dec>%2.2f</dec> alt=<alt>%2.2f</alt> az=<az>%3.2f</az> X=<x>%d</x> Y=<y>%d</y> C1-b=<C1-B>%d</C1-B> C5-B=<C5-B>%d</C5-B> C9-B=<C9-B>%d</C9-B> C16-B=<C16-B>%d</C16-B> C25-B=<C25-B>%d</C25-B><canonical=%s><br>",image_name,image_name,current_date,ra,dec,alt,az,x+1,height-y,C1_B,C5_B,C9_B,C16_B,C25_B,baseline_image_name);
	          printf(" <a href=\"%s/%s\">%s</a> <a href=\"%s/%s\">jpg</a> jd=<jd>%.6f</jd> ra=~<ra>%2.2f</ra> dec=~<dec>%2.2f</dec> alt=<alt>%2.2f</alt> az=<az>%3.2f</az> X=<x>%d</x> Y=<y>%d</y> Jx=%d Jy=%d<canonical=%s><br>",dir,image_name,image_name,dir,jpg_name,current_date,ra,dec,alt,az,x+1,height-y,(long)((x+1)/2),(long)(y/2),baseline_image_name);
                  if (star_points[star_point_index].bright_star==TRANSIENT) printf("</transient>\n");
                  else printf("</luminance_change>\n");
		}
	     }
	  }
end_loop:
          star_point_index++;
       }
     }
     /* average the two frames and save for farther useage */
     if (IsClear(star_points,current_location.clear_threshold)) /* only if the image is a clear image */
     {  baseline_image->AverageMatrix(this,0.5);
        baseline_image->SaveFits(baseline_image_name);
     }
     baseline_image->FreeImage();
     delete sorted_by_x;
     delete baseline_image;
     delete baseline_peaks;
     baseline_peaks=NULL;
     /* draw the scale */
     if (!(font_file=fopen("fonts/courier18.bfnt","r"))) printf("could not open file 'fonts/courier18.bfnt'\n");
     font=cd_loadfont(font_file);
     cd_font_height[font]=15;
     fclose(font_file);
     cd_setstyle(font_buffer,200,200,200,255,255,255,1,0);
     cd_drawtext(bitmap,font,36,70,font_buffer,"dM",NO,NO,NO);
     y=100;
     for (scale_item=0;scale_item<1;scale_item+=0.2)
     { char buffer[20];
       unsigned short opoq;
       opoq=(unsigned short)(255*pow((0.99-scale_item)/0.99,1.00));
       if (opoq>255) opoq=255;
       cd_setstyle(cloud_style,0,0,255,opoq,opoq,opoq,1,0);
       dist=-1*(2.5*log10(scale_item));
       sprintf(buffer,"%.2f",dist);
       cd_drawtext(bitmap,font,20,y,font_buffer,buffer,NO,NO,NO);
       cd_rect(bitmap,80,y,110,y+cd_font_height[font],1,font_buffer,font_buffer);
       cd_rect(bitmap,80,y,110,y+cd_font_height[font],1,cloud_style,cloud_style);
       y+=cd_font_height[font]+20;
     }
     return(1);
}


void ImageMatrix::FindClouds2(BITMAP bitmap,star_point *star_point,star_point *baseline_peaks)
{ /* get clouds areas */
  unsigned char cloud_style[10];
  long baseline_point_index,star_point_index,matched_point_index,x,y;
  double dist,scale_item,max_dist;
  FONT font;
  unsigned char font_buffer[10];
  FILE *font_file;
  /* get all the less visible star points */
  baseline_point_index=0;
  while (baseline_peaks[baseline_point_index].x>0)
  {  double min_dist;
     star_point_index=0;
     min_dist=10000.0;
     while (star_points[star_point_index].x>0)
     {  /* check if this is the same star point */
        dist=sqrt(pow(star_points[star_point_index].x-baseline_peaks[baseline_point_index].x,2)+pow(star_points[star_point_index].y-baseline_peaks[baseline_point_index].y,2));
	if (dist<min_dist)
	{  min_dist=dist;
	   matched_point_index=star_point_index;
	}
	star_point_index++;
     }
     if (min_dist>2) 
     {
       if (min_dist>25) baseline_peaks[baseline_point_index].cloud=255;  /* no stars around at all */
     }
     else
     {  
	/* skyglow_current  = skyglow_canonical  x cloud opacity + cloud glow 
	   starglow_current = starglow_canonical x cloud opacity + cloud glow
           equations results:
           cloud_opacity = (skyglow_current - starglow_current) / (skyglow_canonical - starglow_canonical)
           cloud_glow = starglow_current - starglow_canonical x cloud_opacity
        */
       baseline_peaks[baseline_point_index].star=star_points[matched_point_index].star;
       baseline_peaks[baseline_point_index].cloud_opacity=(double)(star_points[matched_point_index].background-star_points[matched_point_index].color)/(double)(baseline_peaks[baseline_point_index].background-baseline_peaks[baseline_point_index].color);
       baseline_peaks[baseline_point_index].cloud_glow=star_points[matched_point_index].color - (unsigned short)(baseline_peaks[baseline_point_index].color * baseline_peaks[baseline_point_index].cloud_opacity);
       star_points[matched_point_index].cloud_opacity=baseline_peaks[baseline_point_index].cloud_opacity;
       star_points[matched_point_index].cloud_glow=baseline_peaks[baseline_point_index].cloud_glow;
	 //        if ((star_points[matched_point_index].color-star_points[matched_point_index].background)*1.1<baseline_peaks[baseline_point_index].color-baseline_peaks[baseline_point_index].background)
	if (star_points[matched_point_index].cloud_opacity<0.99 && star_points[matched_point_index].cloud_glow<AvgColor*9.75 && star_points[matched_point_index].background>0.75*AvgColor)
	  {  /* the parameter of the pow function determines the effect of the visibility of the clouds */
           star_point[matched_point_index].cloud=(unsigned char)(255*pow((0.99-star_points[matched_point_index].cloud_opacity)/0.99,1.00));        /* the star is covered by a cloud */
	   baseline_peaks[baseline_point_index].cloud=star_points[matched_point_index].cloud;
        }
     }
     baseline_point_index++;
   }
  
  /* initialize the stars pix array */
  for (y=0;y<height;y++)
    for (x=0;x<width;x++) starspix[x][y]=0;   
  /* draw the scale for each pixel;
  max_dist=350; /* !!!!! */
  for (y=0;y<height;y++)
    for (x=0;x<width;x++)
      if (y>230 && (x<460 || y<690-(x-460))) /* !!! */
//     if (data[x][y]>AvgColor*0.75)
      {
	if (sqrt(pow(x-MID_X,2)+pow(y-MID_Y,2))<350) /* !!! */
        {   /* find the closest baseline peaks to that pixel */
	    long closests[10],close_index;
	    double dists[10],close_sum,close_num;
            long closest,second_closest,third_closest;
            double min_dist1=10000,min_dist2=10000,min_dist3=10000;
            long opoq;
	    //    unsigned char cloud_style[10];
	    for (close_index=0;close_index<10;close_index++)   /* initialize the closests star's array */
	      (dists[close_index]=100000.0);
	    baseline_point_index=0;
            while(baseline_peaks[baseline_point_index].x>0)
	      {  /* find the closest visible star
//              if (((long)(baseline_peaks[baseline_point_index].star)<1000) || (baseline_peaks[baseline_point_index].star->type[0]!='A' && baseline_peaks[baseline_point_index].star->type[0]!='B' && baseline_peaks[baseline_point_index].star->type[0]!='O')) /* ignore the blue stars */
     //	      if (baseline_peaks[baseline_point_index].cloud<255)
               {  dist=sqrt(pow(x-baseline_peaks[baseline_point_index].x,2)+pow(y-baseline_peaks[baseline_point_index].y,2));
               if (dist<50.0)  /* save time by not searching in every distance */
	       for (close_index=0;close_index<10;close_index++)
                 if (dist<dists[close_index]) 
		 {  memcpy(&(dists[close_index+1]),&(dists[close_index]),sizeof(double)*(9-close_index));
		    memcpy(&(closests[close_index+1]),&(closests[close_index]),sizeof(long)*(9-close_index));
		    dists[close_index]=dist;
		    closests[close_index]=baseline_point_index;
		    break;
		 }
	       }
	       baseline_point_index++;
            }
            opoq=0;
            close_sum=0.0;
            close_num=0.0;
            /* count the non detected stars */
	    for (close_index=0;close_index<10;close_index++)
              if(baseline_peaks[closests[close_index]].cloud==255) close_num=close_num+1;
	        if (close_num<4)
	          for (close_index=0;close_index<10;close_index++)
	     	    if(baseline_peaks[closests[close_index]].cloud==255) baseline_peaks[closests[close_index]].cloud=128;
            close_num=0.0;
            for (close_index=0;close_index<5;close_index++)
              if (dists[close_index]<50)
	      {  close_sum+=baseline_peaks[closests[close_index]].cloud*(dists[close_index]/10);
	      //        printf("%d: %d #",close_index,baseline_peaks[closests[close_index]].cloud);
	         close_num=(close_num+dists[close_index]/10);
              }
            if (close_num==0) opoq=255; /* no stars where detected at this area at all */
            else opoq=(unsigned short)(close_sum/close_num);
            if (opoq>255) opoq=255;  /* make sure the value does not excceed 255 */
            if (opoq<20 && opoq>3) opoq=opoq+20;
            starspix[x][y]=opoq;
         }
      }
   /* draw the pixels */
   for (y=30;y<height-30;y++)
     for(x=30;x<width-30;x++)
  //     if (starspix[x][y]>0)
     {  unsigned char opoq;
        unsigned char cloud_style[10];
        long xc,yc,sum=0,num=0;
        for (yc=y-20;yc<y+20;yc++)
          for(xc=x-20;xc<x+20;xc++)
    //      if (starspix[xc][yc]>0) 
	    {  sum+=starspix[xc][yc];
	       num++;
	    }
        if (num>0) opoq=(unsigned char)(sum/num);
        else opoq=0;
        if (opoq>0) 
        {  cd_setstyle(cloud_style,0,0,255,opoq,opoq,opoq,1,0);
	   cd_plot(bitmap,x,y,cloud_style);
        }
     }
   /* draw the scale */
   if (!(font_file=fopen("fonts/courier18.bfnt","r"))) printf("could not open file 'fonts/courier18.bfnt'\n");
   font=cd_loadfont(font_file);
   cd_font_height[font]=15;
   fclose(font_file);
   cd_setstyle(font_buffer,200,200,200,255,255,255,1,0);
   cd_drawtext(bitmap,font,36,70,font_buffer,"dM",NO,NO,NO);
   y=100;
   for (scale_item=0;scale_item<1;scale_item+=0.2)
   { char buffer[20];
     unsigned short opoq;
     opoq=(unsigned short)(255*pow((0.99-scale_item)/0.99,1.00));
     if (opoq>255) opoq=255;
     cd_setstyle(cloud_style,0,0,255,opoq,opoq,opoq,1,0);
     dist=-1*(2.5*log10(scale_item));
     sprintf(buffer,"%.2f",dist);
     cd_drawtext(bitmap,font,20,y,font_buffer,buffer,NO,NO,NO);
     cd_rect(bitmap,80,y,110,y+cd_font_height[font],1,font_buffer,font_buffer);
     cd_rect(bitmap,80,y,110,y+cd_font_height[font],1,cloud_style,cloud_style);
     y+=cd_font_height[font]+20;
   } 
return;
   /* draw the clouds areas */
   /* first clear the starpix data */
#define max_cloud_dist 100
   for (y=0;y<height;y++)
     for(x=0;x<width;x++) starspix[x][y]=0;
   baseline_point_index=0;
   while (baseline_peaks[baseline_point_index].x>0)
   {  unsigned short opoq;
      if (baseline_peaks[baseline_point_index].cloud>0)
      {  star_point_index=baseline_point_index+1;
	 while (baseline_peaks[star_point_index].x>0)
	 { dist=sqrt(pow(baseline_peaks[baseline_point_index].x-baseline_peaks[star_point_index].x,2)+pow(baseline_peaks[baseline_point_index].y-baseline_peaks[star_point_index].y,2));
	   if (dist<max_cloud_dist)
	   {  /* paint the cloud */
	      unsigned char area_cloud;
              long x1,x2,y1,y2;
	      area_cloud=(unsigned char)((baseline_peaks[baseline_point_index].cloud+baseline_peaks[baseline_point_index].cloud)/2); /* the area cloud value is the avarage of the two cloud point */
               /* set the rectangle between the two stars */
              if (baseline_peaks[baseline_point_index].x<baseline_peaks[star_point_index].x)
	      { x1=baseline_peaks[baseline_point_index].x;
	        x2=baseline_peaks[star_point_index].x;
	      }
              else
              {  x1=baseline_peaks[star_point_index].x;
                 x2=baseline_peaks[baseline_point_index].x;
              }
              if (baseline_peaks[baseline_point_index].y<baseline_peaks[star_point_index].y)
              {  y1=baseline_peaks[baseline_point_index].y;
      	         y2=baseline_peaks[star_point_index].y;
              }
              else
	      {  y1=baseline_peaks[star_point_index].y;
	         y2=baseline_peaks[baseline_point_index].y;
	      }
	      /* paint the rectangle */
	      for (x=x1-(long)(0.3*(x2-x1));x<x2+(long)(0.3*(x2-x1));x++)
		for (y=y1-(long)(0.3*(y2-y1));y<y2+(long)(0.3*(y2-y1));y++)
		{ double radius,dist_from_center;  /* for drwaing circles instead of rectangles */
                  unsigned short cent_x,cent_y;
	          radius=(sqrt(pow(x2-x1,2)+pow(y2-y1,2))/2)*1.2;
	          cent_x=(unsigned short)((float)(x2+x1)/2); /* coordinate of the center of the square */
	          cent_y=(unsigned short)((float)(y2+y1)/2); /* coordinate of the center of the square */
                  dist_from_center=sqrt(pow(x-cent_x,2)+pow(y-cent_y,2));
	          if (dist_from_center<radius && x>0 && x<width-1 && y>0 && y<height-1)
    //            { if (starspix[x][y]==0)
		    {  opoq=(unsigned char)((area_cloud*1+((data[x][y]-AvgColor)/AvgColor)*255*0)/(1+0));
                       if (data[x][y]<1*AvgColor) opoq=0;   /* no cloud at this point */
                       if (opoq>255) opoq=255;
	               if (starspix[x][y]>0 && (float)opoq/(float)starspix[x][y]>1.3) /* smoth the edges */
		       {  opoq=(unsigned short)((1.0*opoq*(1-dist_from_center/radius)+starspix[x][y])/(1+1.0*(1-dist_from_center/radius))-starspix[x][y]);
	               }
	               else if (starspix[x][y]>0) opoq=0;
		       //   if (starspix[x][y]==0 && opoq>starspix[x][y])
		       //  {  opoq=opoq-(unsigned short)(opoq*(x<x1)*(((x1-x)/(0.3*(x2-x1))))-opoq*(x>x2)*(((x-x2)/(0.3*(x2-x1))))-opoq*(y<y1)*(((y1-y)/(0.3*(y2-y1))))-opoq*(y>y2)*(((y-y2)/(0.3*(y2-y1)))));
		       //   }
		       // else if (starspix[x][y]>0) opoq=0;
	               if (opoq>0)
	               {
	   	         cd_setstyle(cloud_style,0,0,255,opoq,opoq,opoq,1,0);
	                 cd_plot(bitmap,x,y,cloud_style);
	                 starspix[x][y]+=opoq;  /* mark this location so it won't be colored again */
	               }
		    }
		  
                } 
	   }
	   star_point_index++;
       }   
    }
    baseline_point_index++;
   }
   return;
}

/* ************************************************************************** */
/* **************************** Average Matrix ******************************
purpose: makes an average of two fits frames
input: matrix - a pointer to an ImageMatrix structures to compare with.
       weight - the weight of the new matrix in the total average.
output: 1 if succeeded, 0 if failed
*/
long ImageMatrix::AverageMatrix(ImageMatrix *matrix, double weight)
{  long x,y;
   if (width!=matrix->width) return(0);
   if (height!=matrix->height) return(0);
   for (x=0;x<width;x++)
     for (y=0;y<height;y++)
       if (matrix->data[x][y]!=0)
//         if ((double)data[x][y]/(double)matrix->data[x][y]>0.4) /* make sure that bad pixels or clouds will not be included */
           data[x][y]=(COLOR_TYPE)(((double)(data[x][y]+matrix->data[x][y]*weight))/(1.0+weight));
   return(1); /* function succeeded */
}

/* ************************************************************************** */
/* **************************** IsClear ******************************
purpose: check if a frame is clear and not cloudy
input: star_points - the list of PSFs.
       clear_threshold - a number (0-1) indicates the minimum number of detected stars
output: 1 if clear, 0 if cloudy or wet (or any other reason)
comments: the functions determines if the frame is clear by the number of
visible stars and it also checks that stars are visible in all parts of the image
paper reference: "All-sky Relative Opacity Mapping Using Night Time Panoramic Images", PASP
*/
long ImageMatrix::IsClear(star_point *star_points, double clear_threshold)
{  int star_index=0,checked=0,not_checked=0,double_star_index; 
   while (stars[star_index].starname[0])
   {  if (stars[star_index].altitude>20 && stars[star_index].x>0 
      && !(stars[star_index].x>current_location.bad.x1 && stars[star_index].x<current_location.bad.x2
           && stars[star_index].y>current_location.bad.y1 && stars[star_index].y<current_location.bad.y2)
      && stars[star_index].magnitude<4.5) /* check all the stars brighter than 4.5 */
      {  if (stars[star_index].checked) checked++;
         else 
	 {  for (double_star_index=0;double_star_index<star_index;double_star_index++)
	      if (stars[star_index].x==stars[double_star_index].x && stars[star_index].y==stars[double_star_index].y) 
		{  not_checked--;   /* don't count a double star if one has already been counted */
		   break;
                }
            not_checked++;
         }
//printf("%s %d %d\n",stars[star_index].starname,stars[star_index].checked,star_index);
      }
      star_index++;
   }
//printf("%d %d %d\n",star_index,checked,not_checked);
   if (checked+not_checked==0) return(0);  /* no stars were found. also prevent divison by zero error */
   if ((double)checked/(double)(checked+not_checked)>clear_threshold) return(1);
   else return(0);

  //   short stars_spread[36][2];
  //   long distance,angle;
  //   long star_point_index=0;
  //   if (exposure_time<160 || exposure_time>200) return(0);  /* verify that the exposure time is right */
  //   /* initialize the stars spread array */
  //   for (distance=0;distance<2;distance++)
  //     for (angle=0;angle<36;angle++)
  //       stars_spread[angle][distance]=0;
  //   while (star_points[star_point_index].x>-1)
  //   {  if (star_points[star_point_index].x!=MID_X && star_points[star_point_index].y!=MID_Y)
  //      {  if (sqrt(pow(star_points[star_point_index].x-MID_X,2)+pow(star_points[star_point_index].y-MID_Y,2))<0.5*(0.5*width)) distance=0; /* short distance */
  //         else distance=1; /* long distance from the center of the image*/
  //         angle=180+(long)(RAD2DEG*atan2(star_points[star_point_index].x-MID_X,star_points[star_point_index].y-MID_Y));
  //         angle=angle/10;  /* set angle to a value between 0-35 */
  //         stars_spread[angle][distance]=stars_spread[angle][distance]+1;  /* increment the stars found there */
  //      }
  //      star_point_index++;
  //   }
  //   if (star_point_index<600) return(0);   /* check if there are enough stars visible */
  //   /* check that the frame is covered with stars */
  //   for (distance=0;distance<2;distance++)
  //     for (angle=0;angle<36;angle++)
  //       if (stars_spread[angle][distance]<4) return(0);   /* an area with very few stars was found */
  //   return(1); /* function succeeded */
}

 
/* ************************************************************************** */
/* *************************** ClearBackground ******************************
purpose: remove all the stars from the image and save it.
parameters: star_points - a list of all the stars detected in the image.
            input_image - the name of the fits image file.
comments: this function generates a fits file
          It also changes the data in the "data" array
*/
void ImageMatrix::ClearBackground(char *input_image,BITMAP bitmap,double max_dist)
{  unsigned short x,y,x1,x2,y1,y2,bg_color,xc,yc;
   float left_color,right_color,up_color,down_color;
   long sum,num,darkest_pixel_value=100000;
   char background_image_name[60],*p_background_image_name,*p;
   /* set the filename of the background image */
   if (input_image)              /* save the background image as a file */
   {  strcpy(background_image_name,input_image);
      p=strrchr(background_image_name,'.');
      if (!p) p=&(background_image_name[strlen(background_image_name)]); /* no suffix in the image name */
      if (p[-1]=='p') p=p-1; /* replace the letter "p" with the letter "b" (only if the letter p appears in the name) */
      strcpy(p,"b.fits");    /* now background_image_name contains the output file name of the background image */
      p_background_image_name=background_image_name;
   }
   else p_background_image_name=NULL; /* don't save the background as file */
   /* detect the stars and cover them with the background's value */
   for(y=5;y<height-5;y++)
     for(x=5;x<width-5;x++)
       {  if (data[x][y]<darkest_pixel_value && sqrt(pow(x-MID_X,2)+pow(y-MID_Y,2))<max_dist) darkest_pixel_value=data[x][y];
        if (data[x][y]>data[x-1][y]+data[x][y]>data[x+1][y]+data[x][y]>data[x][y-1]+data[x][y]>data[x][y+1]+data[x][y]>data[x-1][y+1]+data[x][y]>data[x-1][y-1]+data[x][y]>data[x+1][y-1]+data[x][y]>data[x+1][y+1]<=1)
        {  if (data[x][y]<data[x-2][y] || data[x][y]<data[x+2][y] || data[x][y]<data[x][y-2] || data[x][y]<data[x][y+2]) continue;
           if (data[x][y]<AvgColor*0.8) continue;
           x1=x2=x;
           y1=y2=y;
           /* find the broders of the star (the area of the star) */
           while(x1>3 && data[x1-2][y]>0.6*(float)AvgColor && data[x1-3][y]>0.4*(float)AvgColor && ((float)data[x1][y]/(float)data[x1-2][y]>1.02 || (float)data[x1][y]/(float)data[x1-3][y]>1.04 || data[x1][y]>10*AvgColor)) x1--;
           while(x2<width-3 && data[x2+2][y]>0.6*(float)AvgColor && data[x2+3][y]>0.4*(float)AvgColor && ((float)data[x2][y]/(float)data[x2+2][y]>1.02 || (float)data[x2][y]/(float)data[x2+3][y]>1.04 || data[x2][y]>10*AvgColor)) x2++;
           while(y1>3 && data[x][y1-2]>0.6*(float)AvgColor && data[x][y1-3]>0.4*(float)AvgColor && ((float)data[x][y1]/(float)data[x][y1-2]>1.02 || (float)data[x][y1]/(float)data[x][y1-3]>1.04 || data[x][y1]>10*AvgColor)) y1--;
           while(y2<height-5 && data[x][y2+2]>0.6*(float)AvgColor && data[x][y2+3]>0.4*(float)AvgColor && ((float)data[x][y2]/(float)data[x][y2+2]>1.02 || (float)data[x][y2]/(float)data[x][y2+3]>1.04 || data[x][y2]>10*AvgColor)) y2++;
           /* now x1,y1,x2,y2 is the smallest rectangle that covers the star */
           /* check if it is only one dot */
           if (x1==x2 && y1==y2)
           {  data[x][y]=(unsigned short)((float)(data[x-1][y]+data[x+1][y]+data[x][y-1]+data[x][y+1])/4);
              continue;
           }
           /* check if it is a line */
           if (x1==x2)
           {  for(yc=y1;yc<y2;yc++) data[x][yc]=(unsigned short)((float)(data[x][y-1]+data[x][y+1])/2);
              continue;
           }
           if (y1==y2)
           {  for(xc=x1;xc<x2;xc++) data[xc][y]=(unsigned short)((float)(data[xc][y]+data[xc][y])/2);
              continue;
           }
           /* get the colors of the sides of the rectangle */
           sum=num=0;
           for(yc=y1;yc<=y2;yc++) if (data[x1][yc]>AvgColor/5)
           {  sum+=data[x1][yc];
              num++;
           }
           if (num) left_color=sum/num;
           else left_color=0;
           sum=num=0;
           for(yc=y1;yc<=y2;yc++) if (data[x2][yc]>AvgColor/5)
           {  sum+=data[x2][yc];
              num++;
           }
           if (num) right_color=sum/num;
           else right_color=0;
           sum=num=0;
           for(xc=x1;xc<=x2;xc++) if (data[xc][y1]>AvgColor/5)
           {  sum+=data[xc][y1];
              num++;
           }
           if (num) up_color=sum/(x2-x1+1);
           else up_color=0;
           sum=num=0;
           for(xc=x1;xc<=x2;xc++) if (data[xc][y2]>AvgColor/5)
           {  sum+=data[xc][y2];
              num++;
           }
           if (num) down_color=sum/num;
           else down_color=0;
           /* make sure no color is black */
           if (left_color<AvgColor/3) left_color=right_color;
           if (right_color<AvgColor/3) right_color=left_color;
           if (up_color<AvgColor/3) up_color=down_color;
           if (down_color<AvgColor/3) down_color=up_color;
           /* color the rectangle */
           for (yc=y1;yc<=y2;yc++)
           for (xc=x1;xc<=x2;xc++)
	   {  float x_color,y_color;
	      unsigned short pix_color;
	      x_color=(left_color*(x2-xc)+right_color*(xc-x1))/(x2-xc+xc-x1);
	      y_color=(up_color*(y2-yc)+down_color*(yc-y1))/(y2-yc+yc-y1);
	      pix_color=(unsigned short)((x_color*(x2-x1)+y_color*(y2-y1))/(x2-x1+y2-y1));
	      if (pix_color<data[xc][yc]) data[xc][yc]=pix_color;
	   }
       }
     }
   if (p_background_image_name) SaveFits(p_background_image_name);
   if (bitmap>=0) GetPicture(bitmap);
}

/* ************************************************************************** */
/* ***************************** DarkestPixel *******************************
purpose: find the lowest pixel's value.
parameters: dist_from_center - maximum distance from center to search in.
comments: the returned value is the midean of the 10th lowest pixels counts
*/
unsigned short ImageMatrix::DarkestPixel(double dist_from_center)
{  long x,y,count_index;
   unsigned short min_counts[10];
   /* initalize the counts array */
   for (count_index=0;count_index<10;count_index++) min_counts[count_index]=0xFFFF;
     for (y=(long)(MID_Y-dist_from_center);y<MID_Y+(long)dist_from_center;y++)
       for (x=(long)(MID_X-dist_from_center);x<MID_X+(long)dist_from_center;x++)
         if (y>0 && y<height && x>0 && x<width)
           if (sqrt(pow(x-MID_X,2)+pow(y-MID_Y,2))<dist_from_center)
	     if (data[x][y]>0.1*AvgColor) /* make sure it is not an empty pixel */
             {  for (count_index=0;count_index<10;count_index++)
	        if (data[x][y]<min_counts[count_index])
	        {  memcpy(&(min_counts[count_index+1]),&(min_counts[count_index]),sizeof(unsigned short)*(9-count_index));
	           min_counts[count_index]=data[x][y];
	           break;
	        }
             }
    return(min_counts[4]);   /* return the midean value */
}

/* ************************************************************************* */
/* *************************** CompareBackground ******************************
purpose: remove all the stars from the image and save it.
parameters: input_image - the name of the fits image file.
                          if this parameter is null then the image is compared to
                          its lowest pixel's value.
            bitmap - the bitmap to draw the output on.
comments: this function generates a fits file
          before calling this function. It is better that the function
          "ClearBackground" will be called.
*/
void ImageMatrix::CompareBackground(BITMAP bitmap,char *baseline_image_name,double max_dist)
{  ImageMatrix *baseline_image=NULL;
   long x,y;
   double image_jd,scale,scale_start,scale_step;
   char *p;
   unsigned short darkest_pixel,y_top;
   /* get the name of the right baseline image */
 
   if (baseline_image_name)
   {  /* set the filename of the compared image */
      baseline_image=new ImageMatrix;
      baseline_image->LoadFits(baseline_image_name);
      baseline_image->ClearBackground(NULL,-1,max_dist);
//    strcpy(compared_image_name,input_image);
//    p=strrchr(compared_image_name,'.');
//    if (!p) p=&(compared_image_name[strlen(compared_image_name)]); /* no suffix in the image name */
//   if (p[-1]=='p') p=p-1; /* replace the letter "p" with the letter "b" (only if the letter p appears in the name) */
//   strcpy(p,"c.fits");    /* now compared_image_name contains the output file name of the compared image */
//     baseline_image=new ImageMatrix;
//     if (!baseline_image->LoadFits(baseline_image_name))
//  {  printf("Can not open file '%s'\n",baseline_image_name);
//    return;
//  }
   }
   darkest_pixel=DarkestPixel(200.0);
   for (y=0;y<height-1;y++)
     for (x=0;x<width-1;x++)
       if (data[x][y]>darkest_pixel)
       {  if (baseline_image && data[x][y]<=baseline_image->data[x][y]) data[x][y]=baseline_image->data[x][y]+1;
	  if (baseline_image && baseline_image->data[x][y]==0) data[x][y]=0;
          else
	  {  double delta;
	     if (baseline_image) delta=1000*2.5*log10((double)data[x][y]/(double)baseline_image->data[x][y]);
             else delta=1000*2.5*log10((double)data[x][y]/(double)darkest_pixel);
	     if (delta<10) delta=0;
	     if (delta>0xFFFF) delta=0;     /* this pixel is noise */
             if (baseline_image) delta*=2;  /* when comparing to another frame */
	     data[x][y]=(unsigned short)delta;
	 }
       }
       else data[x][y]=0;
   /* darw the scale colored rectangles */
   y_top=100;
   if (baseline_image) 
   {  scale_start=2.75;
      scale_step=0.25;
   }
   else 
   {  scale_start=7.0;
      scale_step=0.75;
   }
   for (scale=scale_start;scale>=1.0;scale-=scale_step)
   {  for (y=0;y<30;y++)
        for(x=50;x<100;x++)
	  data[x][y_top+y]=((baseline_image!=NULL)+1)*(long)(1000*2.5*log10(scale));//(1000*(scale));
      y_top+=50;
   }
   GetPicture(bitmap);  /* get the fits data to a jpeg image */
   /* draw the labels of the scale */
   { FILE *font_file;
     FONT font;
     unsigned char font_buffer[10];
     char buffer[20];
     if (!(font_file=fopen("fonts/courier18.bfnt","r"))) printf("could not open file 'fonts/courier18.bfnt'\n");
     font=cd_loadfont(font_file);
     cd_font_height[font]=16;
     fclose(font_file);
     cd_setstyle(font_buffer,200,200,200,255,255,255,1,0);
     cd_drawtext(bitmap,font,20,70,font_buffer,"dM",NO,NO,NO);
     y_top=110;
     for (scale=scale_start;scale>=1.0;scale-=scale_step)
     { if (baseline_image) sprintf(buffer,"%.2f",2.5*log10(scale));
       else sprintf(buffer,"%.1f",21-2.5*log10(scale));
       cd_drawtext(bitmap,font,6,y_top,font_buffer,buffer,NO,NO,NO);
       y_top+=50;
     }
   }     
   //   SaveFits("temp.fits");
   if (baseline_image) baseline_image->FreeImage();
   if (baseline_image) delete baseline_image;
}


/* **************** MarkUnknownStars *************** 
   mark the non-catalogued objects.
   input: star_points - a list of star points of the frame.
   returned value: the number of unkown expected objects.
   comments: this function should be called after FindUnkownStars or FIndUnkownStars2.
             these functions leave a null value in the 'star' field in case of a 
             non-catalogued point index.
*/
long ImageMatrix::MarkUnknownStars(star_point *star_points,stripe_point *stripes,char *filename,BITMAP bitmap)
{  long x,y,close_stars_num,b,star_point_index=0,unexpected_num,label_index=0,stripe_index=0;
   double alt,az,ra,dec;
   long C1_B,C5_B,C9_B,C16_B,C25_B;
   char jpgfile[100];

   /* count the number of transients - if there are too much then it is because the picture is bad */
   unexpected_num=0;
   while (star_points[star_point_index].x>=0) 
     if (star_points[star_point_index++].star==NULL) unexpected_num++;
   if (unexpected_num>=20) return(0);  /* probably bad image */

   /* find the index of the last label */
   while (labels[label_index].x>0) label_index++;
   if (strchr(filename,'/')) filename=&(strrchr(filename,'/')[1]);
   strcpy(jpgfile,filename);
   if (strchr(jpgfile,'.')) strcpy(strchr(jpgfile,'.'),".jpg"); /* change the extention to "jpg" */
   /* write the stripes transients*/
   while (stripes[stripe_index].x>0 && stripe_index<STRIPES_NUM)
   {  if (!IsCloudArea((stripes[stripe_index].min_x+stripes[stripe_index].max_x)/2,(stripes[stripe_index].min_y+stripes[stripe_index].max_y)/2)) /* found a stripe */
      {  x=stripes[stripe_index].x;
         y=stripes[stripe_index].y;
         labels[label_index].x=x;
	 labels[label_index].y=y;
	 labels[label_index].type=STAR_TYPE_UNEXPECTED;
	 labels[label_index].text[0]='\0'; /* no label for the unexpected object */
	 labels[label_index+1].x=-1;
	 label_index++;
	 /* write the information to the std output stream */
	 xy2altaz(x,y,&alt,&az);
	 altaz2radec(&current_location,alt,az,&ra,&dec);  
         printf("<stripe>");
	 printf(" stripe:<a href=\"%s\">%s</a>|<a href=\"%s\">%s</a> jd=<jd>%.6f</jd> ra=~<ra>%2.2f</ra> dec=~<dec>%2.2f</dec> alt=<alt>%2.2f</alt> az=<az>%3.2f</az> X=<x>%d</x> Y=<y>%d</y><br>",filename,filename,jpgfile,jpgfile,current_date,ra,dec,alt,az,x+1,height-y);
         printf("</stripe>\n");
      }
      stripe_index++;
   }
   /* write the transients */
   unexpected_num=0;
   star_point_index=0;
   while (star_points[star_point_index].x>=0)
   { if (star_points[star_point_index].star==NULL)
     {  x=star_points[star_point_index].x;
        y=star_points[star_point_index].y;
#ifndef FITS
        b=star_point_index+1;
        if (star_points[star_point_index].bright_star==1) close_stars_num=10;  /* a bright star should be marked */
        else close_stars_num=0;
        while (star_points[b].x>=0)
        {  if (star_points[b].star==NULL)
           if (sqrt(pow(x-star_points[b].x,2)+pow(y-star_points[b].y,2))<25) close_stars_num++;
           b++;
        }
        if (close_stars_num>=1)  /* "unexpected" star(s) were found */
#endif
        {  if (!IsCloudArea(x,y))    /* the area is not a cloud */
	   /* assign a value to all the stars around so they won't be marked also */
	   {  b=star_point_index+1;    /* make sure no square will appear close to this square */
	      while(star_points[b].x>=0)
	      {  if (star_points[b].star==NULL)
	         if (sqrt(pow(x-star_points[b].x,2)+pow(y-star_points[b].y,2))<25)
	         star_points[b].star=(cat_entry *)1;   /* so it won't be circles again as an "unexpected star" */
	         b++;
              }
              /* add the square to the labels */
	      labels[label_index].x=x;
	      labels[label_index].y=y;
	      labels[label_index].type=STAR_TYPE_UNEXPECTED;
	      labels[label_index].text[0]='\0'; /* no label for the unexpected object */
	      labels[label_index+1].x=-1;
	      label_index++;
              /* write the information to the std output stream */
	      xy2altaz(x,y,&alt,&az);
              altaz2radec(&current_location,alt,az,&ra,&dec);
              C1_B=star_points[star_point_index].C1-star_points[star_point_index].background;
              C5_B=star_points[star_point_index].C5-star_points[star_point_index].background;
              C9_B=star_points[star_point_index].C9-star_points[star_point_index].background;
              C16_B=star_points[star_point_index].C16-star_points[star_point_index].background;
              C25_B=star_points[star_point_index].C25-star_points[star_point_index].background;
              printf("<transient>");
	      printf(" <a href=\"%s\">%s</a>|<a href=\"%s\">%s</a> jd=<jd>%.6f</jd> ra=~<ra>%2.2f</ra> dec=~<dec>%2.2f</dec> alt=<alt>%2.2f</alt> az=<az>%3.2f</az> X=<x>%d</x> Y=<y>%d</y> C1-b=<C1-B>%d</C1-B> C5-B=<C5-B>%d</C5-B> C9-B=<C9-B>%d</C9-B> C16-B=<C16-B>%d</C16-B> C25-B=<C25-B>%d</C25-B><br>",filename,filename,jpgfile,jpgfile,current_date,ra,dec,alt,az,x+1,height-y,C1_B,C5_B,C9_B,C16_B,C25_B);
	      printf("</transient>\n");
	      unexpected_num++;
           }
        }
      }
      star_point_index++;
   }
   return(unexpected_num);
}

/* ****************************************** */
/* **************** unwrap ****************** */
/* transforms the image into a panoramic view
   parameters - bitmap - the bitmap object on which
    the image is painted.
    scale - double: I.e., scale=2.0 means the output
    bitmap is 720x180
/* ****************************************** */

int ImageMatrix::unwrap(BITMAP h_bitmap, double scale)
{  double az,alt;
   long x,y;
   long val;
   for (az=0;az<360;az=az+1/scale)
     for (alt=0;alt<90;alt=alt+1/scale)
     {
        altaz2xy(alt,az,&x,&y);
        val=data[x][y];
        if (val<AvgColor*0.5) val=0;
        else
	if (val>2.0*AvgColor) val=255;  /* 3 */
        else
	val=(long)(((double)val/(2.0*AvgColor))*255);
        x=(long)(scale*az);
        y=(long)(scale*(90-alt));
        cd_rawassign(h_bitmap,x,y,0,(unsigned char)val);
        cd_rawassign(h_bitmap,x,y,1,(unsigned char)val);
        cd_rawassign(h_bitmap,x,y,2,(unsigned char)val);
     }
}