#include <stdio.h>
#include <math.h>

/* nearneighbor.c -- Computes the nearest neighbor interpolation of MOLA data 

   Oleg Abramov
   Planetary Image Research Laboratory
   abramovo@pirl.lpl.arizona.edu
 
   version 0.2 

*/
 
#define GRIDX   600
#define GRIDY   200
#define MIN_X   161
#define MAX_X   167
#define MIN_Y   72
#define MAX_Y   74

#define N       50
#define WEIGHT(x) 1/(x*x)
#define NUMIN   80732

#define MAX_DIST 1000000000  

main(int argc, char *argv[])
{
  
  FILE *in;

  float delta_x, delta_y, distance, *largest_element;
  float lon_array[NUMIN], lat_array[NUMIN], elev_array[NUMIN], closest_distance[N], *closest_elevation[N];
  float x_array[GRIDX+1], y_array[GRIDY+1] ;
  int nitems, raw_point, raw_count = 0, pointcount = 0, count_x, count_y, grid_point_x, grid_point_y, temp, temp2;
  float weighted_average, average_numerator, average_denominator;
  double x, y, increment_x, increment_y;


  delta_x = MAX_X - MIN_X; 
  delta_y = MAX_Y - MIN_Y;
  increment_x = delta_x / GRIDX; 
  increment_y = delta_y / GRIDY;

  /* Check if the file argument is present */
  if (!argv[1]) {
     printf( "\nNo input file specified\n" );
     exit(0);
  }

  /* Check if the file exists */
  in=fopen(argv[1],"r");
  if (in == NULL) {
      printf( "\nUnable to open file\n" );
      exit(0);
  } 
  
  /* Read in longitude, latitude, and elevation from file */
  while(1)
  {
    nitems = fscanf(in, "%f\t%f\t%f\n", &lon_array[raw_count], &lat_array[raw_count], &elev_array[raw_count]);
    if (nitems == EOF)
      break;
    raw_count++;

  }   

  /* Set up a grid */

  count_x = 0;
  count_y = 0;

  for (y=MIN_Y; y <= MAX_Y; y = MIN_Y + count_y * increment_y ) {
    y_array[count_y] = y;
    count_y++;
  }
 
  for (x = MIN_X; x <= MAX_X; x = MIN_X + count_x * increment_x ) {
    x_array[count_x] = x;
    count_x++;
  }


  /* Find the nearest neighbor */

  for (grid_point_y = 0; grid_point_y <= GRIDY; grid_point_y++) {
    for (grid_point_x = 0; grid_point_x <= GRIDX; grid_point_x++) {

      /* Set every element in the dist_array to MAX_DIST */
    for (temp=0; temp < N; temp++) {
      closest_distance[temp] = MAX_DIST;
    }
 
    for (raw_point = 0; raw_point < raw_count; raw_point++ )  {
      distance = sqrt((x_array[grid_point_x] - lon_array[raw_point]) * (x_array[grid_point_x] - lon_array[raw_point]) + (y_array[grid_point_y] - lat_array[raw_point]) * (y_array[grid_point_y] - lat_array[raw_point]));
      
      /*find the largest element in the closest_distance array */
      largest_element = &closest_distance[0];
      temp2=0;
      for (temp=0; temp < N; temp++) {
        if ( *largest_element < closest_distance[temp] )
           {  largest_element = &closest_distance[temp]; temp2 = temp; 
        }
      }
      if (distance < *largest_element) {
       closest_distance[temp2] = distance; closest_elevation[temp2] = &elev_array[raw_point];  
      }
    }

    /* Compute weighted average */

    average_numerator = 0;
    average_denominator = 0;

    for (temp=0; temp < N; temp++) {
      /* Avoiding division by zero */
      if (closest_distance[temp] == 0) { closest_distance[temp] = 0.001; } 
      average_numerator =  average_numerator + *closest_elevation[temp] * WEIGHT(closest_distance[temp]);  
      average_denominator = average_denominator + WEIGHT(closest_distance[temp]);  
    }

    weighted_average = average_numerator / average_denominator;
    printf ("%f\t%f\t%f\n", x_array[grid_point_x], y_array[grid_point_y], weighted_average);   
  }
  }
}