// DLL source file pascal.cpp
#include "pascal.h"
#include <stdio.h>
#include <string.h>
#include <stdlib.h>     /* used for EXIT_SUCCESS */
#include <time.h>
#pragma warning(disable : 4996)
//LPSAFEARRAY lpsa1;
//LPSAFEARRAY lpsa2;

#define FALSE   0
#define TRUE    1


double __stdcall  MyFunc_at2(double a, double b)
{
    return atan2(a,b);
}

double __stdcall  MyFunc_at(double a)

{
    return atan(a);
}


double __stdcall  MyFunc_fl(double a)
{
    return floor(a);
}


int __stdcall XYRectify(char *infilename, char *outfilename)
{
	// This routine takes a Riegl RGB image and rectifies it into a pinhole image
	// Finally the image is interpolated for black pixels inside the image area

   char *buffer; int ImWidth, ImHeight, AddDim;
   int i, j, l;
   double c = 0.0516906;
   double xp = -0.0749;
   double yp = -0.0491;
   double k1 = 0.0000152902;
   double k2 = -0.00000000534204;
   double k3 = 8.63902E-13;
   double p1 = -0.00000311216;
   double p2 = 0.00000322987;
   double pixW = 0.0000046;
   ImWidth = 11608; // Dimensions of the input (RAW) image with distortions
   ImHeight = 8708;
   AddDim = 200; // We add marginal to the rectified image since the pinhole image has round edges
   double xm, ym;


   buffer = (char *) malloc(200 * sizeof(char)); // Needed for the msgbox() routine
   FILE *fin, *fout;

   unsigned char *InImage[12000]; // A row in the input image has a pointer stored 
   
   fin = fopen(infilename,"rb");  // Open in raw image
   if (fin==NULL)
   {
	   l  = sprintf(buffer,"Tiedoston in avaus ei onnistunut \n");
       MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); 
	   free (buffer);
	   return 0;
   }
   // We are ok to continue, start allocating memory for the input image
   for (i=0;i<=ImHeight-1;i++) // Loop the rows
    {   
	 InImage[i] =  (unsigned char *) calloc(ImWidth*3,sizeof(unsigned char));
     if (InImage[i]==NULL)
	 {
      j  = sprintf(buffer,"XYrectify: Could not allocate enough memory, row: InImage %d \n",i);
      MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK);
	  fclose(fin);
      for (j=0;j<=i;j++)   { free(InImage[j]); }
  	  free(buffer);
	  return 0;
	 }
    }
   // We read the inimage
    int numcount, paikka;
    numcount = 0; paikka = 0;
    for (i=0;i<=ImHeight-1;i++)
	 {
  	  paikka=i*ImWidth*3;
      fseek(fin, paikka, SEEK_SET);
      numcount += fread(&InImage[i][0], sizeof(unsigned char)*ImWidth*3, 1, fin);
	  }

   i  = sprintf(buffer,"InImage READ rows: %d    \n", (numcount) );
   MessageBox ( NULL, buffer,"rajat",MB_OK); 
   fclose(fin);
//	return 1;

	// Now compute the rectified image
     unsigned char *OutImage;
	 OutImage =  (unsigned char *)calloc((ImHeight+AddDim)*(ImWidth+AddDim)*3,sizeof(unsigned char));
     if (OutImage==NULL)
	 {
      j  = sprintf(buffer,"XYrectify: Could not allocate enough memory for output image, row: InImage %d \n",i);
      MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK);
      for (j=0;j<=ImHeight-1;j++)   { free(InImage[j]); }
  	  free(buffer);
	  return 0;
	 }


   int row, col, rowi, coli; double x, y, r2, R, dr;
   unsigned char red, grn, blu;

   fout = fopen(outfilename,"wb");  // Open in raw image
   if (fout==NULL)
   {
	   l  = sprintf(buffer,"Tiedoston in avaus ei onnistunut \n");
       MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); 
 	   free (buffer);
       for (j=0;j<=ImHeight-1;j++)   { free(InImage[j]); }
	   free(OutImage); 
       free(buffer);
	   return 0;
   }
   double xc, yc;
   for (row=0;row<=ImHeight-1;row++) // Loop the rows
   {
	   for (col = 0; col <= ImWidth-1;col++)
	   {
         red = InImage[row][col*3];
         grn = InImage[row][col*3+1];
		 blu = InImage[row][col*3+2];
         xm = (col - (ImWidth - 1) / 2.0) * pixW * 1000.0;
         ym = ((ImHeight - 1) / 2.0 - row) * pixW * 1000.0;
         x = xm - xp;
         y = ym - yp;
         r2 = pow(x,2) + pow(y,2);
         R = pow(r2,0.5);
         dr = k1 * pow(R,3) + k2 * pow(R,5) + k3 * pow(R,7); 
         xc = xm - xp + x * dr / R + p1 * (r2 + 2 * pow(x,2)) + 2 * p2 * x * y;
         yc = ym - yp + y * dr / R + p2 * (r2 + 2 * pow(y,2)) + 2 * p1 * x * y;
         coli = int((ImWidth + AddDim - 1) / 2 + xc / (pixW * 1000.0));
         rowi = int((ImHeight + AddDim - 1) / 2 - yc / (pixW * 1000.0));
         // If CInt(coli) < ImWidth + 200 And CInt(coli) > 0 And CInt(rowi) < ImHeight + 200 And CInt(rowi) > 0 Then
		 // fseek(fout, long(rowi) * long(ImWidth+AddDim)*3 + coli * 3, SEEK_SET);
		 paikka = long(rowi) * long(ImWidth+AddDim)*3 + coli * 3;
		 OutImage[paikka]=red;
		 OutImage[paikka+1]=grn;
		 OutImage[paikka+2]=blu;
		 //Put #2, 1 + CLng(coli) * 3 + (CLng(rowi)) * (ImWidth + 200) * 3, pixel
         //fwrite(&red, sizeof(unsigned char), 1, fout);
         //fwrite(&grn, sizeof(unsigned char), 1, fout);         
		 //fwrite(&blu, sizeof(unsigned char), 1, fout);
	   }
   }

   fwrite(&OutImage[0], sizeof(unsigned char)*(ImHeight+AddDim)*(ImWidth+AddDim)*3, 1, fout);
   i  = sprintf(buffer,"InImage WROTE rows: %d    \n", (row) );
   MessageBox ( NULL, buffer,"rajat",MB_OK); 

   fclose(fout);
   for (i=0;i<=ImHeight-1;i++){free (InImage[i]);}
   free(OutImage); 
   free(buffer);
   return 1;

   }






int __stdcall MyFuncLargeFile(char *infilename, int startrec, int NumofRecs, char *outfilename)
{  
   int l, k;
   char *buffer;
   buffer = (char *) malloc(200 * sizeof(char));
   k  = sprintf(buffer,"Start at %s record\n",infilename);
   MessageBox (NULL, buffer,"Hello 64-bit world",MB_OK);

   k  = sprintf(buffer,"Start at %d\n record",startrec);
   MessageBox (NULL, buffer,"Hello 64-bit world",MB_OK);
   k  = sprintf(buffer,"Num of  %d\n records",NumofRecs);
   MessageBox (NULL, buffer,"Hello 64-bit world",MB_OK);

   FILE *fin, *fout;
   __int64 i; __int64 offset;

   unsigned char Hdr[227];
   unsigned char Point[28];
   fin = fopen(infilename,"rb");
   if (fin==NULL)
   {
	   l  = sprintf(buffer,"Tiedoston in avaus ei onnistunut \n");
       MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); 
	   return 0;
   }
   fout = fopen(outfilename,"wb");
   if (fout==NULL)
   {
	   l  = sprintf(buffer,"Tiedoston out avaus ei onnistunut \n");
       MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); 
	   return 0;
   }

   l = _fseeki64(fin,0,SEEK_SET);
   i = _ftelli64(fin);
   k = fread(&Hdr[0],sizeof(Hdr),1,fin); // Read the public header block of 227 bytes
   l = fwrite(&Hdr[0],sizeof(Hdr),1,fout); // Write the public header block of 227 bytes
   i = _ftelli64(fout);
   int m;
   m = 0;
   for (l = startrec; l <= (NumofRecs+startrec);l++) // read the 28 -byte records
   {
	m++;
    offset = __int64(sizeof(Hdr)) + __int64(l-1) * __int64(28);
	k = _fseeki64(fin,offset,SEEK_SET);
    k = fread(&Point[0],28,1,fin); // Read the 28-byte information
	offset = __int64(sizeof(Hdr)) + __int64(m-1) * __int64(28);
	k = _fseeki64(fout,offset,SEEK_SET);
    k = fwrite(&Point[0],28,1,fout); // Write the 28-byte information
    //k  = sprintf(buffer,"Wrote at %I64u\n record",offset);
    //MessageBox (NULL, buffer,"Hello 64-bit world",MB_OK);
   }
   fclose(fin);
   fclose(fout);
   k  = sprintf(buffer,"Wrote at %I64u\n Bytes",offset);
   MessageBox (NULL, buffer,"Hello 64-bit world",MB_OK);
   free(Hdr); 
   free(Point);
   
   return 0;
}

int __stdcall MyFunc_Get_LiDAR(Point3D *LiDAR, Point3D *P_top, int Np, double Zbase, double maxrad, double TanSun, double Base_elev, double cf, double pw, double cont, double rmse)

{

// Reads LiDAR data files and returns the points that are used for testing occlusion/shading, and relative height estimates
	
FILE *fIn1, *fIn;
int i, *Nha,Npulses, j, k;
char *buffer;
char *tiedostonimi;

struct LidarRecord
 { double GPSTime;
   char PulseCount;
   struct Point3D Pulse[4];
   short  Intensity[4];
   double Range[4];
   double Angle;
   double Roll;
   double Pitch;
   double Heading;
   struct Point3D Poslidar;
   short StripNum;
   char SyncBit;
   char Res1;
   char Res2;
   char Res3;
} LidR;

Nha = (int *) malloc(10 * sizeof(int));
buffer = (char *) malloc(200 * sizeof(char));
tiedostonimi = (char *) malloc(100 * sizeof(char));
double Htree; double appu; double radiustest;double radius08;double testrad;double baseheight;int N;double candrad;
int l; int m;

    N=0;
 	fIn1=fopen("C:\\data\\Filenames.txt","r");
    if (fIn1==NULL)
   	{
     i  = sprintf(buffer,"Filenames.txt:n avaus ei onnistunut \n");
     MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); 
	 return 0;
	}
	j = 0; i = 0; Npulses = 0; m=0;
	while (i!=EOF)
	{
		i = fscanf(fIn1,"%s", tiedostonimi);
	    if (i==EOF){break;}
        fIn=fopen(tiedostonimi,"rb");
        if (fIn==NULL)
		{
         l  = sprintf(buffer,"LiDAR-bin tiedoston avaus ei onnistunut \n");
         MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 
		 fclose(fIn);
	     return 0;
		}
		m = m+1;
		l = fread(&Nha[m],sizeof(int),1,fIn);

        Htree = P_top->Z - Zbase;  baseheight = P_top->Z - Htree;
//		k  = sprintf(buffer," Pulses %d %d\n",Nha[m],m);
//        MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 

		for (k = 0; k<=Nha[m]-1;k++) 
		{
  		l = fread(&LidR,sizeof(LidarRecord),1,fIn);

  	    if (LidR.Pulse[3].Z > Base_elev) 
	    {
		 //candrad = pow(pow((LidR.Pulse[3].X - P_top->X),2)+pow((LidR.Pulse[3].X - P_top->Y),2),0.5);

         for (j = 3; j >= (4-LidR.PulseCount);j--)
		 {
		    appu = (P_top->Z - LidR.Pulse[j].Z) / Htree;
            radiustest = (cf * Htree * (2.5 * pow(appu,pw)) + pow(cont,2));
            radius08 = (cf * Htree * (2.5 * pow(0.2,pw)) + pow(cont,2));
            testrad = radiustest + pow(rmse,1.5);
            candrad = pow(pow((LidR.Pulse[j].X - P_top->X),2) + pow((LidR.Pulse[j].Y - P_top->Y),2),0.5);
			
			if ((candrad > testrad) && (candrad < maxrad) )
			{

				if (LidR.Pulse[j].Z > Base_elev && N < Np)  // restrict to the upper boud of the array
				{
                  LiDAR[N].X = LidR.Pulse[j].X;
				  LiDAR[N].Y = LidR.Pulse[j].Y;
				  LiDAR[N].Z = LidR.Pulse[j].Z;
                  N++;
				}
			}
		 }
		}
		}
		fclose(fIn);
	}
     
 //   while loop
free(Nha) ;  
fclose(fIn1);
N--;
//	k  = sprintf(buffer," N %d\n",N);
//	MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 

return N;
}

int _stdcall MyFunc_Cast_Image_ray_and_Find_Closest_Riegl_LiDAR(Point3D *X1, Point3D *Ray, Point3D *X2, Point3D *PMin, Point3D *PMax, double distlimit, int *MaxI, Point3D *TopPoint)
{

	// May 19, 2016 this routines opens LiDAR point files and finds the closest LiDAR to a ray
	// infile points to file with filenames of BINary LiDAR files. Adjusted for Riegl format.

struct Point3D 
{
double X, Y, Z;
} ;

FILE *fIn1, *fIn;
int i, *Nha,Npulses, j, l, Pointer, *WFPointer, k, m, jj, Nechoes, allocated;
char *buffer, *tiedostonimi, echocount;
double GpsTime;
char ByteLAS, ScanAngle, StripNum;
struct Point3D PosLiDAR;
struct Point3D *pos;
short Intensity, WaveSampleSize;
float RetPointWFLoc;
Nha = (int *) malloc(10 * sizeof(int));
buffer = (char *) malloc(200 * sizeof(char));
tiedostonimi = (char *) malloc(100 * sizeof(char));

    // The file c:\data\Filenames.txt has the filenames for the 1-ha binary LiDAR data
 	fIn1=fopen("C:\\data\\RieglFilenames.txt","r");
    if (fIn1==NULL)	{
     jj  = sprintf(buffer,"Filenames.txt:n avaus ei onnistunut \n");
     MessageBox ( NULL, buffer,"MyFunc_Cast_Image_ray_and_Find_Closest_Riegl_LiDAR",MB_OK); 
	 return 0;}
   // if (fIn1!=NULL)	{
     //jj  = sprintf(buffer,"Filenames.txt:n avaus onnistui \n");
     //MessageBox ( NULL, buffer,"MyFunc_Cast_Image_ray_and_Find_Closest_Riegl_LiDAR",MB_OK); }

	j = 0;
	i = 0;
	Npulses = 0;
	l = 0;
    Nechoes = 0;  // pos[] stores echoes, not points
	allocated = 0;
	while (i!=EOF)
	{

		i = fscanf(fIn1,"%s", tiedostonimi);
        if (i==EOF){break;}
        fIn=fopen(tiedostonimi,"rb");
	//	MessageBox ( NULL, tiedostonimi,"MyFunc_Cast_Image_ray_and_Find_Closest_Riegl_LiDAR",MB_OK);
        if (fIn==NULL)
		{
         jj  = sprintf(buffer,"LiDAR-bin tiedoston avaus ei onnistunut \n");
         MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 
		 fclose(fIn);
	     return 0;
		}
		j = j+1;
		i = fread(&Nha[j],sizeof(int),1,fIn);
	// jj = sprintf( buffer, "%d", Nha[j] );
	//	MessageBox ( NULL,buffer ,"MyFunc_Cast_Image_ray_and_Find_Closest_Riegl_LiDAR",MB_OK);
		if (i!=NULL) 
		{
		Npulses = Npulses + Nha[j];
		// On 1st time call malloc, declare storage for average five echoes per pulse

    	if (j == 1) 
		{
			pos = (Point3D *) malloc(4 * Npulses * sizeof(Point3D));
			allocated = 4 * Npulses;
		}

		}
		if (pos == NULL || Npulses==0)
		{
         jj  = sprintf(buffer,"Muistinvaraus ei onnistunut\n");
         MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 
		 fclose(fIn);
		 return 0;
		}
		
		//jj = sprintf( buffer, "%d", ftell(fIn));
        //MessageBox ( NULL,buffer ,"MFile position",MB_OK);
        if (Nha[j]>0)
		{
			k = 0;
			Pointer = 4;
			while (k < (Nha[j]-1)) // Lets read the pulses
			{
				k = k + 1; // In this file
				if (k == 1) {fseek(fIn,4,SEEK_SET);} // first pulse starts at file position 4
                if (k != 1) {fseek(fIn, Pointer,SEEK_SET);} // Else read from Position Pointer
                
		   //     jj = sprintf( buffer, "%d %d", ftell(fIn), k);
           //     MessageBox ( NULL,buffer ,"Pointer",MB_OK);

			    // general pulse info
				i = fread(&echocount,sizeof(char),1,fIn);
				Nechoes = Nechoes + (int)echocount;
   		        if (Nechoes > allocated)  {
		          jj = sprintf( buffer, "Cannot assign memory" );
                  MessageBox ( NULL,buffer ,"Cast",MB_OK);
				  free(pos);
				  return 0;
				}

				i = fread(&GpsTime,sizeof(double),1,fIn);
				if (GpsTime < 1000)
				{
				 jj = sprintf( buffer, "%d %f", ftell(fIn), GpsTime);
                 MessageBox ( NULL,buffer ,"Pointer GPS_Time",MB_OK);

				}
				i = fread(&ByteLAS,sizeof(char),1,fIn);

                for (m = 1; m<=(int)echocount;m++)
				{
				  l = l + 1; // Total number of echoes
			      i = fread(&pos[l].X,sizeof(double),1,fIn);
				  i = fread(&pos[l].Y,sizeof(double),1,fIn);
				  i = fread(&pos[l].Z,sizeof(double),1,fIn);
			      i = fread(&Intensity,sizeof(short),1,fIn);
				  i = fread(&WFPointer,sizeof(int),1,fIn);
				  i = fread(&WaveSampleSize,sizeof(short),1,fIn);
				  i = fread(&RetPointWFLoc,sizeof(float),1,fIn);
				}
				i = fread(&ScanAngle,sizeof(char),1,fIn);
                i = fread(&PosLiDAR,sizeof(Point3D),1,fIn);
                i = fread(&StripNum,sizeof(char),1,fIn);
				Pointer = Pointer + 1 + 8 + 1 + (int)echocount * 36 + 1 + 24 + 1;
			}
		}
   
		 if (i==NULL)
		{
         i  = sprintf(buffer,"LiDAR datan luku ei onnistunut\n");
         MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 
		 fclose(fIn1);
		 return 0;
		}

		fclose(fIn);
	}

 // At this point we should have the data in LidR[] struct-array
    
  Point3D X2mX1, X1mX0, XxY;
  double Norm2_X2X1;
  X2mX1.X = X2->X - X1->X;
  X2mX1.Y = X2->Y - X1->Y;
  X2mX1.Z = X2->Z - X1->Z;
  Norm2_X2X1 = pow(X2mX1.X*X2mX1.X + X2mX1.Y*X2mX1.Y + X2mX1.Z*X2mX1.Z,0.5);
  double dist, Hmax, H;
  int Nin;
  Nin = 0; Hmax = 0.0;

  //i  = sprintf(buffer,"Starting to loop %f %f\n", pos[1].X, PMin->X );
  //MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 
  //i  = sprintf(buffer,"Npulses %d\n", Nechoes);
  //MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 
  //i  = sprintf(buffer,"Pmin %f %f %f\n", PMin->X, PMin->Y, PMin->Z);
  //MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 
  //i  = sprintf(buffer,"Pmax %f %f %f\n", PMax->X, PMax->Y, PMax->Z);
  //MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 


   for (k = 1; k<=Nechoes-1;k++)
   {
	  // if (k % 8000 == 0) {
  		 //    i  = sprintf(buffer,"Ekaluuppi %f %\n", pos[k].X);
         //    MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 
	   //} 

	   if ((pos[k].X > PMin->X) && (pos[k].X < PMax->X)) 
	   {
		   if ((pos[k].Y > PMin->Y) && (pos[k].Y < PMax->Y)) 
		   {
			   if ((pos[k].Z > PMin->Z) && (pos[k].Z < PMax->Z))
			   {
				     // This is a LiDAR point for which Line-To-Point-3D-distance is needed
  		             X1mX0.X = X1->X - pos[k].X;
					 X1mX0.Y = X1->Y - pos[k].Y;
					 X1mX0.Z = X1->Z - pos[k].Z;
					 // We need cross product of X2mX1 x X1mX0
                     XxY.X = X2mX1.Y * X1mX0.Z - X2mX1.Z * X1mX0.Y;
                     XxY.Y = X2mX1.Z * X1mX0.X - X2mX1.X * X1mX0.Z;
                     XxY.Z = X2mX1.X * X1mX0.Y - X2mX1.Y * X1mX0.X;
                     dist = pow(XxY.X*XxY.X + XxY.Y*XxY.Y + XxY.Z*XxY.Z,0.5) / Norm2_X2X1;
		
					 if (dist < distlimit) 
					 {						   
  		//			  i  = sprintf(buffer,"dist %f\n", dist);
        //                MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 
					   Nin++;
					   H=pos[k].Z;
					   if (H > Hmax)
					   {
						   Hmax = H;
						   *MaxI = k;
					   }
					 }

			   }
		   }
	   }

   }


   TopPoint->X = pos[*MaxI].X;
   TopPoint->Y = pos[*MaxI].Y;
   TopPoint->Z = pos[*MaxI].Z;

   //i  = sprintf(buffer,"Finished %f %f %f\n", pos[*MaxI].X, pos[*MaxI].Y, pos[*MaxI].Z);
   //MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 

   free(pos);
   free(Nha);
   free(buffer);
   free(tiedostonimi);

   return 1;

}

int _stdcall MyFunc_Cast_Image_ray_and_Find_Closest_43Byte_LiDAR(Point3D *X1, Point3D *Ray, Point3D *X2, Point3D *PMin, Point3D *PMax, double distlimit, int *MaxI, Point3D *TopPoint)
{

	// May 19, 2016 this routines opens LiDAR point files and finds the closest LiDAR to a ray
	// infile points to file with filenames of BINary LiDAR files. Adjusted for Riegl format.
	// April 16, a new 43 byte per point format in use

struct fPoint3D 
{
double X;
double Y;
float Z;
} ;
struct dirvector
{
	float i;
	float j;
	float k;
};

struct LidarRecord
 { struct fPoint3D echo;
   short  Intensity;
   unsigned char ScanAngle;
   unsigned char echocount;
   unsigned char echonum;
   unsigned char wavelen;
   unsigned char agc;
   float range;
   struct dirvector dirvect; 
} *LidR;

struct Point3D {
	double X, Y, Z;
} ;

FILE *fIn1, *fIn;
int i, *Nha,Npulses, j, l, jj, Nechoes;
char *buffer, *tiedostonimi; 

Nha = (int *) malloc(10 * sizeof(int));         // array if several bin files comprise the point cloud
buffer = (char *) malloc(200 * sizeof(char));
tiedostonimi = (char *) malloc(100 * sizeof(char));

    // The file c:\data\Filenames.txt has the filenames for the 1-ha binary LiDAR data
 	fIn1=fopen("C:\\data\\43ByteFilenames.txt","r");
    if (fIn1==NULL)	
	 {jj  = sprintf(buffer,"Filenames.txt:n avaus ei onnistunut \n");
     MessageBox ( NULL, buffer,"MyFunc_Cast_Image_ray_and_Find_Closest_43Byte_LiDAR",MB_OK); 
	 return 0;}
   // if (fIn1!=NULL)	{
     //jj  = sprintf(buffer,"Filenames.txt:n avaus onnistui \n");
     //MessageBox ( NULL, buffer,"MyFunc_Cast_Image_ray_and_Find_Closest_Riegl_LiDAR",MB_OK); }

	j = 0;	i = 0; 	Npulses = 0; 	l = 0;   Nechoes = 0;  // pos[] stores echoes, not points 

	while (i!=EOF)
	 {
		i = fscanf(fIn1,"%s", tiedostonimi);
        if (i==EOF){break;}
        fIn=fopen(tiedostonimi,"rb");
        if (fIn==NULL)
		  {jj  = sprintf(buffer,"LiDAR-bin tiedoston avaus ei onnistunut \n");
          MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 
		  fclose(fIn);
	      return 0;}
		j = j+1;
		i = fread(&Nha[j],sizeof(int),1,fIn);
	    //jj = sprintf( buffer, "%d", Nha[j] );
		//MessageBox ( NULL,buffer ,"MyFunc_Cast_Image_ray_and_Find_Closest_Riegl_LiDAR",MB_OK);
		if (i!=NULL) 
		  { Npulses = Npulses + Nha[j];
		  // On 1st time call malloc, declare storage for average five echoes per pulse, now echo records!
		  // On 1st time call malloc, later call realloc
		  if (j == 1) {LidR = (LidarRecord *) malloc(Npulses * sizeof(LidarRecord));}
		  if (j > 1)  {LidR = (LidarRecord *) realloc(LidR, Npulses * sizeof(LidarRecord));}	}

		if (LidR == NULL || Npulses==0)
		  {jj  = sprintf(buffer,"Muistinvaraus ei onnistunut\n");
          MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 
		  fclose(fIn);
		  return 0;}
		
        if (Nha[j]>0) 
		   {i = fread(&LidR[Npulses-Nha[j]],sizeof(LidarRecord),Nha[j],fIn);}

		 if (i==NULL)
   		   {i  = sprintf(buffer,"LiDAR datan luku ei onnistunut\n");
           MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 
		   fclose(fIn1);
		   return 0; }
		
		} // While
   

		fclose(fIn);
	



   Point3D X2mX1, X1mX0, XxY;
   double Norm2_X2X1;
   X2mX1.X = X2->X - X1->X;
   X2mX1.Y = X2->Y - X1->Y;
   X2mX1.Z = X2->Z - X1->Z;
   Norm2_X2X1 = pow(X2mX1.X*X2mX1.X + X2mX1.Y*X2mX1.Y + X2mX1.Z*X2mX1.Z,0.5);
   double dist, Hmax, H;
   int Nin, k;
   Nin = 0; Hmax = 0.0;

//   i  = sprintf(buffer,"Starting to loop %f %f\n",LidR[0].FirstPulse.X, PMin->X );
//   MessageBox ( NULL, buffer,"MyFuncCast",MB_OK); 


   for (k = 0; k<=Npulses-1;k++)
   {
	   if ((LidR[k].echo.X > PMin->X) && (LidR[k].echo.X < PMax->X) ) // require only first echo data && (LidR[k].echonum==1)

	   {
		   if ((LidR[k].echo.Y > PMin->Y) && (LidR[k].echo.Y < PMax->Y)) 
		   {
			   if ((LidR[k].echo.Z > PMin->Z) && (LidR[k].echo.Z < PMax->Z))
			   {
				     // This is a LiDAR point for which Line-To-Point-3D-distance is needed
  		             X1mX0.X = X1->X - LidR[k].echo.X;
					 X1mX0.Y = X1->Y - LidR[k].echo.Y;
					 X1mX0.Z = X1->Z - LidR[k].echo.Z;
					 // We need cross product of X2mX1 x X1mX0
                     XxY.X = X2mX1.Y * X1mX0.Z - X2mX1.Z * X1mX0.Y;
                     XxY.Y = X2mX1.Z * X1mX0.X - X2mX1.X * X1mX0.Z;
                     XxY.Z = X2mX1.X * X1mX0.Y - X2mX1.Y * X1mX0.X;
                     dist = pow(XxY.X*XxY.X + XxY.Y*XxY.Y + XxY.Z*XxY.Z,0.5) / Norm2_X2X1;
		
					 if (dist < distlimit) 
					 {
					   Nin++;
					   H=LidR[k].echo.Z;
					   if (H > Hmax)
					   {
						   Hmax = H;
						   *MaxI = k;
					   }
					 }

			   }
		   }
	   }

   }


   TopPoint->X = LidR[*MaxI].echo.X;
   TopPoint->Y = LidR[*MaxI].echo.Y;
   TopPoint->Z = LidR[*MaxI].echo.Z;

   free(LidR);
   free(Nha);
   free(buffer);
   free(tiedostonimi);

   return 1;

}