/***********************************************************************************************************
* prominence -f /path/demfile [-d 0] [-r 100000]                                                           *
*                                                                                                          *
* Calculates prominence of peaks based on elevation data                                                   *
*                                                                                                          *
* For installation you need the GDAL-libraries (debian: package libgdal-dev)                               *
* Compiling with "cc -Wall -o prominence prominence.c -lgdal -lm"                                          *
*                                                                                                          *
************************************************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include "gdal/gdal.h"
#include "gdal/cpl_conv.h"
#include "gdal/gdal_frmts.h"

struct list_peak{
 long long int    id;
 long long int    prom_id;
 double           lon,lat,ele,prominence;
 struct list_peak *next;
};    
  
   

#define PI             3.1415927
int debuglevel=0;


void printhelp(){
 printf("Usage:   prominence -f demfile [options]\n");
 printf("         -f demfile\n         -d debuglevel (0-4, default=0)\n         -r maximal search radius (100..500000,default 100000)\n         -o output format csv or sql (default csv)\n");
 printf("demfile is a geotiff (1 band int16, SRS=EPSG 4326)\n");
 printf("stdin   is a csv with \"id;lon;lat;ele\" where \n");
 printf("         ele is integer, float or empty (not mixed like \"1234m\")\n");
 printf("stdout  is a csv with \"id;lon;lat;prominence\", where\n");
 printf("         prominence is integer [0..radius]\n");
 printf("More documentation you will find in the OSM-Wiki:\n\n");
}


double get_height(double lon,double lat,double *adfGeoTransform,GDALRasterBandH hBand,long int xsize,long int ysize){
/* get height at lon/lat with elevation data from hBand (parameters in adfGeoTransform) */

 long int   xpx,ypx;
 double     h;
 int16_t    *area2x2;

 /* calculate upper left pixel next to lon/lat  */

 xpx=(long int)floor((lon-adfGeoTransform[0])/adfGeoTransform[1]);
 ypx=(long int)floor((lat-adfGeoTransform[3])/adfGeoTransform[5]);
    
/* read a 2x2 area from demfile at this pixel, so getting all 4 pixel next to lon/lat  */
/* take the max. height from these 4 values                                            */
 
 area2x2=malloc(4*sizeof(int16_t));
 GDALRasterIO(hBand,GF_Read,xpx,ypx,2,2,area2x2,2,2,GDT_Int16,0,0);
 if(debuglevel>3){printf("  Getting x/y:%ld %ld at %lf %lf h=%d %d %d %d\n",xpx,ypx,adfGeoTransform[0]+xpx*adfGeoTransform[1],adfGeoTransform[3]+ypx*adfGeoTransform[5],area2x2[0],area2x2[1],area2x2[2],area2x2[3]);}
 h=area2x2[0];
 if(area2x2[1]>h){h=area2x2[1];}
 if(area2x2[2]>h){h=area2x2[2];}
 if(area2x2[3]>h){h=area2x2[3];}
 free(area2x2);
 if(h<-32000){h=0.0;}
 return h;
}


void get_prominence_by_ele(struct list_peak *firstpeak,double radius){
/* gets a chained list of peaks, calculate prominence of each peak, returns nothing */

 struct list_peak *peak=NULL,*otherpeak=NULL;
 double           ccos,d,dx,dy;
    
 peak=firstpeak;
 while(peak!=NULL){
  ccos=cos(peak->lat*PI/180);
  otherpeak=peak->next;
  if(peak->prominence==100){peak->prominence=radius;}
  if(debuglevel>3){printf("calc dist for peak %lld\n",peak->id);}
  while(otherpeak!=NULL){
   if(otherpeak->prominence==100){otherpeak->prominence=radius;}
   dx=(otherpeak->lon-peak->lon)*ccos*40000000/360;
   dy=(otherpeak->lat-peak->lat)*40000000/360;
   d=sqrt(dx*dx+dy*dy);
   if(d>radius){d=radius;}
   if(d>100){
    if((d < otherpeak->prominence)&&(peak->ele >= otherpeak->ele)){
     otherpeak->prominence=d;
     otherpeak->prom_id=peak->id;
     if(debuglevel>3){printf("   assign prominence to other peak %lld, d=%f\n",otherpeak->id,d);}
    }
    if((d < peak->prominence)     &&(peak->ele <= otherpeak->ele)){
     peak->prominence=d;
     peak->prom_id=otherpeak->id;
     if(debuglevel>3){printf("   found higher peak %lld for peak %lld, d=%f\n",otherpeak->id,peak->id,d);}
    }
   }
   otherpeak=otherpeak->next;
  }
  peak=peak->next;
 }
}
                                    


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

 GDALDatasetH     hDataset;
 GDALRasterBandH  hBand;
 double           adfGeoTransform[6];
 char             *demfile=NULL;
 char             *outputformat=NULL;
 char             *textrest;
 double           radius=100000.0;
 char             line[100],elestring[100];
 long long int    id;
 long int         i=0,st_s=0,st_g=0;
 double           lon,lat,ele;
 struct list_peak *firstpeak=NULL,*peak=NULL,*lastpeak=NULL;

/* Parse command line parameters */

 if(argc<3){printhelp();exit(1);}
 for(i=1;i<=argc-2;i+=2){
  if     (strcmp(argv[i],"-d")==0){debuglevel=atoi(argv[i+1]);}
  else if(strcmp(argv[i],"-f")==0){demfile=argv[i+1];}          
  else if(strcmp(argv[i],"-r")==0){radius=strtod(argv[i+1],&textrest);}          
  else if(strcmp(argv[i],"-o")==0){outputformat=argv[i+1];}          
  else if(strcmp(argv[i],"-h")==0){printhelp();exit(0);}          
  else if(strcmp(argv[i],"-?")==0){printhelp();exit(0);}
 }   
 if(!demfile){printhelp();exit(1);}
 if(radius<100){radius=100;}
 if(radius>500000){radius=500000;}

/* open demfile, get srs, origin, pixel size and pointrer to the band 1 */

 GDALRegister_GTiff();
 hDataset = GDALOpen(demfile,GA_ReadOnly);
 if( hDataset == NULL ){
  fprintf(stderr,"Cannot open File %s\n",demfile);
  exit(1);
 }
 if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None ){
  if(debuglevel>0){
   printf( "Reading %s %dpx x %dpx x %d band\n",demfile,GDALGetRasterXSize(hDataset),GDALGetRasterYSize(hDataset),GDALGetRasterCount(hDataset));
   printf( "Origin =(%.6f,%.6f) ",adfGeoTransform[0],adfGeoTransform[3]);
   printf( "Pixel Size=(%.6f,%.6f)\n",adfGeoTransform[1],adfGeoTransform[5]);
  }
  hBand=GDALGetRasterBand(hDataset,1);
 }else{
  fprintf(stderr,"Could not fetch transformation parameters from %s\n",demfile);
  exit(1);
 }

/* Read stdin, split each line in id,lon,lat,direction */

 while(fgets(line,100,stdin)){
  i=sscanf(line,"%lld;%lf;%lf;%[^\n]",&id,&lon,&lat,elestring);
  if(i==3){elestring[0]='\0';}
  strtok(elestring,";\n");
  ele=strtod(elestring,&textrest);

/* If ele is not a float, get it from DEMfile... */

  if(((ele<0.1)&&(ele>-0.1))||(textrest[0]!='\0')){
   ele=get_height(lon,lat,adfGeoTransform,hBand,GDALGetRasterXSize(hDataset),GDALGetRasterYSize(hDataset));
   elestring[0]='\0';
   st_s++;
   if(debuglevel>2){
    printf("\nInput: %s",line);
    printf(" -> ele: id=%lld lon=%f lat=%f ele=%f\n",id,lon,lat,ele);
   }
  }
  if(i<3){id=0;}
 
/* Append peak to linked list of peaks */  
  
  if(id){
   peak=malloc(sizeof(struct list_peak));
   if(!peak){
    fprintf(stderr,"Out of memory at peak no %ld\n",st_s);
    exit(1);
   }
   peak->id=id;peak->lon=lon;peak->lat=lat;peak->ele=ele;peak->prominence=100;peak->next=NULL;
   if(firstpeak==NULL){
    firstpeak=peak;
   }else{
    lastpeak->next=peak;
   }
   lastpeak=peak;
  }
  st_g++;
 }
 if(debuglevel>0){printf("Got %ld peaks, %ld without parseable ele\n",st_g,st_s);}
 get_prominence_by_ele(firstpeak,radius);
 if(debuglevel>0){printf("Calculated peaks\n");}
 peak=firstpeak;
 while(peak!=NULL){
  if((!outputformat)||(strcmp(outputformat,"csv")==0)){printf("%lld;%.7lf;%.7lf;%.0lf;%.0lf;%lld\n",peak->id,peak->lon,peak->lat,peak->ele,peak->prominence,peak->prom_id);}
  peak=peak->next;
 } 
 exit(0);
}