Attribute VB_Name = "Module1"
'Option Explicit

Rem Public variables declared here
Rem CONSTANTS

Rem Maximum number of images to be displayed (0 to MAXIMA-1), e.g. 8 makes 9 images at max

Public Const MAXIMA = 9
Public Const TO_RADIANS = 1.74532925199433E-02
Public Const TO_DEGREES = 57.2957795130823
Public Const PI = 3.14159265358979
Public Const CAM_TO_IMA = 0
Public Const IMA_TO_CAM = 1
Public Const N_CAND_MAX = 1000
Public Const MAXITE = 100
Rem Number of images in a SET
Public NumOfImages As Integer
Rem The descriptive name for the SET
Public ImagesetName As String
Rem An index
Public pointnumber As Integer

Rem FileOut holds "C:\TEMP\PIC#.BMP" 16 character - string, Image takes palce I/O via BMP-files and
Rem Picture -objects LOADIMAGE() -function

Dim FileOut As String * 16
Public BatchFileName As String
Public DemFileName As String

Rem ***************************************************************************
Rem Structs (Type) & VARS that have to do with saving / retrieving a measurement

Type PlotInfo
  Number As Long          ' Plot has an Id
  Radius As Double        ' Radius of the plot
  Buffer As Double        ' Additional radius, gives a buffer around the plot
  X As Double             ' X of center, on ground
  Y As Double             ' Y "
  z As Double             ' Z "
  H_lidar As Double       ' Maximum height of Lidar-points inside plot
  N_trees As Long         ' Counter for tree observations
  X_shift As Double       ' if shifted along X
  Y_shift As Double       ' if shifted along Y
  Ctext As String * 60    ' Description of the plot
End Type
Public Plot_Info As PlotInfo
Type MeasuredPoint
  plot As Long                      ' Point belongs to plot, key
  Num As Integer                    ' Point Number (Id for the 3D-measurement)
  X As Double                       ' X of space intersection
  Y As Double                       ' Y of space intersection
  z As Double                       ' Z of space intersection
  z_DTM As Single                   ' Z of GND
  RMSE As Single                    ' RMSE in um of the space intersection
  Corr_Match As Double              ' Correlation of the image matching
  images As Byte                    ' Number of images involved
  TreeSpecies As Byte               ' Tree species as observed by the user: 1=Pine, 2=Norway spruce, 3 = Birch, 5 = Aspen, 6 Other bl.
  cf As Double                      ' Parameters of the crown model
  pw As Double
  cont As Double
  Crown_RMSE As Single              ' RMSE of the lidar fit
  Lidar_Dcrm As Single
  Foto_Dcrm As Single               ' Fotogrammetric Crown width
  Foto_h As Single
  Foto_Corr As Single               ' Correlation of the Template MAtcihng for Crown width
  dbh_foto_Dcrm As Single           ' dbh computed based on height, sp and foto Dcrm
  dbh_lidar_Dcrm As Single          ' dbh   "                              lidar Dcr,
  H_lidar As Single                 ' Height of the highest lidar hit
  X_lidar As Double                 ' X of the highest lidar hit
  Y_lidar As Double                 ' Y of the highest lidar hit
  Z_lidar As Double                 ' Z of the highest lidar hit
  FieldNumber As Integer            ' Tree Number in field Data that the tree is exected to correspond to
  RefImage As Byte                  ' Index of the image used as reference in image matching
  Ima(0 To MAXIMA - 1) As Byte      ' Indeces of Images involved in space intersection
  ImageCode(0 To MAXIMA - 1) As Long  ' Indeces of Images involved in space intersection
  ima_x(0 To MAXIMA - 1) As Single  ' Image (COL) coordinates
  ima_y(0 To MAXIMA - 1) As Single  ' Image (ROW) coordinates
  TreeString As String * 30         ' Any notes
End Type

Public LidarFitDone As Boolean
Public FotoTMdone As Boolean

Type MeasuredDist
  Num As Integer                    ' Point Number (Id for the width-measurement)
  X As Double                       ' X of space intersection
  Y As Double                       ' Y of space intersection
  z As Double                       ' Z of space intersection
  FieldNumber As Integer            ' Tree Number in field Data that the tree is exected to correspond to
  FreeCode As Integer               ' Free coding
  Ima As Byte                       ' Index of Image (0 to MAXIMA)
  ImageCode As Long                 ' Image code number (Label, identifies)
  ima_x(0 To 1) As Double           ' Image x-coordinates, start-end
  ima_y(0 To 1) As Double           ' Image y-coordinates, start-end
  Scalea As Single                   ' Image scale at given Z
  Dist As Single                    ' Measured distance
  TreeString As String * 30         ' Any notes, 30 characters
  SunElev As Double
  SunAzi  As Double
  SunAziDiff  As Double
  NadirAngle As Double
  ObjectToCameraAzi As Double
End Type
Rem Measuring crown widths
Public WidthMeasurements(0 To MAXIMA - 1) As MeasuredDist


Rem Size of an element of MeasuredPoint is: 2+4+4+4+4+1+1+2+30+9+4*9+4*9 = 133 Bytes / element (or more, there's no sizeof)

Public MeasurementCounter As Integer         ' a counter that is updated automatically after saving a record
Public Measurement As MeasuredPoint          ' for writing (saving) a record
Public SavedMeasurement As MeasuredPoint     ' for reading a record
Public TreeData(1 To 100) As MeasuredPoint  ' storage to hold data (needed?)
Public TreedataFilename As String            ' Path & Filename where records are written / read from
Public TreeDataFilesaveOk As Boolean         ' A boolean that tells if the user is satisfied with the contents of the record (for SAVE)
' Public RMSE_sol As Double                    ' RMSE in um of current space intersection solution is stored here (copied to struct measurement.RMSE)


Rem *****************************************
Rem X,Y,Z Different versions of 3D-point

Public X_cen As Double, Y_cen As Double, Z_cen As Double  ' if the user wishes to give a coordinate to which all images are centered
Public X_sol As Double, Y_sol As Double, Z_sol As Double  ' current solution of space intersection is stored in these vars (copied to measurement.x, .y, .z)
Public X_ini As Double, Y_ini As Double, Z_ini As Double  ' Initial approx. for space resection (previous solution, or in case of failure nadir point of 2nd image | Z 2000)
Public X_start As Double, Y_start As Double, Z_start As Double ' When program initiates, nadir point of 2nd image, used for initial approx. if others fail

Rem **************************************************************************************************
Rem vars that have to do with warning the user of possible error in solving the correspondence problem

Public initialize As Boolean
Public RMSE_Beep_Level_One As Double                      ' A Beep -sound (alert) is set after space intersection to warn the user, specifies a level _1_ in um of RMSE)
Public RMSE_Beep_Level_Two As Double                      ' A Beep -sound (alert) is set after space intersection to warn the user, specifies a level _2_ in um of RMSE)
Public Beep_One_Frequency As Long                         ' Frequency in Hertz of level _1_ Beep (Beep API-function requires LONG)
Public Beep_Two_Frequency As Long                         ' Frequency in Hertz of level _2_ Beep
Public Beep_One_Duration As Long                          ' Duration in ms of level _1_ Beep
Public Beep_Two_Duration As Long                          ' Duration in ms of level _1_ Beep


Rem Other control vars

Public Epi_Depth As Double                                ' Depth of epipolar search in + - meters
Public StepValue As Double                                ' 3D-point (starting from solution)

Public itestop_X As Double, itestop_Y As Double, itestop_Z As Double

Public Heigth_Limit As Double
Public Numbers_Plotted As Boolean
Public Call_Clear_And_Plot_Measurements As Boolean
Public After_Mouse_down_Call_Epipolar_Line As Boolean


Rem ************************************************************
Rem Declarations for Elliptic kernel Cut  / TEMPLATE ACQUISITION

Type point3D
 X As Double
 Y As Double
 z As Double
End Type

Type Vector3D
 X As Double
 Y As Double
 z As Double
End Type

Type Point2D
 X As Double
 Y As Double
End Type
 
Type POINTAPI
 X As Long
 Y As Long
End Type

 
Type Ellipse
 A As Double  ' length in Z-direction (um)
 B As Double  ' length in perpendicular direction (um)
 X As Double  ' midpoint, m, camera coords
 Y As Double  ' midpoint, m,
 dx As Double ' translation because of Zasymmetry, m
 dy As Double ' translation because of Zasymmetry, m
 dx_col As Integer ' " pixels
 dx_row As Integer ' " pixels
 alpha As Double ' rotation angle (with respect to camera x-axis (positive counterclockwise)
End Type

Public Ellipse(0 To MAXIMA - 1) As Ellipse

Rem **************************************************************************************
Rem Declarations needed in creating correlation images & computing backprojections etc.

Rem Metric parameters of Matching algorithm
Rem Elliptic kernel
Public Zdiff As Double, TemplateWidth As Double, EllipseZasymmetry As Double
Rem Plot's radius, location, and XY-extent of search area
Public Plotradius As Double
Public square_side As Double

Type PlotCenter
 Code As String * 4
 XYangle As Double
 X As Double
 Y As Double
 z As Double
End Type


Public PlotCenter As PlotCenter
Rem Clustering parameters
Public r_limit As Double, xy_thin As Double
Rem XYZ-point set
Public Meanheight As Double, Zasymmetry As Double
Public gridXextent As Double, gridYextent As Double, gridZextent As Double
Public gridXtess As Double, gridYtess As Double, gridZtess As Double
Rem Time used for computations
Public MatchDate As String * 10   ' stores the date
Public CorTime(0 To MAXIMA - 1) As Double
Public BackProjTime As Double
Public MeshTime As Double
Rem Image Codes, LONG integers
Public ImageCodes(0 To MAXIMA - 1) As Long


Public Const maximasize = 1200
Public cor_arr(0 To MAXIMA - 1, 0 To maximasize, 0 To maximasize) As Double
Public R_tau() As Double
Public a_tau() As Single
Public n_tau() As Long
Public mc As Integer  ' Number of clusters found
Public n_candidate_points As Long
Public cluscoords() As point3D
Public clus() As Single
Public SearchSpaceData() As point3D
Public atausum() As Double


Type ClusMatchStruct
 Index As Integer            ' Index to arrays holding clus-data
 IsMatch As Boolean          ' True for a match with a tree
 IsInside As Boolean
 IsInBorderZone As Boolean   ' True for a cluster near (within half matchistance) border
 MatchTreeNum As Integer     ' Field number of tree matching with
 MatchTreeSpec As Integer
 MatchTreeStat As Integer
 MatchTreed13 As Double
 MatchTreeh As Double
 MatchTreeZbuttTach As Double
 MatchTreeZbuttDem As Double
 Xtree As Double             ' Coordinates of the tree top matched with
 YTree As Double             '     "
 Ztree As Double             '     "
 IsCommission As Boolean     ' True for a commision error (no match found for a candidate)
 Dist3D As Double            ' 3D distance to the matched tree
 NTreesInCylinder As Integer ' How many tree tops where in the match-cylinder
 IndecesTreesInCyl(1 To 10) As Integer
 Dist2dToTreesInCyl(1 To 10) As Double
 X As Double                 ' Coordinates of the cluster
 Y As Double                 '
 z As Double                 '
 Rvalue As Double            ' 3D-correlation value
 Npoints As Integer          ' Number of 3D-points that formed the cluster
End Type

Type TreeMatchStruct
 Num As Integer              ' Tree field number
 Spec As Integer             ' Species
 Status As Integer           ' Status
 d13 As Double               ' dbh
 h As Double                 ' height
 X As Double                 ' Treetop's coordinates (using correct z-model)
 Y As Double                 '
 z As Double                 '
 ZButtTach As Double           ' Elevation for butt, obtained with the
 ZbuttDEM As Double
 NclusInCylinder As Integer  ' Number of clusters in the match-cylinedr of the tree
 IndecesClusInCyl(1 To 10) As Integer
 Dist2dToClusInCyl(1 To 10) As Double
 IndexMatchClus As Integer   ' Index of the matched cluster
 Dist3D As Double            ' 3D distance betweem the top and the cluster
 IsMatch As Boolean          ' True for a match
 IsOmission As Boolean       ' True for a missed tree
 IsInBorderZone As Boolean   ' True for a cluster near (within half matchistance) border
 IsInside As Boolean
 Xclus As Double             ' clusters coordinates
 Yclus As Double             '
 Zclus As Double             '
End Type


Type TreeVect
 X As Double
 Y As Double
 Zbutt As Double
 Zdem As Double
 Ztop As Double
 d13 As Double
 Height As Double
 Num As Long
 Species As Long
 Status As Long
 Vol As Double
End Type

Rem Accuracy checking

Public FTrees(1 To 600) As TreeVect
Public ClusMatchStruct() As ClusMatchStruct
Public TreeMatchStruct() As TreeMatchStruct

Type ImageInfoForOutput
 ImCode As Long
 CorrTime As Double
 SunA As Double
 SunE As Double
 x0 As Double
 y0 As Double
 Z0 As Double
End Type

Type OutputVector
  PlotCode As String * 4
  Date As String * 10
  Time As String * 8
  NImagesInMatch As Integer
  ImInfo(0 To MAXIMA - 1) As ImageInfoForOutput
  BkProjTime As Double
  ModelTreeNum As Integer
  ModelTreeX As Double
  ModelTreeY As Double
  ModelTreeZ As Double
  Rlimit As Double
  XYthin As Double
  XYMatchDist As Double
  ZMatchDist As Double
  TWidth As Double
  EllZasym As Double
  Zdiff As Double
  ZDepth As Double
  Zasym As Double
  Meshdist As Double
  Meanheight As Double
  Plotradius As Double
  PlotCenterX As Double
  PlotCenterY As Double
  PlotCenterZ As Double
  DigString As String * 8
  NtreesInside As Integer
  NClusInside As Integer
  Nmatched As Integer
  Nomission As Integer
  NCommission As Integer
  MatchP As Double
  Mrate As Double
  ZbiasMatched As Double
  ZbiasAll As Double
  Xbiasmatched As Double
  Ybiasmatched As Double
End Type
  
  
Public OV As OutputVector

Rem Each item of MeasuredPoint is 2+4+4+4+4+1+1+2+(1*6)+(4*6)+(4*6) = 76 Bytes (assuming max 6 images)
Rem So a Table of 500 items consumes 38000 Bytes

Public pointset3D(1 To 300, 1 To 4) As Double


Type Orientation_information
   ImageCode As Long              ' Long integer for each image
   FileName As String * 60        ' Description of the image
   Color As Byte                  ' 1 for 24-bit color, 0 for 8-bit greyscale
   o_row As Double               ' origin of subimage
   o_col As Double               ' origin of subimage
   Width As Integer               ' width in pixels of aerial photo (main image)
   Height As Integer              ' height in pixels of aerial photo
   sub_c_row As Integer           ' sub image pan center row
   sub_c_col As Integer           ' sub image pan center col
   sub_width As Double           ' width of subimage
   sub_height As Double          ' height of subimage
   C As Double                    ' Focal lenght of camera
   x_ps As Double                 ' PPA in FC-x
   y_ps As Double                 ' PPA in FC-y
   lambda As Double               ' Helmert scale factor
   alpha As Double                ' Helmert angle
   mean_x As Double               ' Mean of fiducial marks x's in PS-coords
   mean_y As Double
   X_mean As Double               ' Mean of fid's in IMA-coords
   Y_mean As Double
   a_ As Double                   ' Affine coefficients
   b_ As Double
   c_ As Double
   d_ As Double
   e_ As Double
   f_ As Double
   omega As Double                ' exterior orientation parameters
   phi As Double
   kappa As Double
   Xo As Double
   Yo As Double
   Zo As Double
   Sun_azimuth As Double
   Sun_elevation As Double
      StartOf_string As String * 60      ' Starting row for additional exposures
   Num_of_addit_expos As Long  ' Number of additional exposures
   AdditFileName(1 To 4) As String * 60    ' Image location
   AdditType(1 To 4) As Long         ' 0 = BW 8-byte,1 = RGB 8-byte,2 = BW 16-byte, 3 = RGBIR 16-byte
   AdditWidth(1 To 4) As Long        ' width in pixels of aerial photo (main image)
   AdditHeight(1 To 4) As Long       ' height in pixels of aerial photo
   EndOf_string As String * 60      ' Ending row for additional exposures

End Type

Type Window_information
   win_o_col As Integer          ' the col value(in main image coord system) of the origo in the picture-box
   win_o_row As Integer
   win_width As Integer          ' width of the picture-box
   win_height As Integer
   pan_x As Double               ' Zoom factor pan_x = pan_y , 1 = No Zooming
   pan_y As Double
End Type

Rem Each pixel in the RAW-image files consists of R-,G- and B-value. All scaled [0,255].
Type RGBtriplet
 r As Byte
 g As Byte
 B As Byte
End Type
Rem Each Each pixels in a Windows BMP-file consists of a B,G R-value (in this order)
Type RGBQ
 B As Byte
 g As Byte
 r As Byte
End Type

Rem win_info() -array stores pan-information
Public win_info(0 To MAXIMA - 1) As Window_information
Public cor_win_info(0 To MAXIMA - 1) As Window_information
Rem image_info() -array stores image information
' Public image_info(0 To MAXIMA - 1) As Orientation_information
Public image_info(0 To 700) As Orientation_information
Public cor_ima_info(0 To MAXIMA - 1) As Orientation_information
Rem 3D-rotation Matrices for (0 to MAXIMA-1) images
' Public A(1 To 3, 1 To 3, 0 To MAXIMA - 1)
Public A(1 To 3, 1 To 3, 0 To 700)

Rem Image-data (RGB) is stored for (0 to MAXIMA-1) images in a RGBtriplet array kuva() (kuva is Finnish for image)
' Public kuva(0 To maxcol, 0 To maxrow, 0 To MAXIMA - 1) As RGBtriplet
Public kuva() As RGBtriplet
Rem Stores index of last mousedown-event in a pictureBox (used with VLL)
Public LastImageClicked As Integer
Rem Once there's a solution (X,Y,Z) of space intersection, this VAR i set to True, and we may use the solution
Rem of the previous measured 3D-point as the startpoint (approximation) of the next iteration.
Public SolutionExists As Boolean

Public Imagesdisplayed(0 To 700) As Integer
Public Const AER_IMA = 1
Public Const CORR_IMA = 2

Rem Global variables to VARS that hold min and max indeces of pixel-locations that make up the picture-box's area.

Rem All picture Boxes are of equal size (an even number). win_w is width and win_h height.
Public Win_w As Integer
Public win_h As Integer

Rem CCOR is used for computing VLL crosscorrelation
Declare Function MYFUNC_TEMPMEANSTD Lib "c:\data\pascaldll.dll" (ByRef v1 As Byte, ByVal w As Long, ByVal h As Long, ByRef meanofimage As Double, ByRef STdOfImage As Double) As Double
Declare Function MYFUNC_CCORA Lib "c:\data\pascaldll.dll" (ByRef im1 As Byte, ByRef im2 As Byte, ByRef r As Double, ByRef meanA As Double, ByRef TempS As Long) As Double
Declare Function MYFUNC_CCORB Lib "c:\data\pascaldll.dll" (ByRef tempim1 As Byte, ByRef refim2 As Byte, ByRef r As Double, ByRef meantempim As Double, ByRef TempWidth As Long, ByRef TempHeight As Long) As Double


Declare Function MYFUNC_CCOR Lib "c:\data\pascaldll.dll" (ByRef v1 As Byte, ByVal n1 As Integer, ByVal m1 As Integer, ByRef v2 As Byte, ByVal n2 As Integer, ByVal m2 As Integer, ByVal ima1_x As Integer, ByVal ima1_y As Integer, ByVal ima2_x As Integer, ByVal ima2_y As Integer, ByVal w_temp As Integer, ByVal h_temp As Integer) As Double

Rem FLOOR is just a call to C-type floor() -function (which does not exist in VB)
Declare Function MYFUNC_FLOOR Lib "c:\data\pascaldll.dll" (ByVal v1 As Double) As Double
Declare Function MYFUNC_TEST Lib "c:\data\pascaldll.dll" (ByRef v1 As String) As Double
Declare Function MYFUNC_ASIN Lib "c:\data\pascaldll.dll" (ByVal v1 As Double) As Double
Declare Function MYFUNC_ACOS Lib "c:\data\pascaldll.dll" (ByVal v1 As Double) As Double
Declare Function MYFUNC_ATAN2 Lib "c:\data\pascaldll.dll" (ByVal v1 As Double, ByVal v1 As Double) As Double
Declare Function MYFUNC_ATAN Lib "c:\data\pascaldll.dll" (ByVal v1 As Double) As Double

Rem These functions may be used if one wants to send the DIB (BMP) device-independant bitmap directly (fast) to the picture-box
Declare Function MYFUNC_CREATEBMP Lib "c:\data\pascaldll.dll" (ByRef lfile As Byte, ByRef rfile As Byte, ByRef ofile As Byte, ByVal Index As Long, ByVal lstartcol As Long, ByVal rstartcol As Long, ByVal lstartrow As Long, ByVal rstartrow As Long, ByVal lendcol As Long, ByVal lendrow As Long, ByVal rendcol As Long, ByVal rendrow As Long, ByVal pan_x As Double, ByVal pan_y As Double, ByVal nwidth As Long, ByVal nheight As Long, ByVal stereoheight As Long, ByVal stereowidth As Long, ByVal colormodel As Long) As Double
Declare Function MYFUNC_CORIMA Lib "c:\data\pascaldll.dll" (ByRef v1 As Byte, ByVal na As Integer, ByVal ma As Integer, ByRef t1 As Byte, ByVal wt As Integer, ByVal ht As Integer, ByVal c_col As Integer, ByVal c_row As Integer, ByRef rk As Double) As Double
Declare Function MYFUNC_CREATENORMALIZEDIMAGE Lib "c:\data\pascaldll.dll" (ByVal C As Double, ByVal nwidth As Long, ByVal nheight As Long, ByVal Owidth As Long, ByVal Oheight As Long, ByRef filename_in As Byte, ByRef filename_out As Byte, ByRef pixelsize As Double, ByRef af As Double, ByRef r As Double) As Double
Declare Function MYFUNC_READBINARYFILES Lib "c:\data\pascaldll.dll" (ByRef v1 As Byte, ByVal lstartrow As Long, ByVal rstartrow As Long, ByVal lstartcol As Long, ByVal rstartcol As Long, ByVal lendcol As Long, ByVal nwidth As Long) As Double
Declare Function MYFUNC_MAKE3DPOINTS Lib "c:\data\pascaldll.dll" (ByVal origoX As Double, ByVal origoY As Double, ByVal origoZ As Double, ByVal gridXextent As Double, ByVal gridYextent As Double, ByVal gridZextent As Double, ByVal gridXtess As Double, ByVal gridYtess As Double, ByVal gridZtess As Double, ByVal XYrotangle As Double, ByRef FileName As Byte, ByVal strlen As Long) As Long
Public Declare Function GetObject Lib "gdi32" Alias "GetObjectA" (ByVal hObject As Long, ByVal ncount As Long, lpObject As Any) As Long
Public Declare Function SetBitmapBits Lib "gdi32" (ByVal hBitmap As Long, ByVal dwCount As Long, lpBits As Any) As Long
Public Declare Function GetBitmapBits Lib "gdi32" (ByVal hBitmap As Long, ByVal dwCount As Long, lpBits As Any) As Long
Public Declare Function Beep Lib "kernel32" (ByVal dwFreq As Long, ByVal dwDuration As Long) As Long
Public Declare Function MYFUNC_COMPUTEBKPROJ Lib "c:\data\pascaldll.dll" (ByVal Nimage As Long, ByVal Npoints As Long, ByRef af As Double, ByRef ori As Double, ByRef cimainfo As Double, ByRef atau As Single, ByRef SSpace As Double, ByRef cor_arr As Double) As Long
Declare Function MYFUNC_MATLABCALLNOMSGS Lib "c:\data\pascaldll.dll" (ByVal v1 As Double) As Double

Public Declare Function Polygon Lib "gdi32" (ByVal hdc As Long, lpPoint As POINTAPI, ByVal ncount As Long) As Long
Public Declare Function LineTo Lib "gdi32" (ByVal hdc As Long, ByVal X As Long, ByVal Y As Long) As Long
Public Declare Function MoveToEx Lib "gdi32" (ByVal hdc As Long, ByVal X As Long, ByVal Y As Long, lpPoint As POINTAPI) As Long
Public Declare Function rectangle Lib "gdi32" Alias "Rectangle" (ByVal hdc As Long, ByVal X1 As Long, ByVal Y1 As Long, ByVal X2 As Long, ByVal Y2 As Long) As Long
Public Declare Function Ellipsi Lib "gdi32" Alias "Ellipse" (ByVal hdc As Long, ByVal X1 As Long, ByVal Y1 As Long, ByVal X2 As Long, ByVal Y2 As Long) As Long
Public Declare Function TextOut Lib "gdi32" Alias "TextOutA" (ByVal hdc As Long, ByVal X As Long, ByVal Y As Long, ByVal lpstring As String, ByVal ncount As Long) As Long
Public Type Size
        cx As Long
        cy As Long
End Type
Public Declare Function GetWindowExtEx Lib "gdi32" (ByVal hdc As Long, lpSize As Size) As Long


Const BI_RGB = 0&
Rem RGBQUAD type is needed for DIB
Type RGBQUAD1
        B As Byte
        g As Byte
        r As Byte
        N As Byte
End Type
Type RGBQUAD
        rgbblue As Byte
        rgbGreen As Byte
        rgbred As Byte
        rgbreserved As Byte
End Type

Rem The 14 byte file header that starts the BMP-file
Type BITMAPFILEHEADER
        bfType As Integer         ' always BM
        bfSize As Long            ' size of file in bytes
        bfReserved1 As Integer    ' always zero
        bfReserved2 As Integer    ' always zero
        bfOffBits As Long         ' the offset to first pixel (55 for true color images with no color-table)
End Type

Rem The 40-Byte BITMAPINFOHEADER follows the file header
Type BITMAPINFOHEADER
        biSize As Long            ' 40 for 40 Bytes (this could be sizeof(BITMAPINFOHEADER), but there's no sizeof() in VB)
        biWidth As Long           ' image width in pixels
        biHeight As Long          ' image height in pixels, if a negative number, the image is stored with low-left corner first
        biPlanes As Integer       ' 1
        biBitCount As Integer     ' 24 for true-color
        biCompression As Long     ' 0        "
        biSizeImage As Long       ' 0        "
        biXPelsPerMeter As Long   ' 0        "
        biYPelsPerMeter As Long   ' 0        "
        biClrUsed As Long         ' 0        "
        biClrImportant As Long    ' 0        "
End Type

Type BITMAP
 bmType As Long
 bmWidth As Long
 bmHeight As Long
 bmWidthBytes As Long
 bmPlanes As Integer
 bmBitsPixel As Integer
 bmBits As Long
End Type
 
Type BITMAPINFO
        bmiHeader As BITMAPINFOHEADER
        bmiColors As RGBQUAD
End Type

Type ALS_Point
 GPS_time As Double  ' 8
 Xl As Double        ' 8
 Yl As Double        ' 8
 Zl As Single        ' 4
 omega As Single     ' 4
 phi As Single       ' 4
 kappa As Single     ' 4
 Scan_angle As Single ' 4
 Xf As Double         ' 8
 Yf As Double         ' 8
 Zf As Single         ' 4
 Intf As Single       ' 4
 Rangef As Single     ' 4
 Xs As Double         ' 8
 Ys As Double         ' 8
 Zs As Single         ' 4
 Ints As Single       ' 4
 Rangel As Single     ' 4
End Type ' 100 bytes per pulse

Public LP As ALS_Point
Public LP_Arr() As ALS_Point

Type LP_Espoo
 Echo As Byte
 X As Double
 Y As Double
 z As Single
End Type

Public LPE_AR() As LP_Espoo

Type Ltree
 X As Double
 Y As Double
 z As Single
 h As Single
 NLP As Long
 cf As Double       ' coefficient of crown model
 pw As Double       ' coefficient of crown model
 cont2 As Double    ' coefficient of crown model
End Type

Public Trees() As Ltree


Rem *****************************************
Rem Declarations for height models (DTM / CHM)

Type height_model
 ncols As Double
 nrows As Double
 xllcorner As Double
 yllcorner As Double
 CellSize As Double
 nodata_value As Double
 filepath As String * 80
End Type

Public Zmodel As height_model
Public ZfXY() As Single
Public CHMmodel As height_model
Public CHMZfXY() As Single

Public RasterModelReady As Boolean

Rem TINs

Public TINModelReady As Boolean
Public TINCHMModelReady As Boolean
Rem Declarations for measuring crown diameter
Public Pixel_width As Double
Public FirstPointMeasured As Boolean
Public SecondPointMeasured As Boolean
Public Firstpixel As Point2D
Public SecondPixel As Point2D
Public crown_width As Double

Rem DATA type for checking points in massive sets of stereo-pairs
Type Ima_obs
 Point As Long
 Im As Long
 X As Double
 Y As Double
End Type
Public IM_OBS() As Ima_obs
Type Ground_OBS
 Point As Long
 Type As Long
 X As Double
 Y As Double
 z As Double
End Type
Public GND_OBS() As Ground_OBS
 
Rem Declarations for TIN-dem
Type D3point
 Id As Long
 X As Double
 Y As Double
 z As Double
End Type
Public Nodes() As D3point

Type Point
 X As Double
 Y As Double
End Type

Type triangle
 i As D3point
 j As D3point
 k As D3point
 poly(0 To 2) As Point
 A As Double
 B As Double
 C As Double
 D As Double  ' Surface equation Ax + By + Cz = D
 xmin As Double
 Ymin As Double
End Type

Public N_CHM_tri As Long
Public N_CHM_Nodes As Long
Public CHM_Nodes() As D3point
Public CHM_Tri() As triangle
Public East_IND() As Long
Public CHM_TRI_east_IND() As Long
Public N_Nodes As Long
Public N_tri As Long ' number of triangles in the Tri-NODES -TIN model
Public MinX1 As Double, MinY1 As Double, MaxY1 As Double, MaxX1 ' 1-Hectare addresses with respect to these


'Public Nodes() As D3point
Public Tri() As triangle


Public KuvioRaj() As Point


Rem Declarations needed in computing ZONAL_calc images

Type ROwColIdTrip
 col As Long
 row As Long
 Id As Long
End Type

Type SaplingStandObs
 X As Double
 Y As Double
 z As Double
 Id As Long
 Radius As Double
End Type

Rem Laserpoints (ALTM 2033 data)
Type Lpoint
 Type As Byte
 X As Double
 Y As Double
 z As Single
 Ints As Single
End Type

Rem 2006 ALTM 3100 data
Type LidarRecord
 GPSTime As Double
 PulseCount As Byte
 Returns(1 To 4) As point3D
 Intensity(1 To 4) As Integer
 Range(1 To 4) As Double
 angle As Double
 Roll As Double
 Pitch As Double
 Heading As Double
 Poslidar As point3D
 StripNum As Integer
 SyncBit As Byte
 Res1 As Byte
 Res2 As Byte
 Res3 As Byte
End Type
Public LidR() As LidarRecord

Public NCounter As Long

Public Const Inside = 1
Public Const OUTSIDE = 0

Public Sub Mark_And_Plot(RMSE As Double, Hlimit As Double, Hlowest As Double, maxA As Long, Ntrees As Long, LPE_SP() As Byte, Hftau() As Single, Hfind() As Long, sade() As Single, xydist() As Single, zdist() As Single, Trees() As Ltree, LPTreeHit() As Boolean)


Rem Now should we mark again those lidar points that have been used
Rem

Dim Colorpix As Long
Dim p_x As Double, p_y As Double, h As Double

Select Case LPE_SP(Ntrees)
 Case 1
  Colorpix = RGB(255, 0, 0)
 Case 2
  Colorpix = RGB(0, 255, 0)
 Case 3
  Colorpix = RGB(0, 0, 255)
 Case 4
  Colorpix = RGB(255, 255, 0)
 Case 0
  Colorpix = RGB(128, 128, 128)
End Select

aa = Trees(Ntrees).NLP ' = 0  ' why is zero?
h = Trees(Ntrees).h

Select Case LPE_SP(Ntrees)
 Case 0, 1, 3, 4
 Hlow = 0.7 * Trees(Ntrees).h
 Case 2
 Hlow = 0.5 * Trees(Ntrees).h
End Select

'If Hlow < Hlimit Then
' Hlow = 0.7 * Hlimit
'End If
'Aw = Hlimit
 For i = UBound(Hftau) To 1 Step -1
   If Hftau(Hfind(i)) > (Hlow) And LPTreeHit(Hfind(i)) = False Then
    xydist(maxA) = Sqr((LPE_AR(Hfind(i)).X - Trees(Ntrees).X) ^ 2 + (LPE_AR(Hfind(i)).Y - Trees(Ntrees).Y) ^ 2)
    zdist(maxA) = (Trees(Ntrees).h - Hftau(Hfind(i)))
    Rem Crown radius, given height and species and current three parameters cf, pw, cont2
     'sade(maxA) = ((Trees(Ntrees).cf * h * (zdist(maxA) / h) ^ Trees(Ntrees).pw + Trees(1).cont2)) * 1.15
     If zdist(maxA) < 0 Then zdist(maxA) = 0.01
     sade(maxA) = (Trees(Ntrees).cf * h * Sin(4 * ((zdist(maxA) + 0.001) / h)) ^ Trees(Ntrees).pw + Trees(Ntrees).cont2) * 1.2
     Rem Increase sade by 20 %
     If sade(maxA) > xydist(maxA) Then  ' Or xydist(maxA) < (h / 20) Then
     ' Call r_transform_ground_to_pixel(1, LPE_AR(Hfind(i)).X, LPE_AR(Hfind(i)).Y, LPE_AR(Hfind(i)).z, p_x, p_y)
     ' Form1.Picture1(1).DrawWidth = 2
     ' Form1.Picture1(1).PSet ((p_x - (image_info(1).o_col + win_info(1).win_o_col)) * win_info(1).pan_x - 1, ((image_info(1).Height - 1) - p_y - (image_info(1).o_row + win_info(1).win_o_row)) * win_info(1).pan_y - 1), Colorpix
      DoEvents
      Trees(Ntrees).NLP = Trees(Ntrees).NLP + 1
      LPTreeHit(Hfind(i)) = True
      End If
     End If
TakeNext:  Next i

'Form1.Label10.Caption = "Done Fitting, iterations: " & Nite & " RMSE " & Format$(RMSE, "0.00 m")
Dim Xt As Double, Zt As Double, Yt As Double, Xf As Double, Yf As Double
Dim Zdiff As Double, p_x_beg As Double, p_y_beg As Double, p_x_end As Double, p_y_end As Double
'Colorpix = RGB(255, 255, 255)

'If RMSE > 0.7 Then GoTo OhitaPiirto

If RMSE > 0.5 Then Colorpix = RGB(255, 255, 255)
 For i = 0 To NumOfImages - 1
   Call r_transform_ground_to_pixel(i, Trees(Ntrees).X, Trees(Ntrees).Y, Trees(Ntrees).z, p_x, p_y)
   Form1.Picture1(i).DrawWidth = 5
   Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1), Colorpix
 Next i
  
 'If RMSE > 1 Then GoTo OhitaPiirto
' GoTo OhitaPiirto
Select Case LPE_SP(Ntrees)
 Case 1
  Colorpix = RGB(255, 0, 0)
 Case 2
  Colorpix = RGB(0, 255, 0)
 Case 3
  Colorpix = RGB(0, 100, 255)
 Case 4
  Colorpix = RGB(255, 255, 0)
 Case 0
  Colorpix = RGB(128, 128, 128)
End Select

'apu = MIN(Hlowest, 0.3 * Trees(Ntrees).h)
apu = 0.4 * Trees(Ntrees).h
For Zdiff = 0 To apu Step 1.5
  Zt = Trees(Ntrees).z - Zdiff
  Radius = (Trees(Ntrees).cf * Trees(Ntrees).h * Sin(4 * ((Zdiff) / h)) ^ Trees(Ntrees).pw + Trees(Ntrees).cont2) * 1
  mx = 0
  Shift = 3.14 / 10
  For phi = -3.14 To 3.15 Step 3.14 / 8
   Xt = Radius * Cos(phi) + Trees(Ntrees).X
   Yt = Radius * Sin(phi) + Trees(Ntrees).Y
   Xf = Radius * Cos(phi + Shift) + Trees(Ntrees).X
   Yf = Radius * Sin(phi + Shift) + Trees(Ntrees).Y
   Rem Using Line-method draw the crown
    For i = 0 To NumOfImages - 1
     'DistA = ((image_info(i).Xo - Xt) ^ 2 + (image_info(i).Yo - Yt) ^ 2)
     'DistB = ((image_info(i).Xo - Trees(Ntrees).X) ^ 2 + (image_info(i).Yo - Trees(Ntrees).Y) ^ 2)
     'If (DistA + 1) > DistB Then
      Call r_transform_ground_to_pixel(i, Xt, Yt, Zt, p_x_beg, p_y_beg)
      Call r_transform_ground_to_pixel(i, Xf, Yf, Zt, p_x_end, p_y_end)
      Form1.Picture1(i).DrawWidth = 1
      Form1.Picture1(i).ForeColor = RGB(255, 255, 255)
      Form1.Picture1(i).Line ((p_x_beg - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, ((image_info(i).Height - 1) - p_y_beg - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1)-((p_x_end - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, ((image_info(i).Height - 1) - p_y_end - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1), Colorpix
     'End If
    Next i
   Next phi
  Next Zdiff
'  NCounter = NCounter + 1
'  SavePicture Form1.Picture1(1).Image, "c:\temp\AVI\" & Format$(NCounter, "0000") & ".bmp"
OhitaPiirto:
DoEvents


End Sub
Public Sub Find_Next_Tree(RMSE As Double, maxA As Long, Hlowest As Double, Ntrees As Long, Hlimit As Double, Hftau() As Single, Hfind() As Long, Trees() As Ltree, LPE_AR() As LP_Espoo, LPTreeHit() As Boolean, xydist() As Single, zdist() As Single, sade() As Single, Inpoints() As point3D, Z_gnd() As Single, LPE_SP() As Byte, cf() As Double, pw() As Double, cont2() As Double, A_MAT() As Double, LA() As Double)
 
Rem (c) Ilkka Korpela May 2006
Rem
Rem This subroutine operates with
Rem  (1) a set of already mapped tree crowns in Trees() - struct
Rem  (2) a set of lidar points, some of which are already marked as "used" in LPTreeHit() -array
Rem
Rem  Array LPE_AR() contains the raw lidar data
Rem  Array Z_gnd() has ground elevations
Rem  Indec array Hfind() has the sorted order of ground normalized heights in Hftau()

Dim Xsum As Double, Ysum As Double, h As Double
Dim i As Long, k As Long, j As Long
Dim red As Double, green As Double, blue As Double, testpoint As point3D
Dim p_x As Double, p_y As Double

Rem The search for treetop positions is limited above Hlimit

 For i = UBound(Hftau) - 1 To 1 Step -1
  If Hftau(Hfind(i)) > Hlimit And LPTreeHit(Hfind(i)) = False Then
      Rem We found the next highest, non-used point, it makes a tree top candidate
      Ntrees = Ntrees + 1
      Trees(Ntrees).X = LPE_AR(Hfind(i)).X
      Trees(Ntrees).Y = LPE_AR(Hfind(i)).Y
      Trees(Ntrees).z = LPE_AR(Hfind(i)).z
      Trees(Ntrees).h = LPE_AR(Hfind(i)).z - Z_gnd(Hfind(i))
      h = Trees(Ntrees).h
      testpoint.X = Trees(Ntrees).X
      testpoint.Y = Trees(Ntrees).Y
      testpoint.z = Trees(Ntrees).z
      LPTreeHit(Hfind(i)) = True
'      Call r_transform_ground_to_pixel(1, Trees(Ntrees).X, Trees(Ntrees).Y, Trees(Ntrees).z, p_x, p_y)
'      Form1.Picture1(1).DrawWidth = 5
'      Form1.Picture1(1).PSet ((p_x - (image_info(1).o_col + win_info(1).win_o_col)) * win_info(1).pan_x - 1, ((image_info(1).Height - 1) - p_y - (image_info(1).o_row + win_info(1).win_o_row)) * win_info(1).pan_y - 1), RGB(255, 255, 255)
 '     Form1.Picture1(1).DrawWidth = 1
      DoEvents
      Rem Solve species from image j, assign basic crown parameters
      j = 1
      Call RGB_Vector_For_Point(red, green, blue, testpoint, CLng(j))
      LPE_SP(Ntrees) = FindSpecies(CDbl(red), CDbl(green), CDbl(blue), CDbl(1))
      Trees(Ntrees).cf = cf(LPE_SP(Ntrees))
      Trees(Ntrees).pw = pw(LPE_SP(Ntrees))
      Trees(Ntrees).cont2 = cont2(LPE_SP(Ntrees))
      k = i
      GoTo Found
    End If
  Next i
Found:
Rem The search for crown points is restricted. Assume relative crown lenght according to species

Select Case LPE_SP(Ntrees)
 Case 0, 1, 3, 4
 Hlow = 0.7 * Trees(Ntrees).h
 Case 2
 Hlow = 0.5 * Trees(Ntrees).h
End Select

Hlowest = 0

Rem Find lidar points that possibly make the crown points
 For i = k To 1 Step -1
  If Hftau(Hfind(i)) > Hlow And LPTreeHit(Hfind(i)) = False Then
     Rem Distances to candidate
     xydist(maxA) = Sqr((LPE_AR(Hfind(i)).X - Trees(Ntrees).X) ^ 2 + (LPE_AR(Hfind(i)).Y - Trees(Ntrees).Y) ^ 2)
     zdist(maxA) = (Trees(Ntrees).h - Hftau(Hfind(i)))
     Rem Crown radius, given height and species and current three parameters cf, pw, cont2
     sade(maxA) = Trees(Ntrees).cf * Trees(Ntrees).h * Sin(4 * (zdist(maxA) / Trees(Ntrees).h)) ^ Trees(Ntrees).pw + Trees(Ntrees).cont2
     
     If sade(maxA) > xydist(maxA) And zdist(maxA) < (h - Hlow) Then
      Rem Distance to top (downwards)
       LPTreeHit(Hfind(i)) = True
       Trees(Ntrees).NLP = Trees(Ntrees).NLP + 1
       xydist(Trees(Ntrees).NLP) = xydist(maxA)
      ' Exit Sub
       zdist(Trees(Ntrees).NLP) = zdist(maxA)
       If zdist(maxA) > Hlowest Then Hlowest = zdist(maxA)
       Rem Collect observations: Points XYZ
       Inpoints(Trees(Ntrees).NLP).X = LPE_AR(Hfind(i)).X
       Xsum = Xsum + LPE_AR(Hfind(i)).X
       Inpoints(Trees(Ntrees).NLP).Y = LPE_AR(Hfind(i)).Y
       Ysum = Ysum + LPE_AR(Hfind(i)).Y
       Inpoints(Trees(Ntrees).NLP).z = LPE_AR(Hfind(i)).z
     '  LPTreeHit(Hfind(i)) = True
       Call r_transform_ground_to_pixel(1, LPE_AR(Hfind(i)).X, LPE_AR(Hfind(i)).Y, LPE_AR(Hfind(i)).z, p_x, p_y)
       Form1.Picture1(1).DrawWidth = 1
       Form1.Picture1(1).PSet ((p_x - (image_info(1).o_col + win_info(1).win_o_col)) * win_info(1).pan_x - 1, ((image_info(1).Height - 1) - p_y - (image_info(1).o_row + win_info(1).win_o_row)) * win_info(1).pan_y - 1), RGB(255, 156, 60)
 
      End If
     End If
  Next i
 Close (6)
 
 
If Trees(Ntrees).NLP < (25) Then
 Rem It is not a tree, leave the routine
 Form1.Label10.Caption = "Fitting FAIL, too few points: " & Trees(Ntrees).NLP
 Ntrees = Ntrees - 1
 Exit Sub
End If
 
 Rem Check the distribution of the point cloud, omit if not a center
 
 Xsum = Xsum / CDbl((Trees(Ntrees).NLP))
 If Abs((Trees(Ntrees).X - Xsum)) < 0.6 Then Trees(Ntrees).X = Xsum
 Ysum = Ysum / CDbl((Trees(Ntrees).NLP))
 If Abs((Trees(Ntrees).Y - Ysum)) < 0.6 Then Trees(Ntrees).Y = Ysum
 
 NITE = 0
StartOfFitModel:
'Dim RMSE As Double
RMSE = 0

Call FitModel(RMSE, Hlowest, Ntrees, xydist(), sade(), zdist(), Trees(), LA(), A_MAT())
apu = MYFUNC_MATLABCALLNOMSGS(CDbl(1))

NITE = NITE + 1
'Exit Sub
Open "c:\data\corr_vectN.txt" For Input As 1
Input #1, Dcf
Input #1, Dpw
Input #1, Dcont2
Close (1)
If Abs(Dcf) < 0.001 And Abs(Dpw) < 0.001 And Abs(Dcont2) < 0.001 Then GoTo DoneFitting
If NITE > 20 Then GoTo DoneFitting
 
 Trees(Ntrees).cf = Trees(Ntrees).cf + Dcf
 Trees(Ntrees).pw = Trees(Ntrees).pw + Dpw
 Trees(Ntrees).cont2 = Trees(Ntrees).cont2 + Dcont2
 If Trees(Ntrees).pw < 0 Or RMSE > 3 Then
   Trees(Ntrees).pw = 0.7
   Trees(Ntrees).cont2 = 0.6
   Trees(Ntrees).cf = 0.15
 End If
GoTo StartOfFitModel

Rem Let's draw the crown in the images
DoneFitting:

If RMSE > 1.2 Then
 Rem It is not a tree
 Form1.Label10.Caption = "Fitting FAIL, iterations: " & NITE & " RMSE " & Format$(RMSE, "0.00 m") & " in " & Trees(Ntrees).NLP & " points"
 Ntrees = Ntrees - 1
 Exit Sub
End If

Rem Check if the h and CW(h-0.5h) make sense
 sade(maxA) = Trees(Ntrees).cf * Trees(Ntrees).h * Sin(4 * ((Trees(Ntrees).h * 0.5) / Trees(Ntrees).h)) ^ Trees(Ntrees).pw + Trees(Ntrees).cont2
 If sade(maxA) / Trees(Ntrees).h < 0.05 Or sade(maxA) / Trees(Ntrees).h > 0.15 Then
 Rem It is not a tree!
  Form1.Label10.Caption = "Fitting FAIL/CROWN RATIO: " & NITE & " RMSE " & Format$(RMSE, "0.00 m") & " in " & Trees(Ntrees).NLP & " points"
  Ntrees = Ntrees - 1
  Exit Sub
 End If

Form1.Label10.FontSize = 8
Form1.Label10.Caption = "Done Fitting, iterations: " & NITE & " RMSE " & Format$(RMSE, "0.00 m") & " in " & Trees(Ntrees).NLP & " points"

End Sub
Public Sub FitModel(RMSE As Double, Hlowest As Double, Ntrees As Long, xydist() As Single, sade() As Single, zdist() As Single, Trees() As Ltree, LA() As Double, A_MAT() As Double)
    
    Rem We have observations of crown radius in xydist() -array
    Rem         estimated         "      "   in sade() -array
    Dim h As Double
    Open "c:\data\A_txt.txt" For Output As 1
    Open "c:\data\L_txt.txt" For Output As 2
    Close (3)
    Open "c:\data\weights.txt" For Output As 3
     h = Trees(Ntrees).h
     RMSE = 0
     For i = 1 To Trees(Ntrees).NLP
       If zdist(i) > Hlowest Then Hlowest = zdist(i)
       j = j + 1
       sade(i) = Trees(Ntrees).cf * h * Sin(4 * (zdist(i) / h)) ^ Trees(Ntrees).pw + Trees(Ntrees).cont2
       LA(j) = xydist(i) - sade(i)
       RMSE = RMSE + LA(j) ^ 2
       A_MAT(j, 1) = h * Sin(4 * (zdist(i) / h)) ^ Trees(Ntrees).pw                    ' d(Sade)/d(cf)
       A_MAT(j, 2) = Trees(Ntrees).cf * h * Sin(4 * (zdist(i) / h)) ^ Trees(Ntrees).pw * Log(Sin(4 * ((zdist(i) + 0.01) / h)))
       A_MAT(j, 3) = 1
       For kX = 1 To 3
        Print #1, j & " " & kX & " " & A_MAT(j, kX)
       Next kX
       Print #2, LA(j)
       Print #3, Format$(1#, "0.00")  '  + 1 / ((zdist(i) / h + 0.01) ^ 2)

     Next i
    Open "c:\Data\matlab.txt" For Output As 4
    Print #4, j & " " & 3 & " " & 3 * j
    Close (4)
  
 Close (1)
 Close (2)
 Close (3)
RMSE = Sqr(RMSE / Trees(Ntrees).NLP)
Rem Solve the system of equations

End Sub



Public Function InsidePolygon(Polygon() As Point, ByVal N As Long, ByRef p As Point) As Long
  Dim counter As Long
  Dim i As Long
  Dim xinters As Double
  Dim p1 As Point, p2 As Point
  p1 = Polygon(0)
  For i = 1 To N Step 1
    p2 = Polygon(i Mod N)
    If (p.Y > MIN(p1.Y, p2.Y)) Then
      If (p.Y <= MAX(p1.Y, p2.Y)) Then
        If (p.X <= MAX(p1.X, p2.X)) Then
          If (p1.Y <> p2.Y) Then
            xinters = (p.Y - p1.Y) * (p2.X - p1.X) / (p2.Y - p1.Y) + p1.X
            If ((p1.X = p2.X) Or (p.X <= xinters)) Then
              counter = counter + 1
            End If
          End If
        End If
      End If
    End If
    p1 = p2
  Next
  If (counter Mod 2 = 0) Then
    InsidePolygon = OUTSIDE
  Else
    InsidePolygon = Inside
  End If
End Function

Private Function MIN(X, Y)
   If X < Y Then MIN = X Else MIN = Y
End Function

Private Function MAX(X, Y)
   If X > Y Then MAX = X Else MAX = Y
End Function

 
Public Sub Main()
Load Form1
Form1.Visible = True

BAPPPATH = "C:\data"

End Sub

Public Function getCHMheight(ByRef X As Double, Y As Double) As Double
' If CHMRasterModelReady = False Then
'  MsgBox ("Read CHM model first!")
'  getCHMheight = 0#
'  Exit Function
' End If
 
Dim row As Integer, col As Integer
 'GoTo skip
 Rem solve row, col values in array zfxy()

 col = CInt((X - CHMmodel.xllcorner) / CHMmodel.CellSize)
' Exit Function
 'row = (CHMmodel.nrows - 1) - CInt((y - CHMmodel.yllcorner) / CHMmodel.CellSize)
 row = -Int((CHMmodel.yllcorner - Y) / CHMmodel.CellSize)
 'row = -CInt((y - (Zmodel.yllcorner + Zmodel.CellSize * (Zmodel.nrows - 1))) / Zmodel.CellSize)
  
 ' Xcell = xllcorner + x * cellsize
 ' (YCell - ((yllcorner + (nrows - 1) * cellsize)))/cellsize =  y * cellsize
 
 If col < 0 Or col > CHMmodel.ncols - 1 Or row < 0 Or row > CHMmodel.nrows - 1 Then
'  MsgBox ("Point falls outside Zmodel!")
  getCHMheight = 0#
  Exit Function
 End If
 getCHMheight = CHMZfXY(col, row)

 Exit Function

End Function




Public Function getheight(ByRef X As Double, Y As Double) As Double

On Error GoTo Point_Falls_Out

 If RasterModelReady = False Then
  'MsgBox ("Read height model first!")
  getheight = 0#
  Exit Function
 End If
 
Dim row As Integer, col As Integer
 'GoTo skip
 Rem solve row, col values in array zfxy()

 col = CInt((X - Zmodel.xllcorner) / Zmodel.CellSize)
' Exit Function
'row = (Zmodel.nrows - 1) - CInt((Y - Zmodel.yllcorner) / Zmodel.CellSize)
 row = -Int((Zmodel.yllcorner - Y) / Zmodel.CellSize)
 ' row = -CInt((y - (Zmodel.yllcorner + Zmodel.CellSize * (Zmodel.nrows - 1))) / Zmodel.CellSize)
  
 ' Xcell = xllcorner + x * cellsize
 ' (YCell - ((yllcorner + (nrows - 1) * cellsize)))/cellsize =  y * cellsize
 
 If col < 0 Or col > Zmodel.ncols - 1 Or row < 0 Or row > Zmodel.nrows - 1 Then
'  MsgBox ("Point falls outside Zmodel!")
  getheight = 0#
  Exit Function
 End If
 getheight = ZfXY(col, row)

 Exit Function
skip:
 Rem *******************
 Rem UP
 
 ' Dim row As Integer, col As Integer
 Dim decimal_x As Double, decimal_y As Double
 
 Rem Truncate
 col = Int((X - Zmodel.xllcorner) / Zmodel.CellSize)
 'Exit Function
 row = Int((Zmodel.yllcorner - Y) / Zmodel.CellSize)
  
 Rem Calculate the middle point coords
 
 x_p = Zmodel.xllcorner + col * Zmodel.CellSize + Zmodel.CellSize / 2
 y_p = Zmodel.yllcorner - row * Zmodel.CellSize - Zmodel.CellSize / 2
 
 ReDim k(1 To 3, 1 To 3) As Double
 ReDim Zw(1 To 3, 1 To 3) As Double
 Rem Coefficients
 
 k(1, 1) = 0.25 * (Abs(x_p - (X - Zmodel.CellSize)) * Abs(y_p - (Y - Zmodel.CellSize))) / Zmodel.CellSize ^ 2
 k(1, 2) = 0.25 * (Abs(x_p - (X - Zmodel.CellSize)) * Abs(y_p - (Y + Zmodel.CellSize))) / Zmodel.CellSize ^ 2
 k(2, 1) = 0.25 * (Abs(x_p - (X + Zmodel.CellSize)) * Abs(y_p - (Y + Zmodel.CellSize))) / Zmodel.CellSize ^ 2
 k(2, 2) = 0.25 * (Abs(x_p - (X + Zmodel.CellSize)) * Abs(y_p - (Y - Zmodel.CellSize))) / Zmodel.CellSize ^ 2
 
 If Abs(1 - (k(1, 1) + k(1, 2) + k(2, 1) + k(2, 2))) > 0.02 Then MsgBox ("Summa?")
 col = Int(((X - Zmodel.CellSize) - Zmodel.xllcorner) / Zmodel.CellSize)
 row = Int((Zmodel.yllcorner - (Y - Zmodel.CellSize)) / Zmodel.CellSize)
 If col < 0 Or col > Zmodel.ncols - 1 Or row < 0 Or row > Zmodel.nrows - 1 Then GoTo Point_Falls_Out
 Zw(1, 1) = k(1, 1) * ZfXY(col, row)
 
 col = Int(((X - Zmodel.CellSize) - Zmodel.xllcorner) / Zmodel.CellSize)
 row = Int((Zmodel.yllcorner - (Y + Zmodel.CellSize)) / Zmodel.CellSize)
 If col < 0 Or col > Zmodel.ncols - 1 Or row < 0 Or row > Zmodel.nrows - 1 Then GoTo Point_Falls_Out
 Zw(1, 2) = k(1, 2) * ZfXY(col, row)
 
 col = Int(((X + Zmodel.CellSize) - Zmodel.xllcorner) / Zmodel.CellSize)
 row = Int((Zmodel.yllcorner - (Y + Zmodel.CellSize)) / Zmodel.CellSize)
 If col < 0 Or col > Zmodel.ncols - 1 Or row < 0 Or row > Zmodel.nrows - 1 Then GoTo Point_Falls_Out
 Zw(2, 1) = k(2, 1) * ZfXY(col, row)
 
 col = Int(((X + Zmodel.CellSize) - Zmodel.xllcorner) / Zmodel.CellSize)
 row = Int((Zmodel.yllcorner - (Y - Zmodel.CellSize)) / Zmodel.CellSize)
 If col < 0 Or col > Zmodel.ncols - 1 Or row < 0 Or row > Zmodel.nrows - 1 Then GoTo Point_Falls_Out
 Zw(2, 2) = k(2, 2) * ZfXY(col, row)
 
 
 
 col = Int((X - Zmodel.xllcorner) / Zmodel.CellSize)
 row = Int((Zmodel.yllcorner - Y) / Zmodel.CellSize)
 If col < 0 Or col > Zmodel.ncols - 1 Or row < 0 Or row > Zmodel.nrows - 1 Then GoTo Point_Falls_Out
 getheight = 1# / 3# * (Zw(1, 1) + Zw(1, 2) + Zw(2, 1) + Zw(2, 2)) + (2# / 3#) * ZfXY(col, row)
 
 
 
 Exit Function
 
Point_Falls_Out:
' MsgBox ("Error in retrieving elevation")
 getheight = 0#
 Exit Function
  
 ' row = -CInt((Y - (Zmodel.yllcorner + Zmodel.CellSize * (Zmodel.nrows - 1))) / Zmodel.CellSize)
 ' Xcell = xllcorner + x * cellsize
 ' (YCell - ((yllcorner + (nrows - 1) * cellsize)))/cellsize =  y * cellsize
 
 Rem *************************
 Rem Interpolation
 decimal_x = (((X - Zmodel.xllcorner) / Zmodel.CellSize) - col) * Zmodel.CellSize + Zmodel.CellSize / 2
 decimal_y = (((Y - Zmodel.yllcorner) / Zmodel.CellSize) - row) * Zmodel.CellSize + Zmodel.CellSize / 2
 
 
 
End Function
Public Sub Get_Angular_Info(ByVal i As Integer, ByRef NadAng As Double, ByRef SolAziDiff As Double, ByRef ObjToCamAzi As Double)

Dim Fii As Double, Theta As Double
Dim sunray As Vector3D
Dim cam_vec_x As Vector3D, cam_vec_y As Vector3D, cam_vec_z As Vector3D
Dim raystart As Vector3D
Dim p_x As Double, p_y As Double
Dim sun As Double

If SolutionExists = False Then Exit Sub

Dim cout As String


   Rem Fii e [-Pi,Pi], (270W,270W), 90E=0
   Fii = (90 - image_info(i).Sun_azimuth * TO_DEGREES) * TO_RADIANS
   Rem Theta: 0 = zenith, +Pi/2 = horizon
   Theta = (90 - image_info(i).Sun_elevation * TO_DEGREES) * TO_RADIANS
        
   Dim camera_x As Double, camera_y As Double
   Dim nadir_x As Double, nadir_y As Double
   Dim sun_x As Double, sun_y As Double
   Dim cam_x As Double, cam_y As Double
   Dim alpha As Double, z As Double, apu As Double
   
   Call r_transform_3D(i, X_sol, Y_sol, Z_sol, camera_x, camera_y)
   Call r_transform_3D(i, image_info(i).Xo, image_info(i).Yo, 100, nadir_x, nadir_y)
   
   Call create_vector(cam_vec_x, A(1, 1, i), A(2, 1, i), A(3, 1, i))
   Call create_vector(cam_vec_y, A(1, 2, i), A(2, 2, i), A(3, 2, i))
   Call create_vector(cam_vec_z, A(1, 3, i), A(2, 3, i), A(3, 3, i))

   z = -image_info(i).C
 
   cam_x = -camera_x
   cam_y = -camera_y
 
   Rem Direction from tree to image center
 
   alpha = MYFUNC_ATAN2(cam_y, cam_x)
   alpha = alpha + PI
   
   sunray.X = Sin(Theta) * Cos(Fii)
   sunray.Y = Sin(Theta) * Sin(Fii)
   sunray.z = Cos(Theta)
   
   raystart.X = (camera_x) * cam_vec_x.X + (camera_y) * cam_vec_y.X + z * cam_vec_z.X
   raystart.Y = (camera_x) * cam_vec_x.Y + (camera_y) * cam_vec_y.Y + z * cam_vec_z.Y
   raystart.z = (camera_x) * cam_vec_x.z + (camera_y) * cam_vec_y.z + z * cam_vec_z.z
   
   sun = MYFUNC_ATAN2(100 * sunray.Y, 100 * sunray.X)
   ' MsgBox ("Sun's azimuth: " & -(-PI / 2 + sun) * TO_DEGREES)
   cout = cout & Format$(-(-PI / 2 + sun) * TO_DEGREES, "#.0") & Chr$(9)
   ' MsgBox ("Sun's elevation: " & CStr(image_info(i).Sun_elevation * TO_DEGREES))
   'Cout = Cout & Format$(image_info(i).Sun_elevation * TO_DEGREES, "#.0") & Chr$(9)
   
   'sun = sun + PI
   
   Rem Alpha & beeta -Pi..Pi with respect to East
   alpha = MYFUNC_ATAN2(100 * raystart.Y, 100 * raystart.X)
   If alpha < -PI / 2 Then
    'MsgBox ("Camera ray's azimuth: " & -(-PI / 2 + alpha) * TO_DEGREES)
   apu = -(-PI / 2 + alpha) * TO_DEGREES
   If apu > 180 Then
    apu = apu - 180
    GoTo cout1
   End If
   If apu < 180 Then apu = apu + 180
cout1:
   cout = cout & Format$(apu, "#.0") & Chr$(9)
   End If
   If alpha >= -PI / 2 Then
   'MsgBox ("Camera ray's azimuth: " & -(-PI / 2 + alpha) * TO_DEGREES)
   apu = -(-PI / 2 + alpha) * TO_DEGREES
   If apu > 180 Then
    apu = apu - 180
    GoTo cout2
   End If
   If apu < 180 Then apu = apu + 180
cout2:
   cout = cout & Format$(apu, "#.0") & Chr$(9)
   End If
   
   Rem ************************
   Rem CAMERA RAY's azimuth APU
   ObjToCamAzi = apu
   Rem ************************
   Rem apu holds azimuth
   alpha = alpha + PI
   sun = -(-PI / 2 + sun)
   sun = sun * TO_DEGREES
   Rem Sun holds azimuth
   
   Call normalize(raystart)
   Call normalize(sunray)
   Call r_transform_3D(i, X_sol + 1000 * sunray.X, Y_sol + 1000 * sunray.Y, Z_sol + 1000 * sunray.z, sun_x, sun_y)

   Dim diff
      diff = sun - apu
       If diff < -180 Then
        diff = ((180 + diff) + 180)
        GoTo skip3
       End If
      If diff > 180 Then
       diff = (180 - (diff - 180)) * -1
      End If
skip3:
   
    
   Rem *****************************************
   Rem Diffrence in AZImuth -180,...,180 degrees
   SolAziDiff = diff
   
   'If (Abs(diff) >= 0 And Abs(diff) < 45) Then cout = cout & "B" & Chr$(9)
   'If (Abs(diff) >= 45 And Abs(diff) < 90) Then cout = cout & "BS" & Chr$(9)
   'If (Abs(diff) >= 90 And Abs(diff) < 135) Then cout = cout & "FS" & Chr$(9)
   'If (Abs(diff) >= 135 And Abs(diff) <= 180) Then cout = cout & "F" & Chr$(9)

   cam_vec_z.z = -cam_vec_z.z
   cam_vec_z.Y = -cam_vec_z.Y
   cam_vec_z.X = -cam_vec_z.X
   
   Call normalize(cam_vec_z)
   Rem ***************
   Rem Nadir angle (off-angle)

   NadAng = vector_angle(cam_vec_z, raystart) * TO_DEGREES



End Sub
' Linear function with various extent
Public Function Bartlett_Linear(ByVal X#, kernel_size As Double) As Double
    X = Abs(X)
    If X < kernel_size Then Bartlett_Linear = (1 - Abs(X) / kernel_size) / kernel_size
End Function
Public Sub Bilinear(k As Long, i As Long, w_temp As Long, h_temp As Long, SizeArr() As Byte, tx() As Byte, t1() As Byte)
Dim X&, Y&, X1&, Y1&, m&, N&, kX#, kY#, fX#, fY#
Dim IR#, iG#, iB#, r1#, r2#
If i = 2 And k = 6 Then
aa = 1
End If
srcWidth = w_temp
srcHeight = h_temp
dstWidth = SizeArr(k, 1)
dstHeight = SizeArr(k, 2)

 kX = (srcWidth - 1) / (dstWidth - 1)
 kY = (srcHeight - 1) / (dstHeight - 1)
    For Y = dstHeight - 1 To 0 Step -1
        fY = Y * kY ' Exact position (floating-point number)
        Y1 = Int(fY) ' Integer position (integer part of number)
        fY = fY - Y1 ' Fraction part of number (integer+fraction=exact)
        For X = 0 To dstWidth - 1
            fX = X * kX
            X1 = Int(fX)
            fX = fX - X1
            X1 = X1
            IR = 0: iG = 0: iB = 0
            ' Uses various kernel size
            kernel_size = 1
            For m = -kernel_size + 1 To kernel_size
                r1 = Bartlett_Linear(m - fY, CDbl(kernel_size))
                For N = -kernel_size + 1 To kernel_size
                    r2 = Bartlett_Linear(fX - N, CDbl(kernel_size))
                    iB = iB + tx(X1 + N, Y1 + m) * r1 * r2
                Next
            Next
           t1(X, Y) = iB
        Next
    Next

Open "c:\data\test_" & CStr(i) & "_" & CStr(k) & ".img" For Output As 10
Print #10, dstWidth, dstHeight, SizeArr(k, 3), SizeArr(k, 4)
Close (10)

Open "c:\data\test_" & CStr(i) & "_" & CStr(k) & ".raw" For Binary As 10
Put #10, , t1
Close (10)


End Sub

Public Sub Measure_tree_tops_By_Matching_sub(i_ref As Long)
Rem This routine uses a reference image, that it clicked for its tree top image position
Rem The reference image ray is discretized, and along it image correlation is calculated
Rem for templates of varying size (the routine calls resize templates ())
Rem
Rem Correlation is aggregated for points in a 2D space: Size x position (XYZ), and
Rem the maximum value is given. The reference image template has an associated crown
Rem width, which is re-scaled by the Size-parameter (multiplier)
Rem MUFUNC_corima() routine is modified and used in VB for the cross-correlation between image patch

Dim p_x As Double, p_y As Double
Dim Z_gnd As Double, Y_Gnd As Double, X_Gnd As Double

Rem First we solve the XY-position of the ray, given Z as Z_sol. We get X_gnd, Y_gnd_Z_gnd and the pixel-coordinates
Z_gnd = Z_sol

Call solve_3D_X_Y_from_Z_x_y(i_ref, Z_sol, CDbl(Form1.Text1(i_ref * 2).Text), CDbl(Form1.Text1(i_ref * 2 + 1).Text), X_Gnd, Y_Gnd)
Rem Pixel values
Call r_transform_ground_to_pixel(i_ref, X_Gnd, Y_Gnd, Z_gnd, p_x, p_y)

Rem Loop the Z-values, discretize the reference ray in i_ref
Dim Z_range As Double, Z_step As Double, X_test As Double, Y_test As Double, Z_test As Double
Dim l As Long
ReDim testpoints(1 To 200) As point3D
Z_range = 9
Z_step = 0.2

For Z_test = (Z_gnd - Z_range / 2) To (Z_gnd + Z_range / 2) Step Z_step
 l = l + 1
 Call solve_3D_X_Y_from_Z_x_y(i_ref, Z_test, CDbl(Form1.Text1(i_ref * 2).Text), CDbl(Form1.Text1(i_ref * 2 + 1).Text), X_test, Y_test)
 testpoints(l).X = X_test
 testpoints(l).Y = Y_test
 testpoints(l).z = Z_test
Next Z_test
Dim N_Of_TestPoints As Long
Dim R_Max As Double, R0_Max As Double
Dim ScaleMax As Long, MaxPoint As Long
R0_Max = -1
R_Max = -1
N_Of_TestPoints = l
  
Rem This array will hold the aggregated correlation, for the 11 cases
ReDim R_array(1 To N_Of_TestPoints, 1 To 11) As Double
ReDim R0_Array(1 To 11) As Double
Dim kX As Long, kY As Long
Rem The discretized reference ray is now stored as points in the point3D array testpoints()
Rem Which images are to be used
Dim w_temp As Long, h_temp As Long, c_col As Long, c_row As Long, address As Long
For i = 0 To NumOfImages - 1
  
  If Form1.Check1(i).Value = 1 Then
   ReDim Color_photo_Patch(0 To image_info(i).sub_width - 1, 0 To 159) As RGBtriplet
 '  ReDim Color_Photo_temp(0 To 99, 0 To 139) As RGBtriplet
 '  ReDim Color_Photo_Row(0 To 99, 0 To 139) As RGBtriplet
'  ReDim BW_Photo_Patch(0 To image_info(i).sub_width - 1, 0 To 139) As Byte
'  ReDim BW_Photo_temp(0 To 99, 0 To 139) As Byte
'  ReDim BW_Photo_Row(0 To 139) As Byte
  ReDim startRows(0 To MAXIMA - 1) As Long
  ReDim startCols(0 To MAXIMA - 1) As Long
 
 Rem Let's read the image patches for the first testpoint, store StartRows()
 Call r_transform_ground_to_pixel(i, X_Gnd, Y_Gnd, Z_gnd, p_x, p_y)
 address = CLng((image_info(i).sub_height - 1) - (p_y + 80)) * image_info(i).sub_width * 3 ' + p_x * 3
 startRows(i) = CLng((image_info(i).sub_height - 1) - (p_y + 80))
 startCols(i) = 0
 
 'Open image_info(i).FileName For Binary As 1
 Close (1)
 Open image_info(i).FileName For Binary As 1
 
 'Exit Sub
 Get #1, address + 1, Color_photo_Patch
 'Close (1)
 'Exit Sub
    
    For k = 1 To 11
     Open "c:\data\test_" & CStr(i) & "_" & CStr(k) & ".img" For Input As 12
'Exit Sub
     Input #12, w_temp, h_temp, c_col, c_row
     Close (12)
     ReDim template(0 To w_temp - 1, 0 To h_temp - 1) As Byte
     Open "c:\data\test_" & CStr(i) & "_" & CStr(k) & ".raw" For Binary As 12
     Get #12, , template
     Close (12)
     Rem Now we should take the values, from the image, let's take green
     ReDim Image(0 To w_temp - 1, 0 To h_temp - 1) As Byte
     For j = 1 To N_Of_TestPoints
      Call r_transform_ground_to_pixel(i, testpoints(j).X, testpoints(j).Y, testpoints(j).z, p_x, p_y)
     Rem we read a 100 by 100 image patch to be used with every template
     Rem Compute the row in the 100 row high image
     Rem p_x and p_y are in the hot_spot; we must shift
     apu = CLng((image_info(i).sub_height - 1) - p_y) - startRows(i)
     Rem apu is the shift in row direction, of the hot_spot row
     apu_y = apu - c_row
     apu_x = p_x - c_col
     For kX = 0 To w_temp - 1
       For kY = 0 To h_temp - 1
        ' Image(kX, kY) = CByte((CLng(Color_photo_Patch(apu_x + kX, apu_y + kY).r) + CLng(Color_photo_Patch(apu_x + kX, apu_y + kY).g) + CLng(Color_photo_Patch(apu_x + kX, apu_y + kY).B)) / 3#)
        Image(kX, kY) = CByte((Color_photo_Patch(apu_x + kX, apu_y + kY).g))
       ' Exit Sub
       Next kY
     Next kX
     Dim meaIm As Double, stdIm As Double
     apu = MYFUNC_TEMPMEANSTD(template(0, 0), CLng(UBound(template, 1)), CLng(UBound(template, 2)), meaIm, stdIm)
     Dim rx As Double
     Rem Compute cross-correlation between the scaled template and the image patch
     apu = MYFUNC_CCORB(template(0, 0), Image(0, 0), rx, meaIm, w_temp, h_temp)
     R_array(j, k) = R_array(j, k) + rx
    
     If R_array(j, k) > R_Max Then
      R_Max = R_array(j, k)
      ScaleMax = k
      MaxPoint = j
     End If
   Next j
  Next k
  End If
 
 Next i
 
 Rem *********** TRY IMAGE 0 AGAIN AROUND SOLUTION "
 i = 0
 Dim testpoint As point3D
 testpoint.X = testpoints(MaxPoint).X
' Exit Sub
 testpoint.Y = testpoints(MaxPoint).Y
 testpoint.z = testpoints(MaxPoint).z
 Call Match_template(CLng(i), testpoint)
 
' Call r_transform_ground_to_pixel(i, testpoints(MaxPoint).X, testpoints(MaxPoint).Y, testpoints(MaxPoint).z, p_x, p_y)
'Exit Sub
' address = CLng((image_info(i).sub_height - 1) - (p_y + 70)) * image_info(i).sub_width * 3 ' + p_x * 3
' startRows(i) = CLng((image_info(i).sub_height - 1) - (p_y + 70))
' startCols(i) = 0
 
' Open image_info(i).FileName For Binary As 1
' Get #1, address + 1, Color_photo_Patch
' Close (1)
' r_maxc = -1
'    For k = 1 To 11
'     Open "d:\test_" & CStr(i) & "_" & CStr(k) & ".img" For Input As 12
'     Input #12, w_temp, h_temp, c_col, c_row
'     Close (12)
'     ReDim template(0 To w_temp - 1, 0 To h_temp - 1) As Byte
'     Open "d:\test_" & CStr(i) & "_" & CStr(k) & ".raw" For Binary As 12
'     Get #12, , template
'     Close (12)
'     ReDim Image(0 To w_temp - 1, 0 To h_temp - 1) As Byte
'     For d_p_x = -2 To 2
'      For d_p_y = -2 To 2
'
'      apu = CLng((image_info(i).sub_height - 1) - (p_y + d_p_y)) - startRows(i)
'      apu_y = apu - c_row
'      apu_x = (p_x + d_p_x) - c_col
'      For kX = 0 To w_temp - 1
'        For kY = 0 To h_temp - 1
'         Image(kX, kY) = CByte((Color_photo_Patch(apu_x + kX, apu_y + kY).g))
'        Next kY
 '     Next kX
'      apu = MYFUNC_TEMPMEANSTD(template(0, 0), CLng(UBound(template, 1)), CLng(UBound(template, 2)), meaIm, stdIm)
'      apu = MYFUNC_CCORB(template(0, 0), Image(0, 0), rx, meaIm, w_temp, h_temp)
'      If rx > r_maxc Then
'        r_maxc = rx
'        ScaleMaxi = k
'       End If
'   Next d_p_y
'  Next d_p_x
' Next k
 
 '  For i = 0 To NumOfImages - 1
'    Form1.Picture1(i).Cls
 '   Call r_transform_ground_to_pixel(i, testpoints(MaxPoint).X, testpoints(MaxPoint).Y, testpoints(MaxPoint).z, p_x, p_y)
    'Exit Sub
 '   Form1.Picture1(i).DrawWidth = 3
 '   Form1.Picture1(i).FontSize = 8
 '   Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y), RGB(255, 255, 0)
    'If i = 1 Then
    'Form1.Picture1(i).CurrentX = Form1.Picture1(i).CurrentX + 5
    'Form1.Picture1(i).CurrentY = Form1.Picture1(i).CurrentY + 5
    'Form1.Picture1(i).Print Format$(R_Max / CDbl(NumOfImages + 1), "0.00")
    'End If
 '  Next i

Form1.Label10.Caption = "Height: " & Format$(testpoints(MaxPoint).z - getheight(testpoints(MaxPoint).X, testpoints(MaxPoint).Y), "0.00 m") & " r = " & Format$(R_Max / CDbl(NumOfImages + 1), "0.00")
    

Rem draw a circle
'Ulim = 11
'ReDim ScalingArr(1 To Ulim) As Double
'ScalingArr(1) = 0.4: ScalingArr(2) = 0.5: ScalingArr(3) = 0.6
'ScalingArr(4) = 0.7: ScalingArr(5) = 0.75: ScalingArr(6) = 0.8
'ScalingArr(7) = 0.85: ScalingArr(8) = 0.9: ScalingArr(9) = 1#: ScalingArr(10) = 1.1
'ScalingArr(11) = 1.2:

'i = 0
'Form1.Picture1(i).DrawWidth = 2
'Dim Xt As Double, Yt As Double, Xf As Double, Yf As Double
'Dim Shift As Double, p_x_beg  As Double, p_y_beg As Double, p_x_end As Double, p_y_end As Double
'Shift = 3.14 / 10
'Radius = CDbl(Form1.Text5.Text) * ScalingArr(R0_ind)
'Radius = CDbl(Form1.Text5.Text) * ScalingArr(ScaleMaxi)
'Exit Sub
'  For phi = -3.14 To 3.15 Step 3.14 / 8
'   Colorpix = RGB(120, 255, 0)   '  255 - Zdiff * 10)
'   Xt = Radius * Cos(phi) + testpoints(MaxPoint).X
'   Yt = Radius * Sin(phi) + testpoints(MaxPoint).Y
'   Xf = Radius * Cos(phi + Shift) + testpoints(MaxPoint).X
'   Yf = Radius * Sin(phi + Shift) + testpoints(MaxPoint).Y
'   Zt = testpoints(MaxPoint).z - 1
'   Call r_transform_ground_to_pixel(i, Xt, Yt, Zt, p_x_beg, p_y_beg)
'   Call r_transform_ground_to_pixel(i, Xf, Yf, Zt, p_x_end, p_y_end)
'   Form1.Picture1(i).Line ((p_x_beg - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, ((image_info(i).Height - 1) - p_y_beg - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1)-((p_x_end - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, ((image_info(i).Height - 1) - p_y_end - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1), Colorpix
'Next phi
 
 X_sol = testpoints(MaxPoint).X
 Y_sol = testpoints(MaxPoint).Y
 Z_sol = testpoints(MaxPoint).z
 Form1.Label4(0).Caption = Format$(X_sol, "0.00")
 Form1.Label4(1).Caption = Format$(Y_sol, "0.00")
 Form1.Label4(2).Caption = Format$(Z_sol, "0.00")
 DoEvents

 Rem ********'species
 
For Index = 0 To NumOfImages - 1
Dim RedVal As Integer, GreenVal As Integer, BlueVal As Integer, IrVal As Integer
If image_info(Index).Num_of_addit_expos = 2 Then
additnum = 2
Call r_transform_ground_to_pixel(Index, testpoints(MaxPoint).X, testpoints(MaxPoint).Y, testpoints(MaxPoint).z, p_x, p_y)
p_y = (image_info(Index).Height - 1) - p_y
p_x = Int(p_x * (image_info(Index).AdditWidth(additnum) / image_info(Index).Width))
p_y = Int(p_y * (image_info(Index).AdditHeight(additnum) / image_info(Index).Height))
 If image_info(Index).AdditType(additnum) = 3 Then
'    Open image_info(Index).AdditFileName(additnum) For Binary As 20
      'Get #20, 1 + ((p_y) * 4) * image_info(Index).AdditWidth(additnum) * 2 + (p_x) * 4 * 2, RedVal
'Close (20)
'Exit Sub
'      Get #20, 1 + ((p_y) * 4) * image_info(Index).AdditWidth(additnum) * 2 + (p_x) * 4 * 2 + 2, GreenVal
'      Get #20, 1 + ((p_y) * 4) * image_info(Index).AdditWidth(additnum) * 2 + (p_x) * 4 * 2 + 4, BlueVal
'      Get #20, 1 + ((p_y) * 4) * image_info(Index).AdditWidth(additnum) * 2 + (p_x) * 4 * 2 + 6, IrVal
    Close (20)
    sp = FindSpecies(CDbl(RedVal), CDbl(GreenVal), CDbl(BlueVal), CDbl(IrVal))
End If
End If
Next Index

Form1.Label10.Caption = sp & " H: " & Format$(testpoints(MaxPoint).z - getheight(testpoints(MaxPoint).X, testpoints(MaxPoint).Y), "0.00 m") & " r = " & Format$(R_Max / CDbl(NumOfImages + 1), "0.00")

Rem ************* species
'Measurement.Corr_Match = R_Max
 
'Open "c:\temp\mk_trees.txt" For Input As 6
 Close (6)
' Open "c:\temp\mk_trees.txt" For Input As 6
' Do Until EOF(6)
' Input #6, Xm, Ym, z_butt, TreeNum, TreeSpecies, Status, TreeHeight, TreeDiam
'  If Sqr((X_sol - Xm) ^ 2 + (Y_sol - Ym) ^ 2) < 1 Then
'  Rem The tree is found
'  GoTo foundTRee
'  End If
' Loop
 
 Rem If here, then no tree was found
 Close (6)
 Exit Sub
foundTRee:
Close (6)

'Open "d:\ko3_koe.txt" For Append As 12
'Print #12, X_sol, Y_sol, Z_sol, Z_sol - getheight(X_sol, Y_sol), Radius * 2, Xm, Ym, TreeNum, TreeHeight, TreeDiam
'Close (12)
 
'Open "d:\r_s.txt" For Append As 10
'For j = 1 To N_Of_TestPoints
' For k = 1 To 11
' Print #10, R_array(j, k) & "," & j & "," & k & "," & testpoints(j).X & "," & testpoints(j).Y & "," & testpoints(j).z
' Next k
'Next j
'Close (10)
 

End Sub
Public Sub Match_template(i As Long, testpoint As point3D)
'Exit Sub
ReDim Color_photo_Patch(0 To image_info(i).sub_width - 1, 0 To 139) As RGBtriplet

Rem *********** TRY IMAGE 0 AGAIN AROUND SOLUTION "

Dim p_x As Double, p_y As Double, address As Long
ReDim startRows(0 To MAXIMA - 1) As Long

Call r_transform_ground_to_pixel(i, testpoint.X, testpoint.Y, testpoint.z, p_x, p_y)
address = CLng((image_info(i).sub_height - 1) - (p_y + 70)) * image_info(i).sub_width * 3 ' + p_x * 3
startRows(i) = CLng((image_info(i).sub_height - 1) - (p_y + 70))

Close (1)
 Open image_info(i).FileName For Binary As 1
 
 Get #1, address + 1, Color_photo_Patch
 Close (1)
 r_maxc = -1
    For k = 1 To 11
     apu = Dir("c:\data\test_" & CStr(i) & "_" & CStr(k) & ".img")
     If apu = "" Then Exit Sub
     Open "c:\data\test_" & CStr(i) & "_" & CStr(k) & ".img" For Input As 12
     'Exit Sub
     Input #12, w_temp, h_temp, c_col, c_row
     Close (12)
     ReDim template(0 To w_temp - 1, 0 To h_temp - 1) As Byte
     Open "c:\data\test_" & CStr(i) & "_" & CStr(k) & ".raw" For Binary As 12
     Get #12, , template
     Close (12)
     ReDim Image(0 To w_temp - 1, 0 To h_temp - 1) As Byte
     For d_p_x = -4 To 4
      For d_p_y = -4 To 4
     
      apu = CLng((image_info(i).sub_height - 1) - (p_y + d_p_y)) - startRows(i)
      apu_y = apu - c_row
      apu_x = (p_x + d_p_x) - c_col
      For kX = 0 To w_temp - 1
        For kY = 0 To h_temp - 1
         Image(kX, kY) = CByte((Color_photo_Patch(apu_x + kX, apu_y + kY).g))
'         Exit Sub
        Next kY
      Next kX
      apu = MYFUNC_TEMPMEANSTD(template(0, 0), CLng(UBound(template, 1)), CLng(UBound(template, 2)), meaIm, stdIm)
      apu = MYFUNC_CCORB(template(0, 0), Image(0, 0), rx, meaIm, w_temp, h_temp)
      If rx > r_maxc Then
        r_maxc = rx
        ScaleMaxi = k
       End If
   Next d_p_y
  Next d_p_x
 Next k
 
   For i = 0 To NumOfImages - 1
    Call r_transform_ground_to_pixel(i, testpoint.X, testpoint.Y, testpoint.z, p_x, p_y)
    Form1.Picture1(i).DrawWidth = 3
    Form1.Picture1(i).FontSize = 8
    Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1), RGB(255, 255, 0)
   Next i
Ulim = 11
ReDim ScalingArr(1 To Ulim) As Double
ScalingArr(1) = 0.4: ScalingArr(2) = 0.5: ScalingArr(3) = 0.6
ScalingArr(4) = 0.7: ScalingArr(5) = 0.75: ScalingArr(6) = 0.8
ScalingArr(7) = 0.85: ScalingArr(8) = 0.9: ScalingArr(9) = 1#: ScalingArr(10) = 1.1
ScalingArr(11) = 1.2:

i = 0
Form1.Picture1(i).DrawWidth = 2
Dim Xt As Double, Yt As Double, Xf As Double, Yf As Double
Dim Shift As Double, p_x_beg  As Double, p_y_beg As Double, p_x_end As Double, p_y_end As Double
Shift = 3.14 / 10
'Radius = CDbl(Form1.Text5.Text) * ScalingArr(R0_ind)
'Radius = CDbl(Form1.Text5.Text) * ScalingArr(ScaleMaxi)
Exit Sub
  For phi = -3.14 To 3.15 Step 3.14 / 8
   Colorpix = RGB(120, 255, 0)   '  255 - Zdiff * 10)
   Xt = Radius * Cos(phi) + testpoint.X
   Yt = Radius * Sin(phi) + testpoint.Y
   Xf = Radius * Cos(phi + Shift) + testpoint.X
   Yf = Radius * Sin(phi + Shift) + testpoint.Y
   Zt = testpoint.z - 1
   Call r_transform_ground_to_pixel(i, Xt, Yt, Zt, p_x_beg, p_y_beg)
   Call r_transform_ground_to_pixel(i, Xf, Yf, Zt, p_x_end, p_y_end)
   Form1.Picture1(i).Line ((p_x_beg - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, ((image_info(i).Height - 1) - p_y_beg - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1)-((p_x_end - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, ((image_info(i).Height - 1) - p_y_end - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1), Colorpix
Next phi

FotoTMdone = True
Measurement.Foto_Dcrm = Radius * 2
Measurement.Corr_Match = r_maxc

End Sub
Public Sub savemeasurement()
Rem This routine adds a measurement

If LidarFitDone = True And FotoTMdone = True Then
 Open "C:\data\marv1_trees.txt" For Append As 1
 Dim capu As String
capu = Format$(Measurement.plot, "0") & ","
capu = capu & Format$(Measurement.Num, "0") & ","
capu = capu & Format$(Measurement.TreeSpecies, "0") & ","
capu = capu & Format$(Measurement.X, "0.00") & ","
capu = capu & Format$(Measurement.Y, "0.00") & ","
capu = capu & Format$(Measurement.z, "0.00") & ","
capu = capu & Format$(Measurement.z_DTM, "0.00") & ","
capu = capu & Format$(Measurement.Foto_h, "0.00") & ","
capu = capu & Format$(Measurement.H_lidar, "0.00") & ","
capu = capu & Format$(Measurement.Foto_Dcrm, "0.00") & ","
capu = capu & Format$(Measurement.Lidar_Dcrm, "0.00") & ","
capu = capu & Format$(Measurement.Crown_RMSE, "0.00") & ","
capu = capu & Format$(Measurement.cf, "0.0000") & ","
capu = capu & Format$(Measurement.pw, "0.0000") & ","
capu = capu & Format$(Measurement.cont, "0.0000")
Print #1, capu
Close (1)

Measurement.TreeSpecies = 0
Measurement.X = 0
Measurement.Y = 0
Measurement.z = 0
Measurement.z_DTM = 0
Measurement.Foto_h = 0
Measurement.H_lidar = 0
Measurement.Foto_Dcrm = 0
Measurement.Lidar_Dcrm = 0
Measurement.Crown_RMSE = 0
Measurement.cf = 0
Measurement.pw = 0
Measurement.cont = 0
LidarFitDone = False
FotoTMdone = False
Call PlotMeasurements
Else
 MsgBox ("You forgot either Multi-Size TM or crown modelling with lidar")
End If


End Sub


Public Function CCOR1() As Double

' This routine computes [-1,1] correlation image using cross-correlation. Correlation is computed
' for image akuva[] and template[] for the common area.
' BW-image BYTE (uchar) array *akuva[], declared as akuva(0 to width-1, 0 to height-1) in VB
' Template is the BYTE (uchar) array *temp[], declared also as temp(0 to wt-1, 0 to ht-1) in VB
' The center of the template (tree top's location) is given in (c_col, c_row)
' The output goes to the correlation double-array, declared same size as akuva-image-array.
'
     
'     Dim i As Long, j As Long
'     Dim m As Long, startrow As Long, startcol As Long, StartIndex As Long
'     Dim lkm_a As Long, lkm_b As Long
'     Dim nominator As Double, denominator1 As Double, denominator2 As Double, meanA As Double, meanB As Double
'     lkm_a = 0
'     meanA = 0#
    
'  Loop the template area, calculate mean of template, omit all zero values (outside  !
'    For l = 0 To wt - 1 ' cols
'      For m = 0 To ht - 1 ' rows
'            If temp_a(l, m) <> 0 Then
'              meanA = meanA + CDbl(temp_a(l, m))
'             lkm_a = lkm_a + 1
'            End If
'       Next m
'     Next l
'    meanA = meanA / CDbl(lkm_a)
' '  Too be sure we do not travel outside the image area, main loop is restricted to the area
' '  with template width & height removed from top&bottom, left & right.
' '  Constants
'    denominator1 = 0#
'    For l = 0 To wt - 1 ' loop template columns
'      For m = 0 To ht - 1 ' loop template rows
'        If temp_a(l, m) <> 0 Then '   it's a valid element
'         denominator1 = denominator1 + (CDbl(temp_a(l, m)) - meanA) * (CDbl(temp_a(l, m)) - meanA)
'         nom_term1(l, m) = CDbl(temp_a(l, m)) - meanA
'        End If
'      Next m
'    Next l


'    for (i= 1 ;i<=(width-1)-wt;i++) {          // i loops main image columns
'        for (j= 1;j<=(height-1)-ht;j++) {     // j loops main image rows

'    // Cross-correlation is now computed for point (i,j), => rkuva[j*width+i]
'    // c_col and c_row give the location of the HOTSPOT-point of the template.


            startrow = j ' ; // -(c_row);  // the location of the upper left corner of the template on the image
            startcol = i ' ; // -(c_col);  //                   "
'            startindex = startrow*width+startcol;  // " index in akuva[]
'
'    lkm_b=0;
'    meanB=0.0;
'    for (l=0;l<=wt-1;l++) {             // loop template columns
'        for (m=0;m<=ht-1;m++)   {       // loop template rows
'            if (temp_a(l,m)!=0)        // it's a valid element
'            {
'                lkm_b++;
'                temp_b(l,m)=akuva[startindex+(m*width)+l];
'                meanB += double(temp_b(l,m));
'            }
'
'        }
'    }
'    meanB /= double(lkm_b);

'    nominator=0.0;
'    denominator2=0.0;

    

'    for (l=0;l<=wt-1;l++) {            // loop template cols
'        for (m=0;m<=ht-1;m++)   {      // loop template rows
'            if (temp_a(l,m)!=0)
'            {
'             nominator +=  nom_term1(l,m) * (double(temp_b(l,m)) - meanB);
'             denominator2 +=  (double(temp_b(l,m)) - meanB) * (double(temp_b(l,m)) - meanB);
'            }
'        }
'    }
'   if ((denominator1 * denominator2) < 0.000002 )
'   {
'       j  = sprintf(buffer,"Trying to divide by zero in Myfunc_corima");
'       MessageBox ( NULL, buffer,"MyFunc_corima",MB_OK);
'       return -1.0;
'   }

'   // Here we put the results to the location of the hotspot
'   rkuva[(j+(c_row))*width+(i+(c_col))] = nominator/pow((denominator1 * denominator2),0.5);
'        }  // main loop rows
'    }      // main loop cols
'    return 1.0;
' }

End Function

Public Function ASIN(X As Double) As Double
ASIN = Atn(X / (-X * X + 1) ^ 0.5)
End Function
Public Sub CaptureEllipticTemplate(i As Integer, Ellipse As Ellipse)

Rem This Routine receives the middlepoint photo-coordinates (metric) of the
Rem template image, and the parameters of the ellipse ( a, b, rot angle, x, y).

Rem Routine scans a 101 x 101 pixel-area around the midpoint. If a scanline does not contain
Rem a single pixel belonging to the ellipse, scanline is dropped out, same applies for scan-
Rem columns.

Rem Parameter i refers to image (0 to numofimages-1)

Rem Type (struct) ellipse holds data:
Rem .alpha (rotation angle, rotation about mid-point, in photo-coordinate system)
Rem .a     (axis length in the direction of radial displacement (affected by Zdiff-parameter))
Rem .b     (axis length in crown width direction, affected by parameter TemplateWidth)
Rem .x     (x coord of center point in camera coordinates, affected by EllipseZasymmetry)
Rem .y     (y coord of center point in camera coordinates,  "       "        ")
Rem
Rem I.e. We let the ellipse be translated along the direction of the radial displacement, as
Rem defined by the metric EllipseZasymmetry -parameter, that 'tells' were the template center
Rem is located along the tree stem / ot may also be above the top, if parameter is set a positive
Rem value.

Dim p_x As Double, p_y As Double, apu_x As Double, apu_y As Double
Dim CenCol As Integer, CenRow As Integer     ' The image coordinates of the ellipse center (tree top)
Dim startrow As Integer, startcol As Integer
Dim endrow As Integer, endcol As Integer
Dim j As Integer, k As Integer, Ninside As Integer, jx As Integer, kX As Integer, apu
Dim Colorpix As RGBtriplet
Dim BWpix As Byte
Dim Inside As Boolean
Dim firstcol As Integer, firstrow As Integer, lastcol As Integer, lastrow As Integer
Dim firstrowdet As Boolean
Dim firstcoldet As Boolean

Dim HALFTEMPSIZE As Integer

HALFTEMPSIZE = 50

Rem Get the image pixel-coordinates of template center
Call a_transform_affine(i, 0, Ellipse.X, Ellipse.Y, p_x, p_y)
Rem Re-assigning (necessary?) and round to nearest integer value
CenCol = CInt(p_x - image_info(i).o_col): CenRow = CInt(((image_info(i).Height - 1) - p_y) - image_info(i).o_row)

Rem Define the outer coords of the scan area, which will be 101 x 101 in size
startrow = CenRow - HALFTEMPSIZE
endrow = CenRow + HALFTEMPSIZE
startcol = CenCol - HALFTEMPSIZE
endcol = CenCol + HALFTEMPSIZE

Rem Declare an array to hold the the template image (grayscale, Byte-array)
ReDim template(0 To endcol - startcol, 0 To endrow - startrow) As Byte
Rem Declare an array of same size. Boolean, each pixel value is TRUE for 'a member of the template' or FALSE for 'outside'
ReDim insidetemplate(0 To endcol - startcol, 0 To endrow - startrow) As Boolean

Rem The Template-image filenames are c:\temp\temp#.raw
Rem It's first deleted before opening in binary mode
Close (2)
Close (1)
Kill "c:\data\temp" & CStr(i) & ".raw"
Open "c:\data\temp" & CStr(i) & ".raw" For Binary As 2
'Exit Sub
Rem Open the Image file for input also  (for retrieving the pixels)
Open image_info(i).FileName For Binary As 1
'Exit Sub
Rem start scanning columns and rows, counters j and k are assigned pixel coordinates of the main image!
Rem counter jx and kx act...
Rem Counter Ninside
Rem Boolean inside
 Dim row As Double, col As Double
 Inside = False
 jx = -1
 For j = startcol To endcol     ' loop columns (x), there are 101, j is in the sub-image coords!
  Ninside = 0
  jx = jx + 1
  kX = -1
  
  For k = startrow To endrow    ' loop rows  (y) 101
    kX = kX + 1
    
    Rem Obtain photo-coordinates for this pixel location
    Call a_transform_affine(i, 1, CDbl(j + image_info(i).o_col), CDbl((image_info(i).Height - image_info(i).o_row - 1) - (k)), p_x, p_y)
    
    Rem Transfer this point's photo-coordinates (p_x, p_y) to ellipse center
    apu_x = p_x - Ellipse.X
    apu_y = p_y - Ellipse.Y
    
    Rem Rotate it to normal position, by alpha, alpha is adjusted by 90 degrees WHY???
    Rem There's a problem here.
    p_x = Cos(Ellipse.alpha) * apu_x + Sin(Ellipse.alpha) * apu_y
    p_y = -Sin(Ellipse.alpha) * apu_x + Cos(Ellipse.alpha) * apu_y
    
    Rem Check if the point is inside normalized ellipse
     Inside = False
     
     If ((p_x ^ 2) / (Ellipse.A ^ 2) + (p_y ^ 2) / (Ellipse.B ^ 2)) < 1# Then
      Ninside = Ninside + 1            ' Add Counter
      insidetemplate(jx, kX) = True    ' Set Boolean
     End If
     Rem According to image's color model perform color -> greyscale transformation
      
      Select Case image_info(i).Color
      Case 0 ' It is a BW image, we put pixel value in all three rgb-componenets
             ' Image begins at address 1, 2nd line at width + 1, k = row, j = column (0,..w-1)
      Get #1, CLng(k) * image_info(i).sub_width + CLng(j) + 1, BWpix
      Colorpix.r = BWpix
      Colorpix.g = BWpix
      Colorpix.B = BWpix
      Case 1 ' It is a COLOR image
      Get #1, CLng(k) * 3 * image_info(i).sub_width + CLng(j) * 3 + 1, Colorpix
      'Exit Sub
      End Select
      
      If insidetemplate(jx, kX) = False Then
       Colorpix.r = 0
       Colorpix.g = 0
       Colorpix.B = 0
      End If
      
      Rem template-array was a byte array, let's combine here. NOTE! There are more elegant
      Rem 3->1 transformations than average!
      ' apu = ((CInt(Colorpix.R) + CInt(Colorpix.g) + CInt(Colorpix.B)) / 3)
      'apu = ((CInt(0.333 * Colorpix.r) + CInt(0.333 * Colorpix.g) + CInt(0.333 * Colorpix.B)))
      apu = Colorpix.g
      If apu < 0 Then apu = 0
      If apu > 255 Then apu = 255
      template(jx, kX) = CByte(apu)
       
     Next k  ' Next row  in the 101x101 search area
    Next j ' Next column in the 101x101 search area
  
  Rem Now we have a 101x101 template, that needs to be polished from lines and rows
  Rem containing only zero-values (insidetemplate(col,row) = FALSE)
  
  
  firstrowdet = False
  firstcoldet = False
  firstcol = 0
  
  For j = 0 To endcol - startcol     ' Loop 101 columns
   Ninside = 0
   For k = 0 To endrow - startrow    ' Loop 101 rows
    If insidetemplate(j, k) = True Then
     Ninside = 1 ' there's a value in this column
    End If
   Next k
    If firstcoldet = False And Ninside = 1 Then
     Rem We've found the first column where there's an inside pixel
     firstcol = j
     firstcoldet = True
     GoTo nextcol   ' Hop to find the last column
    End If
    If firstcoldet = True And Ninside = 0 Then
     Rem We've found the last column
     lastcol = j - 1
     GoTo testrows ' Hop for finding the row-bounds
    End If
  
nextcol:  Next j

Rem We should never end up executing the next line!
MsgBox ("Image template n:o " & CStr(i) & " cannot define col borders!")

Rem Start to look for row-bounds
testrows:
firstrow = 0
  For k = 0 To endrow - startrow  ' Loop thru 101 rows
  Ninside = 0
   For j = 0 To endcol - startcol ' Loop thru 101 columns
    If insidetemplate(j, k) = True Then
     Ninside = 1 ' there's a value in this row
    End If
   Next j ' Next column
   If Ninside = 1 And firstrowdet = False Then
     firstrow = k
     firstrowdet = True
     GoTo nextrow
   End If
   If Ninside = 0 And firstrowdet = True Then
     lastrow = k - 1
     GoTo templatesizedetermined
   End If

nextrow:   Next k    ' next row

Rem We should never end up executing the next line!
MsgBox ("Image template n:o " & CStr(i) & " cannot define row borders!")

templatesizedetermined:
Rem The size of the template image is determined
Rem Declare the array

  ReDim temp(0 To (lastcol - firstcol), 0 To (lastrow - firstrow)) As Byte
  
  jx = -1
  For j = firstcol To lastcol   ' Loop columns in the 101x101 area
  jx = jx + 1
  kX = -1
   For k = firstrow To lastrow  ' Loop rows in the 101x101 area
   kX = kX + 1
   temp(jx, kX) = template(j, k)
   Next k
  Next j
 
 Rem write the header and the template
 Open "c:\data\temp" & CStr(i) & ".img" For Output As 3
 Print #3, lastcol - firstcol + 1    ' width
 Print #3, lastrow - firstrow + 1    ' height
 Print #3, HALFTEMPSIZE - firstcol + Ellipse.dx_col  ' c_col HOTSPOT, it's not here, it's off by
 Print #3, HALFTEMPSIZE - firstrow - Ellipse.dx_row  ' c_row HOTSPOT
 Close (3)
 Put #2, 1, temp
 Close (2)
 Close (1)


End Sub

Public Sub RGB_Vector_For_Point(ByRef red As Double, green As Double, blue As Double, testpoint As point3D, Image As Long)
 Dim p_x As Double, p_y As Double
 Dim j As Long
 j = Image
 Call r_transform_ground_to_pixel(j, testpoint.X, testpoint.Y, testpoint.z, p_x, p_y)
 Close (2)
 Open image_info(j).FileName For Binary As 2
 Dim paikka As Long, startrow As Long, startcol As Long
 ReDim Winrow(1 To 3) As RGBtriplet
 startcol = p_x - 1
 startrow = (image_info(j).Height - 1) - p_y - 1
 'Dim red As Double, blue As Double, green As Double
   For lx = -1 To 1
   paikka = CLng(startrow + lx) * (image_info(j).sub_width) * CLng(3) + CLng(startcol) * 3
   Get #2, paikka + 1, Winrow
   red = red + Winrow(1).r + Winrow(2).r + Winrow(3).r
   green = green + Winrow(1).g + Winrow(2).g + Winrow(3).g
   blue = blue + Winrow(1).B + Winrow(2).B + Winrow(3).B
  Next lx
 red = red / 9#
 green = green / 9#
 blue = blue / 9#
Close (2)
End Sub


Public Sub a_transform_affine(ByVal i As Integer, ByVal direction As Integer, X As Double, Y As Double, p_x As Double, p_y As Double)
Dim a_ As Double, b_ As Double, c_ As Double, d_ As Double, e_ As Double, f_ As Double

 a_ = image_info(i).a_
 b_ = image_info(i).b_
 c_ = image_info(i).c_
 d_ = image_info(i).d_
 e_ = image_info(i).e_
 f_ = image_info(i).f_
  
 Select Case direction
 Case 0   ' From camera coordinates to image cordinates (pixels)
 p_x = a_ * X + b_ * Y + c_
 p_y = d_ * X + e_ * Y + f_
 Case 1  ' From image coordinates (pixels) to camera cordinates (mm)
 p_x = (-e_ * c_ + e_ * X + f_ * b_ - Y * b_) / (a_ * e_ - b_ * d_)
 p_y = -(a_ * f_ - a_ * Y - c_ * d_ + X * d_) / (a_ * e_ - b_ * d_)
' p_x = p_x + camera.ps_x
' p_y = p_y + camera.ps_y
 End Select
End Sub

Public Sub PlotMeasurements()
Rem For plotting MARV4-data
Open "C:\data\marv1_trees.txt" For Input As 1
'Exit Sub
Dim Mplot As Long, MNum As Long, MTreeSpecies As Long, mx As Double, MY As Double, Mz As Double
Dim Mz_DTM As Double, MFoto_h As Double, MH_lidar As Double, MFoto_Dcrm As Double, MLidar_Dcrm As Double
Dim MCrown_RMSE As Double, Mcf As Double, Mpw As Double, Mcont As Double
Dim p_x As Double, p_y As Double
 For i = 0 To NumOfImages - 1
  Form1.Picture1(i).Cls
 Next i
Newplot = 0
Do Until EOF(1)
Input #1, Mplot, MNum, MTreeSpecies, mx, MY, Mz, Mz_DTM, MFoto_h, MH_lidar, MFoto_Dcrm, MLidar_Dcrm, MCrown_RMSE, Mcf, Mpw, Mcont, dummy, dummy, dummy
'Close (1)
'Exit Sub
If Mplot <> Newplot Then
 i = 0
 Call r_transform_ground_to_pixel(i, mx + 35, MY, Mz, p_x, p_y)
 Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, -1 + ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y), Colorpix
 Form1.Picture1(i).FontSize = 15
 Form1.Picture1(i).Print Mplot
End If
Newplot = Mplot
'Exit Sub
'If Mplot = Measurement.plot Then
 For i = 0 To NumOfImages - 1
 Call r_transform_ground_to_pixel(i, mx, MY, Mz, p_x, p_y)
 Select Case MTreeSpecies
 Case 1
 Colorpix = RGB(255, 0, 0)
 Case 2
 Colorpix = RGB(0, 255, 0)
 Case 3
 Colorpix = RGB(100, 120, 255)
 Case 4
 Colorpix = RGB(225, 225, 25)
 End Select
 Form1.Picture1(i).DrawWidth = 4
 Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, -1 + ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y), Colorpix
 Form1.Picture1(i).FontSize = 10
 'Form1.Picture1(i).Print MNum
 Form1.Picture1(i).DrawWidth = 1
 Next i
'End If
Loop
Close (1)
End Sub

Public Sub Plot_measurements()

Dim p_x As Double, p_y As Double
Dim Pnumber As Double
Dim inputstring As String
Dim s_start As Integer, s_end As Integer
Dim s_case As Integer, i As Integer
'On Error GoTo error_in_reading_and_plotting
Close (11)

ReDim ColorMap(1 To 64, 1 To 3) As Byte
l = 0
Open "C:\data\JetColorMAP.csv" For Input As 2
Do Until EOF(2)
l = l + 1
Input #2, ColorMap(l, 1), ColorMap(l, 2), ColorMap(l, 3)
Loop
Close (2)

MaxDTM = 35
MinDTM = 0



'Open "f:\temp\kuiva\Kuiva_polyline.txt" For Input As 1
ReDim PlotPoly(0 To 200) As Point
'i = 0
'Do Until EOF(1)
' Input #1, PlotPoly(i).X, PlotPoly(i).Y, dummy
' i = i + 1
'Loop
Close (1)
Rem Make last point meet first
PlotPoly(i).X = PlotPoly(0).X
PlotPoly(i).Y = PlotPoly(0).Y
ReDim Preserve PlotPoly(0 To i) As Point
Close (1)

Dim teksti As String

Open TreedataFilename For Input As 11
'Exit Sub
Do Until EOF(11)
'Exit Sub


' Input #11, SavedMeasurement.Num, teksti, SavedMeasurement.X, SavedMeasurement.Y, SavedMeasurement.z ', teksti ' , dummy
' Input #11, SavedMeasurement.X, SavedMeasurement.Y, SavedMeasurement.z, teksti ' , dummy
' Input #11, SavedMeasurement.X, SavedMeasurement.Y, SavedMeasurement.z, dummy, dummy, dummy, dummy, dummy, dummy, dummy, dummy, dummy, dummy
'Close (11)
'Exit Sub
' Input #11, SavedMeasurement.X, SavedMeasurement.Y, SavedMeasurement.z ', hf
'Close (11)
'Exit Sub
 Line Input #11, inputstring
'GoTo SkipParse
Rem Parse it
s_case = 0
s_start = 0
 
 For i = 2 To 100
 If Mid$(inputstring, i, 1) = "," Then
  s_end = i
  s_case = s_case + 1
  Select Case s_case
  Case 1 ' Num
  SavedMeasurement.Num = CInt(Mid$(inputstring, s_start + 1, s_end - s_start - 1))
  s_start = s_end
  Case 2 ' X
  SavedMeasurement.X = CDbl(Mid$(inputstring, s_start + 1, s_end - s_start - 1))
  s_start = s_end
  Case 3 ' Y
  SavedMeasurement.Y = CDbl(Mid$(inputstring, s_start + 1, s_end - s_start - 1))
  s_start = s_end
  Case 4 ' Z
  SavedMeasurement.z = CDbl(Mid$(inputstring, s_start + 1, s_end - s_start - 1))
  s_start = s_end
  End Select
 End If
Next i
SkipParse:
' MsgBox (Pnumber & " " & X_ini & " " & Y_ini & " " & Z_ini)
If SavedMeasurement.X > 1 Then
Z_maa = getheight(SavedMeasurement.X, SavedMeasurement.Y)

'Z_maa = SavedMeasurement.z - Z_maa
'Open "c:\Data\ilkka.txt" For Append As 200
'Print #200, SavedMeasurement.x, SavedMeasurement.y, SavedMeasurement.z, Z_maa
'Close (200)
Dim testpoint As Point
testpoint.X = SavedMeasurement.X
testpoint.Y = SavedMeasurement.Y

'INOUT = InsidePolygon(PlotPoly, UBound(PlotPoly), testpoint)
'If INOUT = Inside Then


     Apu_Z = 255 - (255 - (Z_maa - MinDTM) * 255 / (MaxDTM - MinDTM)) ' DTM
     Rem HSV-code
     If Apu_Z < 4 Then Apu_Z = 4
     If Apu_Z > 254 Then Apu_Z = 253
     Apu_Z = Apu_Z / 4
     Colorpix = RGB(ColorMap(Apu_Z, 1), ColorMap(Apu_Z, 2), ColorMap(Apu_Z, 3))

Colorpix = RGB(255, 255, 255)
Dim maa_x As Double, maa_y As Double
For i = 0 To NumOfImages - 1
Call r_transform_ground_to_pixel(i, SavedMeasurement.X, SavedMeasurement.Y, SavedMeasurement.z, p_x, p_y)
alku_x = (p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 0
alku_y = ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 0
Call r_transform_ground_to_pixel(i, SavedMeasurement.X, SavedMeasurement.Y, Z_maa, p_x, p_y)
p_x = (p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 0
p_y = ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 0



Form1.Picture1(i).DrawWidth = 1
'Form1.Picture1(i).Line (alku_x, alku_y)-(p_x, p_y), RGB(225, 225, 225)

Call r_transform_ground_to_pixel(i, SavedMeasurement.X, SavedMeasurement.Y, SavedMeasurement.z, p_x, p_y)


Form1.Picture1(i).DrawWidth = 3
Form1.Picture1(i).FontBold = False
Form1.Picture1(i).Font = Arial
Form1.Picture1(i).FontSize = 10
Form1.Picture1(i).ForeColor = RGB(255, 255, 0)
Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, -1 + ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y), Colorpix
col_a = (p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1
row_a = -1 + ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y
Form1.Picture1(i).FillStyle = 1
'Form1.Picture1(i).Circle (col_a, row_a), 5


'Exit Sub
'Form1.Picture1(i).Print Format$(SavedMeasurement.Num, "0")
'Form1.Picture1(i).Print " " & teksti

Form1.Picture1(i).ForeColor = RGB(255, 255, 255)
' Form1.Picture1(i).Print Pnumber & " " & Format$(Z_ini, "#.0")
   If Numbers_Plotted = True Then
   ' Form1.Picture1(i).Print SavedMeasurement.Num ' Format$(SavedMeasurement.z - Z_maa, "0.00 m") '
   ' Form1.Picture1(i).Print Format$(SavedMeasurement.z - Z_maa, "0.0") '
   End If
 Form1.Picture1(i).DrawWidth = 1
 Next i
End If
'End If

Loop
Close (11)

Exit Sub
error_in_reading_and_plotting:
MsgBox ("An error occurred in reading the file: " & TreedataFilename & " , check that it has no blank lines at the end or eq.")
Exit Sub

End Sub
Public Sub calculate_region(ByVal i As Integer, CenterCol As Integer, CenterRow As Integer, ByVal Width As Integer, ByVal Height As Integer, ByRef Start_Col As Integer, Start_Row As Integer, End_Col As Integer, End_Row As Integer, pan_x As Double, pan_y As Double)

Rem For image i, check if for region width, height, centered at CenterCol, CenterRow
Rem may be defined, return values for StartCol, StartRow, EndCol, EndRow
Rem If Overflow in any direction (Col, Row), positive, negative
Rem    then adjust CenterCol or CenterRow to fit main image

Rem i , CenterCol, CenterRow  width,  height
Rem StartCol , Start_Row EndCol , EndRow
Rem calculate MinCol, MaxCol, MinRow, MaxRow
Dim MinCol, maxcol, MinRow, maxrow As Integer
Rem check first that window isn't larger than image

'MsgBox (width)
'MsgBox (height)

If Width > image_info(i).sub_width Or Height > image_info(i).sub_height Then
  MsgBox ("Pan window is too large")
  Start_Col = 0
  Start_Row = 0
  End_Col = image_info(i).sub_width - 1
  End_Row = image_info(i).sub_height - 1
 Exit Sub
End If

Rem startCol & end Col

Recalculate:
If Width Mod 2 = 1 Then           ' odd width
MinCol = CInt(CenterCol - (Width - 1) / 2)
If MinCol < 0 Then
 MinCol = 0
 maxcol = MinCol + (Width) - 1
 GoTo checkheight
End If
maxcol = CInt(CenterCol + (Width - 1) / 2)
 If (maxcol > (image_info(i).sub_width - 1)) Then
  maxcol = image_info(i).sub_width - 1
  MinCol = maxcol - Width
  GoTo checkheight
 End If
End If
If Width Mod 2 = 0 Then
MinCol = CenterCol - (Width / 2) + 1
 If MinCol < 0 Then
  MinCol = 0
  maxcol = MinCol + Width - 1
  GoTo checkheight
 End If
maxcol = CenterCol + Width / 2
 If maxcol > (image_info(i).sub_width - 1) Then
  maxcol = image_info(i).sub_width - 1
  MinCol = maxcol - Width + 1
  GoTo checkheight
 End If
End If

checkheight:
If Height Mod 2 = 1 Then           ' odd height
MinRow = CInt(CenterRow - (Height - 1) / 2)
If MinRow < 0 Then
 MinRow = 0
 maxrow = MinRow + Height - 1
 GoTo loppu
End If
maxrow = CInt(CenterRow + (Height - 1) / 2)
If maxrow > (image_info(i).sub_height - 1) Then
 maxrow = image_info(i).sub_height - 1
 MinRow = maxrow - Height + 1
 GoTo loppu
End If
End If
If Height Mod 2 = 0 Then
MinRow = CenterRow - (Height / 2) + 1
 If MinRow < 0 Then
  MinRow = 0
  maxrow = MinRow + Height - 1
  GoTo loppu
 End If
maxrow = CenterRow + Height / 2
 If maxrow > (image_info(i).sub_height - 1) Then
  maxrow = (image_info(i).sub_height - 1)
  MinRow = maxrow - Height + 1
 End If
End If

loppu:
Start_Col = MinCol
Start_Row = MinRow
End_Col = maxcol
End_Row = maxrow

win_info(i).win_o_col = MinCol
win_info(i).win_o_row = MinRow

End Sub

Public Sub Move_By_Step(key As String)
Dim p_x As Double, p_y As Double
Dim Xr As Double, Yr As Double, Zr As Double

' get the current solution

Xr = X_sol
Yr = Y_sol
Zr = Z_sol

' Update Yr
Rem Select the case
Select Case key
Case "E" ' Move to East
Xr = Xr + StepValue
Case "W" ' Move to West
Xr = Xr - StepValue
Case "N" ' Move to North
Yr = Yr + StepValue
Case "S" ' Move to South
Yr = Yr - StepValue
Case "U" ' Move Up
Zr = Zr + StepValue
Case "D" ' Move Down
Zr = Zr - StepValue
End Select

 X_sol = Xr
 Y_sol = Yr
 Z_sol = Zr
 
 Form1.Label4(0).Caption = Format$(X_sol, "#.00")
 Form1.Label4(1).Caption = Format$(Y_sol, "#.00")
 Form1.Label4(2).Caption = Format$(Z_sol, "#.00")

' Update Images 1, 2 and 3
Dim i As Integer
For i = 0 To NumOfImages - 1
 Form1.Picture1(i).DrawWidth = 3
 Call r_transform_ground_to_pixel(i, Xr, Yr, Zr, p_x, p_y)
 Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1), RGB(255, 255, 255)
 Form1.Picture1(i).DrawWidth = 1
Next i

'Call center_to_xyz

Form1.Label10.Caption = "Solution moved as " & key & " pressed by " & StepValue & " meters"

End Sub
Public Sub MATVECMULT_D(Y() As Double, Dim_Y As Integer, A() As Double, na As Integer, ma As Integer, YOUT() As Double)

Rem Multiply matrix A(i,j) with vector Y(j) to get vector X(i)

If ma <> Dim_Y Then
  MsgBox ("In MATVECMULT_D dimensions of Matrix and Vector conflict!")
  Exit Sub
End If
Dim ix As Integer, jx As Integer
For ix = 1 To na
  For jx = 1 To ma
    YOUT(ix) = YOUT(ix) + A(ix, jx) * Y(jx)
  Next jx
 Next ix
 

End Sub

Public Sub indexx(ByVal N As Long, ra() As Single, indx() As Long)

Rem Heapsort from Numerical recipes
Rem indx(1 to N) is an index table, returned
Rem ra(1 to N) is the data to be sorted, left as is.

Dim i As Long, j As Long, l As Long, IR As Long, indxt As Long
Dim q As Single
Dim apu As Integer

On Error GoTo virhe_sort

For j = 1 To N
 indx(j) = j
Next j

l = N / 2 + 1
IR = N

10:
 If (l > 1) Then
  l = l - 1
  indxt = indx(l)
  q = ra(indxt)
 Else
  indxt = indx(IR)
  q = ra(indxt)
  indx(IR) = indx(1)
  IR = IR - 1
   If (IR = 1) Then
    indx(1) = indxt
   Exit Sub
   End If
 End If
i = l
j = l + l
20: If j <= IR Then
     If j < IR Then
      If (ra(indx(j)) < ra(indx(j + 1))) Then j = j + 1
     End If
     If (q < ra(indx(j))) Then
      indx(i) = indx(j)
      i = j
      j = j + j
     Else
     j = IR + 1
     End If
     GoTo 20
    End If
   indx(i) = indxt
   GoTo 10
virhe_sort:
MsgBox ("Error in Heapsort, l is " & l)
End Sub
Public Function FindSpecies(red As Double, green As Double, blue As Double, IR As Double) As Byte

'Exit Function
Rem Minimum distance classifier
'IR = 2424
'blue = 246


For i = 0 To 3
 If Form1.Check2(i).Value = 1 Then
  FindSpecies = i + 1
 End If
Next i
Exit Function

ReDim meanvect(1 To 4) As Point2D
ReDim VarCov(1 To 4, 1 To 2, 1 To 2) As Double
ReDim Detc(1 To 4) As Double
meanvect(1).X = 262.990566
meanvect(1).Y = 1930.501048
meanvect(2).X = 213.0305944
meanvect(2).Y = 1707.183566
meanvect(3).X = 246.0503686
meanvect(3).Y = 2424.57371
meanvect(4).X = 285.8220339
meanvect(4).Y = 870.4576271
VarCov(1, 1, 1) = 0.00163134
VarCov(1, 1, 2) = -0.00014061
VarCov(1, 2, 1) = -0.00014061
VarCov(1, 2, 2) = 0.00002091
VarCov(2, 1, 1) = 0.00186288
VarCov(2, 1, 2) = -0.00011466
VarCov(2, 2, 1) = -0.00011466
VarCov(2, 2, 2) = 0.00001327
VarCov(3, 1, 1) = 0.00127503
VarCov(3, 1, 2) = -0.00006913
VarCov(3, 2, 1) = -0.00006913
VarCov(3, 2, 2) = 0.00000797
VarCov(4, 1, 1) = 0.00023066
VarCov(4, 1, 2) = -0.00004818
VarCov(4, 2, 1) = -0.00004818
VarCov(4, 2, 2) = 0.00001856
Detc(1) = 7.843681841
Detc(2) = 7.936855713
Detc(3) = 8.26850951
Detc(4) = 8.707540423

ReDim Dist(1 To 4) As Double
ReDim vect1(1 To 2) As Double
ReDim Vect2(1 To 2) As Double
Dim maxdist As Double
maxdist = -100000

For i = 1 To 4
  vect1(1) = (meanvect(i).X - blue) * VarCov(i, 1, 1) + (meanvect(i).Y - IR) * VarCov(i, 2, 1)
  vect1(2) = (meanvect(i).X - blue) * VarCov(i, 1, 2) + (meanvect(i).Y - IR) * VarCov(i, 2, 2)
  Dist(i) = -0.5 * Detc(i) - 0.5 * (vect1(1) * (meanvect(i).X - blue) + vect1(2) * (meanvect(i).Y - IR))
  If Dist(i) > maxdist Then
   maxdist = Dist(i)
   sp = i
  End If
Next i
'Form1.Picture1(0).Cls
Form1.Picture1(0).Print sp
Form1.Picture1(0).CurrentY = 10
FindSpecies = sp
Exit Function


ReDim Spvect(1 To 4) As Vector3D
Rem Assign 1 = Pine, 2 = Spruce , 3 = Birch , 4 = Aspen
Spvect(1).X = 194.5: Spvect(1).Y = 100.2: Spvect(1).z = 98.1
Spvect(2).X = 191.2: Spvect(2).Y = 63.8:  Spvect(2).z = 69.4
Spvect(3).X = 212.2: Spvect(3).Y = 185.4: Spvect(3).z = 141.6
Spvect(4).X = 148.2: Spvect(4).Y = 151.4: Spvect(4).z = 131.9
ReDim SpDist(1 To 4) As Double
SpDist(1) = 0: SpDist(2) = 0: SpDist(3) = 0: SpDist(4) = 0:
Dim Species As Long, Mindist As Double
Mindist = 255#
For j = 1 To 4
 SpDist(j) = Sqr((Spvect(j).X - red) ^ 2 + (Spvect(j).Y - green) ^ 2 + (Spvect(j).z - blue) ^ 2)
 If SpDist(j) < Mindist Then
  Mindist = SpDist(j)
  If Mindist < 100 Then
   Species = j
  Else
  Species = 0
 End If
End If
Next j
'FindSpecies = Species
End Function
Public Sub MMULT_D(A() As Double, na As Integer, ma As Integer, B() As Double, NB As Integer, MB As Integer, C() As Double)

Rem multiply matrix B*A !! to get C

If na <> MB Then
  MsgBox ("Error in MMUL_D, cannot multiply matrices with dimension conflict!")
  Exit Sub
End If
Dim ix As Integer, jx As Integer, kX As Integer

For ix = 1 To NB  ' loop B's rows
 For jx = 1 To ma  ' loop A's cols
   C(ix, jx) = 0
   For kX = 1 To na ' loop A's rows
    C(ix, jx) = C(ix, jx) + B(ix, kX) * A(kX, jx)
    ' NA x MB    NA x MB     NB x MB     NA x MA
   Next kX
  Next jx
 Next ix

End Sub

Public Sub read_set_file_for_an_image(i As Integer)
Input #1, image_info(i).ImageCode   ' Image long integer code
Input #1, image_info(i).sub_width   ' sub-image (to be stored in an array) width
Input #1, image_info(i).sub_height  ' sub-image height
Input #1, image_info(i).FileName    ' filename to be opened, containing raw-data (RGB)
Input #1, image_info(i).Color       ' 1 for 24-bit color, 0 for greyscale
Input #1, image_info(i).sub_c_col   ' sub-image will be centered at this col when program initiates
Input #1, image_info(i).sub_c_row   ' sub-image will be centered at this row when program initiates
Input #1, image_info(i).o_col       ' the col value for the sub-image origo, in the main image coord. system with Y(row)-axis pointing down
Input #1, image_info(i).o_row       ' the row value for the sub-image origo, in the main image coord. system
Input #1, image_info(i).Width       ' width of the main image
Input #1, image_info(i).Height      ' heigth of the main image
Input #1, image_info(i).C           ' camera constant
Input #1, image_info(i).x_ps        ' in FC-coordinates the x-coordinate of the PPA-point (principal point)
Input #1, image_info(i).y_ps        ' in FC-coordinates the y-coordinate of the PPA-point
Input #1, image_info(i).lambda      ' Helmert rotation-angle, [rad]
Input #1, image_info(i).alpha       ' Helmert scale factor
Input #1, image_info(i).mean_x      ' Helmert mean x of camera coords
Input #1, image_info(i).mean_y      ' Helmert mean y of camera coords
Input #1, image_info(i).X_mean      ' Helmert mean ROW of image coords
Input #1, image_info(i).Y_mean      ' Helmert mean COL of image coords
Input #1, image_info(i).a_          ' Affine a
Input #1, image_info(i).b_          ' Affine b
Input #1, image_info(i).c_          ' Affine c
Input #1, image_info(i).d_          ' Affine d
Input #1, image_info(i).e_          ' Affine e
Input #1, image_info(i).f_          ' Affine f
Input #1, image_info(i).omega       ' ext. orientation angle omega (X)
Input #1, image_info(i).phi         ' ext. orientation angle phi (Y)
Input #1, image_info(i).kappa       ' ext. orientation angle kappa (Z)
Input #1, image_info(i).Xo          ' proj. center X-coord.
Input #1, image_info(i).Yo          ' proj. center Y-coord.
Input #1, image_info(i).Zo          ' proj. center Z-coord.
Input #1, image_info(i).Sun_azimuth    ' azimuth (direction clockwise from KKJ-North in radians)
Input #1, image_info(i).Sun_elevation  ' sun's elevation in radians over horizon

Input #1, image_info(i).StartOf_string              ' Starting row for additional exposures
Input #1, image_info(i).Num_of_addit_expos           ' Number of additional exposures
Rem If we have more exposures, we read types, sizes and locations of them
If image_info(i).Num_of_addit_expos > 0 Then
    For j = 1 To image_info(i).Num_of_addit_expos
        Input #1, image_info(i).AdditType(j)         ' 0 = BW 8-byte,1 = RGB 8-byte,2 = BW 16-byte, 3 = RGBIR 16-byte
        Input #1, image_info(i).AdditFileName(j)     ' Image location
        Input #1, image_info(i).AdditWidth(j)        ' width in pixels of aerial photo (main image)
        Input #1, image_info(i).AdditHeight(j)       ' height in pixels of aerial photo
    Next j
End If
 Imagesdisplayed(i) = AER_IMA
Call r_transform_matrix(i)

End Sub
Public Sub GetEightDigitLong(ByRef DigitString As String)

Rem DigitString has 8 chars

Dim NameLong As Long
NameLong = Int(Rnd() * 100000000#)
DigitString = Format$(NameLong, "00000000")

End Sub
Public Sub set_window_sizes()
If NumOfImages > 0 And NumOfImages < 4 Then
Rem There'll be just one row of images (picture-boxes), make their width an even number and height equal to width
 If (((Form1.ScaleWidth - 80) / NumOfImages)) < (Form1.ScaleHeight - 80) Then
  win_h = (Form1.ScaleWidth - 80) / NumOfImages
  win_h = win_h + (win_h * 3) Mod 4
  Win_w = win_h
 Else
 win_h = (Form1.ScaleHeight - 80)
 win_h = win_h + (win_h * 3) Mod 4
 Win_w = win_h
 End If
End If

If NumOfImages > 3 And NumOfImages < 7 Then
Rem There'll be two rows of images on the main program window
 If (((Form1.ScaleWidth - 140) / 3)) < ((Form1.ScaleHeight - 140) / 2) Then
  win_h = (Form1.ScaleWidth - 140) / 3
  win_h = win_h + (win_h * 3) Mod 4
  Win_w = win_h
 Else
 win_h = ((Form1.ScaleHeight - 140) / 2)
 win_h = win_h + (win_h * 3) Mod 4
 Win_w = win_h
 End If
End If

If NumOfImages > 6 And NumOfImages < 10 Then
Rem There'll be three rows of images
 If (((Form1.ScaleWidth - 140) / 3)) < ((Form1.ScaleHeight - 140) / 3) Then
  win_h = (Form1.ScaleWidth - 140) / 3
  win_h = win_h + (win_h * 3) Mod 4
  Win_w = win_h
 Else
 win_h = ((Form1.ScaleHeight - 200) / 3)
 win_h = win_h + (win_h * 3) Mod 4
 Win_w = win_h
 End If
End If

End Sub
Public Sub center_to_Sub_c()

Dim i As Integer, j As Integer, length As Integer
Dim startcol As Integer, startrow As Integer, endcol As Integer, endrow As Integer

On Error GoTo Error_in_center_to_sub_c

Form1.MousePointer = 11
DoEvents
For i = 0 To NumOfImages - 1
 Call calculate_region(i, image_info(i).sub_c_col, image_info(i).sub_c_row, CInt(1 / win_info(i).pan_x * Win_w), CInt(1 / win_info(i).pan_y * win_h), startcol, startrow, endcol, endrow, win_info(i).pan_x, win_info(i).pan_y)
 win_info(i).win_o_col = startcol
 win_info(i).win_o_row = startrow
 win_info(i).win_width = Win_w
 win_info(i).win_height = win_h

 If image_info(i).Color = 1 Then

 length = Stringlength(image_info(i).FileName)
 ReDim filename_in(0 To length) As Byte
 For j = 0 To length - 1
  filename_in(j) = CByte(Asc(Mid$(image_info(i).FileName, j + 1, 1)))
 Next j
 filename_in(length) = 0
 
 FileOut = "c:\data\pic" & CStr(i) & ".bmp"
 length = Len(FileOut)
 ReDim filename_out(0 To length) As Byte
 For j = 0 To length - 1
  filename_out(j) = CByte(Asc(Mid$(FileOut, j + 1, 1)))
 Next j
 filename_out(length) = 0
  
Call create_bmp(i, startcol, startrow, endcol, endrow, "c:\data\pic" & CStr(i) & ".bmp", win_info(i).pan_x, win_info(i).pan_y)
Call create_bmp_header(FileOut, CLng(Win_w), CLng(win_h))
' apu = MYFUNC_CREATEBMP(CLng(i), CLng(startcol), CLng(startrow), CLng(endcol), CLng(endrow), filename_in(0), filename_out(0), CDbl(win_info(i).pan_x), CDbl(win_info(i).pan_y), CLng(win_h), CLng(Win_w), CLng(image_info(i).sub_width))
ElseIf image_info(i).Color = 0 Then
 Call create_BW_bmp(i, startcol, startrow, endcol, endrow, "c:\data\pic" & CStr(i) & ".bmp", win_info(i).pan_x, win_info(i).pan_y)
End If

' Label10.Caption = CStr(startcol) & "," & CStr(startrow) & "," & CStr(endcol) & "," & CStr(endrow)
Form1.Picture1(i).Picture = LoadPicture("c:\data\pic" & CStr(i) & ".bmp")
Form1.Picture1(i).DrawWidth = 1
DoEvents

Next i
Form1.MousePointer = 1
Exit Sub
Error_in_center_to_sub_c:
MsgBox ("An error occurred in sub Error_in_center_to_sub_c_Click() ")
Close (10)

End Sub

Public Sub center_to_pan_center()
Dim i As Integer, j As Integer, k As Integer, N As Integer, l As Integer, m As Integer
Dim ix As Integer, jx As Integer, kX As Integer, lx As Integer, mx As Integer, nx As Integer, Num As Integer
Dim Num_of_ima As Integer, c_col As Integer, c_row As Integer
Dim cam(0 To MAXIMA) As Integer
Dim startcol As Integer, startrow As Integer, endcol As Integer, endrow As Integer, length  As Integer
Dim p_x As Double, p_y As Double

On Error GoTo Error_in_center_to_pan_center


X_ini = X_cen
Y_ini = Y_cen
Z_ini = Z_cen

' Check from which images there are observations


Form1.MousePointer = 11
DoEvents
For i = 0 To NumOfImages - 1
 Call r_transform_ground_to_pixel(i, X_ini, Y_ini, Z_ini, p_x, p_y)
 c_col = p_x - (image_info(i).o_col)
 c_row = (image_info(i).Height - 1) - p_y - (image_info(i).o_row)
 Call calculate_region(i, c_col, c_row, 1 / win_info(i).pan_x * Win_w, 1 / win_info(i).pan_y * win_h, startcol, startrow, endcol, endrow, win_info(i).pan_x, win_info(i).pan_y)
 win_info(i).win_o_col = startcol
 win_info(i).win_o_row = startrow
 win_info(i).win_width = Win_w
 win_info(i).win_height = win_h

 If image_info(i).Color = 1 Then

 length = Stringlength(image_info(i).FileName)
 ReDim filename_in(0 To length) As Byte
 For j = 0 To length - 1
  filename_in(j) = CByte(Asc(Mid$(image_info(i).FileName, j + 1, 1)))
 Next j
 filename_in(length) = 0
 
 FileOut = "c:\data\pic" & CStr(i) & ".bmp"
 length = Len(FileOut)
 ReDim filename_out(0 To length) As Byte
 For j = 0 To length - 1
  filename_out(j) = CByte(Asc(Mid$(FileOut, j + 1, 1)))
 Next j
 filename_out(length) = 0
  
 Call create_bmp(i, startcol, startrow, endcol, endrow, "c:\data\pic" & CStr(i) & ".bmp", win_info(i).pan_x, win_info(i).pan_y)
'Call create_bmp_header(FileOut, CLng(Win_w), CLng(win_h))
'apu = MYFUNC_CREATEBMP(CLng(i), CLng(startcol), CLng(startrow), CLng(endcol), CLng(endrow), filename_in(0), filename_out(0), CDbl(win_info(i).pan_x), CDbl(win_info(i).pan_y), CLng(win_h), CLng(Win_w), CLng(image_info(i).sub_width))
'ElseIf image_info(i).Color = 0 Then
' Call create_BW_bmp(i, startcol, startrow, endcol, endrow, "c:\temp\pic" & CStr(i) & ".bmp", win_info(i).pan_x, win_info(i).pan_y)
End If

' Label10.Caption = CStr(startcol) & "," & CStr(startrow) & "," & CStr(endcol) & "," & CStr(endrow)
Form1.Picture1(i).Picture = LoadPicture("c:\data\pic" & CStr(i) & ".bmp")
Form1.Picture1(i).DrawWidth = 1
For l = -1 To 1
For m = -1 To 1
 Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - m, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - l), RGB(0, 255, 0)
Next m
Next l
DoEvents

Next i
Form1.MousePointer = 1
Exit Sub
Error_in_center_to_pan_center:
MsgBox ("An error occurred in sub Center_to_pan_center_Click() ")
Close (10)


End Sub

Public Sub center_to_xyz()
Dim i As Integer, j As Integer, l As Integer, m As Integer
Dim startcol As Integer, startrow As Integer, endcol As Integer, endrow As Integer
Dim c_col As Integer, c_row As Integer, length As Integer
Dim p_x As Double, p_y As Double

For i = 0 To NumOfImages - 1
 Imagesdisplayed(i) = AER_IMA
Next i


'If SolutionExists = False Then GoTo Error_in_GCP_space_intersection2

X_ini = X_sol
Y_ini = Y_sol
Z_ini = Z_sol

' Check from which images there are observations


Form1.MousePointer = 11
DoEvents
For i = 0 To NumOfImages - 1
 
 
 
 Call r_transform_ground_to_pixel(i, X_ini, Y_ini, Z_ini, p_x, p_y)
 If p_x > 32000 Then p_x = 32000
 If p_x < -32000 Then p_x = -32000
 c_col = p_x - (image_info(i).o_col)
 If Abs((image_info(i).Height - 1) - p_y - (image_info(i).o_row)) > 32000 Then
  c_row = 0
  Else
  c_row = (image_info(i).Height - 1) - p_y - (image_info(i).o_row)
 End If
 Call calculate_region(i, c_col, c_row, 1 / win_info(i).pan_x * Win_w, 1 / win_info(i).pan_y * win_h, startcol, startrow, endcol, endrow, win_info(i).pan_x, win_info(i).pan_y)
 win_info(i).win_o_col = startcol
 win_info(i).win_o_row = startrow
 win_info(i).win_width = Win_w
 win_info(i).win_height = win_h

 If image_info(i).Color = 1 Then

 length = Stringlength(image_info(i).FileName)
 ReDim filename_in(0 To length) As Byte
 For j = 0 To length - 1
  filename_in(j) = CByte(Asc(Mid$(image_info(i).FileName, j + 1, 1)))
 Next j
 filename_in(length) = 0
 
 FileOut = "c:\data\pic" & CStr(i) & ".bmp"
 length = Len(FileOut)
 ReDim filename_out(0 To length) As Byte
 For j = 0 To length - 1
  filename_out(j) = CByte(Asc(Mid$(FileOut, j + 1, 1)))
 Next j
 filename_out(length) = 0
  
Call create_bmp(i, startcol, startrow, endcol, endrow, "c:\data\pic" & CStr(i) & ".bmp", win_info(i).pan_x, win_info(i).pan_y)
' Call create_bmp_header(FileOut, CLng(Win_w), CLng(win_h))
' apu = MYFUNC_CREATEBMP(CLng(i), CLng(startcol), CLng(startrow), CLng(endcol), CLng(endrow), filename_in(0), filename_out(0), CDbl(win_info(i).pan_x), CDbl(win_info(i).pan_y), CLng(win_h), CLng(Win_w), CLng(image_info(i).sub_width))
ElseIf image_info(i).Color = 0 Then
 Call create_BW_bmp(i, startcol, startrow, endcol, endrow, "c:\data\pic" & CStr(i) & ".bmp", win_info(i).pan_x, win_info(i).pan_y)
End If

' Label10.Caption = CStr(startcol) & "," & CStr(startrow) & "," & CStr(endcol) & "," & CStr(endrow)
Form1.Picture1(i).Picture = LoadPicture("c:\data\pic" & CStr(i) & ".bmp")
Form1.Picture1(i).DrawWidth = 3
Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y), RGB(255, 255, 255)

GoTo SkipCross:

For l = -25 To 25
 m = 0
 If Abs(l) > 7 Then
 Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x + m, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y + l), RGB(255, 255, 255)
 End If
Next l
For m = -25 To 25
 l = 0
 If Abs(m) > 7 Then
 Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x + m, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y + l), RGB(255, 255, 255)
 End If
Next m

SkipCross:
Form1.Picture1(i).CurrentX = Form1.Picture1(i).CurrentX + 20
Form1.Picture1(i).CurrentY = Form1.Picture1(i).CurrentY + 20

Form1.Picture1(i).FontSize = 18
Form1.Picture1(i).FontBold = True
Form1.Picture1(i).Font = "Arial"
Form1.Picture1(i).ForeColor = RGB(255, 255, 255)
'Form1.Picture1(i).Print Measurement.FieldNumber
Form1.Picture1(i).FontSize = 20
Form1.Picture1(i).FontBold = True
Form1.Picture1(i).Font = "Arial"

Form1.Picture1(i).CurrentX = 5
Form1.Picture1(i).CurrentY = 5
Form1.Picture1(i).Print image_info(i).ImageCode & "    1:" & Format$((image_info(i).Zo - 160) / image_info(i).C, "00000")


DoEvents

Next i
Form1.MousePointer = 1
Form1.Label10.Caption = "Images centered to given XYZ"
Dim key As String
'key = "N"
'For lm = 1 To 10
'Move_By_Step (key)
'Next lm
'key = "S"
'For lm = 1 To 10
'Move_By_Step (key)
'Next lm
' NCounter = NCounter + 1
' SavePicture Form1.Picture1(0).Image, "c:\temp\AVI\" & Format$(NCounter, "0000") & ".bmp"


Exit Sub
Error_in_GCP_space_intersection2:
MsgBox ("An error occurred in sub Center_to_xyz_Click() ")
Form1.Label10.Caption = "Image centering failed"
Close (10)

End Sub
Public Sub Place_Objects_On_The_Main_Form()

Rem This routine puts objects: picture1(i) -boxes that will hold image patches (image observations, superimposing)
Rem                            text1(i)    -boxes assiciated with a picture to hold photocoords (used)
Rem                            Label1(i)   for writing pixel coordinates (debugging)
Rem                            Check(i)    -boxes to indicate if an image is used (space intersection, matching..)
                              

Dim addition As Integer
Dim i As Integer, j As Integer

For i = 0 To NumOfImages - 1

 If i >= 0 And i <= 2 Then j = 0
 If i >= 3 And i <= 5 Then j = 1
 If i >= 6 And i <= 8 Then j = 2
 If j > 0 Then addition = 40

 Rem Pixture-boxes that hold the sub-images
 Form1.Picture1(i).Visible = True
 Form1.Picture1(i).Enabled = True
 Form1.Picture1(i).Left = 3 + (i Mod 3) * Win_w + (i Mod 3) * 3
 Form1.Picture1(i).Top = 96 + j * (win_h) + j * addition
 Form1.Picture1(i).Width = Win_w + 4
 Form1.Picture1(i).Height = win_h + 4
 Rem Labels that hold the image / camera coords for each sub-image
 Form1.Label1(i).Left = 3 + (i Mod 3) * Win_w + (i Mod 3) * 3
 Form1.Label1(i).Top = Form1.Picture1(i).Top - 20
 
 Rem Check-boxes that are used the define if an image is used or not in calculating space intersection
 Form1.Check1(i).Left = 3 + (i Mod 3) * Win_w + (i Mod 3) * 3 + 200
 Form1.Check1(i).Top = Form1.Picture1(i).Top - 25
 Form1.Check1(i).Width = 200
 Form1.Check1(i).Caption = "Use: " & image_info(i).ImageCode
 Rem Text-boxes that will hold the camera coords of the mouse-pointed image points
 Form1.Text1(i * 2).Top = Form1.Picture1(i).Top - 40
 Form1.Text1((i * 2) + 1).Top = Form1.Picture1(i).Top - 40
 Form1.Text1(i * 2).Left = 3 + (i Mod 3) * Win_w + (i Mod 3) * 3
 Form1.Text1((i * 2) + 1).Left = 3 + (i Mod 3) * Win_w + (i Mod 3) * 3 + 60
 Form1.Label1(i).Width = Form1.Text1(i).Width * 3
 Rem Text-boxes that will hold the widths
 Form1.Text12(i).Top = Form1.Picture1(i).Top - 40
 Form1.Text12(i).Left = 3 + (i Mod 3) * Win_w + (i Mod 3) * 3 + 2 * Form1.Text1(i * 2).Width
 Form1.Text12(i).Height = Form1.Text1(i).Height
 Form1.Text12(i).Width = Form1.Text1(i).Width
 Form1.Text12(i).Visible = True
 
 Rem Enable & make visible (by default these objects are not in use, because we don't know their number or location untill
 Rem we know the number of images
 Form1.Label1(i).Enabled = True
 Form1.Label1(i).Visible = True
 Form1.Text1(i * 2).Enabled = True
 Form1.Text1(i * 2).Visible = True
 Form1.Text1((i * 2) + 1).Enabled = True
 Form1.Text1((i * 2) + 1).Visible = True
 Form1.Check1(i).Enabled = True
 Form1.Check1(i).Visible = True
 Form1.Label2.Caption = "Measurement number: " & CStr(MeasurementCounter)
 DoEvents
Next i

Rem Make the visibility of unnecessary objects false (just to be sure it happens)
For i = NumOfImages To MAXIMA - 1
 Form1.Picture1(i).Visible = False
 Form1.Text1(i * 2).Visible = False
 Form1.Text1(i * 2 + 1).Visible = False
 Form1.Text12(i).Visible = False
 Form1.Check1(i).Visible = False
 Form1.Label1(i).Visible = False
Next i

End Sub
Public Sub Check_parameter_sanity(ByVal i As Integer)

Rem This routine checks that the parameters given in the set file are OK
Rem (!! To some degree only)

Dim A As Integer
If image_info(i).sub_width < 0 Or image_info(i).sub_width > 2 ^ 15 Then MsgBox ("Image (0-(N-1)): " & CStr(i) & "Parameter: Sub-image width:" & CStr(image_info(i).sub_width) & " in set file?")  ' sub-image (to be stored in an array) width
If image_info(i).sub_height < 0 Or image_info(i).sub_height > 2 ^ 15 Then MsgBox ("Image (0-(N-1)): " & CStr(i) & "Parameter: Sub-image height:" & CStr(image_info(i).sub_height) & " in set file?")  ' sub-image height
Rem image_info(i).filename    ' filename to be opened, containing raw-data (RGB)
If image_info(i).Color < 0 Or image_info(i).Color > 1 Then MsgBox ("Image (0-(N-1)): " & CStr(i) & "Parameter: image color:" & CStr(image_info(i).Color) & " in set file should be 1 (color) or 0 (grayscale)?")     ' 1 for 24-bit color, 0 for greyscale
'Input #1, image_info(i).sub_c_col   ' sub-image will be centered at this col when program initiates
'Input #1, image_info(i).sub_c_row   ' sub-image will be centered at this row when program initiates
'Input #1, image_info(i).o_col       ' the col value for the sub-image origo, in the main image coord. system with Y(row)-axis pointing down
'Input #1, image_info(i).o_row       ' the row value for the sub-image origo, in the main image coord. system
'Input #1, image_info(i).width       ' width of the main image
'Input #1, image_info(i).height      ' heigth of the main image
'Input #1, image_info(i).c           ' camera constant
'Input #1, image_info(i).x_ps        ' in FC-coordinates the x-coordinate of the PPA-point (principal point)
'Input #1, image_info(i).y_ps        ' in FC-coordinates the y-coordinate of the PPA-point
'Input #1, image_info(i).lambda      ' Helmert rotation-angle, [rad]
'Input #1, image_info(i).alpha       ' Helmert scale factor
'Input #1, image_info(i).mean_x      ' Helmert mean x of camera coords
'Input #1, image_info(i).mean_y      ' Helmert mean y of camera coords
'Input #1, image_info(i).X_mean      ' Helmert mean ROW of image coords
'Input #1, image_info(i).Y_mean      ' Helmert mean COL of image coords
'Input #1, image_info(i).a_          ' Affine a
'Input #1, image_info(i).b_          ' Affine b
'Input #1, image_info(i).c_          ' Affine c
'Input #1, image_info(i).d_          ' Affine d
'Input #1, image_info(i).e_          ' Affine e
'Input #1, image_info(i).f_          ' Affine f
'Input #1, image_info(i).omega       ' ext. orientation angle omega (X)
'Input #1, image_info(i).phi         ' ext. orientation angle phi (Y)
'Input #1, image_info(i).kappa       ' ext. orientation angle kappa (Z)
'Input #1, image_info(i).Xo          ' proj. center X-coord.
'Input #1, image_info(i).Yo          ' proj. center Y-coord.
'Input #1, image_info(i).Zo          ' proj. center Z-coord.
'Input #1, image_info(i).Sun_azimuth
'Input #1, image_info(i).Sun_elevation



End Sub

Public Sub CheckAccuracy(DigitString As String, TreeFile As String)

 Exit Sub

Rem This Routine finds
Rem
Rem (1) Number of found trees; for a cluster there exists a match in the Data
Rem     - Trees are listed, and statistics computed
Rem
Rem (2) Number of commission errors; a cluster is not found a match, within the circular plot + buffer
Rem     - Clusters are listed, and statistics computed
Rem
Rem (3) Number of omission errors; a tree falling within the circular plot is not found a match
Rem     - Trees are listed, and statistics computed
Rem
Rem Arguments:
Rem
Rem DigitString     has the ######### -eight numbers
Rem TreeFile        contains measured tree-data (filename)
Rem Height-diameter curve parameters h = 1.3 + d^2/(a+b*d)^2
Rem  Open the file  with the ground truth
Rem  Clusters are in global storage: clus() -array



Rem *************************************************************************************************
Rem Find MATCHES -section; here we take a cluster, and try to find a match for it in the ground truth

Dim MatchZdist As Double, MatchZrmse As Double
Dim N_tree_inside  As Integer
Dim N_trees_in_plot As Integer
Dim MeanHeightField As Double
Dim N_tree_inside_border As Integer
Dim Match2Ddist As Double, MatchXYrmse As Double
ReDim ClusMatchStruct(1 To mc) As ClusMatchStruct

Dim i As Integer

Rem Maximum allowed XY-distance in meters between field measured tree top and candidate tree top
Match2Ddist = 1.3
Rem Maximum allowed XY-distance in meters between field measured tree top and candidate tree top
MatchZdist = 3

OV.XYMatchDist = Match2Ddist
OV.ZMatchDist = MatchZdist

Rem ***********************************************************************************
Rem *** First, lets read all trees, and count the number of trees within the plot area,
Rem *** and in the ring defined by the Match2Dist
Rem

TreeFile = "c:\data\maku_vis_trees.txt"
TreeFile = "c:\data\ma4_vis_trees.txt"
TreeFile = "c:\data\mk_vis_trees.txt"
N_tree_inside_border = 0
N_tree_inside = 0
MeanHeightField = 0#
ReDim TreeMatchStruct(1 To 350) As TreeMatchStruct

Open TreeFile For Input As 44
'Exit Sub
Dim angle As Double, X_origo As Double, Y_origo As Double, Z_origo As Double
Dim X_shift As Double, Y_shift As Double, Z_shift As Double, X As Double, Y As Double, z As Double
Dim sina As Double, cosa As Double, lx As Integer, jx As Integer, kX As Integer
Dim d2_dist As Double
Input #44, angle: Input #44, X_origo: Input #44, Y_origo: Input #44, Z_origo
Input #44, X_shift: Input #44, Y_shift: Input #44, Z_shift
angle = angle * TO_RADIANS: cosa = Cos(angle): sina = Sin(angle)

 Do Until EOF(44)
  lx = lx + 1
  Input #44, X, Y, z, FTrees(lx).d13, FTrees(lx).Height, FTrees(lx).Num, FTrees(lx).Species, FTrees(lx).Status
  FTrees(lx).X = cosa * X - sina * Y + X_origo + X_shift
  FTrees(lx).Y = sina * X + cosa * Y + Y_origo + Y_shift
  FTrees(lx).Zbutt = z + Z_origo + Z_shift
  FTrees(lx).Ztop = FTrees(lx).Zbutt + FTrees(lx).Height
  'FTrees(lx).Zdem = getTINheight(FTrees(lx).x, FTrees(lx).y, CLng(0))
  FTrees(lx).Zdem = getheight(FTrees(lx).X, FTrees(lx).Y)
  Rem Check if the tree falls inside circular plot's XY-extent
  d2_dist = sqrt((FTrees(lx).X - PlotCenter.X) ^ 2 + (FTrees(lx).Y - PlotCenter.Y) ^ 2)
  
  If d2_dist <= (Plotradius + Match2Ddist) Then
   N_tree_inside_border = N_tree_inside_border + 1
   kX = N_tree_inside_border                             ' Number of trees inside the plot
   TreeMatchStruct(kX).Num = FTrees(lx).Num: TreeMatchStruct(kX).Spec = FTrees(lx).Species
   TreeMatchStruct(kX).Status = FTrees(lx).Status: TreeMatchStruct(kX).d13 = FTrees(lx).d13
   TreeMatchStruct(kX).h = FTrees(lx).Height: TreeMatchStruct(kX).X = FTrees(lx).X
   TreeMatchStruct(kX).Y = FTrees(lx).Y: TreeMatchStruct(kX).z = FTrees(lx).Ztop
   TreeMatchStruct(kX).ZbuttDEM = FTrees(lx).Zdem:  TreeMatchStruct(kX).ZButtTach = FTrees(lx).Zbutt
   TreeMatchStruct(kX).NclusInCylinder = 0
   TreeMatchStruct(kX).IndexMatchClus = 0: TreeMatchStruct(kX).Dist3D = -1#
   TreeMatchStruct(kX).IsMatch = False: TreeMatchStruct(kX).IsOmission = False
   TreeMatchStruct(kX).Xclus = 0: TreeMatchStruct(kX).Yclus = 0
   TreeMatchStruct(kX).Zclus = 0
   TreeMatchStruct(kX).IsInside = False
  End If
  Rem Check if the tree falls inside the plot border
  If d2_dist < Plotradius Then
   N_tree_inside = N_tree_inside + 1      ' N of field trees within plotradius
   MeanHeightField = MeanHeightField + FTrees(lx).Ztop
   TreeMatchStruct(kX).IsInside = True
  End If
 Loop ' Get the next ground truth tree observation
Close (44)
N_trees_in_plot = N_tree_inside

Form1.Label10.Caption = "Field plot trees read " & lx & " trees"

Dim MeanHeightFoto As Double
MatchXYrmse = 0#: MeanHeightFoto = 0#: MatchZrmse = 0#
MeanHeightField = MeanHeightField / CDbl(N_tree_inside)
'Exit Sub
Rem Looping thru Clusters

Dim N_cand_inside As Integer
Dim N_cand_inside_border As Integer
Dim j As Integer
Dim xp As Double, yp As Double, zp As Double
 
 For i = 1 To mc
   xp = cluscoords(i).X / (atausum(i))
   yp = cluscoords(i).Y / (atausum(i))
   zp = cluscoords(i).z / (atausum(i))
   d2_dist = ((xp - PlotCenter.X) ^ 2 + (yp - PlotCenter.Y) ^ 2) ^ 0.5
   If d2_dist <= (Plotradius + Match2Ddist) Then
    N_cand_inside_border = N_cand_inside_border + 1
    j = N_cand_inside_border
    MeanHeightFoto = MeanHeightFoto + zp
    ClusMatchStruct(j).IsInside = False
    ClusMatchStruct(j).X = xp
    ClusMatchStruct(j).Y = yp
    ClusMatchStruct(j).z = zp
    ClusMatchStruct(j).Index = j
    ClusMatchStruct(j).IsMatch = False
    ClusMatchStruct(j).IsCommission = True
   End If
   If d2_dist < (Plotradius) Then
    N_cand_inside = N_cand_inside + 1
    ClusMatchStruct(j).IsInside = True
   End If
 Next i ' Next cluster found
 
 Rem The mean height (elevation) of the candidates within the plot
 MeanHeightFoto = MeanHeightFoto / CDbl(N_cand_inside)
 
 
 Rem ************************************************************************
 Rem *** FINDING MATCHES FOR THE FIELD TREES --> CANDIDATE WITHIN THE PLOT***
 Rem ************************************************************************
 
 Rem 1) For each tree and cluster, compute the number of clusters and Treetops in their
 Rem    Match-Cylinders (loop clusters and loop trees, O(N^2)), compute closest pairs (minimum distances)
 Rem
 Rem    Loop through Trees, On the first round Mark the with 0 c / 0 t and 1 c / 1 tree cases
 Rem    to be SURE omission-errors and Matches
 Rem
 Rem    At this point we are left with tree tops that have cases (1 c / > 1 t) and (>1 c / >1 t)
 Rem                                   clusters that have 0 trees in their Match-cylinder
 Rem
 Rem    Loop thru clusters (within boundary) and Mark 0 t -cases as commission errors
 Rem
 Rem    For the (> 1 t / 1 c) cases check for that tree, that
 Rem                                it is matched with a cluster closest to it, and mark other trees
 Rem                                missed.
 Rem    Loop thru the clusters:
 Rem    For the (>1 t / > 1c) we have clusters with multiple trees in their match-cylinder, trees
 Rem    which have multiple clusters in their matching cylinder.
 Rem
 Rem    Processed or NOT ?
 Rem
 
 Dim minXYZdist As Double, NinMatchCyl As Integer, NTreeInClusCyl As Integer, d3_dist As Double
 Rem Loop thru candidates within radius + border (N_cand_inside)
 Dim NCluswithzerotrees As Integer, NCluswithonetree As Integer, NCluswithmanytrees As Integer
 
  For j = 1 To N_cand_inside_border
    minXYZdist = 50: NTreeInClusCyl = 0
    xp = ClusMatchStruct(j).X
    yp = ClusMatchStruct(j).Y
    zp = ClusMatchStruct(j).z
     For jx = 1 To N_tree_inside_border
      
      d2_dist = sqrt((xp - TreeMatchStruct(jx).X) ^ 2 + (yp - TreeMatchStruct(jx).Y) ^ 2)
      d3_dist = sqrt((xp - TreeMatchStruct(jx).X) ^ 2 + (yp - TreeMatchStruct(jx).Y) ^ 2 + (zp - TreeMatchStruct(jx).z) ^ 2)
      If d2_dist < Match2Ddist And Abs(zp - TreeMatchStruct(jx).z) < MatchZdist Then
        NTreeInClusCyl = NTreeInClusCyl + 1
        ClusMatchStruct(j).Dist2dToTreesInCyl(NTreeInClusCyl) = d3_dist
        ClusMatchStruct(j).IndecesTreesInCyl(NTreeInClusCyl) = jx
        TreeMatchStruct(jx).NclusInCylinder = TreeMatchStruct(jx).NclusInCylinder + 1
        TreeMatchStruct(jx).IndecesClusInCyl(TreeMatchStruct(jx).NclusInCylinder) = j
        TreeMatchStruct(jx).Dist2dToClusInCyl(TreeMatchStruct(jx).NclusInCylinder) = d3_dist
       End If
       If d3_dist <= minXYZdist Then
        minXYZdist = d3_dist ' new minimum
        Rem At this point we store the minimum tree number and distance
        ClusMatchStruct(j).MatchTreeNum = TreeMatchStruct(jx).Num
        ClusMatchStruct(j).Dist3D = minXYZdist
       End If
     Rem We have also use for the number of trees in a clusters match-cylinder
     ClusMatchStruct(j).NTreesInCylinder = NTreeInClusCyl
    Next jx
    If NTreeInClusCyl = 0 Then NCluswithzerotrees = NCluswithzerotrees + 1
    If NTreeInClusCyl = 1 Then NCluswithonetree = NCluswithonetree + 1
    If NTreeInClusCyl > 1 Then NCluswithmanytrees = NCluswithmanytrees + 1
  Next j
  
 Dim Ntreeswithzeroclus As Integer, Ntreeswithoneclus As Integer, Ntreeswithmanyclus As Integer
 
  For jx = 1 To N_tree_inside_border
   If TreeMatchStruct(jx).IsInside = True Then
    If TreeMatchStruct(jx).NclusInCylinder = 0 Then Ntreeswithzeroclus = Ntreeswithzeroclus + 1
    If TreeMatchStruct(jx).NclusInCylinder = 1 Then Ntreeswithoneclus = Ntreeswithoneclus + 1
    If TreeMatchStruct(jx).NclusInCylinder > 1 Then Ntreeswithmanyclus = Ntreeswithmanyclus + 1
   End If
  Next jx
 Dim min_n As Integer, min_d As Double, min_h As Double, min_stat As Integer, min_spec As Integer
 Dim NTreesindorderZone As Integer, N_match_field_cand As Integer
 kX = 0: N_match_field_cand = 0: jx = 0
 

 
 Dim minXYdist As Double, min_x As Double, min_y As Double, min_z As Double
 Dim Nin As Integer, NClusWith1to1 As Integer
 Rem LOOP THRU TREES TO FIND 0/0 and 1/1 CASES
 Dim N_omis As Integer, N_match As Integer
  For kX = 1 To N_tree_inside_border
   If TreeMatchStruct(kX).IsInside = True Then
     If TreeMatchStruct(kX).NclusInCylinder = 0 Then
     Rem A SURE OMISSION ERROR
     TreeMatchStruct(kX).IsOmission = True
     N_omis = N_omis + 1
     GoTo NextTreeForCheck
    End If
    If TreeMatchStruct(kX).NclusInCylinder = 1 And ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).NTreesInCylinder = 1 Then
    Rem A SURE MATCH
     N_match = N_match + 1
     GoTo markamatch
    End If
    GoTo NextTreeForCheck
markamatch:
    TreeMatchStruct(kX).IndexMatchClus = TreeMatchStruct(kX).IndecesClusInCyl(1)
    TreeMatchStruct(kX).Dist3D = TreeMatchStruct(kX).Dist2dToClusInCyl(1)
    TreeMatchStruct(kX).IsMatch = True
    TreeMatchStruct(kX).Xclus = ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).X
    TreeMatchStruct(kX).Yclus = ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).Y
    TreeMatchStruct(kX).Zclus = ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).z
    Rem Recall to mark the cluster
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).IsMatch = True
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).IsCommission = False
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreeNum = TreeMatchStruct(kX).Num
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreeSpec = TreeMatchStruct(kX).Spec
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreeStat = TreeMatchStruct(kX).Status
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreed13 = TreeMatchStruct(kX).d13
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreeh = TreeMatchStruct(kX).h
    
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreeZbuttDem = TreeMatchStruct(kX).ZbuttDEM
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreeZbuttTach = TreeMatchStruct(kX).ZButtTach

    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).Xtree = TreeMatchStruct(kX).X
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).YTree = TreeMatchStruct(kX).Y
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).Ztree = TreeMatchStruct(kX).z
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).Dist3D = minXYZdist
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).Npoints = clus(TreeMatchStruct(kX).IndecesClusInCyl(1), 1)
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).Rvalue = atausum(TreeMatchStruct(kX).IndecesClusInCyl(1)) / clus(TreeMatchStruct(kX).IndecesClusInCyl(1), 1)
    NClusWith1to1 = NClusWith1to1 + 1
   End If
NextTreeForCheck:
Next kX
Dim op As Integer, apu As Long, ind As Long
op = N_cand_inside_border

Dim m As Integer, l As Integer, NT As Integer, NC As Integer, add As Integer
  For jx = 1 To N_tree_inside_border  ' Loop Trees inside, allow a match from border!
    
    If TreeMatchStruct(jx).NclusInCylinder > 1 And TreeMatchStruct(jx).IsMatch = False Then
      ' This is a tree that has several clusters in it's cylinder
      ' If jx = 103 Then MsgBox ("103")
      minXYZdist = 50
      apu = 0
      For j = 1 To TreeMatchStruct(jx).NclusInCylinder
       ind = TreeMatchStruct(jx).IndecesClusInCyl(j)
       apu = apu + ClusMatchStruct(ind).NTreesInCylinder
       ' If ClusMatchStruct(ind).Index = 50 Then MsgBox ("50")
       If TreeMatchStruct(jx).Dist2dToClusInCyl(j) < minXYZdist Then
        minXYZdist = TreeMatchStruct(jx).Dist2dToClusInCyl(j)
        min_n = TreeMatchStruct(jx).IndecesClusInCyl(j)
       End If
      Next j
      If apu = TreeMatchStruct(jx).NclusInCylinder Then
      Rem This tree has only clusters which on their part have just one tree, choose
      Rem the tree to attach the cluster with, it is the (min_n) -tree, make other trees missed!
      Rem There are N=apu clusters to choose from
      TreeMatchStruct(jx).IndexMatchClus = min_n
      TreeMatchStruct(jx).Dist3D = minXYZdist
      TreeMatchStruct(jx).IsMatch = True
      TreeMatchStruct(jx).Xclus = ClusMatchStruct(min_n).X
      TreeMatchStruct(jx).Yclus = ClusMatchStruct(min_n).Y
      TreeMatchStruct(jx).Zclus = ClusMatchStruct(min_n).z
      Rem Recall to mark the cluster
      ClusMatchStruct(min_n).IsMatch = True
      ClusMatchStruct(min_n).IsCommission = False
      ClusMatchStruct(min_n).MatchTreeNum = TreeMatchStruct(jx).Num
      ClusMatchStruct(min_n).MatchTreeSpec = TreeMatchStruct(jx).Spec
      ClusMatchStruct(min_n).MatchTreeStat = TreeMatchStruct(jx).Status
      ClusMatchStruct(min_n).MatchTreed13 = TreeMatchStruct(jx).d13
      ClusMatchStruct(min_n).MatchTreeh = TreeMatchStruct(jx).h
      ClusMatchStruct(min_n).MatchTreeZbuttDem = TreeMatchStruct(jx).ZbuttDEM
      ClusMatchStruct(min_n).MatchTreeZbuttTach = TreeMatchStruct(jx).ZButtTach
      ClusMatchStruct(min_n).Xtree = TreeMatchStruct(jx).X
      ClusMatchStruct(min_n).YTree = TreeMatchStruct(jx).Y
      ClusMatchStruct(min_n).Ztree = TreeMatchStruct(jx).z
      ClusMatchStruct(min_n).Dist3D = minXYZdist
      ClusMatchStruct(min_n).Npoints = clus(min_n, 1)
      ClusMatchStruct(min_n).Rvalue = atausum(min_n) / clus(min_n, 1)
      GoTo NextMTreeToCheck
      End If
      Rem Check if This cluster is used
      Rem We have now NT trees with > 2 clus
      If ClusMatchStruct(min_n).IsMatch = True Then MsgBox ("Duplicate!")
     
      TreeMatchStruct(jx).IndexMatchClus = min_n
      TreeMatchStruct(jx).Dist3D = minXYZdist
      TreeMatchStruct(jx).IsMatch = True
      TreeMatchStruct(jx).Xclus = ClusMatchStruct(min_n).X
      TreeMatchStruct(jx).Yclus = ClusMatchStruct(min_n).Y
      TreeMatchStruct(jx).Zclus = ClusMatchStruct(min_n).z
      Rem Recall to mark the cluster
      ClusMatchStruct(min_n).IsMatch = True
      ClusMatchStruct(min_n).IsCommission = False
      ClusMatchStruct(min_n).MatchTreeNum = TreeMatchStruct(jx).Num
      ClusMatchStruct(min_n).MatchTreeSpec = TreeMatchStruct(jx).Spec
      ClusMatchStruct(min_n).MatchTreeStat = TreeMatchStruct(jx).Status
      ClusMatchStruct(min_n).MatchTreed13 = TreeMatchStruct(jx).d13
      ClusMatchStruct(min_n).MatchTreeh = TreeMatchStruct(jx).h
      ClusMatchStruct(min_n).MatchTreeZbuttDem = TreeMatchStruct(jx).ZbuttDEM
      ClusMatchStruct(min_n).MatchTreeZbuttTach = TreeMatchStruct(jx).ZButtTach
      ClusMatchStruct(min_n).Xtree = TreeMatchStruct(jx).X
      ClusMatchStruct(min_n).YTree = TreeMatchStruct(jx).Y
      ClusMatchStruct(min_n).Ztree = TreeMatchStruct(jx).z
      ClusMatchStruct(min_n).Dist3D = minXYZdist
    End If
NextMTreeToCheck:
Next jx

 

  
Rem Now remains to polish the Clusters-commission errors for which there lie's
Rem a tree outside the circular plot area. These clusters must lie in the bor-
Rem dering zone: (radius-match2ddist) - (radius)

 Dim NborderMatches As Integer
  For j = 1 To N_cand_inside_border   ' loop clusters
    If ClusMatchStruct(j).NTreesInCylinder > 0 And ClusMatchStruct(j).IsMatch = False Then
    Rem We have a case with NON-matched cluster
     If TreeMatchStruct(ClusMatchStruct(j).IndecesTreesInCyl(1)).NclusInCylinder > 0 Then
     Rem A tree has been found for the border cluster
     ind = ClusMatchStruct(j).IndecesTreesInCyl(1)
     If TreeMatchStruct(ind).IsMatch = False Then GoTo matchthisclus
     If TreeMatchStruct(ind).IsMatch = True Then
     m = 1
incresem:
     If m < ClusMatchStruct(j).NTreesInCylinder Then
      m = m + 1
      ind = ClusMatchStruct(j).IndecesTreesInCyl(m)
      If TreeMatchStruct(ind).IsMatch = False Then GoTo matchthisclus
      GoTo incresem
     End If
     End If
     GoTo NextCand
matchthisclus:
     TreeMatchStruct(ind).IndexMatchClus = j
     TreeMatchStruct(ind).Dist3D = ClusMatchStruct(j).Dist3D
     TreeMatchStruct(ind).IsMatch = True
     TreeMatchStruct(ind).Xclus = ClusMatchStruct(j).X
     TreeMatchStruct(ind).Yclus = ClusMatchStruct(j).Y
     TreeMatchStruct(ind).Zclus = ClusMatchStruct(j).z
     ClusMatchStruct(j).IsMatch = True
     ClusMatchStruct(j).IsCommission = False
     ClusMatchStruct(j).MatchTreeNum = TreeMatchStruct(ind).Num
     ClusMatchStruct(j).MatchTreeSpec = TreeMatchStruct(ind).Spec
     ClusMatchStruct(j).MatchTreeStat = TreeMatchStruct(ind).Status
     ClusMatchStruct(j).MatchTreed13 = TreeMatchStruct(ind).d13
     ClusMatchStruct(j).MatchTreeh = TreeMatchStruct(ind).h
     ClusMatchStruct(j).MatchTreeZbuttDem = TreeMatchStruct(ind).ZbuttDEM
     ClusMatchStruct(j).MatchTreeZbuttTach = TreeMatchStruct(ind).ZButtTach
     ClusMatchStruct(j).Xtree = TreeMatchStruct(ind).X
     ClusMatchStruct(j).YTree = TreeMatchStruct(ind).Y
     ClusMatchStruct(j).Ztree = TreeMatchStruct(ind).z
     ClusMatchStruct(j).Dist3D = TreeMatchStruct(ind).Dist3D
     ClusMatchStruct(j).Npoints = clus(ClusMatchStruct(j).Index, 1)
     ClusMatchStruct(j).Rvalue = atausum(ClusMatchStruct(j).Index) / ClusMatchStruct(j).Npoints
     NborderMatches = NborderMatches + 1
     
    End If
   End If
   
NextCand:      Next j


DoEvents
   For i = 0 To NumOfImages - 1 ' Loop images
    Form1.Picture1(i).Cls
   Next i
DoEvents
Rem Let's plot matches, omissions and comissions
Dim p_x As Double, p_y As Double
Dim dXave As Double, dYave As Double, dZave As Double
Dim CountMatches As Integer, Countomissions As Integer
Dim ZmatchedmeanTrees As Double, ZmatchedmeanClus As Double
Dim Colorpix As Long
 CountMatches = 0: Countomissions = 0: ZmatchedmeanTrees = 0: ZmatchedmeanClus = 0
 dXave = 0: dYave = 0: dZave = 0:
    Colorpix = RGB(255, 255, 255)
  For j = 1 To N_tree_inside_border ' Loop thru trees
   If TreeMatchStruct(j).IsInside = True Then
   ZmatchedmeanTrees = ZmatchedmeanTrees + TreeMatchStruct(j).z
   Select Case TreeMatchStruct(j).IsMatch
   Case True
    CountMatches = CountMatches + 1
    dXave = dXave + (TreeMatchStruct(j).X - TreeMatchStruct(j).Xclus)
    dYave = dYave + (TreeMatchStruct(j).Y - TreeMatchStruct(j).Yclus)
    dZave = dZave + (TreeMatchStruct(j).z - TreeMatchStruct(j).Zclus)
   Case False
    Countomissions = Countomissions + 1
   End Select
   For i = 0 To NumOfImages - 1 ' Loop images
    Form1.Picture1(i).DrawWidth = 1
    
    Select Case TreeMatchStruct(j).IsMatch
    Case True

     Colorpix = RGB(0, 255, 41)
    Call r_transform_ground_to_pixel(i, TreeMatchStruct(j).Xclus, TreeMatchStruct(j).Yclus, TreeMatchStruct(j).Zclus, p_x, p_y)
    For l = -6 To 6
     For m = -6 To 6
     If Abs(m) > 4 Or Abs(l) > 4 Then
      Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1 + l, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1 + m), Colorpix
     End If
     Next m
    Next l
    Case False
    Colorpix = RGB(0, 10, 255)
    Call r_transform_ground_to_pixel(i, TreeMatchStruct(j).X, TreeMatchStruct(j).Y, TreeMatchStruct(j).z, p_x, p_y)
    For l = -1 To 1
     For m = -7 To 7
     ' If Abs(m) >= 0 And Abs(l) = 0 Then
      Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1 + l, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1 + m), Colorpix
     'End If
     Next m
    Next l
    For l = -7 To 7
     For m = -1 To 1
     ' If Abs(m) >= 0 And Abs(l) = 0 Then
      Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1 + l, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1 + m), Colorpix
     'End If
     Next m
    Next l
    For l = -7 To 7
     For m = -7 To 7
     If Abs(m) > 5 Or Abs(l) > 5 Then
      Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1 + l, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1 + m), Colorpix
     End If
     Next m
    Next l

    End Select
   Next i
   End If
  Next j

 dXave = dXave / CountMatches:  dYave = dYave / CountMatches:   dZave = dZave / CountMatches
 'Close (All)
 'Exit Sub
 Dim CountCommissions As Integer
 CountCommissions = 0
  For j = 1 To N_cand_inside_border ' Loop thru clusters
     If ClusMatchStruct(j).IsInside = True Then ZmatchedmeanClus = ZmatchedmeanClus + ClusMatchStruct(j).z
     If ClusMatchStruct(j).IsCommission = True And ClusMatchStruct(j).IsInside = True Then
      CountCommissions = CountCommissions + 1
     End If
   For i = 0 To NumOfImages - 1 ' Loop images
     If ClusMatchStruct(j).IsCommission = True And ClusMatchStruct(j).IsInside = True Then
     Colorpix = RGB(255, 0, 0)
     Form1.Picture1(i).DrawWidth = 1
     Call r_transform_ground_to_pixel(i, ClusMatchStruct(j).X, ClusMatchStruct(j).Y, ClusMatchStruct(j).z, p_x, p_y)
    For l = -7 To 7
     For m = -1 To 1
     ' If Abs(m) >= 0 And Abs(l) = 0 Then
      Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1 + l, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1 + m), Colorpix
     'End If
     Next m
    Next l
        For l = -1 To 1
         For m = -7 To 7
     ' If Abs(m) >= 0 And Abs(l) = 0 Then
      Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1 + l, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1 + m), Colorpix
     'End If
     Next m
    Next l

     Form1.Picture1(i).DrawWidth = 1
    End If
    Next i
  Next j
    
ZmatchedmeanTrees = ZmatchedmeanTrees / CDbl(N_tree_inside)
ZmatchedmeanClus = ZmatchedmeanClus / CDbl(N_cand_inside)
Dim Zbias As Double
Zbias = ZmatchedmeanTrees - ZmatchedmeanClus
Form1.Label10.Caption = "Cs: " & N_cand_inside & " Ts: " & N_tree_inside & " Hs: " & CountMatches & " Os: " & Countomissions & " Cs: " & CountCommissions & " AI-%: " & Format$(((N_tree_inside - CountCommissions - Countomissions) / N_tree_inside) * 100, "#.0") & "% Zbias: " & Format$(dZave, "#.00 m")
' MsgBox (Format$(dXave, "#.00 m") & " " & Format$(dYave, "#.00 m") & " " & Format$(dZave, "#.00 m"))

'OV.NtreesInside = N_tree_inside_border
'OV.NClusInside = N_cand_inside_border

OV.NtreesInside = N_tree_inside
OV.NClusInside = N_cand_inside
OV.Nmatched = CountMatches
OV.Nomission = Countomissions
OV.NCommission = CountCommissions
OV.MatchP = (CountMatches / N_tree_inside) * 100#
OV.Mrate = ((N_tree_inside - CountCommissions - Countomissions) / N_tree_inside) * 100
OV.ZbiasMatched = dZave
OV.ZbiasAll = Zbias
OV.Xbiasmatched = dXave
OV.Ybiasmatched = dYave

Rem ******************************************************************************
Rem RESULTS FILE-PRINTING SECTION; DigitString-holds a unique 8-digit for strorage
Rem Parameters are stored in global vars, or as locals of this procedure
Rem
Rem Rem Open in Append mode


Call PrintResults

Exit Sub
End Sub
 
 
Public Sub CheckAccuracyForPoly(DigitString As String, TreeFile As String)

Rem This Routine finds
Rem
Rem (1) Number of found trees; for a cluster there exists a match in the Data
Rem     - Trees are listed, and statistics computed
Rem
Rem (2) Number of commission errors; a cluster is not found a match, within the circular plot + buffer
Rem     - Clusters are listed, and statistics computed
Rem
Rem (3) Number of omission errors; a tree falling within the circular plot is not found a match
Rem     - Trees are listed, and statistics computed
Rem
Rem Arguments:
Rem
Rem DigitString     has the ######### -eight numbers
Rem TreeFile        contains measured tree-data (filename)
Rem Height-diameter curve parameters h = 1.3 + d^2/(a+b*d)^2
Rem  Open the file  with the ground truth
Rem  Clusters are in global storage: clus() -array



Rem *************************************************************************************************
Rem Find MATCHES -section; here we take a cluster, and try to find a match for it in the ground truth

Dim MatchZdist As Double, MatchZrmse As Double
Dim N_tree_inside  As Integer
Dim N_trees_in_plot As Integer
Dim MeanHeightField As Double
Dim N_tree_inside_border As Integer
Dim Match2Ddist As Double, MatchXYrmse As Double
ReDim ClusMatchStruct(1 To mc) As ClusMatchStruct

Dim i As Integer

Rem Maximum allowed XY-distance in meters between field measured tree top and candidate tree top
Match2Ddist = 1.3
Rem Maximum allowed XY-distance in meters between field measured tree top and candidate tree top
MatchZdist = 3

OV.XYMatchDist = Match2Ddist
OV.ZMatchDist = MatchZdist

Rem ***********************************************************************************
Rem *** First, lets read all trees, and count the number of trees within the plot area,
Rem *** and in the ring defined by the Match2Dist
Rem

TreeFile = "c:\data\maku_vis_trees.txt"
N_tree_inside_border = 0
N_tree_inside = 0
MeanHeightField = 0#
ReDim TreeMatchStruct(1 To 450) As TreeMatchStruct

Open TreeFile For Input As 44
'Close (44)
'Exit Sub
Dim angle As Double, X_origo As Double, Y_origo As Double, Z_origo As Double
Dim X_shift As Double, Y_shift As Double, Z_shift As Double, X As Double, Y As Double, z As Double
Dim sina As Double, cosa As Double, lx As Integer, jx As Integer, kX As Integer
Dim d2_dist As Double
Input #44, angle: Input #44, X_origo: Input #44, Y_origo: Input #44, Z_origo
Input #44, X_shift: Input #44, Y_shift: Input #44, Z_shift
angle = angle * TO_RADIANS: cosa = Cos(angle): sina = Sin(angle)
Dim piste As Point

 Do Until EOF(44)
  lx = lx + 1
  Input #44, X, Y, z, FTrees(lx).d13, FTrees(lx).Height, FTrees(lx).Num, FTrees(lx).Species, FTrees(lx).Status
  FTrees(lx).X = cosa * X - sina * Y + X_origo + X_shift
  FTrees(lx).Y = sina * X + cosa * Y + Y_origo + Y_shift
  FTrees(lx).Zbutt = z + Z_origo + Z_shift
  FTrees(lx).Ztop = FTrees(lx).Zbutt + FTrees(lx).Height
  'FTrees(lx).Zdem = getTINheight(FTrees(lx).x, FTrees(lx).y, CLng(0))
  FTrees(lx).Zdem = getheight(FTrees(lx).X, FTrees(lx).Y)
  Rem Check if the tree falls inside circular plot's XY-extent
  'd2_dist = sqrt((FTrees(lx).x - Plotcenter.x) ^ 2 + (FTrees(lx).y - Plotcenter.y) ^ 2)
  piste.X = FTrees(lx).X
  piste.Y = FTrees(lx).Y
  If InsidePolygon(KuvioRaj, UBound(KuvioRaj) - 1, piste) = Inside Then
  
  'If d2_dist <= (Plotradius + Match2Ddist) Then
   N_tree_inside_border = N_tree_inside_border + 1
   kX = N_tree_inside_border                             ' Number of trees inside the plot
   TreeMatchStruct(kX).Num = FTrees(lx).Num: TreeMatchStruct(kX).Spec = FTrees(lx).Species
'Exit Sub
   TreeMatchStruct(kX).Status = FTrees(lx).Status: TreeMatchStruct(kX).d13 = FTrees(lx).d13
   TreeMatchStruct(kX).h = FTrees(lx).Height: TreeMatchStruct(kX).X = FTrees(lx).X
   TreeMatchStruct(kX).Y = FTrees(lx).Y: TreeMatchStruct(kX).z = FTrees(lx).Ztop
   TreeMatchStruct(kX).ZbuttDEM = FTrees(lx).Zdem:  TreeMatchStruct(kX).ZButtTach = FTrees(lx).Zbutt
   TreeMatchStruct(kX).NclusInCylinder = 0
   TreeMatchStruct(kX).IndexMatchClus = 0: TreeMatchStruct(kX).Dist3D = -1#
   TreeMatchStruct(kX).IsMatch = False: TreeMatchStruct(kX).IsOmission = False
   TreeMatchStruct(kX).Xclus = 0: TreeMatchStruct(kX).Yclus = 0
   TreeMatchStruct(kX).Zclus = 0
   TreeMatchStruct(kX).IsInside = False
  End If
  Rem Check if the tree falls inside the plot border
'  If d2_dist < Plotradius Then
   N_tree_inside = N_tree_inside + 1      ' N of field trees within plotradius
   MeanHeightField = MeanHeightField + FTrees(lx).Ztop
   TreeMatchStruct(kX).IsInside = True
 ' End If
 Loop ' Get the next ground truth tree observation
Close (44)
N_trees_in_plot = N_tree_inside

Form1.Label10.Caption = "Field plot trees read " & lx & " trees"

Dim MeanHeightFoto As Double
MatchXYrmse = 0#: MeanHeightFoto = 0#: MatchZrmse = 0#
MeanHeightField = MeanHeightField / CDbl(N_tree_inside)
'Exit Sub
Rem Looping thru Clusters

Dim N_cand_inside As Integer
Dim N_cand_inside_border As Integer
Dim j As Integer
Dim xp As Double, yp As Double, zp As Double
 
 For i = 1 To mc
   xp = cluscoords(i).X / (atausum(i))
   yp = cluscoords(i).Y / (atausum(i))
   zp = cluscoords(i).z / (atausum(i))
   piste.X = xp
   piste.Y = yp
   ' d2_dist = ((xp - Plotcenter.x) ^ 2 + (yp - Plotcenter.y) ^ 2) ^ 0.5
     If InsidePolygon(KuvioRaj, UBound(KuvioRaj) - 1, piste) = Inside Then
   ' If d2_dist <= (Plotradius + Match2Ddist) Then
    N_cand_inside_border = N_cand_inside_border + 1
    j = N_cand_inside_border
    MeanHeightFoto = MeanHeightFoto + zp
    ClusMatchStruct(j).IsInside = False
    ClusMatchStruct(j).X = xp
    ClusMatchStruct(j).Y = yp
    ClusMatchStruct(j).z = zp
    ClusMatchStruct(j).Index = j
    ClusMatchStruct(j).IsMatch = False
    ClusMatchStruct(j).IsCommission = True
   End If
   'If d2_dist < (Plotradius) Then
    N_cand_inside = N_cand_inside + 1
    ClusMatchStruct(j).IsInside = True
   'End If
 Next i ' Next cluster found
 
 Rem The mean height (elevation) of the candidates within the plot
 MeanHeightFoto = MeanHeightFoto / CDbl(N_cand_inside)
 
 
 Rem ************************************************************************
 Rem *** FINDING MATCHES FOR THE FIELD TREES --> CANDIDATE WITHIN THE PLOT***
 Rem ************************************************************************
 
 Rem 1) For each tree and cluster, compute the number of clusters and Treetops in their
 Rem    Match-Cylinders (loop clusters and loop trees, O(N^2)), compute closest pairs (minimum distances)
 Rem
 Rem    Loop through Trees, On the first round Mark the with 0 c / 0 t and 1 c / 1 tree cases
 Rem    to be SURE omission-errors and Matches
 Rem
 Rem    At this point we are left with tree tops that have cases (1 c / > 1 t) and (>1 c / >1 t)
 Rem                                   clusters that have 0 trees in their Match-cylinder
 Rem
 Rem    Loop thru clusters (within boundary) and Mark 0 t -cases as commission errors
 Rem
 Rem    For the (> 1 t / 1 c) cases check for that tree, that
 Rem                                it is matched with a cluster closest to it, and mark other trees
 Rem                                missed.
 Rem    Loop thru the clusters:
 Rem    For the (>1 t / > 1c) we have clusters with multiple trees in their match-cylinder, trees
 Rem    which have multiple clusters in their matching cylinder.
 Rem
 Rem    Processed or NOT ?
 Rem
 
 Dim minXYZdist As Double, NinMatchCyl As Integer, NTreeInClusCyl As Integer, d3_dist As Double
 Rem Loop thru candidates within radius + border (N_cand_inside)
 Dim NCluswithzerotrees As Integer, NCluswithonetree As Integer, NCluswithmanytrees As Integer
 
  For j = 1 To N_cand_inside_border
    minXYZdist = 50: NTreeInClusCyl = 0
    xp = ClusMatchStruct(j).X
    yp = ClusMatchStruct(j).Y
    zp = ClusMatchStruct(j).z
     For jx = 1 To N_tree_inside_border
      
      d2_dist = sqrt((xp - TreeMatchStruct(jx).X) ^ 2 + (yp - TreeMatchStruct(jx).Y) ^ 2)
      d3_dist = sqrt((xp - TreeMatchStruct(jx).X) ^ 2 + (yp - TreeMatchStruct(jx).Y) ^ 2 + (zp - TreeMatchStruct(jx).z) ^ 2)
      If d2_dist < Match2Ddist And Abs(zp - TreeMatchStruct(jx).z) < MatchZdist Then
        NTreeInClusCyl = NTreeInClusCyl + 1
        ClusMatchStruct(j).Dist2dToTreesInCyl(NTreeInClusCyl) = d3_dist
        ClusMatchStruct(j).IndecesTreesInCyl(NTreeInClusCyl) = jx
        TreeMatchStruct(jx).NclusInCylinder = TreeMatchStruct(jx).NclusInCylinder + 1
        TreeMatchStruct(jx).IndecesClusInCyl(TreeMatchStruct(jx).NclusInCylinder) = j
        TreeMatchStruct(jx).Dist2dToClusInCyl(TreeMatchStruct(jx).NclusInCylinder) = d3_dist
       End If
       If d3_dist <= minXYZdist Then
        minXYZdist = d3_dist ' new minimum
        Rem At this point we store the minimum tree number and distance
        ClusMatchStruct(j).MatchTreeNum = TreeMatchStruct(jx).Num
        ClusMatchStruct(j).Dist3D = minXYZdist
       End If
     Rem We have also use for the number of trees in a clusters match-cylinder
     ClusMatchStruct(j).NTreesInCylinder = NTreeInClusCyl
    Next jx
    If NTreeInClusCyl = 0 Then NCluswithzerotrees = NCluswithzerotrees + 1
    If NTreeInClusCyl = 1 Then NCluswithonetree = NCluswithonetree + 1
    If NTreeInClusCyl > 1 Then NCluswithmanytrees = NCluswithmanytrees + 1
  Next j
  
 Dim Ntreeswithzeroclus As Integer, Ntreeswithoneclus As Integer, Ntreeswithmanyclus As Integer
 
  For jx = 1 To N_tree_inside_border
   If TreeMatchStruct(jx).IsInside = True Then
    If TreeMatchStruct(jx).NclusInCylinder = 0 Then Ntreeswithzeroclus = Ntreeswithzeroclus + 1
    If TreeMatchStruct(jx).NclusInCylinder = 1 Then Ntreeswithoneclus = Ntreeswithoneclus + 1
    If TreeMatchStruct(jx).NclusInCylinder > 1 Then Ntreeswithmanyclus = Ntreeswithmanyclus + 1
   End If
  Next jx
 Dim min_n As Integer, min_d As Double, min_h As Double, min_stat As Integer, min_spec As Integer
 Dim NTreesindorderZone As Integer, N_match_field_cand As Integer
 kX = 0: N_match_field_cand = 0: jx = 0
 

 
 Dim minXYdist As Double, min_x As Double, min_y As Double, min_z As Double
 Dim Nin As Integer, NClusWith1to1 As Integer
 Rem LOOP THRU TREES TO FIND 0/0 and 1/1 CASES
 Dim N_omis As Integer, N_match As Integer
  For kX = 1 To N_tree_inside_border
   If TreeMatchStruct(kX).IsInside = True Then
     If TreeMatchStruct(kX).NclusInCylinder = 0 Then
     Rem A SURE OMISSION ERROR
     TreeMatchStruct(kX).IsOmission = True
     N_omis = N_omis + 1
     GoTo NextTreeForCheck
    End If
    If TreeMatchStruct(kX).NclusInCylinder = 1 And ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).NTreesInCylinder = 1 Then
    Rem A SURE MATCH
     N_match = N_match + 1
     GoTo markamatch
    End If
    GoTo NextTreeForCheck
markamatch:
    TreeMatchStruct(kX).IndexMatchClus = TreeMatchStruct(kX).IndecesClusInCyl(1)
    TreeMatchStruct(kX).Dist3D = TreeMatchStruct(kX).Dist2dToClusInCyl(1)
    TreeMatchStruct(kX).IsMatch = True
    TreeMatchStruct(kX).Xclus = ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).X
    TreeMatchStruct(kX).Yclus = ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).Y
    TreeMatchStruct(kX).Zclus = ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).z
    Rem Recall to mark the cluster
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).IsMatch = True
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).IsCommission = False
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreeNum = TreeMatchStruct(kX).Num
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreeSpec = TreeMatchStruct(kX).Spec
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreeStat = TreeMatchStruct(kX).Status
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreed13 = TreeMatchStruct(kX).d13
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreeh = TreeMatchStruct(kX).h
    
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreeZbuttDem = TreeMatchStruct(kX).ZbuttDEM
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).MatchTreeZbuttTach = TreeMatchStruct(kX).ZButtTach

    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).Xtree = TreeMatchStruct(kX).X
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).YTree = TreeMatchStruct(kX).Y
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).Ztree = TreeMatchStruct(kX).z
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).Dist3D = minXYZdist
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).Npoints = clus(TreeMatchStruct(kX).IndecesClusInCyl(1), 1)
    ClusMatchStruct(TreeMatchStruct(kX).IndecesClusInCyl(1)).Rvalue = atausum(TreeMatchStruct(kX).IndecesClusInCyl(1)) / clus(TreeMatchStruct(kX).IndecesClusInCyl(1), 1)
    NClusWith1to1 = NClusWith1to1 + 1
   End If
NextTreeForCheck:
Next kX
Dim op As Integer, apu As Long, ind As Long
op = N_cand_inside_border

Dim m As Integer, l As Integer, NT As Integer, NC As Integer, add As Integer
  For jx = 1 To N_tree_inside_border  ' Loop Trees inside, allow a match from border!
    
    If TreeMatchStruct(jx).NclusInCylinder > 1 And TreeMatchStruct(jx).IsMatch = False Then
      ' This is a tree that has several clusters in it's cylinder
      ' If jx = 103 Then MsgBox ("103")
      minXYZdist = 50
      apu = 0
      For j = 1 To TreeMatchStruct(jx).NclusInCylinder
       ind = TreeMatchStruct(jx).IndecesClusInCyl(j)
       apu = apu + ClusMatchStruct(ind).NTreesInCylinder
       ' If ClusMatchStruct(ind).Index = 50 Then MsgBox ("50")
       If TreeMatchStruct(jx).Dist2dToClusInCyl(j) < minXYZdist Then
        minXYZdist = TreeMatchStruct(jx).Dist2dToClusInCyl(j)
        min_n = TreeMatchStruct(jx).IndecesClusInCyl(j)
       End If
      Next j
      If apu = TreeMatchStruct(jx).NclusInCylinder Then
      Rem This tree has only clusters which on their part have just one tree, choose
      Rem the tree to attach the cluster with, it is the (min_n) -tree, make other trees missed!
      Rem There are N=apu clusters to choose from
      TreeMatchStruct(jx).IndexMatchClus = min_n
      TreeMatchStruct(jx).Dist3D = minXYZdist
      TreeMatchStruct(jx).IsMatch = True
      TreeMatchStruct(jx).Xclus = ClusMatchStruct(min_n).X
      TreeMatchStruct(jx).Yclus = ClusMatchStruct(min_n).Y
      TreeMatchStruct(jx).Zclus = ClusMatchStruct(min_n).z
      Rem Recall to mark the cluster
      ClusMatchStruct(min_n).IsMatch = True
      ClusMatchStruct(min_n).IsCommission = False
      ClusMatchStruct(min_n).MatchTreeNum = TreeMatchStruct(jx).Num
      ClusMatchStruct(min_n).MatchTreeSpec = TreeMatchStruct(jx).Spec
      ClusMatchStruct(min_n).MatchTreeStat = TreeMatchStruct(jx).Status
      ClusMatchStruct(min_n).MatchTreed13 = TreeMatchStruct(jx).d13
      ClusMatchStruct(min_n).MatchTreeh = TreeMatchStruct(jx).h
      ClusMatchStruct(min_n).MatchTreeZbuttDem = TreeMatchStruct(jx).ZbuttDEM
      ClusMatchStruct(min_n).MatchTreeZbuttTach = TreeMatchStruct(jx).ZButtTach
      ClusMatchStruct(min_n).Xtree = TreeMatchStruct(jx).X
      ClusMatchStruct(min_n).YTree = TreeMatchStruct(jx).Y
      ClusMatchStruct(min_n).Ztree = TreeMatchStruct(jx).z
      ClusMatchStruct(min_n).Dist3D = minXYZdist
      ClusMatchStruct(min_n).Npoints = clus(min_n, 1)
      ClusMatchStruct(min_n).Rvalue = atausum(min_n) / clus(min_n, 1)
      GoTo NextMTreeToCheck
      End If
      Rem Check if This cluster is used
      Rem We have now NT trees with > 2 clus
      If ClusMatchStruct(min_n).IsMatch = True Then MsgBox ("Duplicate!")
     
      TreeMatchStruct(jx).IndexMatchClus = min_n
      TreeMatchStruct(jx).Dist3D = minXYZdist
      TreeMatchStruct(jx).IsMatch = True
      TreeMatchStruct(jx).Xclus = ClusMatchStruct(min_n).X
      TreeMatchStruct(jx).Yclus = ClusMatchStruct(min_n).Y
      TreeMatchStruct(jx).Zclus = ClusMatchStruct(min_n).z
      Rem Recall to mark the cluster
      ClusMatchStruct(min_n).IsMatch = True
      ClusMatchStruct(min_n).IsCommission = False
      ClusMatchStruct(min_n).MatchTreeNum = TreeMatchStruct(jx).Num
      ClusMatchStruct(min_n).MatchTreeSpec = TreeMatchStruct(jx).Spec
      ClusMatchStruct(min_n).MatchTreeStat = TreeMatchStruct(jx).Status
      ClusMatchStruct(min_n).MatchTreed13 = TreeMatchStruct(jx).d13
      ClusMatchStruct(min_n).MatchTreeh = TreeMatchStruct(jx).h
      ClusMatchStruct(min_n).MatchTreeZbuttDem = TreeMatchStruct(jx).ZbuttDEM
      ClusMatchStruct(min_n).MatchTreeZbuttTach = TreeMatchStruct(jx).ZButtTach
      ClusMatchStruct(min_n).Xtree = TreeMatchStruct(jx).X
      ClusMatchStruct(min_n).YTree = TreeMatchStruct(jx).Y
      ClusMatchStruct(min_n).Ztree = TreeMatchStruct(jx).z
      ClusMatchStruct(min_n).Dist3D = minXYZdist
    End If
NextMTreeToCheck:
Next jx

 

  
Rem Now remains to polish the Clusters-commission errors for which there lie's
Rem a tree outside the circular plot area. These clusters must lie in the bor-
Rem dering zone: (radius-match2ddist) - (radius)

 Dim NborderMatches As Integer
  For j = 1 To N_cand_inside_border   ' loop clusters
    If ClusMatchStruct(j).NTreesInCylinder > 0 And ClusMatchStruct(j).IsMatch = False Then
    Rem We have a case with NON-matched cluster
     If TreeMatchStruct(ClusMatchStruct(j).IndecesTreesInCyl(1)).NclusInCylinder > 0 Then
     Rem A tree has been found for the border cluster
     ind = ClusMatchStruct(j).IndecesTreesInCyl(1)
     If TreeMatchStruct(ind).IsMatch = False Then GoTo matchthisclus
     If TreeMatchStruct(ind).IsMatch = True Then
     m = 1
incresem:
     If m < ClusMatchStruct(j).NTreesInCylinder Then
      m = m + 1
      ind = ClusMatchStruct(j).IndecesTreesInCyl(m)
      If TreeMatchStruct(ind).IsMatch = False Then GoTo matchthisclus
      GoTo incresem
     End If
     End If
     GoTo NextCand
matchthisclus:
     TreeMatchStruct(ind).IndexMatchClus = j
     TreeMatchStruct(ind).Dist3D = ClusMatchStruct(j).Dist3D
     TreeMatchStruct(ind).IsMatch = True
     TreeMatchStruct(ind).Xclus = ClusMatchStruct(j).X
     TreeMatchStruct(ind).Yclus = ClusMatchStruct(j).Y
     TreeMatchStruct(ind).Zclus = ClusMatchStruct(j).z
     ClusMatchStruct(j).IsMatch = True
     ClusMatchStruct(j).IsCommission = False
     ClusMatchStruct(j).MatchTreeNum = TreeMatchStruct(ind).Num
     ClusMatchStruct(j).MatchTreeSpec = TreeMatchStruct(ind).Spec
     ClusMatchStruct(j).MatchTreeStat = TreeMatchStruct(ind).Status
     ClusMatchStruct(j).MatchTreed13 = TreeMatchStruct(ind).d13
     ClusMatchStruct(j).MatchTreeh = TreeMatchStruct(ind).h
     ClusMatchStruct(j).MatchTreeZbuttDem = TreeMatchStruct(ind).ZbuttDEM
     ClusMatchStruct(j).MatchTreeZbuttTach = TreeMatchStruct(ind).ZButtTach
     ClusMatchStruct(j).Xtree = TreeMatchStruct(ind).X
     ClusMatchStruct(j).YTree = TreeMatchStruct(ind).Y
     ClusMatchStruct(j).Ztree = TreeMatchStruct(ind).z
     ClusMatchStruct(j).Dist3D = TreeMatchStruct(ind).Dist3D
     ClusMatchStruct(j).Npoints = clus(ClusMatchStruct(j).Index, 1)
     ClusMatchStruct(j).Rvalue = atausum(ClusMatchStruct(j).Index) / ClusMatchStruct(j).Npoints
     NborderMatches = NborderMatches + 1
     
    End If
   End If
   
NextCand:      Next j


DoEvents
   For i = 0 To NumOfImages - 1 ' Loop images
    Form1.Picture1(i).Cls
   Next i
DoEvents
Rem Let's plot matches, omissions and comissions
Dim p_x As Double, p_y As Double
Dim dXave As Double, dYave As Double, dZave As Double
Dim CountMatches As Integer, Countomissions As Integer
Dim ZmatchedmeanTrees As Double, ZmatchedmeanClus As Double
Dim Colorpix As Long
 CountMatches = 0: Countomissions = 0: ZmatchedmeanTrees = 0: ZmatchedmeanClus = 0
 dXave = 0: dYave = 0: dZave = 0:
    Colorpix = RGB(255, 255, 255)
  For j = 1 To N_tree_inside_border ' Loop thru trees
   If TreeMatchStruct(j).IsInside = True Then
   ZmatchedmeanTrees = ZmatchedmeanTrees + TreeMatchStruct(j).z
   Select Case TreeMatchStruct(j).IsMatch
   Case True
    CountMatches = CountMatches + 1
    dXave = dXave + (TreeMatchStruct(j).X - TreeMatchStruct(j).Xclus)
    dYave = dYave + (TreeMatchStruct(j).Y - TreeMatchStruct(j).Yclus)
    dZave = dZave + (TreeMatchStruct(j).z - TreeMatchStruct(j).Zclus)
   Case False
    Countomissions = Countomissions + 1
   End Select
   For i = 0 To NumOfImages - 1 ' Loop images
    Form1.Picture1(i).DrawWidth = 1
    
    Select Case TreeMatchStruct(j).IsMatch
    Case True

   '  Colorpix = RGB(0, 255, 41)
    Call r_transform_ground_to_pixel(i, TreeMatchStruct(j).Xclus, TreeMatchStruct(j).Yclus, TreeMatchStruct(j).Zclus, p_x, p_y)
    For l = -6 To 6
     For m = -6 To 6
     If Abs(m) > 4 Or Abs(l) > 4 Then
      Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1 + l, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1 + m), Colorpix
     End If
     Next m
    Next l
    Case False

    Call r_transform_ground_to_pixel(i, TreeMatchStruct(j).X, TreeMatchStruct(j).Y, TreeMatchStruct(j).z, p_x, p_y)
    For l = -1 To 1
     For m = -7 To 7
     ' If Abs(m) >= 0 And Abs(l) = 0 Then
      Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1 + l, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1 + m), Colorpix
     'End If
     Next m
    Next l
    For l = -7 To 7
     For m = -1 To 1
     ' If Abs(m) >= 0 And Abs(l) = 0 Then
      Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1 + l, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1 + m), Colorpix
     'End If
     Next m
    Next l
    For l = -7 To 7
     For m = -7 To 7
     If Abs(m) > 5 Or Abs(l) > 5 Then
      Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1 + l, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1 + m), Colorpix
     End If
     Next m
    Next l

    End Select
   Next i
   End If
  Next j

 dXave = dXave / CountMatches:  dYave = dYave / CountMatches:   dZave = dZave / CountMatches
 'Close (All)
 'Exit Sub
 Dim CountCommissions As Integer
 CountCommissions = 0
  For j = 1 To N_cand_inside_border ' Loop thru clusters
     If ClusMatchStruct(j).IsInside = True Then ZmatchedmeanClus = ZmatchedmeanClus + ClusMatchStruct(j).z
     If ClusMatchStruct(j).IsCommission = True And ClusMatchStruct(j).IsInside = True Then
      CountCommissions = CountCommissions + 1
     End If
   For i = 0 To NumOfImages - 1 ' Loop images
     If ClusMatchStruct(j).IsCommission = True And ClusMatchStruct(j).IsInside = True Then
     Colorpix = RGB(255, 255, 255)
     Form1.Picture1(i).DrawWidth = 1
     Call r_transform_ground_to_pixel(i, ClusMatchStruct(j).X, ClusMatchStruct(j).Y, ClusMatchStruct(j).z, p_x, p_y)
    For l = -7 To 7
     For m = -1 To 1
     ' If Abs(m) >= 0 And Abs(l) = 0 Then
      Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1 + l, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1 + m), Colorpix
     'End If
     Next m
    Next l
        For l = -1 To 1
         For m = -7 To 7
     ' If Abs(m) >= 0 And Abs(l) = 0 Then
      Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1 + l, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - 1 + m), Colorpix
     'End If
     Next m
    Next l

     Form1.Picture1(i).DrawWidth = 1
    End If
    Next i
  Next j
    
ZmatchedmeanTrees = ZmatchedmeanTrees / CDbl(N_tree_inside)
ZmatchedmeanClus = ZmatchedmeanClus / CDbl(N_cand_inside)
Dim Zbias As Double
Zbias = ZmatchedmeanTrees - ZmatchedmeanClus
Form1.Label10.Caption = "Cs: " & N_cand_inside & " Ts: " & N_tree_inside & " Mcs: " & CountMatches & " msd: " & Countomissions & " Cms: " & CountCommissions & " M-Rate: " & Format$(((N_tree_inside - CountCommissions - Countomissions) / N_tree_inside) * 100, "#.0") & "% Z-bias: " & Format$(dZave, "#.00 m")
' MsgBox (Format$(dXave, "#.00 m") & " " & Format$(dYave, "#.00 m") & " " & Format$(dZave, "#.00 m"))

OV.NtreesInside = N_tree_inside_border
OV.NClusInside = N_cand_inside_border
OV.Nmatched = CountMatches
OV.Nomission = Countomissions
OV.NCommission = CountCommissions
OV.MatchP = (CountMatches / N_tree_inside) * 100#
OV.ZbiasMatched = dZave
OV.ZbiasAll = Zbias
OV.Xbiasmatched = dXave
OV.Ybiasmatched = dYave

Rem ******************************************************************************
Rem RESULTS FILE-PRINTING SECTION; DigitString-holds a unique 8-digit for strorage
Rem Parameters are stored in global vars, or as locals of this procedure
Rem
Rem Rem Open in Append mode


Call PrintResults


Exit Sub
 



End Sub
Public Sub PrintResults()
Rem Vector OV and arrays ClusMatchStruct, TreeMatchStruct are printed

Rem ** Compute here the Zer = a + b * TreeHeight
Dim Zer As Double, Xval As Integer, Yval As Integer
Dim VolMissed As Double, VolMatch As Double
Dim sx As Double, sy As Double, st2 As Double, aa As Double, bb As Double
Dim siga As Double, sigb As Double, ss As Double
Dim Vol As Double, XYer As Double, RMSEXY As Double, RMSEZ As Double
Dim i As Integer, k As Integer
Dim Ymean As Double

For i = 1 To OV.NtreesInside ' Loop thru trees inside the plot area
 
 Vol = 0.011197 * ((TreeMatchStruct(i).d13 / 10) ^ 2.10253) * (0.986 ^ (TreeMatchStruct(i).d13 / 10)) * (TreeMatchStruct(i).h ^ 3.98519) * ((TreeMatchStruct(i).h - 1.3) ^ (-2.659))
 Select Case TreeMatchStruct(i).Spec
 Case 1
 Vol = 0.036089 * ((TreeMatchStruct(i).d13 / 10) ^ 2.01395) * (0.99676 ^ (TreeMatchStruct(i).d13 / 10)) * (TreeMatchStruct(i).h ^ 2.07025) * ((TreeMatchStruct(i).h - 1.3) ^ (-1.07209))
 Case 2
 Vol = 0.022927 * ((TreeMatchStruct(i).d13 / 10) ^ 1.91505) * (0.99146 ^ (TreeMatchStruct(i).d13 / 10)) * (TreeMatchStruct(i).h ^ 2.82541) * ((TreeMatchStruct(i).h - 1.3) ^ (-1.53547))
 End Select
 
 If TreeMatchStruct(i).IsMatch = True And TreeMatchStruct(i).IsInside = True Then
    Zer = (TreeMatchStruct(i).z - TreeMatchStruct(i).Zclus)
    XYer = sqrt((TreeMatchStruct(i).Y - TreeMatchStruct(i).Yclus) ^ 2 + (TreeMatchStruct(i).X - TreeMatchStruct(i).Xclus) ^ 2)
    RMSEXY = RMSEXY + XYer
    RMSEZ = RMSEZ + (Zer ^ 2)
    sx = sx + TreeMatchStruct(i).h
    sy = sy + Zer
    Ymean = Ymean + Zer
    VolMatch = VolMatch + Vol
    k = k + 1
  End If
 If TreeMatchStruct(i).IsMatch = False And TreeMatchStruct(i).IsInside = True Then
    VolMissed = VolMissed + Vol
    k = k + 1
  End If
Next i

RMSEXY = sqrt(RMSEXY / (OV.Nmatched - 1))
'Exit Sub
RMSEZ = sqrt(RMSEZ / (OV.Nmatched - 1))

Dim sx0ss As Double, t As Double, sx2
sx0ss = sx / OV.Nmatched
Ymean = Ymean / OV.Nmatched
For i = 1 To OV.NtreesInside ' Loop thru trees inside the plot area
 If TreeMatchStruct(i).IsMatch = True And TreeMatchStruct(i).IsInside = True Then
  t = (TreeMatchStruct(i).h - sx0ss)
  st2 = st2 + t * t
  bb = bb + t * (TreeMatchStruct(i).z - TreeMatchStruct(i).Zclus)
End If
Next i

bb = bb / st2
aa = (sy - sx * bb) / OV.Nmatched
siga = sqrt((1# + sx * sx / (OV.Nmatched * st2)) / (OV.Nmatched))
sigb = sqrt((1 / st2) / (OV.Nmatched - 2))
Dim ssres As Double, ssreg As Double

For i = 1 To OV.NtreesInside ' Loop thru trees inside the plot area
 If TreeMatchStruct(i).IsMatch = True And TreeMatchStruct(i).IsInside = True Then
  ssres = ssres + ((TreeMatchStruct(i).z - TreeMatchStruct(i).Zclus) - (aa + bb * TreeMatchStruct(i).h)) ^ 2
  ssreg = ssreg + ((aa + bb * TreeMatchStruct(i).h) - Ymean) ^ 2
End If
Next i
Dim r2 As Double, sres As Double, seb As Double
r2 = ssreg / (ssreg + ssres)
sres = sqrt(ssres / (OV.Nmatched - 2))
seb = sres / sqrt(st2)


Open "c:\data\ma4res.txt" For Append As 1
Dim CO As String

CO = OV.PlotCode & ","
CO = CO & OV.Date & ","
CO = CO & OV.Time & ","
CO = CO & Format$(OV.NImagesInMatch, "0") & ","
CO = CO & Format$(OV.ModelTreeNum, "0000") & ","
CO = CO & Format$(OV.ModelTreeX, "#.00") & ","
CO = CO & Format$(OV.ModelTreeY, "#.00") & ","
CO = CO & Format$(OV.ModelTreeZ, "#.00") & ","
CO = CO & Format$(OV.Rlimit, "#.000") & ","
CO = CO & Format$(OV.XYthin, "#.00") & ","
CO = CO & Format$(OV.XYMatchDist, "#.00") & ","
CO = CO & Format$(OV.ZMatchDist, "#.00") & ","
CO = CO & Format$(OV.TWidth, "#.00") & ","
CO = CO & Format$(OV.EllZasym, "#.00") & ","
CO = CO & Format$(OV.Zdiff, "#.00") & ","
CO = CO & Format$(OV.ZDepth, "#.00") & ","
CO = CO & Format$(OV.Zasym, "#.00") & ","
CO = CO & Format$(OV.Meshdist, "#.00") & ","
CO = CO & Format$(OV.Meanheight, "#.00") & ","
CO = CO & Format$(OV.Plotradius, "#.00") & ","
CO = CO & Format$(OV.PlotCenterX, "#.00") & ","
CO = CO & Format$(OV.PlotCenterY, "#.00") & ","
CO = CO & Format$(OV.PlotCenterZ, "#.00") & ","
CO = CO & Format$(OV.DigString, "00000000") & ","
CO = CO & Format$(OV.NtreesInside, "000") & ","
'CO = CO & Format$(k, "000") & ","
CO = CO & Format$(OV.NClusInside, "000") & ","
CO = CO & Format$(OV.Nmatched, "000") & ","
CO = CO & Format$(OV.Nomission, "000") & ","
CO = CO & Format$(OV.NCommission, "000") & ","
CO = CO & Format$(OV.MatchP, "#.00") & ","
CO = CO & Format$(OV.Mrate, "#.00") & ","
CO = CO & Format$(OV.ZbiasMatched, "#.00") & ","
CO = CO & Format$(OV.ZbiasAll, "#.00") & ","
CO = CO & Format$(OV.Xbiasmatched, "#.00") & ","
CO = CO & Format$(OV.Ybiasmatched, "#.00") & ","
CO = CO & Format$(RMSEXY, "#.00") & ","
CO = CO & Format$(RMSEZ, "#.00") & ","
CO = CO & Format$(aa, "#.0000") & ","
CO = CO & Format$(bb, "#.0000") & ","
CO = CO & Format$(seb, "#.0000") & ","
CO = CO & Format$(r2, "#.0000") & ","
CO = CO & Format$(VolMatch, "#.0000") & ","
CO = CO & Format$(VolMissed, "#.0000")


Print #1, CO
Close (1)


Rem Check If the user wants the Cluster-data and tree-data to be stored in the OV.DigString-file
'If Form1.Check2.Value = 1 Then
Open "c:\data\" & OV.DigString & ".tr" For Output As 1


For i = 1 To OV.NtreesInside ' Loop thru trees inside the plot area
 CO = ""
 CO = CO & Format$(TreeMatchStruct(i).Num, "0000") & ","
 CO = CO & Format$(TreeMatchStruct(i).Spec, "00") & ","
 CO = CO & Format$(TreeMatchStruct(i).Status, "00") & ","
 CO = CO & Format$(TreeMatchStruct(i).d13, "000") & ","
 CO = CO & Format$(TreeMatchStruct(i).h, "#.00") & ","
 CO = CO & Format$(TreeMatchStruct(i).X, "#.00") & ","
 CO = CO & Format$(TreeMatchStruct(i).Y, "#.00") & ","
 CO = CO & Format$(TreeMatchStruct(i).z, "#.00") & ","
 CO = CO & Format$(TreeMatchStruct(i).ZbuttDEM, "#.00") & ","
 CO = CO & Format$(TreeMatchStruct(i).ZButtTach, "#.00") & ","
 CO = CO & Format$(TreeMatchStruct(i).NclusInCylinder, "000") & ","
 CO = CO & Format$(TreeMatchStruct(i).IndexMatchClus, "000") & ","
 CO = CO & Format$(TreeMatchStruct(i).Dist3D, "#.00") & ","
 CO = CO & Format$(Abs(TreeMatchStruct(i).IsMatch), "0") & ","
 CO = CO & Format$(Abs(TreeMatchStruct(i).IsOmission), "0") & ","
 CO = CO & Format$(Abs(TreeMatchStruct(i).IsInside), "0") & ","
 CO = CO & Format$(TreeMatchStruct(i).Xclus, "#.00") & ","
 CO = CO & Format$(TreeMatchStruct(i).Yclus, "#.00") & ","
 CO = CO & Format$(TreeMatchStruct(i).Zclus, "#.00") & ","
 Print #1, CO
Next i
Close (1)

Open "c:\data\" & OV.DigString & ".cl" For Output As 1
For i = 1 To mc ' Loop thru all clusters
CO = ""
CO = CO & Format$(ClusMatchStruct(i).Index, "000") & ","
CO = CO & Format$(Abs(ClusMatchStruct(i).IsMatch), "0") & ","
CO = CO & Format$(Abs(ClusMatchStruct(i).IsInside), "0") & ","
CO = CO & Format$(ClusMatchStruct(i).MatchTreeNum, "0000") & ","
CO = CO & Format$(ClusMatchStruct(i).MatchTreeSpec, "00") & ","
CO = CO & Format$(ClusMatchStruct(i).MatchTreeStat, "00") & ","
CO = CO & Format$(ClusMatchStruct(i).MatchTreed13, "000") & ","
CO = CO & Format$(ClusMatchStruct(i).MatchTreeh, "#.00") & ","
CO = CO & Format$(ClusMatchStruct(i).Xtree, "#.00") & ","
CO = CO & Format$(ClusMatchStruct(i).YTree, "#.00") & ","
CO = CO & Format$(ClusMatchStruct(i).Ztree, "#.00") & ","
CO = CO & Format$(ClusMatchStruct(i).MatchTreeZbuttDem, "#.00") & ","
CO = CO & Format$(ClusMatchStruct(i).MatchTreeZbuttTach, "#.00") & ","
CO = CO & Format$(Abs(ClusMatchStruct(i).IsCommission), "0") & ","
CO = CO & Format$(ClusMatchStruct(i).Dist3D, "#.00") & ","
CO = CO & Format$(ClusMatchStruct(i).NTreesInCylinder, "000") & ","
CO = CO & Format$(ClusMatchStruct(i).X, "#.00") & ","
CO = CO & Format$(ClusMatchStruct(i).Y, "#.00") & ","
CO = CO & Format$(ClusMatchStruct(i).z, "#.00") & ","
CO = CO & Format$(ClusMatchStruct(i).Rvalue, "#.000") & ","
CO = CO & Format$(ClusMatchStruct(i).Npoints, "000")
Print #1, CO
Next i
Close (1)

'End If

End Sub

Public Function Stringlength(ByRef A As String) As Integer
Dim i As Integer
For i = 0 To 100
 If (Mid$(A, i + 1, 1)) = "" Then
 Stringlength = i - 1
 Exit Function
 End If
Next i
End Function

Public Function sqrt(ByVal ar As Double) As Double
Rem square root
If ar < 0 Then
MsgBox ("Square root of negative argument!")
sqrt = 0.0001
Exit Function
End If
sqrt = ar ^ 0.5
End Function
Public Function vector_length(ByRef vector As Vector3D) As Double
Rem calculates euclidian lenght of a 3D vector
vector_length = (vector.X ^ 2 + vector.Y ^ 2 + vector.z ^ 2) ^ 0.5
End Function
Public Function vector_angle(vec1 As Vector3D, vec2 As Vector3D) As Double
Rem compute angle between two vectors
Rem angle=acos(inner_product / length*length)
vector_angle = MYFUNC_ACOS((vec1.X * vec2.X + vec1.Y * vec2.Y + vec1.z * vec2.z) / (vector_length(vec1) * vector_length(vec2)))
End Function
Public Function getCHMTINheight(Xm As Double, Ym As Double, ByRef i As Long) As Double
Dim p As Point

   If TINCHMModelReady = False Then
'     MsgBox ("Read the TIN model!")
   End If
   
   p.X = Xm
   p.Y = Ym

   Rem Check if previous Triangle is still valid
   If i > 0 Then
    apu = InsidePolygon(CHM_Tri(CHM_TRI_east_IND(i) - 1).poly(), 3, p)
    If apu = Inside Then
     getCHMTINheight = (1 - CHM_Tri(CHM_TRI_east_IND(i) - 1).A * (CHM_Tri(CHM_TRI_east_IND(i) - 1).xmin - p.X) - CHM_Tri(CHM_TRI_east_IND(i) - 1).B * (CHM_Tri(CHM_TRI_east_IND(i) - 1).Ymin - p.Y)) / CHM_Tri(CHM_TRI_east_IND(i) - 1).C
     GoTo CHMTIN_valmis
    End If
   End If

   Dim j As Long, StartInd As Long

'   For j = 1 To N_CHM_tri
'    If (CHM_Tri(CHM_TRI_east_IND(j) - 1).xmin - p.x) > -21 Then
'     StartInd = CHM_TRI_east_IND(j) - 1
'     GoTo FindCHMTIN
'    End If
'   Next j

FindCHMTIN:

 '  For i = StartInd To N_CHM_tri
   For i = 1 To N_CHM_tri
    apu = InsidePolygon(CHM_Tri(CHM_TRI_east_IND(i) - 1).poly(), 3, p)
    'Exit Function
    If apu = Inside Then
     getCHMTINheight = (1 - CHM_Tri(CHM_TRI_east_IND(i) - 1).A * (CHM_Tri(CHM_TRI_east_IND(i) - 1).xmin - p.X) - CHM_Tri(CHM_TRI_east_IND(i) - 1).B * (CHM_Tri(CHM_TRI_east_IND(i) - 1).Ymin - p.Y)) / CHM_Tri(CHM_TRI_east_IND(i) - 1).C
     GoTo CHMTIN_valmis
    End If
   Next i
CHMTIN_valmis:
If getCHMTINheight = 0 Then
  i = -99
'  MsgBox ("zero elevation")
 End If
End Function

Public Function getTIN_Optimized(Xm As Double, Ym As Double, ByRef i As Long) As Double
Dim p As Point

   getTIN_Optimized = 0
   
   p.X = Xm
   p.Y = Ym
   
   Rem Check bounding box
   If p.X < MinX1 Or p.X > MaxX1 Then
    getTIN_Optimized = -99
   End If
   If p.Y < MinY1 Or p.Y > MaxY1 Then
    getTIN_Optimized = -99
   End If
   
   Rem 1-hectare
   Dim FN1 As String * 3, FN2 As String * 3
   FN1 = Format$(Int((p.X - MinX1) / 100), "000")
   FN2 = Format$(Int((p.Y - MinY1) / 100), "000")

   Rem Check if previous Triangle is still valid
   Rem Recall indexing differences (1 to N_tri) and (0 to N_tri-1)
   Dim NHa As Long, j As Long
   
   If i > 0 And i < N_tri Then
    apu = InsidePolygon(Tri(i).poly(), 3, p)
    If apu = Inside Then
     getTIN_Optimized = (1 - Tri(i).A * (Tri(i).xmin - p.X) - Tri(i).B * (Tri(i).Ymin - p.Y)) / Tri(i).C
     Close (11)
     Exit Function
    End If
   End If
  
   
   capu = Dir("c:\data\Tris\" & FN1 & FN2 & ".bin")
   If capu <> "" Then
    Open "c:\data\Tris\" & FN1 & FN2 & ".bin" For Binary As 11
    Get #11, 1, NHa
     For k = 1 To NHa
     'Exit Function
      Get #11, 1 + k * 4, j
      apu = InsidePolygon(Tri(j).poly(), 3, p)
       If apu = Inside Then
         getTIN_Optimized = (1 - Tri(j).A * (Tri(j).xmin - p.X) - Tri(j).B * (Tri(j).Ymin - p.Y)) / Tri(j).C
         Close (11)
         Exit Function
       End If
      Next k
      Close (11)
     Rem We did not get elevation, check all
      For j = 1 To N_tri
      apu = InsidePolygon(Tri(j).poly(), 3, p)
       If apu = Inside Then
         getTIN_Optimized = (1 - Tri(j).A * (Tri(j).xmin - p.X) - Tri(j).B * (Tri(j).Ymin - p.Y)) / Tri(j).C
         Close (11)
         Exit Function
       End If
      Next j
    End If
      
 getTIN_Optimized = -99

End Function
Public Function getTINheight(Xm As Double, Ym As Double, ByRef i As Long) As Double
Dim p As Point

   getTINheight = 0
   
   If TINModelReady = False Then
'     MsgBox ("Read the TIN model!")
   End If
   
   p.X = Xm
   p.Y = Ym
   
   Rem Check if previous Triangle is still valid
   Rem Recall indexing differences (1 to N_tri) and (0 to N_tri-1)
   
   If i > 0 And i < N_tri Then
    apu = InsidePolygon(Tri(East_IND(i) - 1).poly(), 3, p)
    'Exit Function
    If apu = Inside Then
     getTINheight = (1 - Tri(East_IND(i) - 1).A * (Tri(East_IND(i) - 1).xmin - p.X) - Tri(East_IND(i) - 1).B * (Tri(East_IND(i) - 1).Ymin - p.Y)) / Tri(East_IND(i) - 1).C
     GoTo TIN_valmis
    End If
   End If
     
  StartInd = 1
  GoTo FindTIN
   
   Dim j As Long
   
   For j = 1 To N_tri Step 10
    If (Tri(East_IND(j) - 1).xmin - p.X) > -120 Then
   '  Exit Function
     StartInd = j
     GoTo FindTIN
    End If
   Next j
   
   Rem if we don't find it, start from 1

   
FindTIN:
   
   For i = StartInd To N_tri
  '   For i = 1 To N_tri
     apu = InsidePolygon(Tri(East_IND(i) - 1).poly(), 3, p)
    ' Exit Function
    If apu = Inside Then
     getTINheight = (1 - Tri(East_IND(i) - 1).A * (Tri(East_IND(i) - 1).xmin - p.X) - Tri(East_IND(i) - 1).B * (Tri(East_IND(i) - 1).Ymin - p.Y)) / Tri(East_IND(i) - 1).C
     GoTo TIN_valmis
    End If
   Next i
TIN_valmis:

If getTINheight = 0 Then
 getTINheight = -99
 i = 0
' MsgBox ("Zground -99 m")
End If
End Function


Public Sub MAKE3DPOINTSET(origoX As Double, origoY As Double, origoZ As Double, gridXextent As Double, gridYextent As Double, gridZextent As Double, gridXtess As Double, gridYtess As Double, gridZtess As Double, XYrotangle As Double, Zasymmetry As Double, Meanheight As Double)

Rem This is an alternative routine for creating the 3-D search space
Rem DATA is stored in an array (instead of writing it sequentally on disk file for subsequent read)

Rem Compute storage need, (Xextent/Xtess+1)*(Yextent/Ytess+1)*(Zextent/Ztess+1)

'On Error GoTo Error_In_MAKE3DPOINTSET

Dim Upperlimit As Long
Dim l As Long
Dim X As Double, Y As Double, z As Double
Dim XCELLS As Long, YCELLS As Long, ZCELLS As Long
XCELLS = CInt((gridXextent / gridXtess)) + 1
YCELLS = CInt((gridYextent / gridYtess)) + 1
ZCELLS = CInt((gridZextent / gridZtess)) + 1
Dim i As Long, j As Long

Rem This is the maximum storage

Upperlimit = XCELLS * YCELLS * ZCELLS ' CLng((CLng(gridXextent / gridXtess) + 1) * (CLng(gridYextent / gridYtess) + 1) * (CLng(gridZextent / gridZtess) + 1))
 
' If RasterModelReady = False Then
'  MsgBox ("Read height model first!")
'  Exit Sub
' End If


'Form1.Picture1(1).CurrentX = 5
'Form1.Picture1(1).CurrentY = 5
'Form1.Picture1(1).Print "GND + canopy points"

ReDim SearchSpaceData(0 To Upperlimit - 1) As point3D
Dim XYrotangle_cos As Double, XYrotangle_sin As Double

Dim fHDOM As Double
Dim HDepth As Double

fHDOM = CDbl(Form1.Text8.Text)
HDepth = CDbl(Form1.Text9.Text)
XYrotangle_cos = Cos(XYrotangle)
XYrotangle_sin = Sin(XYrotangle)
l = 0
m = 0
Dim piste As Point
Dim Xm As Double, Ym As Double, Hmean As Double
Dim lask As Long, Nlask As Long
Dim p_x As Double, p_y As Double



'Open "c:\depth.txt" For Append As 10
   For X = -gridXextent / 2# To (gridXextent / 2# + 0.001) Step gridXtess
'    lask = lask + 1
    For Y = -gridYextent / 2# To (gridYextent / 2# + 0.001) Step gridYtess
     Rem Take the MAX of CHM and the DTM, loop
     Xm = XYrotangle_cos * X - XYrotangle_sin * Y + origoX
     Ym = XYrotangle_sin * X + XYrotangle_cos * Y + origoY
     
      
     'Minz = getTINheight(Xm, Ym, i)
     'Maxz = getCHMTINheight(Xm, Ym, j)
 'Exit Sub
     
     Minz = getheight(Xm, Ym)
     maxZ = getCHMheight(Xm, Ym)
     'maxZ = Minz + 22
     maxZ = Minz + fHDOM * (maxZ - Minz)
'

'     If lask Mod 10 = 0 Then
 '    For i = 0 To NumOfImages - 1
  '    Call r_transform_ground_to_pixel(i, Xm, Ym, MinZ, p_x, p_y)
  '    Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y), RGB(255, 255, 255)
  '    DoEvents
  '   Next i
  '  End If
 
'     Hmean = 24
'     MaxZ = Minz + fHDOM * Hmean

     
     
     
    ' If MaxZ - Minz < 13 Then
    '  MaxZ = Minz + 10
    ' End If
     
     'Exit Sub
     'Print #10, Format$(Xm, "#.00") & "," & Format$(Ym, "#.00") & "," & Format$(Maxz, "#.00") & "," & Format$(Minz, "#.00")
     
     For z = Minz + (maxZ - Minz) * HDepth To maxZ Step gridZtess
    '  Nlask = Nlask + 1
      'SearchSpaceData(l).z = getheight(SearchSpaceData(l).x, SearchSpaceData(l).y) + Meanheight + Zasymmetry + z
      Rem
      SearchSpaceData(l).X = Xm
      SearchSpaceData(l).Y = Ym
      SearchSpaceData(l).z = z
    ' If lask Mod 10 = 0 Then
    '  For i = 0 To NumOfImages - 1
    '  Call r_transform_ground_to_pixel(i, Xm, Ym, z, p_x, p_y)
'      Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y), RGB(255, 255, 0)
'      If Nlask Mod 100 = 0 Then
       'SavePicture Form1.Picture1(1).Image, "c:\temp\AVI\" & Format$(lask, "000000") & ".bmp"
'      End If
    '  Next i
    '  DoEvents
     
    ' End If
      
'      If lask Mod 10 = 0 And Nlask Mod 100 = 0 Then
'       SavePicture Form1.Picture1(1).Image, "c:\temp\AVI\" & Format$(Nlask, "000000") & ".bmp"
'      End If

      
      If l = UBound(SearchSpaceData) - 1 Then
       ReDim Preserve SearchSpaceData(0 To UBound(SearchSpaceData) + 10000) As point3D
      End If
      l = l + 1
     Next z
    Next Y
    Next X
 Close (10)
 ReDim Preserve SearchSpaceData(0 To l - 1) As point3D

 
 'If l <> Upperlimit Then MsgBox (l & "  ? = ?   " & Upperlimit)
 Form1.Label10.Caption = "3D-Point mesh creation succesfully done"
 'MsgBox (SearchSpaceData(66372).y)
Exit Sub
Error_In_MAKE3DPOINTSET:
MsgBox ("An error occurred in trying to create 3D search space point mesh!")
Form1.Label10.Caption = "3D-Point mesh creation failed"
End Sub


Public Sub Draw_Epipolar_line()
Dim Xvll As Double, Yvll As Double, Zvll As Double, Zv As Double
Dim p_x As Double, p_y As Double
Dim i As Integer, j As Integer, m As Integer, l As Integer
Dim X1 As Double, X2 As Double, Y1 As Double, Y2 As Double, Xr As Double, Yr As Double, Zr As Double
Dim r11 As Double, r12 As Double, r13 As Double, r21 As Double, r22 As Double, r23 As Double, r31 As Double, r32 As Double, r33 As Double


If LastImageClicked = -1 Or SolutionExists = False Then
 Form1.Label10.Caption = "Cannot draw epipolar line, either solution is lacking or LastImageClicked = -1"
 Exit Sub
End If
j = LastImageClicked

r11 = A(1, 1, j): r12 = A(1, 2, j): r13 = A(1, 3, j)
r21 = A(2, 1, j): r22 = A(2, 2, j): r23 = A(2, 3, j)
r31 = A(3, 1, j): r32 = A(3, 2, j): r33 = A(3, 3, j)
X1 = CDbl(Form1.Text1(j * 2).Text)
Y1 = CDbl(Form1.Text1(j * 2 + 1).Text)
Xvll = X_sol
Yvll = Y_sol
Zvll = Z_sol

For Zv = Zvll - Epi_Depth To Zvll + Epi_Depth Step 0.2
 Xr = image_info(j).Xo + (Zv - image_info(j).Zo) * ((X1 * r11 + Y1 * r12 - r13 * image_info(j).C) / (X1 * r31 + Y1 * r32 - r33 * image_info(j).C))
 Yr = image_info(j).Yo + (Zv - image_info(j).Zo) * ((X1 * r21 + Y1 * r22 - r23 * image_info(j).C) / (X1 * r31 + Y1 * r32 - r33 * image_info(j).C))
 Zr = Zv

For i = 0 To NumOfImages - 1
 Call r_transform_ground_to_pixel(i, Xr, Yr, Zr, p_x, p_y)
 'Form1.Picture1(i).DrawWidth = 3
 Form1.Picture1(i).DrawWidth = 1
 If Imagesdisplayed(i) = CORR_IMA Then
  Form1.Picture1(i).PSet ((p_x - (cor_ima_info(i).o_col + cor_win_info(i).win_o_col)) * cor_win_info(i).pan_x, ((cor_ima_info(i).Height - 1) - p_y - (cor_ima_info(i).o_row + cor_win_info(i).win_o_row)) * cor_win_info(i).pan_y - l), RGB(0, 255, 0)
 End If
 If Imagesdisplayed(i) = AER_IMA Then
  'Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x, -1 + ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y), RGB(0, 255, 255)
  Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y), RGB(0, 255, 255)
  'Form1.Picture1(i).PSet ((p_x - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x, ((image_info(i).Height - 1) - p_y - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y - l), RGB(0, 255, 255)
 End If
Next i
Next Zv

Rem 1.6. koe
 ' NCounter = NCounter + 1
 ' SavePicture Form1.Picture1(3).Image, "c:\temp\AVI\" & Format$(NCounter, "0000") & ".bmp"



End Sub
Public Sub aZoom(ByVal i As Integer, ByVal pan_x As Double, ByVal pan_y As Double)

Dim startcol As Integer
Dim startrow As Integer
Dim endcol As Integer
Dim endrow As Integer
Dim length As Integer, j As Integer
Dim X As Double, Y As Double, x_apu As Double, y_apu As Double
Dim p_x As Double, p_y As Double
On Error GoTo ErrorInaZoom

If Imagesdisplayed(i) = CORR_IMA Then
 MsgBox ("Zoom not alowed for correlation images!")
 Exit Sub
End If

X = CDbl(Form1.Text1(i * 2).Text): Y = CDbl(Form1.Text1(i * 2 + 1).Text)
'Exit Sub
Call a_transform_affine(i, 0, X, Y, p_x, p_y)

x_apu = p_x - image_info(i).o_col
y_apu = (image_info(i).Height - 1) - p_y - (image_info(i).o_row)

win_info(i).pan_x = pan_x
win_info(i).pan_y = pan_y
win_info(i).win_width = Win_w
win_info(i).win_height = win_h

' Call create_bmp(i, image_info(i).sub_width, image_info(i).sub_height, kuva(), "c:\temp\pic" & CStr(i) & ".bmp")
' calculate window area, define start_col, start_row, end_col, end_row
Call calculate_region(i, CInt(x_apu), CInt(y_apu), CInt((1 / win_info(i).pan_x) * Win_w), CInt((1 / win_info(i).pan_y) * win_h), startcol, startrow, endcol, endrow, win_info(i).pan_x, win_info(i).pan_y)

If image_info(i).Color = 1 Then

Rem Wrapper for filenames

 length = Stringlength(image_info(i).FileName)
 ReDim filename_in(0 To length) As Byte
 For j = 0 To length - 1
  filename_in(j) = CByte(Asc(Mid$(image_info(i).FileName, j + 1, 1)))
 Next j
 filename_in(length) = 0
 
 FileOut = "c:\data\pic" & CStr(i) & ".bmp"
 length = Len(FileOut)
 ReDim filename_out(0 To length) As Byte
 For j = 0 To length - 1
  filename_out(j) = CByte(Asc(Mid$(FileOut, j + 1, 1)))
 Next j
 filename_out(length) = 0
  
 Call create_bmp(i, startcol, startrow, endcol, endrow, "c:\data\pic" & CStr(i) & ".bmp", win_info(i).pan_x, win_info(i).pan_y)
' Call create_bmp_header(FileOut, CLng(Win_w), CLng(win_h))
' apu = MYFUNC_CREATEBMP(CLng(i), CLng(startcol), CLng(startrow), CLng(endcol), CLng(endrow), filename_in(0), filename_out(0), CDbl(win_info(i).pan_x), CDbl(win_info(i).pan_y), CLng(win_h), CLng(Win_w), CLng(image_info(i).sub_width))
ElseIf image_info(i).Color = 0 Then
 Call create_BW_bmp(i, startcol, startrow, endcol, endrow, "c:\data\pic" & CStr(i) & ".bmp", win_info(i).pan_x, win_info(i).pan_y)
End If

' Call Create_Bitmap(i, StartCol, StartRow, EndCol, EndRow, kuva(), "c:\temp\pic" & CStr(i) & ".bmp", win_info(i).pan_x, win_info(i).pan_y)
Form1.Picture1(i).Picture = LoadPicture("c:\data\pic" & CStr(i) & ".bmp")
Form1.Label1(i).Caption = Form1.Label1(i).Caption & Format$(win_info(i).pan_x, "#.00") & " " & Format$(win_info(i).win_o_col, "#") & " " & Format$(win_info(i).win_o_row, "#")
DoEvents
Exit Sub
ErrorInaZoom:
MsgBox ("Error in aZoom-function. Image observation missing? (xy-coords. Click a ref. point first)")
Exit Sub

End Sub
Public Sub c_corA(ima1() As Byte, ima2() As Byte, r As Double)
Dim os As Double, nim1 As Double, nim2 As Double, meanA As Double, meanB As Double
os = 0#
nim1 = 0#
nim2 = 0#
meanA = 0#
meanB = 0#
Dim i As Long, j As Long, N As Long
N = 0

 For i = 0 To UBound(ima1, 1)
   For j = 0 To UBound(ima1, 2)
   If ima1(i, j) <> 0 Then
     meanA = meanA + ima1(i, j)
     meanB = meanB + ima2(i, j)
     N = N + 1
   End If
   Next j
 Next i

  meanA = meanA / (N)
  meanB = meanB / (N)
  
 For i = 0 To UBound(ima1, 1)
   For j = 0 To UBound(ima1, 2)
   If ima1(i, j) <> 0 Then
    os = os + (CDbl(ima1(i, j)) - meanA) * (CDbl(ima2(i, j)) - meanB)
    nim1 = nim1 + (CDbl(ima1(i, j)) - meanA) * (CDbl(ima1(i, j)) - meanA)
    nim2 = nim2 + (CDbl(ima2(i, j)) - meanB) * (CDbl(ima2(i, j)) - meanB)
   End If
   Next j
 Next i
   
   If ((nim1 * nim2) < 0.002) Then
     r = -2#
     Exit Sub
   End If
   r = 1 * os / Sqr(nim1 * nim2)

End Sub

Public Sub c_cor(ima1() As RGBtriplet, ima2() As RGBtriplet, r As Double)

Dim os As Double, nim1 As Double, nim2 As Double, meanA As Double, meanB As Double
os = 0#
nim1 = 0#
nim2 = 0#
meanA = 0#
meanB = 0#
Dim i As Long, j As Long

If UBound(ima1, 1) <> UBound(ima2, 1) Then
 MsgBox ("Templates are of unequal size!")
End If

 For i = 0 To UBound(ima1, 1)
   For j = 0 To UBound(ima1, 2)
    meanA = meanA + ima1(i, j).g
    meanB = meanB + ima2(i, j).g
   Next j
 Next i

  meanA = meanA / (UBound(ima1, 1) * UBound(ima1, 2))
  meanB = meanB / (UBound(ima2, 1) * UBound(ima2, 2))
  
 For i = 0 To UBound(ima1, 1)
   For j = 0 To UBound(ima1, 2)
   os = os + (CDbl(ima1(i, j).g) - meanA) * (CDbl(ima2(i, j).g) - meanB)
   nim1 = nim1 + (CDbl(ima1(i, j).g) - meanA) * (CDbl(ima1(i, j).g) - meanA)
   nim2 = nim2 + (CDbl(ima2(i, j).g) - meanB) * (CDbl(ima2(i, j).g) - meanB)
   Next j
 Next i
   
   If ((nim1 * nim2) < 0.002) Then
     r = -2#
     Exit Sub
   End If
   r = 1 * os / Sqr(nim1 * nim2)
 

End Sub

Public Sub create_BW_bmp(ByVal Index As Integer, ByVal startcol As Integer, ByVal startrow As Integer, ByVal endcol As Integer, ByVal endrow As Integer, FileName As String, ByVal pan_x As Double, ByVal pan_y As Double)
Rem Tihs routine will read a greyscale image & make it into a 24-color greyscale BMP-image for display
Dim byteapu As Byte
Dim hederi As BITMAPINFOHEADER ' 40 bytes
Dim varitaulu(0 To 255) As RGBQUAD
Dim bmicolor As RGBQUAD
Dim isohederi As BITMAPINFO
Dim filehederi As BITMAPFILEHEADER ' 14 bytes
Dim i As Integer, j As Integer, k As Integer, l As Integer, m As Integer, N As Integer, lisa As Integer
Dim W_byte_array As Integer, H_byte_array As Integer
Dim i_step As Integer
Dim j_step As Integer
Dim isum As Integer, apusum As Integer, maxh As Integer
Dim paikka As Long

H_byte_array = win_h
W_byte_array = Win_w

'On Error GoTo error_in_creating_bmp_file
Close (1)
Open FileName For Binary As 1

'Exit Sub
' Debugging parameter-values received
' MsgBox ("W_byte_array: " & CStr(W_byte_array))
' MsgBox ("H_byte_array: " & CStr(H_byte_array))
' MsgBox ("Filename: " & CStr(filename))
Dim padding As Integer
lisa = (W_byte_array * 3) Mod 4
If lisa = 0 Then
 padding = 0
 Else
 padding = 4 - lisa
End If


Rem Note the order of BMP-files ! It is not RGB but BGR!

ReDim bkuva(0 To (W_byte_array - 1) + (padding), 0 To H_byte_array - 1) As RGBQ
Rem ckuva is used (temporary) to make sure we don't accidentally go outside OK indeces
ReDim ckuva(-30 To (W_byte_array - 1) + 30 + (padding), -30 To H_byte_array - 1 + 30) As RGBQ


filehederi.bfType = &H4D42
filehederi.bfSize = (CLng(W_byte_array + padding) * CLng(H_byte_array) * CLng(3))
filehederi.bfReserved1 = &H0
filehederi.bfReserved2 = &H0
' filehederi.bfOffBits = &H43A ' (if a color table) 43A=1082
' Offset from the BitMapFileheader to the Bitmap Bits
filehederi.bfOffBits = &H36 ' 55 (54 byte offset, read 55th as first data-byte)

Rem aking Dib - bitmap , header is defined here
hederi.biSize = CLng(40)  ' 4 bytes
hederi.biWidth = CLng(W_byte_array) ' 4
hederi.biHeight = CLng(-H_byte_array) ' 4
Rem By making Height negative we assume (WINDOWS programs assume) it's a BMP with Low-left corner on first scan-line!!
hederi.biPlanes = CInt(1) ' 2
hederi.biBitCount = CInt(24) '2
hederi.biCompression = CLng(0) '4
hederi.biSizeImage = CLng(0) '4
hederi.biXPelsPerMeter = CLng(0) '4
hederi.biYPelsPerMeter = CLng(0) '4
hederi.biClrUsed = CLng(0) '4
hederi.biClrImportant = CLng(0) '4


isohederi.bmiHeader = hederi
' isohederi.bmiColors = bmicolor


' MsgBox (Seek(1)) ' 14 Bytes
Put #1, , filehederi
' MsgBox (Seek(1))
Put #1, , hederi   ' 40 Bytes
' Put #1, , isohederi ' 40 Bytes + size of rgbquad 4

' MsgBox (Seek(1)) ' 1078 (1078+4)= 1082! = H43A
' Put #1, , varitaulu

Rem Skip the Color Table!

' Flip Image upside-down for proper BMP-syntax
' MsgBox ("b_array has  a value: " & CStr(b_array(233, 4)))
' Resize image assuming pan_x = pan_y image N*N
' Flip upside down at the same time

If pan_x > 0.99 And pan_x < 1.1 Then GoTo pan_x_1
If pan_x < 1 Then GoTo ZoomOut
If pan_x > 1 Then GoTo ZoomIn

ZoomOut:
Rem This is the N --> 1 case


ReDim RROW(0 To (endcol - startcol)) As Byte
Open image_info(Index).FileName For Binary As 2

Dim lx As Integer

k = -1
lx = -CInt(1 / pan_y)
For i = 0 To H_byte_array - 1
   lx = lx + CInt(1 / pan_y)
   paikka = CLng(startrow + lx) * (image_info(Index).sub_width) * CLng(1) + startcol * 1
   Get #2, paikka + 1, RROW()
   l = -CInt(1 / pan_x)
    For j = 0 To W_byte_array - 1
     l = l + CInt(1 / pan_x)
     m = startcol + i * CInt(1 / pan_x)
     N = startrow + j * CInt(1 / pan_y)
     If m < 0 Or N < 0 Or m > image_info(Index).sub_width Or N > image_info(Index).sub_height Then
       MsgBox ("In create BMP: exceeding allowed indeces")
       bkuva(i, j).r = 255: bkuva(i, j).g = 255: bkuva(i, j).B = 255:
    GoTo Ohita
    End If
     bkuva(j, i).r = RROW(l)
     bkuva(j, i).g = RROW(l)
     bkuva(j, i).B = RROW(l)
   
Ohita:
  Next j
 Next i
Close (2)
GoTo loppu

ZoomIn:
   
Rem The 1->N Case

ReDim RROW(0 To (endcol - startcol)) As Byte
Open image_info(Index).FileName For Binary As 2
   
i_step = -1
   
For i = startrow To endrow
   i_step = i_step + 1
   j_step = -1
   paikka = CLng(i) * (image_info(Index).sub_width) * 1 + startcol * 1
   Get #2, paikka + 1, RROW()
Rem Parse it
For j = startcol To endcol
  j_step = j_step + 1
    For k = 0 To CInt(pan_x) - 1
     For l = 0 To CInt(pan_y) - 1
      ckuva(j_step * CInt(pan_x) + l, i_step * CInt(pan_y) + k).r = RROW(j_step)
      ckuva(j_step * CInt(pan_x) + l, i_step * CInt(pan_y) + k).g = RROW(j_step)
      ckuva(j_step * CInt(pan_x) + l, i_step * CInt(pan_y) + k).B = RROW(j_step)
     Next l
    Next k
  Next j
 Next i
   
For i = 0 To Win_w - 1
 For j = 0 To win_h - 1
   bkuva(i, j).r = ckuva(i, j).r
   bkuva(i, j).g = ckuva(i, j).g
   bkuva(i, j).B = ckuva(i, j).B
  Next j
Next i
Close (2)
GoTo loppu

pan_x_1:

ReDim RROW(0 To (endcol - startcol)) As Byte
Open image_info(Index).FileName For Binary As 2
k = -1
'Exit Sub
 For i = startrow To endrow
   k = k + 1
   paikka = CLng(i) * (image_info(Index).sub_width) * 1 + startcol * 1
   Rem Read a row
   Get #2, paikka + 1, RROW()
    For m = 0 To (endcol - startcol)
     bkuva(m, k).r = RROW(m)
     bkuva(m, k).g = RROW(m)
     bkuva(m, k).B = RROW(m)
    Next m
  Next i
Close (2)
GoTo loppu


loppu: Put #1, , bkuva
Close (1)

Exit Sub

error_in_creating_bmp_file:
Close (1)
Close (2)
MsgBox ("Error in creating the bmp-file!")
Exit Sub


End Sub


Public Sub create_bmp_header(ByVal File As String, ByVal Win_w As Long, ByVal win_h As Long)
Dim hederi As BITMAPINFOHEADER ' 40 bytes
Dim filehederi As BITMAPFILEHEADER ' 14 bytes
Dim W_byte_array As Integer, H_byte_array As Integer
Dim lisa As Integer, padding As Integer

H_byte_array = win_h
W_byte_array = Win_w

Open File For Binary As 1

lisa = (W_byte_array * 3) Mod 4
If lisa = 0 Then
 padding = 0
 Else
 padding = 4 - lisa
End If
filehederi.bfType = &H4D42
filehederi.bfSize = (CLng(W_byte_array + padding) * CLng(H_byte_array) * CLng(3))
filehederi.bfReserved1 = &H0
filehederi.bfReserved2 = &H0
' filehederi.bfOffBits = &H43A ' (if a color table) 43A=1082
' Offset from the BitMapFileheader to the Bitmap Bits
filehederi.bfOffBits = &H36 ' 55 (54 byte offset, read 55th as first data-byte)

Rem aking Dib - bitmap , header is defined here
hederi.biSize = CLng(40)  ' 4 bytes
hederi.biWidth = CLng(W_byte_array) ' 4
hederi.biHeight = CLng(-H_byte_array) ' 4
Rem By making Height negative we assume (WINDOWS programs assume) it's a BMP with Low-left corner on first scan-line!!
hederi.biPlanes = CInt(1) ' 2
hederi.biBitCount = CInt(24) '2
hederi.biCompression = CLng(0) '4
hederi.biSizeImage = CLng(0) '4
hederi.biXPelsPerMeter = CLng(0) '4
hederi.biYPelsPerMeter = CLng(0) '4
hederi.biClrUsed = CLng(0) '4
hederi.biClrImportant = CLng(0) '4

' MsgBox (Seek(1)) ' 14 Bytes
Put #1, , filehederi
' MsgBox (Seek(1))
Put #1, , hederi   ' 40 Bytes
Close (1)

End Sub
Public Sub normalize(ByRef vector As Vector3D)
Rem Returns a normalized version of the supplied 3D-vector
Dim length As Double
Dim apu
length = vector_length(vector)
vector.X = vector.X / length
vector.Y = vector.Y / length
vector.z = vector.z / length
End Sub

Public Sub old_create_bmp(Shift As Long, W_byte_array As Long, H_byte_array As Long, b_array() As Byte, FileName As String)

Rem This subroutine receives a byte array (intensity image 8-bits / pixel) and calculates then need for byte padding
Rem (4-byte segments / line ) and Creates the BMP-file in filename
Rem

Dim byteapu As Byte
Dim hederi As BITMAPINFOHEADER
Dim varitaulu(0 To 255) As RGBQUAD
Dim isohederi As BITMAPINFO
Dim filehederi As BITMAPFILEHEADER
Dim i As Integer, j As Integer, lisa As Integer

Rem Error Handling
' On Error GoTo error_in_creating_bmp_file

Rem Open the BMP-file

Open FileName For Binary As 1
Dim padding As Integer
Rem Compute need for posiible padding (of zeros at scan line ends)
lisa = (W_byte_array Mod 4)
If lisa = 0 Then
 padding = 0
 Else
 padding = 4 - lisa
End If

Rem Declare storage to hold a copy of the image
ReDim bkuva(1 To W_byte_array + (padding), H_byte_array) As Byte

Rem Prepare the BMP-file's HEADER

For i = 0 To 255
 varitaulu(i).rgbblue = CByte(i)
 varitaulu(i).rgbGreen = CByte(i)
 varitaulu(i).rgbred = CByte(i)
 varitaulu(i).rgbreserved = CByte(0)
Next i

filehederi.bfType = &H4D42           ' 19778 in decimal
filehederi.bfSize = (W_byte_array * H_byte_array)
filehederi.bfReserved1 = &H0
filehederi.bfReserved2 = &H0
filehederi.bfOffBits = &H43A         ' 1082 BYTE offset!

hederi.biSize = CLng(40)             ' 4 bytes
hederi.biWidth = CLng(W_byte_array)  ' 4
hederi.biHeight = CLng(H_byte_array) ' 4
hederi.biPlanes = CInt(1)            ' 2
hederi.biBitCount = CInt(8)          ' 2
hederi.biCompression = CLng(0)       ' 4
hederi.biSizeImage = CLng(0)         ' 4
hederi.biXPelsPerMeter = CLng(0)     ' 4
hederi.biYPelsPerMeter = CLng(0)     ' 4
hederi.biClrUsed = CLng(256)         ' 4
hederi.biClrImportant = CLng(256)    ' 4

isohederi.bmiHeader = hederi

Rem Write the HEADER of the BMP-file

Put #1, , filehederi
Put #1, , isohederi
Put #1, , varitaulu

Rem Flip the image upside-down for. Store this version in bkuva()-array

Rem Since b_array is indexed starting from (0,0), make shift 1

'MsgBox (Seek(1))

Shift = 1

For i = 1 To W_byte_array
 For j = 1 To H_byte_array
  bkuva(i, H_byte_array + 1 - j) = b_array(i - Shift, j - Shift)
 Next j
Next i

Rem Write this array to the BMP-file

loppu: Put #1, , bkuva

Rem Close the file
Close (1)

Exit Sub

error_in_creating_bmp_file:
Close (1)
MsgBox ("Error in creating the bmp-file! ( sub create_bmp() )")
Exit Sub

End Sub
Public Sub create_orto_bmp(orthoimage() As RGBtriplet)

Dim byteapu As Byte
Dim hederi As BITMAPINFOHEADER ' 40 bytes
Dim varitaulu(0 To 255) As RGBQUAD
Dim bmicolor As RGBQUAD
Dim isohederi As BITMAPINFO
Dim filehederi As BITMAPFILEHEADER ' 14 bytes
Dim i As Integer, j As Integer, k As Integer, l As Integer, m As Integer, N As Integer, lisa As Integer
Dim W_byte_array As Integer, H_byte_array As Integer
Dim i_step As Integer
Dim j_step As Integer, padding As Integer, lx As Integer
Dim isum As Integer, apusum As Integer, maxh As Integer
Dim paikka As Long

'orthoimage(0 To N_ortho_columns - 1, 0 To N_ortho_rows - 1)

H_byte_array = UBound(orthoimage, 2) + 1
W_byte_array = UBound(orthoimage, 1) + 1
'Exit Sub
' On Error GoTo error_in_creating_bmp_file
Close (1)
Open "C:\Data\Orto.bmp" For Binary As 1

lisa = (W_byte_array * 3) Mod 4
If lisa = 0 Then
 padding = 0
 Else
 padding = 4 - lisa
End If

Rem Note the order of BMP-files ! It is not RGB but BGR!
ReDim b_array(0 To (W_byte_array - 1) + (padding), 0 To H_byte_array - 1) As RGBtriplet
ReDim bkuva(0 To (W_byte_array - 1) + (padding), 0 To H_byte_array - 1) As RGBQ
Rem ckuva is used (temporary) to make sure we don't accidentally go outside OK indeces
ReDim ckuva(-100 To (W_byte_array - 1) + 100 + (padding), -100 To H_byte_array - 1 + 100) As RGBQ

filehederi.bfType = &H4D42
filehederi.bfSize = (CLng(W_byte_array + padding) * CLng(H_byte_array) * CLng(3))
filehederi.bfReserved1 = &H0
filehederi.bfReserved2 = &H0
' Offset from the BitMapFileheader to the Bitmap Bits
filehederi.bfOffBits = &H36 ' 55 (54 byte offset, read 55th as first data-byte)
Rem aking Dib - bitmap , header is defined here
hederi.biSize = CLng(40)  ' 4 bytes
hederi.biWidth = CLng(W_byte_array) ' 4
hederi.biHeight = CLng(-H_byte_array) ' 4
Rem By making Height negative we assume (WINDOWS programs assume) it's a BMP with Low-left corner on first scan-line!!
hederi.biPlanes = CInt(1) ' 2
hederi.biBitCount = CInt(24) '2
hederi.biCompression = CLng(0) '4
hederi.biSizeImage = CLng(0) '4
hederi.biXPelsPerMeter = CLng(0) '4
hederi.biYPelsPerMeter = CLng(0) '4
hederi.biClrUsed = CLng(0) '4
hederi.biClrImportant = CLng(0) '4

isohederi.bmiHeader = hederi
' isohederi.bmiColors = bmicolor
' MsgBox (Seek(1)) ' 14 Bytes
Put #1, , filehederi
' MsgBox (Seek(1))
Put #1, , hederi   ' 40 Bytes
' Put #1, , isohederi ' 40 Bytes + size of rgbquad 4

' MsgBox (Seek(1)) ' 1078 (1078+4)= 1082! = H43A
' Put #1, , varitaulu

Rem Skip the Color Table!
pan_x_1:

k = -1
 For i = 0 To UBound(orthoimage, 1)
   k = k + 1
   Rem Read a row
    For m = 0 To UBound(orthoimage, 2)
     bkuva(m, k).r = orthoimage(i, m).r
     bkuva(m, k).g = orthoimage(i, m).g
     bkuva(m, k).B = orthoimage(i, m).B
    Next m
  Next i
GoTo loppu


loppu: Put #1, , bkuva
Close (1)

Exit Sub

error_in_creating_bmp_file:
Close (1)
MsgBox ("Error in creating the bmp-file!")
Exit Sub

End Sub
Public Sub create_bmp(ByVal Index As Integer, ByVal startcol As Integer, ByVal startrow As Integer, ByVal endcol As Integer, ByVal endrow As Integer, FileName As String, ByVal pan_x As Double, ByVal pan_y As Double)

Rem Declarations Needed for BMP-printing


Rem This subroutine receives a byte array (intensity image 8-bits / pixel)
Rem Calculates need for byte padding (4-byte segments / line )
Rem Creates the file// E:\ilkka_smoothed.bmp  //
Rem

Dim byteapu As Byte
Dim hederi As BITMAPINFOHEADER ' 40 bytes
Dim varitaulu(0 To 255) As RGBQUAD
Dim bmicolor As RGBQUAD
Dim isohederi As BITMAPINFO
Dim filehederi As BITMAPFILEHEADER ' 14 bytes
Dim i As Integer, j As Integer, k As Integer, l As Integer, m As Integer, N As Integer, lisa As Integer
Dim W_byte_array As Integer, H_byte_array As Integer
Dim i_step As Integer
Dim j_step As Integer, padding As Integer, lx As Integer
Dim isum As Integer, apusum As Integer, maxh As Integer
Dim paikka As Long

H_byte_array = win_h
W_byte_array = Win_w

' On Error GoTo error_in_creating_bmp_file
Close (1)
Open FileName For Binary As 1


lisa = (W_byte_array * 3) Mod 4
If lisa = 0 Then
 padding = 0
 Else
 padding = 4 - lisa
End If

Rem Note the order of BMP-files ! It is not RGB but BGR!
ReDim b_array(0 To (W_byte_array - 1) + (padding), 0 To H_byte_array - 1) As RGBtriplet
ReDim bkuva(0 To (W_byte_array - 1) + (padding), 0 To H_byte_array - 1) As RGBQ
Rem ckuva is used (temporary) to make sure we don't accidentally go outside OK indeces
ReDim ckuva(-100 To (W_byte_array - 1) + 100 + (padding), -100 To H_byte_array - 1 + 100) As RGBQ

filehederi.bfType = &H4D42
filehederi.bfSize = (CLng(W_byte_array + padding) * CLng(H_byte_array) * CLng(3))
filehederi.bfReserved1 = &H0
filehederi.bfReserved2 = &H0
' filehederi.bfOffBits = &H43A ' (if a color table) 43A=1082
' Offset from the BitMapFileheader to the Bitmap Bits
filehederi.bfOffBits = &H36 ' 55 (54 byte offset, read 55th as first data-byte)

Rem aking Dib - bitmap , header is defined here
hederi.biSize = CLng(40)  ' 4 bytes
hederi.biWidth = CLng(W_byte_array) ' 4
hederi.biHeight = CLng(-H_byte_array) ' 4
Rem By making Height negative we assume (WINDOWS programs assume) it's a BMP with Low-left corner on first scan-line!!
hederi.biPlanes = CInt(1) ' 2
hederi.biBitCount = CInt(24) '2
hederi.biCompression = CLng(0) '4
hederi.biSizeImage = CLng(0) '4
hederi.biXPelsPerMeter = CLng(0) '4
hederi.biYPelsPerMeter = CLng(0) '4
hederi.biClrUsed = CLng(0) '4
hederi.biClrImportant = CLng(0) '4


isohederi.bmiHeader = hederi
' isohederi.bmiColors = bmicolor
Dim LOWL As Integer, HIGHL As Integer
LOWL = 12
HIGHL = 21
' MsgBox (Seek(1)) ' 14 Bytes
Put #1, , filehederi
' MsgBox (Seek(1))
Put #1, , hederi   ' 40 Bytes
' Put #1, , isohederi ' 40 Bytes + size of rgbquad 4

' MsgBox (Seek(1)) ' 1078 (1078+4)= 1082! = H43A
' Put #1, , varitaulu

Rem Skip the Color Table!

' Flip Image upside-down for proper BMP-syntax
' MsgBox ("b_array has  a value: " & CStr(b_array(233, 4)))
' Resize image assuming pan_x = pan_y image N*N
' Flip upside down at the same time

'MsgBox (EndCol - StartCol + 1 & " Cols")
'MsgBox (EndRow - StartRow + 1 & " Rows")

Rem For applying contrast enhancement by means of GAmma Correction
Dim gamma As Double
Dim Gcorr As Double
Dim ApplyConstrastEnhancement As Boolean
ApplyConstrastEnhancement = False
 If Form1.Gamma_1.Checked = True Then
  gamma = 0.4
  Gcorr = (255 / (255# ^ gamma))
  ApplyConstrastEnhancement = True
 End If
 If Form1.Gamma_2.Checked = True Then
  gamma = 0.5
  Gcorr = (255 / (255# ^ gamma))
  ApplyConstrastEnhancement = True
 End If
 If Form1.Gamma_3.Checked = True Then
  gamma = 0.6
  Gcorr = (255 / (255# ^ gamma))
  ApplyConstrastEnhancement = True
 End If
 If Form1.Gamma_4.Checked = True Then
  gamma = 1.4
  Gcorr = (255# / (255# ^ gamma))
  ApplyConstrastEnhancement = True
 End If
 If Form1.Gamma_5.Checked = True Then
  gamma = 2
  Gcorr = (255# / (255# ^ gamma))
  ApplyConstrastEnhancement = True
 End If
 If Form1.Gamma_6.Checked = True Then
  gamma = 3
  Gcorr = (255# / (255# ^ gamma))
  ApplyConstrastEnhancement = True
 End If

If pan_x > 0.99 And pan_x < 1.1 Then GoTo pan_x_1
If pan_x < 1 Then GoTo ZoomOut
If pan_x > 1 Then GoTo ZoomIn

ZoomOut:
Rem This is the N --> 1 case

ReDim RROW(0 To (endcol - startcol)) As RGBtriplet
Rem Testing contrast enhancement
'ReDim RROWa(0 To (endcol - startcol), 0 To H_byte_array - 1) As RGBtriplet

Open image_info(Index).FileName For Binary As 2

k = -1
lx = -CInt(1 / pan_y)
For i = 0 To H_byte_array - 1
   lx = lx + CInt(1 / pan_y)
   paikka = CLng(startrow + lx) * (image_info(Index).sub_width) * CLng(3) + CLng(startcol) * 3
   Get #2, paikka + 1, RROW()
   l = -CInt(1 / pan_x)
    For j = 0 To W_byte_array - 1
     l = l + CInt(1 / pan_x)
     m = startcol + i * CInt(1 / pan_x)
     N = startrow + j * CInt(1 / pan_y)
    If m < 0 Or N < 0 Or m > image_info(Index).sub_width Or N > image_info(Index).sub_height Then
'       MsgBox ("In create BMP: exceeding allowed indeces")
       bkuva(i, j).r = 255: bkuva(i, j).g = 255: bkuva(i, j).B = 255:
    GoTo Ohita
    End If
    If ApplyConstrastEnhancement = False Then
     bkuva(j, i).r = RROW(l).r
     bkuva(j, i).g = RROW(l).g
     bkuva(j, i).B = RROW(l).B
    End If
   If ApplyConstrastEnhancement = True Then
     bkuva(j, i).r = CByte(CDbl(RROW(l).r) ^ (gamma) * Gcorr)
     bkuva(j, i).g = CByte(CDbl(RROW(l).g) ^ (gamma) * Gcorr)
     bkuva(j, i).B = CByte(CDbl(RROW(l).B) ^ (gamma) * Gcorr)
   End If
Ohita:
  Next j
 Next i
Close (2)

GoTo loppu

ZoomIn:
   
Rem The 1->N Case

ReDim RROW(0 To (endcol - startcol)) As RGBtriplet

'Open image_info(Index).filename For Binary As 2
Close (2)
Open image_info(Index).FileName For Binary As 2
'Exit Sub
i_step = -1
   
For i = startrow To endrow
   i_step = i_step + 1
   j_step = -1
   paikka = CLng(i) * CLng(image_info(Index).sub_width) * 3 + CLng(startcol) * 3
   Get #2, paikka + 1, RROW()
Rem Parse it
For j = startcol To endcol
  j_step = j_step + 1
    For k = 0 To CInt(pan_x) - 1
     For l = 0 To CInt(pan_y) - 1
      ckuva(j_step * CInt(pan_x) + l, i_step * CInt(pan_y) + k).r = RROW(j_step).r
      ckuva(j_step * CInt(pan_x) + l, i_step * CInt(pan_y) + k).g = RROW(j_step).g
      ckuva(j_step * CInt(pan_x) + l, i_step * CInt(pan_y) + k).B = RROW(j_step).B
     Next l
    Next k
  Next j
 Next i
   
For i = 0 To Win_w - 1
 For j = 0 To win_h - 1
 '   If ckuva(i, j).g > LOWL And ckuva(i, j).g < HIGHL Then
 '        bkuva(i, j).R = (Log(((ckuva(i, j).g - LOWL) / HIGHL + 0.01)) + 4.6) * 55
 '        bkuva(i, j).g = (Log(((ckuva(i, j).g - LOWL) / HIGHL + 0.01)) + 4.6) * 55
 '        bkuva(i, j).b = (Log(((ckuva(i, j).g - LOWL) / HIGHL + 0.01)) + 4.6) * 55
 '   Else
 '        bkuva(i, j).R = 0
 '        bkuva(i, j).g = 0
 '        bkuva(i, j).b = 0
 ' End If
    If ApplyConstrastEnhancement = False Then
     bkuva(i, j).r = ckuva(i, j).r
     bkuva(i, j).g = ckuva(i, j).g
     bkuva(i, j).B = ckuva(i, j).B
    End If
   If ApplyConstrastEnhancement = True Then
     bkuva(i, j).r = CByte(CDbl(ckuva(i, j).r) ^ (gamma) * Gcorr)
     bkuva(i, j).g = CByte(CDbl(ckuva(i, j).g) ^ (gamma) * Gcorr)
     bkuva(i, j).B = CByte(CDbl(ckuva(i, j).B) ^ (gamma) * Gcorr)
   End If
  Next j
Next i
Close (2)
GoTo loppu

pan_x_1:

ReDim RROW(0 To (endcol - startcol)) As RGBtriplet
Close (2)
Open image_info(Index).FileName For Binary As 2
'Exit Sub
k = -1
 For i = startrow To endrow
   k = k + 1
   paikka = CLng(i) * CLng(image_info(Index).sub_width) * 3 + CLng(startcol) * 3
   Rem Read a row
   Get #2, paikka + 1, RROW()
    For m = 0 To (endcol - startcol)
    If ApplyConstrastEnhancement = False Then
     bkuva(m, k).r = RROW(m).r
     bkuva(m, k).g = RROW(m).g
     bkuva(m, k).B = RROW(m).B
   End If
   If ApplyConstrastEnhancement = True Then
     bkuva(m, k).r = CByte(CDbl(RROW(m).r) ^ (gamma) * Gcorr)
     bkuva(m, k).g = CByte(CDbl(RROW(m).g) ^ (gamma) * Gcorr)
    bkuva(m, k).B = CByte(CDbl(RROW(m).B) ^ (gamma) * Gcorr)
   End If
 
 
 '   If RROW(m).g > LOWL And RROW(m).g < HIGHL Then
 '        bkuva(m, k).R = (Log(((RROW(m).g - LOWL) / HIGHL + 0.01)) + 4.6) * 55
 '        bkuva(m, k).g = (Log(((RROW(m).g - LOWL) / HIGHL + 0.01)) + 4.6) * 55
 '        bkuva(m, k).b = (Log(((RROW(m).g - LOWL) / HIGHL + 0.01)) + 4.6) * 55
 '   Else
 '        bkuva(m, k).R = 0
 '        bkuva(m, k).g = 0
 '        bkuva(m, k).b = 0
 ' End If
Rem Code above for measurement of fiducial marks with high contrast enhancement!
    
    
    Next m
  Next i
Close (2)
GoTo loppu


loppu: Put #1, , bkuva
Close (1)

Exit Sub

error_in_creating_bmp_file:
Close (1)
MsgBox ("Error in creating the bmp-file!")
Exit Sub


End Sub


Public Sub create_vector(ByRef vector As Vector3D, ByVal dimx As Double, ByVal dimy As Double, ByVal dimz As Double)
Rem Assigns values & returns the vector
vector.X = dimx
vector.Y = dimy
vector.z = dimz
End Sub

Public Sub Create_Bitmap(ByVal Index As Integer, ByVal startcol As Integer, ByVal startrow As Integer, ByVal endcol As Integer, ByVal endrow As Integer, b_array() As Byte, FileName As String, ByVal pan_x As Double, ByVal pan_y As Double)

Rem Creates a Bitmap and displays it using SetBitMapPixels API-function

Dim bm As BITMAP

Rem Retrieve values
GetObject Form1.Picture1(Index).Image, Len(bm), bm

ReDim bkuva(0 To bm.bmWidth - 1, 0 To bm.bmHeight - 1) As RGBQUAD1
Dim H_byte_array As Integer, W_byte_array As Integer
Dim i As Integer, j As Integer, k As Integer, m As Integer, N As Integer

H_byte_array = win_h
W_byte_array = Win_w



If pan_x > 0.99 And pan_x < 1.1 Then GoTo pan_x_1
If pan_x < 1 Then GoTo ZoomOut
If pan_x > 1 Then GoTo ZoomIn

ZoomOut:
k = -1
For i = 0 To W_byte_array - 1
 For j = 0 To H_byte_array - 1
   m = startcol + i * CInt(1 / pan_x) + (1 / pan_x) / 2 'Step
   N = startrow + j * CInt(1 / pan_y) + CInt(1 / pan_y) / 2
   'isum = 0
   'apusum = 0
   'For k = (-1 / pan_x) / 2 To ((1 / pan_x) - 1) / 2
   '  For l = (-1 / pan_x) / 2 To ((1 / pan_x) - 1) / 2
   '    apusum = apusum + CInt(b_array(StartCol + m + k, EndRow - n + l, Index))
   '    isum = isum + 1
   '  Next l
   ' Next k
   bkuva(i, j).r = b_array(m, N, Index)
   bkuva(i, j).g = b_array(m, N, Index)
   bkuva(i, j).B = b_array(m, N, Index)
   bkuva(i, j).N = 0
  Next j
 Next i
   
GoTo loppu
Dim i_step As Integer, lisaysc As Integer, lisaysr As Integer, j_step As Integer
Dim l As Integer

ZoomIn:
k = -1
i_step = -1
lisaysc = (endcol - startcol) Mod 2 + 1
lisaysr = (endrow - startrow) Mod 2 + 1
For i = startcol To (endcol - lisaysc - 1)
   i_step = i_step + 1
   j_step = -1
  For j = (startrow + lisaysr + 1) To endrow
  j_step = j_step + 1
    For k = 0 To CInt(pan_x) - 1
    For l = 0 To CInt(pan_y) - 1
     bkuva(i_step * CInt(pan_x) + k, j_step * CInt(pan_y) + l).r = b_array(i, j, Index)
     bkuva(i_step * CInt(pan_x) + k, j_step * CInt(pan_y) + l).g = b_array(i, j, Index)
     bkuva(i_step * CInt(pan_x) + k, j_step * CInt(pan_y) + l).B = b_array(i, j, Index)
     bkuva(i_step * CInt(pan_x) + k, j_step * CInt(pan_y) + l).N = 0
    Next l
    Next k
  Next j
 Next i
   
GoTo loppu

pan_x_1:
k = -1
For i = startcol To endcol
 l = -1
 k = k + 1
  For j = startrow To endrow
   l = l + 1
      bkuva(k, l).r = b_array(i, j, Index)
      bkuva(k, l).g = b_array(i, j, Index)
      bkuva(k, l).B = b_array(i, j, Index)
      bkuva(k, l).N = 0
   
Next j
Next i
GoTo loppu


loppu:

SetBitmapBits Form1.Picture1(Index).Image, bm.bmWidthBytes * CLng(win_h), bkuva(0, 0)
Form1.Picture1(Index).Refresh

Exit Sub

error_in_creating_bitmap:
MsgBox ("Error in creating the bitmap for display!")
Exit Sub


End Sub
Public Sub Clear_all_images()
Dim i As Integer
For i = 0 To NumOfImages - 1

 Form1.Check1(i).Value = 0
 Form1.Text1(i * 2).Text = ""
 Form1.Text12(i).Text = ""
 Form1.Text1((i * 2) + 1).Text = ""
 Form1.Picture1(i).Cls
 LastImageClicked = -1
Next i

End Sub


Public Sub Calc_dX_dY_dZ_dx_dy(ByVal cam As Integer, X_ini, Y_ini, Z_ini, dx_dU_K, dy_dU_K, dx_dV_K, dy_dV_K, dx_dW_K, dy_dW_K, dx, dy, X, Y)

' This subroutine calculates the differential quotiens needed in the calculation of space
' intersection and (approximation - obs) residuals
' (elements of normal vector)
' dx_dUK, dy_dU_K,....,dy_dW_K, dx, dy
' cam stores index to image info
' _ini, Y_ini, Z_ini (initial / iteratively improved) XYZ-coordinates

Dim omega, phi, kappa As Double
Dim U, V, w, U0, V0, W0 As Double
Dim z, Z0, x0, y0, x_k, y_k As Double
'Dim x, y as bouble
Dim Zx, Zy, N, C As Double
Dim i, c_index As Integer

Rem Initial approximations
U = X_ini:                  V = Y_ini:              w = Z_ini


omega = image_info(cam).omega: phi = image_info(cam).phi: kappa = image_info(cam).kappa
X = CDbl(Form1.Text1(cam * 2).Text): Y = CDbl(Form1.Text1(cam * 2 + 1).Text)

'Exit Sub

' MsgBox ("x y : " & CStr(x) & " " & CStr(y))
z = 0#
x0 = 0#: y0 = 0#: Z0 = image_info(cam).C
U0 = image_info(cam).Xo: V0 = image_info(cam).Yo: W0 = image_info(cam).Zo

Zx = A(1, 1, cam) * (U - U0) + A(2, 1, cam) * (V - V0) + A(3, 1, cam) * (w - W0)
Zy = A(1, 2, cam) * (U - U0) + A(2, 2, cam) * (V - V0) + A(3, 2, cam) * (w - W0)
N = A(1, 3, cam) * (U - U0) + A(2, 3, cam) * (V - V0) + A(3, 3, cam) * (w - W0)
C = (z - Z0)
 
 
dx_dU_K = (-C / (N ^ 2)) * (N * A(1, 1, cam) - Zx * A(1, 3, cam))
dy_dU_K = (-C / (N ^ 2)) * (N * A(1, 2, cam) - Zy * A(1, 3, cam))
dx_dV_K = (-C / (N ^ 2)) * (N * A(2, 1, cam) - Zx * A(2, 3, cam))
dy_dV_K = (-C / (N ^ 2)) * (N * A(2, 2, cam) - Zy * A(2, 3, cam))
dx_dW_K = (-C / (N ^ 2)) * (N * A(3, 1, cam) - Zx * A(3, 3, cam))
dy_dW_K = (-C / (N ^ 2)) * (N * A(3, 2, cam) - Zy * A(3, 3, cam))


' Calculate the image coordinates from the ground point
x_k = x0 + (z - Z0) * (A(1, 1, cam) * (U - U0) + A(2, 1, cam) * (V - V0) + A(3, 1, cam) * (w - W0)) / (A(1, 3, cam) * (U - U0) + A(2, 3, cam) * (V - V0) + A(3, 3, cam) * (w - W0))
y_k = y0 + (z - Z0) * (A(1, 2, cam) * (U - U0) + A(2, 2, cam) * (V - V0) + A(3, 2, cam) * (w - W0)) / (A(1, 3, cam) * (U - U0) + A(2, 3, cam) * (V - V0) + A(3, 3, cam) * (w - W0))
' calculate error (residual)
dx = x_k - X
dy = y_k - Y
' If cam = 5 Then
' MsgBox (CStr(dx) & "," & CStr(dy))
' End If
End Sub
    Function det_cal2(E() As Double) As Double
    ' This function calculates the determinant of a 2x2  matrix (E)
    Dim term1, term2 As Double
    term1 = E(1, 1) * E(2, 2)
    term2 = E(2, 1) * E(1, 2)
    det_cal2 = term1 - term2
    End Function


Public Sub MATINV3(CC() As Double, DD() As Double)

Rem CC 3x3 input, inverted version DD 3x3

ReDim EE(1 To 2, 1 To 2) As Double
ReDim el(1 To 4) As Double
Rem here starts Cramers rule (with 2x2 sub determinants)
For ix = 1 To 3:  For jx = 1 To 3
  Num = 0
   For kX = 1 To 3:   For lx = 1 To 3
     mx = kX:    nx = lx

Rem Choose elements of subdeterminants
    If lx = jx Then GoTo bypasscolumn
    If kX = ix Then GoTo bypassrow
    Num = Num + 1
    Rem by columns
    el(Num) = CC(kX, lx)
bypasscolumn:
    Next lx
bypassrow:
    Next kX
    EE(1, 1) = el(1):     EE(1, 2) = el(2):   EE(2, 1) = el(3):   EE(2, 2) = el(4)
    DD(jx, ix) = (1 / det_D) * det_cal2(EE())
    If (jx + ix) Mod 2 <> 0 Then
      DD(jx, ix) = -DD(jx, ix)
    End If
  Next jx
 Next ix


End Sub
Function det_cal3(E() As Double) As Double
' This function calculates the determinant of a 3x3
' matrix (E)
Dim term1, term2, term3, term4, term5, term6 As Double
term1 = E(1, 1) * E(2, 2) * E(3, 3)
term2 = E(2, 1) * E(3, 2) * E(1, 3)
term3 = E(3, 1) * E(1, 2) * E(2, 3)
term4 = E(1, 1) * E(3, 2) * E(2, 3)
term5 = E(2, 1) * E(1, 2) * E(3, 3)
term6 = E(3, 1) * E(2, 2) * E(1, 3)
det_cal3 = term1 + term2 + term3 - term4 - term5 - term6
End Function
Public Sub solve_3D_X_Y_from_Z_x_y(ByVal q As Integer, ByVal Z_gnd As Double, ByVal X As Double, ByVal Y As Double, ByRef X_Gnd As Double, ByRef Y_Gnd As Double)

Dim C As Double ' camera constant
Dim k As Double ' denominator in inverse collinear eq.

C = -image_info(q).C

k = A(3, 1, q) * X + A(3, 2, q) * Y + A(3, 3, q) * C

If k <> 0# Then
 X_Gnd = image_info(q).Xo + (Z_gnd - image_info(q).Zo) * ((A(1, 1, q) * (X) + A(1, 2, q) * (Y) + A(1, 3, q) * (C)) / k)
 Y_Gnd = image_info(q).Yo + (Z_gnd - image_info(q).Zo) * ((A(2, 1, q) * (X) + A(2, 2, q) * (Y) + A(2, 3, q) * (C)) / k)
End If

End Sub

Public Sub r_transform_3D(ByVal q As Integer, ByVal X As Double, ByVal Y As Double, ByVal z As Double, p_x, p_y)

Rem Calculates 3D --> 2D transformation, for image q

Dim k As Double    ' inverse of this number is the scale factor

Rem MsgBox ("X,Y,Z in transform_3D: " & Format$(X, "0.00 ") & Format$(Y, "0.00 ") & Format$(Z, "0.00 "))

If (image_info(q).C <> 0#) Then

    k = (A(1, 3, q) * (X - image_info(q).Xo) + A(2, 3, q) * (Y - image_info(q).Yo) + A(3, 3, q) * (z - image_info(q).Zo))
    Rem  MsgBox ("Value for K: " & CStr(K))
   Else
  k = 0#
End If

   If (k <> 0) Then
      p_x = -image_info(q).C * (A(1, 1, q) * (X - image_info(q).Xo) + A(2, 1, q) * (Y - image_info(q).Yo) + A(3, 1, q) * (z - image_info(q).Zo)) / k
      p_y = -image_info(q).C * (A(1, 2, q) * (X - image_info(q).Xo) + A(2, 2, q) * (Y - image_info(q).Yo) + A(3, 2, q) * (z - image_info(q).Zo)) / k
   Else
      p_x = p_y = 0#
   End If
 Rem   MsgBox ("camera coordinates in tf_3D: " & CStr(p_x) & " " & CStr(p_y))
End Sub

Public Sub r_transform_ground_to_pixel(ByVal q As Integer, ByVal X As Double, ByVal Y As Double, ByVal z As Double, p_x As Double, p_y As Double)
Rem from ground coordinates to Pixel values (imagecoordinates)
Dim camera_x As Double, camera_y As Double
Rem MsgBox ("X,Y,Z in transform ground to pixel: " & Format$(X, "0.00 ") & Format$(Y, "0.00 ") & Format$(Z, "0.00 "))
Call r_transform_3D(q, X, Y, z, camera_x, camera_y)
' MsgBox ("camera x & y returned from 3D->2D :" & Format$(camera_x * 100#, "0.000000000 cm") & Format$(camera_y * 100#, "0.00000000 cm"))
Call a_transform_affine(q, 0, camera_x, camera_y, p_x, p_y)
' MsgBox ("image x & y returned from 3D->2D :" & Format$(p_x, "0.00 ") & Format$(p_y, "0.00 "))
End Sub

Public Sub r_transform_matrix(ByVal q As Integer)
    
Rem Rotation matrix for image q

A(1, 1, q) = Cos(image_info(q).phi) * Cos(image_info(q).kappa): A(1, 2, q) = -Cos(image_info(q).phi) * Sin(image_info(q).kappa): A(1, 3, q) = Sin(image_info(q).phi)
A(2, 1, q) = Cos(image_info(q).omega) * Sin(image_info(q).kappa) + Sin(image_info(q).omega) * Sin(image_info(q).phi) * Cos(image_info(q).kappa): A(2, 2, q) = Cos(image_info(q).omega) * Cos(image_info(q).kappa) - Sin(image_info(q).omega) * Sin(image_info(q).phi) * Sin(image_info(q).kappa): A(2, 3, q) = -Sin(image_info(q).omega) * Cos(image_info(q).phi)
A(3, 1, q) = Sin(image_info(q).omega) * Sin(image_info(q).kappa) - Cos(image_info(q).omega) * Sin(image_info(q).phi) * Cos(image_info(q).kappa): A(3, 2, q) = Sin(image_info(q).omega) * Cos(image_info(q).kappa) + Cos(image_info(q).omega) * Sin(image_info(q).phi) * Sin(image_info(q).kappa): A(3, 3, q) = Cos(image_info(q).omega) * Cos(image_info(q).phi)
   
End Sub
Public Sub a_transform_collinear(ByVal x_in As Double, ByVal y_in As Double, ByVal z_in As Double, ByVal C As Double, ByRef Rij() As Double, ByRef x_out As Double, ByRef y_out As Double)

x_out = -C * (Rij(1, 1) * x_in + Rij(1, 2) * y_in - Rij(1, 3) * z_in) / (Rij(3, 1) * x_in + Rij(3, 2) * y_in - Rij(3, 3) * z_in)
y_out = -C * (Rij(2, 1) * x_in + Rij(2, 2) * y_in - Rij(2, 3) * z_in) / (Rij(3, 1) * x_in + Rij(3, 2) * y_in - Rij(3, 3) * z_in)

End Sub
Public Sub a_transform_helmert(ByVal q As Integer, ByVal direction As Integer, X, Y, p_x, p_y)

Rem *****************************************************************************
Rem a_transform_helmert
Rem Transforms according to parameters given
Rem q = image number
Rem direction 0 (x,y) -> (col,row), 1 (row,col)-> (x,y)
Rem *****************************************************************************
Dim x_ori, y_ori As Double

x_ori = X    ' arguments x,y are not allowed to change their value inside this Sub
y_ori = Y
Select Case direction
 Case 0   ' From camera coordinates to image cordinates (pixels)
   X = X - image_info(q).mean_x
   Y = Y - image_info(q).mean_y
   p_x = image_info(q).X_mean + image_info(q).lambda * Cos(image_info(q).alpha) * X - image_info(q).lambda * Sin(image_info(q).alpha) * Y
   p_y = image_info(q).Y_mean + image_info(q).lambda * Sin(image_info(q).alpha) * X + image_info(q).lambda * Cos(image_info(q).alpha) * Y
   X = x_ori
   Y = y_ori
   
 Case 1  ' From image coordinates (pixels) to camera cordinates (mm)
   X = X - image_info(q).X_mean
   Y = Y - image_info(q).Y_mean
   p_x = image_info(q).mean_x + (1# / image_info(q).lambda) * Cos(-image_info(q).alpha) * X - (1# / image_info(q).lambda) * Sin(-image_info(q).alpha) * Y
   p_y = image_info(q).mean_y + (1# / image_info(q).lambda) * Sin(-image_info(q).alpha) * X + (1# / image_info(q).lambda) * Cos(-image_info(q).alpha) * Y
   X = x_ori
   Y = y_ori
End Select
End Sub