Public Sub Match_Local_Crown_A(sp As Byte) If sp = 4 Then sp = 3 If sp = 5 Then sp = 3 If sp = 6 Then sp = 3 If sp = 7 Then sp = 3 If sp = 13 Then sp = 3 If sp = 16 Then sp = 3 If sp = 20 Then sp = 3 If sp = 21 Then sp = 2 If sp < 1 Or sp > 3 Then MsgBox ("Species unresolved") Exit Sub End If Rem November 2006 - April 2014. Algorithm for fitting a crown model (Added from Form1 to BAS) Dim Apu_z As Long, Colorpix As Long, Zmaa As Double, apu As Single, j As Long Dim MaxDTM As Double, MinDTM As Double, i As Long, p_x As Double, p_y As Double Dim FN1 As String * 3, FN2 As String * 3 Dim Xc As Double, Yc As Double Dim LiDARBinPath As String, LiDARBinPath2 As String Dim Npulses As Long, Nha As Long, Nsum As Long, Hlimit As Double, Hscale As Double, DensityFactor As Double Dim Nha2 As Long, P As Double, MPlot As Long Dim Wa As Double, Wc As Double Rem *********** Parameters ******** DensityFactor = 1 On Error Resume Next Hscale = 0.6 ' Down to what relative height do we model the crown DensityFactor = 1 + CDbl(Form1.Text11(0).Text) ' Scale the average height : dcrm -ratio P = 1 ' The weight by which negative residuals (inside crown) are weighted MPlot = 1 Wa = 10 Wc = 30 Close (1) Xs = X_sol: Ys = Y_sol Rem *************** WE NEED TO SAMPLE SEVERAL HECTARES *************** Rem Rem The stored hectares are maintained as long as there is a need to change - re-read Rem There can be 1, 2 or 4 hectares in memory. If the point (x,y) has an x-coordinate Rem close to 100; read both sides; similarly for the y-coordinates FN1 = Format$(Int((Xs - 2510000) / 100), "000") FN2 = Format$(Int((Ys - 6850000) / 100), "000") Rem LiDR -data structure holds data, use it if we are in the same 1-ha cell Dim k As Long Open "c:\data\als2013b_path.hdr" For Input As 1 Input #1, LiDARBinPath Close (1) Close (100) Open LiDARBinPath & FN1 & "_" & FN2 & ".bin" For Binary As 100 Get #100, , Nha Close (100) Npulses = Nha Form1.Caption = "Pulses per m2: " & Format$(Nha / 10000#, "0.0") On Error Resume Next Rem If we are in the same hectare, use current LiDR() - array Ax = -1 Ax = UBound(Lidr) - Nha If Ax = 0 Then GoTo SkipReadingLiDAR ReDim Lidr(1 To Nha) As LidarRecord Open LiDARBinPath & FN1 & "_" & FN2 & ".bin" For Binary As 100 For k = 1 To Nha Get #100, 5 + k * 207, Lidr(k) ' Here we read 207-byte records record by record Next k Close (100) SkipReadingLiDAR: Hlimit = (Z_sol - getheight(X_sol, Y_sol)) * Hscale Htree = (Z_sol - getheight(X_sol, Y_sol)) * 10 ' dm Tree.H = Htree Rem Calculate average, NFI-based DCRM Select Case sp Case 1 ' Pine dcrm = Htree * 0.1636 + 13 ' Case 2, 4 ' Spruce, dead tree dcrm = Htree * 0.15 + 14 Case 3 ' Birch dcrm = Htree * 0.214 + 8 End Select Rem Rem Scale it to crown radiu, consider stand density dcrm = (((dcrm / 10)) / 2) * DensityFactor DistA = dcrm MaxDTM = 25 ' default MinDTM = 0 ' default Rem *********************************************************** Rem Lidar DATA has been read. Here starts the processing of it. Rem *********************************************************** Dim Np As Long ReDim cf(0 To 4) As Double, pw(0 To 4) As Double, cont2(0 To 4) As Double Dim cfa As Double, pwa As Double, cont2a As Double, H As Double Dim red As Double, blue As Double, green As Double, mx As Long, Ntrees As Long, maxA As Long Dim Hlowest As Double maxA = 2000 ReDim Inpoints(1 To maxA) As Point3d ' The points that make the tree ReDim intensities(1 To maxA) As Double ReDim LA(1 To maxA) As Double ' Residuals in crown radius (error vector v, vTv is minimized) ReDim A_MAT(1 To maxA, 1 To 3) As Double ' Partial derivates of the unknwons crown shape parameters, design matrix ReDim XYdist(1 To maxA) As Single ReDim Zdist(1 To maxA) As Single ReDim sade(1 To maxA) As Single Dim testpoint As Point3d Dim rmse As Double Rem Crown model, crown radius as function of tree height (h) and the relative distance down the crown (zdist/h) Rem cf * h * (zdist / h) ^ pw + cont2 Rem Assign basic values per species: pine(1), spruce(2), birch(3), aspen (4) For i = 0 To NumOfImages - 1 Form1.Picture1(i).Cls Next i He_ = Htree / 10# cf(0) = dcrm / He_: cf(1) = dcrm / He_: cf(2) = dcrm / He_: cf(3) = (dcrm / He_): cf(4) = (dcrm / He_) pw(0) = 0.75: pw(1) = 0.5: pw(2) = 0.6: pw(3) = 0.45: pw(4) = 0.7 cont2(0) = 0.4: cont2(1) = 0.6: cont2(2) = 0.6: cont2(3) = 0.9: cont2(4) = 0.6 Tree.cf = cf(sp) Tree.pw = pw(sp) + CDbl(Form1.Text11(1).Text) Tree.cont2 = cont2(sp) + CDbl(Form1.Text11(2).Text) Np = UBound(Lidr()) ReDim LPE_AR(1 To 1) As LP_Espoo ' This array holds the lidar points; no intensity ReDim LPE_AR(1 To maxA) As LP_Espoo ' This array holds the lidar points; no intensity ReDim Z_gnd(1 To maxA) As Single Zmax = 0 l = 0 Hup = 0.3 ' underestiomation of H For k = 1 To Np For j = 4 To 4 Step -1 ' take only first retrun data? If Lidr(k).Returns(j).X > 1 And Lidr(k).Returns(j).z < (Z_sol + Hup) And Lidr(k).PosLiDAR.z < 2500 Then 'And Abs(LidR(k).Returns(1).z - LidR(k).Returns(4).z) < 1 Then Rem Here we collect the observations, we allow .25 m higher testisade = Tree.cf * (Tree.H * 0.1) * (2.5 * ((Z_sol + Hup) - Lidr(k).Returns(j).z) / (Tree.H * 0.1)) ^ Tree.pw + Tree.cont2 ^ 2 Form1.Picture1(0).DrawWidth = 2 ZDIFFTEST = Z_sol - Lidr(k).Returns(j).z disttostem = Sqr((Lidr(k).Returns(j).X - X_sol) ^ 2 + (Lidr(k).Returns(j).Y - Y_sol) ^ 2) If disttostem < 4 Then apu = (Tree.H / 10 - ZDIFFTEST) MinDTM = 0.6 * Tree.H / 10 MaxDTM = Tree.H / 10 If apu < MinDTM Then apu = MinDTM If apu > MaxDTM Then apu = MaxDTM - 0.02 Apu_z = 255 - (255 - (apu - MinDTM) * 255 / (MaxDTM - MinDTM)) ' DTM 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)) mp = Form1.Picture1(0).Width / (Tree.H / 10) - 3 ybot = Form1.Picture1(0).Height - 20 Ytop = 20 Form1.Picture1(0).DrawWidth = 3 If apu / (Tree.H / 10) > 0# Then Form1.Picture1(0).PSet (10 + disttostem * mp, Ytop + ZDIFFTEST * mp), Colorpix ', 'RGB(255, 0, 0) 'Form1.Picture1(0).PSet (10 + ZDIFFTEST * mp, ybot - disttostem * mp), Colorpix ', 'RGB(255, 0, 0) End If End If If Sqr((Lidr(k).Returns(j).X - X_sol) ^ 2 + (Lidr(k).Returns(j).Y - Y_sol) ^ 2) < testisade Then jx = 1 Form1.Picture1(jx).DrawWidth = 3 Call r_transform_ground_to_pixel(CLng(jx), Lidr(k).Returns(j).X, Lidr(k).Returns(j).Y, Lidr(k).Returns(j).z, p_x, p_y) Form1.Picture1(jx).PSet ((p_x - (image_info(jx).o_col + win_info(jx).win_o_col)) * win_info(jx).pan_x - 1, -1 + ((image_info(jx).Height - 1) - p_y - (image_info(jx).o_row + win_info(jx).win_o_row)) * win_info(jx).pan_y), RGB(255, 255, 255) l = l + 1 Z_gnd(l) = getheight(Lidr(k).Returns(j).X, Lidr(k).Returns(j).Y) LPE_AR(l).Echo = j LPE_AR(l).X = Lidr(k).Returns(j).X LPE_AR(l).Y = Lidr(k).Returns(j).Y LPE_AR(l).z = Lidr(k).Returns(j).z If Lidr(k).PosLiDAR.z < 1500 Then ' 1 km point LPE_AR(l).intensity = ((1 - (1201 - Lidr(k).range(j)) / 1201) ^ 2.3) * Lidr(k).intensity(j) + (0.046 * Lidr(k).intensity(j) * (133 - Lidr(k).Res1)) End If If Lidr(k).PosLiDAR.z > 1500 Then ' 2km point LPE_AR(l).intensity = ((1 - (1942 - Lidr(k).range(j)) / 1942) ^ 2.3) * Lidr(k).intensity(j) + (0.0517 * Lidr(k).intensity(j) * (137 - Lidr(k).Res1)) End If End If End If Next j Next k Close (1) NCounter = 0 If l < 2 Then MsgBox ("Less than 2 Points!"): Exit Sub End If ReDim Preserve LPE_AR(1 To l) As LP_Espoo ReDim Preserve Z_gnd(1 To l) As Single ReDim LPTreeHit(1 To l) As Boolean ReDim LPE_SP(1 To l) As Byte Rem Declare storage for ground normalized heights of first pulses (for sorting) i = UBound(LPE_AR) ReDim Hftau(1 To i) As Single Rem Array that hold the indeces for sorted heights ReDim Hfind(1 To i) As Long Rem Copy Hmax = 0 For mx = 1 To i Hftau(mx) = LPE_AR(mx).z - Z_gnd(mx) Hmax = Max(Hmax, Hftau(mx)) Next mx Rem Sort data in order 1st = lowest, last = highest Label10.Caption = "Sorting normalized heights..." DoEvents Call indexx(CLng(UBound(Hftau)), Hftau, Hfind) Rem Declare storage for treetops to be found Label10.Caption = "Fitting crowns..." DoEvents Rem ********** Rem FIRST TREE Rem Make foto-solution as seed Ntrees = 1 Tree.X = X_sol Tree.Y = Y_sol Tree.H = He_ Tree.z = Z_sol Zgnd = getheight(X_sol, Y_sol) Tree.NLP = 0 LPTreeHit(Hfind(i)) = True Rem Determine RGB-vector from image j for point XYZ j = 1 testpoint.X = Tree.X testpoint.Y = Tree.Y testpoint.z = Tree.z Rem OMITTED SAMPLING! Rem Call RGB_Vector_For_Point(red, green, blue, testpoint, CLng(j)) Rem Include all points indide the crown and fit the model Rem Down to rhe height of Hlimit Dim Xsum As Double, Ysum As Double, Hlow As Double Xsum = 0: Ysum = 0 Rem the height of the tree, mark points down to 50% height H = Tree.H Hlow = Hlimit Dim HlidMax As Double HlidMax = 0 MinDTM = 0 MaxDTM = 10 Form1.Picture1(0).DrawWidth = 2 Form1.Picture1(1).DrawWidth = 2 Dim lp As ads40_image_point_struct For i = UBound(Hftau) To 1 Step -1 If Hftau(Hfind(i)) < (H + Hup) And Hftau(Hfind(i)) > Hlow Then XYdist_maxA = Sqr((LPE_AR(Hfind(i)).X - Tree.X) ^ 2 + (LPE_AR(Hfind(i)).Y - Tree.Y) ^ 2) Rem Distance down from top, m Zdist(maxA) = (Tree.H - Hftau(Hfind(i))) If Zdist(maxA) < 0 Then Zdist(maxA) = 0 Rem Crown radius, given height and species and current three parameters cf, pw, cont2 'sade(maxA) = Tree.cf * Tree.H * Sin(4 * (Zdist(maxA) / H)) ^ Tree.pw + Tree.cont2 sade(maxA) = Tree.cf * Tree.H * (2.5 * (Zdist(maxA) / H)) ^ Tree.pw + Tree.cont2 ^ 2 If sade(maxA) > XYdist(maxA) Then ' And zdist(maxA) < Hlow Then Rem Distance to top (downwards) If Hftau(Hfind(i)) > HlidMax Then HlidMax = Hftau(Hfind(i)) Tree.NLP = Tree.NLP + 1 XYdist(Tree.NLP) = XYdist_maxA Zdist(Tree.NLP) = Zdist(maxA) LPTreeHit(Hfind(i)) = True Rem Collect observations: Points XYZ Inpoints(Tree.NLP).X = LPE_AR(Hfind(i)).X Xsum = Xsum + LPE_AR(Hfind(i)).X Inpoints(Tree.NLP).Y = LPE_AR(Hfind(i)).Y Ysum = Ysum + LPE_AR(Hfind(i)).Y Inpoints(Tree.NLP).z = LPE_AR(Hfind(i)).z intensities(Tree.NLP) = LPE_AR(Hfind(i)).intensity MinDTM = 0.6 * Tree.H MaxDTM = Tree.H apu = (MaxDTM - Zdist(maxA)) 'apu = Hftau(Hfind(i)) If apu < MinDTM Then apu = MinDTM If apu > MaxDTM Then apu = MaxDTM - 0.02 Apu_z = 255 - (255 - (apu - MinDTM) * 255 / (MaxDTM - MinDTM)) ' DTM 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)) ' For jx = 0 To 1 ' NumOfImages - 1 ' do we need all images? Form1.Picture1(jx).DrawWidth = 3 ' If image_info(jx).Imagetype = "ADS L0" Then ' Call R_transform_ground_to_ADS40(CLng(jx), LPE_AR(Hfind(i)).X, LPE_AR(Hfind(i)).Y, CDbl(LPE_AR(Hfind(i)).z), lp) ' Form1.Picture1(jx).PSet ((lp.Sample - (image_info(jx).o_col + win_info(jx).win_o_col)) * win_info(jx).pan_x, (lp.Line - (image_info(jx).o_row + win_info(jx).win_o_row)) * win_info(jx).pan_y), Colorpix ' End If ' If image_info(jx).Imagetype = "FRAME" Then Call r_transform_ground_to_pixel(CLng(jx), LPE_AR(Hfind(i)).X, LPE_AR(Hfind(i)).Y, LPE_AR(Hfind(i)).z, p_x, p_y) Form1.Picture1(jx).PSet ((p_x - (image_info(jx).o_col + win_info(jx).win_o_col)) * win_info(jx).pan_x - 1, -1 + ((image_info(jx).Height - 1) - p_y - (image_info(jx).o_row + win_info(jx).win_o_row)) * win_info(jx).pan_y), Colorpix ' End If ' Next jx End If End If Next i Close (6) Close (15) ' Form1.Picture1(1).DrawWidth = 1 Rem Here we shift to mass center LPTreeHit(Hfind(UBound(Hftau))) = True If Tree.NLP < 3 Then MsgBox ("Too few LiDAR points!") GoTo DoneFitting Exit Sub End If Dim A0 As Double, c0 As Double A0 = Tree.cont2 c0 = Tree.pw Rem ******* LOOP ********** NITE = 0 StartOfFitModel: Call FitModelB(rmse, Hlowest, Ntrees, XYdist(), sade(), Zdist(), Trees(), LA(), A_MAT(), P, A0, c0, Wa, Wc) ' Call FitModel(rmse, Hlowest, Ntrees, XYdist(), sade(), Zdist(), Trees(), LA(), A_MAT(), P) apu = MYFUNC_MATLABCALLNOMSGS(CDbl(1)) NITE = NITE + 1 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 > 15 Then GoTo DoneFitting 'Open "c:\data\parameters.txt" For Append As 11 'Print #11, NITE, Tree.NLP, Tree.H, Tree.cf, Tree.pw, Tree.cont2, rmse 'Close (11) Tree.cf = Tree.cf + Dcf Tree.pw = Tree.pw + Dpw Tree.cont2 = Tree.cont2 + Dcont2 If Tree.pw < 0 Or rmse > 3 Then Tree.pw = 0.5 Tree.cont2 = 0.1 Tree.cf = 0.15 End If GoTo StartOfFitModel Rem ******* END GOTO-LOOP ********** Rem Let's draw the crown in the images DoneFitting: Form1.Label10.FontSize = 9 sade(maxA) = 2 * (Tree.cf * Tree.H * (2.5 * (Zdist(maxA) / H)) ^ Tree.pw + Tree.cont2 ^ 2) Form1.Label10.Caption = "Done, iters: " & NITE & " RMS " & Format$(rmse, "0.00 m") & " Dcrm: " & Format$(sade(maxA), "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 Dim LP_BEG As ads40_image_point_struct, lp_end As ads40_image_point_struct For Zdiff = 0.01 To (Tree.H * (1 - Hscale)) Step Tree.H * 0.05 Zt = Tree.z - Zdiff 'Radius = (Tree.cf * Tree.H * Sin(4 * ((Zdiff) / H)) ^ Tree.pw + Tree.cont2) '* 1.2 Radius = Tree.cf * Tree.H * (2.5 * ((Zdiff) / H)) ^ Tree.pw + Tree.cont2 ^ 2 mx = 0 Shift = 3.14 / 10 MaxDTM = (Z_sol - getheight(X_sol, Y_sol)) * 1.05 MinDTM = MaxDTM * 0.55 For phi = -3.14 To 3.15 Step 3.14 / 8 mx = mx + 1 Colorpix = RGB(255 - Zdiff * 10, 255 - Zdiff * 10, 255) ' 255 - Zdiff * 10) Xt = Radius * Cos(phi) + Tree.X Yt = Radius * Sin(phi) + Tree.Y apu = Zt - getheight(Xt, Yt) If apu < MinDTM Then apu = MinDTM If apu > MaxDTM Then apu = MaxDTM - 0.02 Apu_z = 255 - (255 - (apu - MinDTM) * 255 / (MaxDTM - MinDTM)) ' DTM 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)) Xf = Radius * Cos(phi + Shift) + Tree.X Yf = Radius * Sin(phi + Shift) + Tree.Y Rem Using Line-method draw in all images on screen For i = 0 To NumOfImages - 1 Xvec = image_info(i).Xo - Xt: Yvec = image_info(i).Yo - Yt: Zvec = image_info(i).Zo - Zt Dvec = Sqr(Xvec ^ 2 + Yvec ^ 2 + Zvec ^ 2) Xvec = 0.1 * (Xvec / Dvec): Yvec = 0.1 * (Yvec / Dvec): Zvec = 0.1 * (Zvec / Dvec): ZDIFFTEST = Tree.z - (Zt + Zvec) disttostem = Sqr((Tree.X - (Xt + Xvec)) ^ 2 + (Tree.Y - (Yt + Yvec)) ^ 2) radiustest = Tree.cf * Tree.H * (2.5 * ((ZDIFFTEST) / H)) ^ Tree.pw + Tree.cont2 ^ 2 ' Xt, Yt, Zt represent the point, now draw a 10 cm vector towards camera - if inside - not visible Select Case image_info(i).Imagetype Case "FRAME" ' If distB > DistA 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 = 2 Form1.Picture1(i).ForeColor = RGB(255, 255, 255) If disttostem > radiustest Then 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 'RGB(255, 120, 255) End If Case "ADS L0" ' Form1.Picture1(i).DrawWidth = 1 ' Call R_transform_ground_to_ADS40(CLng(i), Xt, Yt, Zt, LP_BEG) ' Call R_transform_ground_to_ADS40(CLng(i), Xf, Yf, CDbl(Zt), lp_end) ' Form1.Picture1(i).Line ((LP_BEG.Sample - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, (LP_BEG.Line - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y)-((lp_end.Sample - (image_info(i).o_col + win_info(i).win_o_col)) * win_info(i).pan_x - 1, (lp_end.Line - (image_info(i).o_row + win_info(i).win_o_row)) * win_info(i).pan_y), RGB(255, 120, 255) End Select Next i DoEvents Next phi Next Zdiff 'Form1.Picture1(1).FontSize = 12 'Form1.Picture1(1).CurrentX = Form1.Picture1(1).CurrentX - 50 'Form1.Picture1(1).CurrentY = Form1.Picture1(1).CurrentY + 20 'Form1.Picture1(1).Print "CW " & Format$((Tree.cf * Tree.H * Sin(4 * ((Zdiff) / H)) ^ Tree.pw + Tree.cont2) * 2, "0.0 m") & " H " & Format$(Tree.H, "0.0 m") Close (6) mp = Form1.Picture1(0).Width / Tree.H - 3 Form1.Picture1(0).DrawWidth = 1 ybot = Form1.Picture1(0).Height - 20 Ytop = 20 'Form1.Picture1(0).Line (10, ybot)-(Tree.H * mp, ybot), RGB(255, 255, 0) Form1.Picture1(0).Line (10, Ytop)-(10, Ytop + Tree.H * mp), RGB(255, 255, 0) Form1.Picture1(0).DrawWidth = 3 For l = 0 To 0.4 Step 0.01 ' Form1.Picture1(0).PSet (10 + l * Tree.H * mp, ybot - (Tree.cf * Tree.H * (2.5 * ((l)) ^ Tree.pw + Tree.cont2 ^ 2) * mp)), RGB(255, 255, 255) Form1.Picture1(0).PSet (10 + (Tree.cf * Tree.H * (2.5 * ((l)) ^ Tree.pw + Tree.cont2 ^ 2) * mp), Ytop + l * Tree.H * mp), RGB(255, 255, 255) Next l Form1.Picture1(0).DrawWidth = 2 Open "c:\temp\" & sp & ".txt" For Append As 8 Noutput = 0 MinDTM = Tree.H * 0.6 MaxDTM = Tree.H For l = 1 To UBound(LPE_AR) Form1.Picture1(0).DrawWidth = 2 ZDIFFTEST = Tree.z - (LPE_AR(l).z) disttostem = Sqr((Tree.X - (LPE_AR(l).X)) ^ 2 + (Tree.Y - LPE_AR(l).Y) ^ 2) radiustest = Tree.cf * Tree.H * (2.5 * ((ZDIFFTEST) / H)) ^ Tree.pw + Tree.cont2 ^ 2 apu = Tree.H - ZDIFFTEST If apu < MinDTM Then apu = MinDTM If apu > MaxDTM Then apu = MaxDTM - 0.02 Apu_z = 255 - (255 - (apu - MinDTM) * 255 / (MaxDTM - MinDTM)) ' DTM 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)) If disttostem < 0.3 And ZDIFFTEST > 0.1 * Tree.H Then Colorpix = RGB(255, 255, 255) Form1.Picture1(0).DrawWidth = 3 End If ' If disttostem < 1.1 * radiustest Then ' Print #8, zdifftest, zdifftest / Tree.H, disttostem, disttostem / radiustest ' Noutput = Noutput + 1 ' End If Form1.Picture1(0).PSet (10 + disttostem * mp, Ytop + ZDIFFTEST * mp), Colorpix ', 'RGB(255, 0, 0) Next l 'Form1.Picture1(5).Print Noutput Close (8) Label10.FontSize = 8 Rem The variables of measurement 'Measurement.plot = Plot_Info.Number Rem Start Collecting Image and LiDAR features Rem Of the images that have 16-bit R,G,B,NIR,PAN values, select one in direct light and Rem One in back-light. Move the point down 10% of tree height, towards sun and opposite dir. Rem This way we have XYZ-points Point.InSun and Point.InShade. Compute the camera-point Rem distances to these, select the closest / furthest Dim imdata As RGBNIR, Fii As Double, Theta As Double Dim Pan_mean As Double, Pan_Max As Double, Pan_Min As Double, Pan_SD As Double Dim objtocam As Vector3D, objtosun As Vector3D Dim PointInLight As Vector3D, PointInShade As Vector3D Dim PlumbLine As Vector3D PlumbLine.X = 0 PlumbLine.Y = 0 PlumbLine.z = 1 Rem Point on surface, in light and shadow 'm Cpu String ' holds output of image features Select Case sp 'LPE_SP(Ntrees) Case 1 Measurement.dbh_lidar_Dcrm = ((-3.155 + 1.323 * sqrt(sade(maxA) * 10) + 0.73 * sqrt(Tree.H * 10)) ^ 2 + 0.811) / 10 d13Foto = ((-3.155 + 1.323 * sqrt(Measurement.Foto_Dcrm * 10) + 0.73 * sqrt(Tree.H * 10)) ^ 2 + 0.811) / 10 Case 2 Measurement.dbh_lidar_Dcrm = ((-3.214 + 1.016 * sqrt(sade(maxA) * 10) + 0.861 * sqrt(Tree.H * 10)) ^ 2 + 0.742) / 10 d13Foto = ((-3.214 + 1.016 * sqrt(Measurement.Foto_Dcrm * 10) + 0.861 * sqrt(Tree.H * 10)) ^ 2 + 0.742) / 10 Case 3 Measurement.dbh_lidar_Dcrm = ((-3.341 + 1.143 * sqrt(sade(maxA) * 10) + 0.7 * sqrt(Tree.H * 10)) ^ 2 + 0.78) / 10 d13Foto = ((-3.341 + 1.143 * sqrt(Measurement.Foto_Dcrm * 10) + 0.7 * sqrt(Tree.H * 10)) ^ 2 + 0.78) / 10 Case 20 Measurement.dbh_lidar_Dcrm = ((-3.214 + 1.016 * sqrt(sade(maxA) * 10) + 0.861 * sqrt(Tree.H * 10)) ^ 2 + 0.811) / 10 d13Foto = ((-3.214 + 1.016 * sqrt(Measurement.Foto_Dcrm * 10) + 0.861 * sqrt(Tree.H * 10)) ^ 2 + 0.811) / 10 End Select Cpu = "" Cpu = Cpu & Format$(MPlot, "0") & "," Cpu = Cpu & Format$(X_sol, "0.00") & "," Cpu = Cpu & Format$(Y_sol, "0.00") & "," Cpu = Cpu & Format$(Z_sol, "0.00") & "," Cpu = Cpu & Format$(LPE_SP(Ntrees), "0") & "," Cpu = Cpu & Format$(Zgnd, "0.00") & "," Cpu = Cpu & Format$(Tree.cf, "0.0000") & "," Cpu = Cpu & Format$(Tree.cont2, "0.0000") & "," Cpu = Cpu & Format$(Tree.pw, "0.0000") & "," Cpu = Cpu & Format$(rmse, "0.00") & "," Cpu = Cpu & Format$(HlidMax, "0.00") & "," Cpu = Cpu & Format$(He_, "0.00") & "," Cpu = Cpu & Format$(sade(maxA), "0.00") & "," Cpu = Cpu & Format$(Measurement.dbh_lidar_Dcrm, "0.0") & "," Cpu = Cpu & Format$(Hscale, "0.00") & "," Cpu = Cpu & Format$(DensityFactor, "0.00") & "," Cpu = Cpu & Format$(P, "0.0") & "," Rem Compute LiDAR-statistics Dim IntSum As Double, IntSum2 As Double, IntN As Double, IntMin As Double, IntMax As Double IntSum = 0: IntSum2 = 0: IntN = 0 IntMin = 5000: IntMax = 0 For i = 1 To Tree.NLP resid = XYdist(i) - sade(i) If Abs(resid) < rmse And intensities(i) > 1 Then IntSum = IntSum + intensities(i) IntSum2 = IntSum2 + intensities(i) ^ 2 IntN = IntN + 1 IntMin = Min(IntMin, intensities(i)) IntMax = Max(IntMax, intensities(i)) End If Next i If IntN > 1 Then IntMean = IntSum / IntN IntSD = Sqr((IntN * IntSum2 - IntSum ^ 2) / (IntN ^ 2)) End If Cpu = Cpu & Format$(IntMean, "0.0") & "," Cpu = Cpu & Format$(IntSD, "0.0") & "," Cpu = Cpu & Format$(IntMin, "0.0") & "," Cpu = Cpu & Format$(IntMax, "0.0") & "," Cpu = Cpu & Format$(IntN, "0") & "," For i = 0 To NumOfImages - 1 Fii = (90 - image_info(i).Sun_azimuth * TO_DEGREES) * TO_RADIANS ' azimuth (direction clockwise from KKJ-North in radians) Theta = (90 - image_info(i).Sun_elevation * TO_DEGREES) * TO_RADIANS ' sun's elevation in radians over horizon objtosun.X = 1 * Sin(Theta) * Cos(Fii) objtosun.Y = 1 * Sin(Theta) * Sin(Fii) objtosun.z = 1 * Cos(Theta) Call normalize(objtosun) Zdiff = 1 ' radius = (Tree.cf * Tree.H * Sin(4 * ((Zdiff) / H)) ^ Tree.pw + Tree.cont2) '* 1.2 PointInLight.X = X_sol + Radius * objtosun.X PointInLight.Y = Y_sol + Radius * objtosun.Y PointInLight.z = Z_sol - Zdiff PointInShade.X = X_sol + Radius * -objtosun.X PointInShade.Y = Y_sol + Radius * -objtosun.Y PointInShade.z = Z_sol - Zdiff Rem Come down 10% towards sun, on the surface of the crown model objtocam.X = image_info(i).Xo - X_sol objtocam.Y = image_info(i).Yo - Y_sol objtocam.z = image_info(i).Zo - Z_sol Call normalize(objtocam) Rem Using vector inner product arg = (objtocam.X * objtosun.X + objtocam.Y * objtosun.Y) / ((objtocam.X ^ 2 + objtocam.Y ^ 2) ^ 0.5 * (objtosun.X ^ 2 + objtosun.Y ^ 2) ^ 0.5) SunObjCamAngleXY = TO_DEGREES * MYFUNC_ACOS(arg) PhaseAngle = TO_DEGREES * vector_angle(objtocam, objtosun) NadirAngle = TO_DEGREES * vector_angle(objtocam, PlumbLine) InLight = -1 Form1.Picture1(i).DrawWidth = 1 If image_info(i).Num_of_addit_expos > 0 And (SunObjCamAngleXY < 45 And SunObjCamAngleXY > -45) Then InLight = 1 Call Get_value_from_16_byte_image(CLng(i), PointInLight.X, PointInLight.Y, PointInLight.z, CLng(2), CInt(1), imdata, Pan_mean, Pan_Max, Pan_Min, Pan_SD) Call r_transform_ground_to_pixel(i, PointInLight.X, PointInLight.Y, PointInLight.z, p_x, p_y) Rem Cyan '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), RGB(255, 0, 255) End If If image_info(i).Num_of_addit_expos > 0 And (SunObjCamAngleXY > 135 Or SunObjCamAngleXY < -135) Then Call Get_value_from_16_byte_image(CLng(i), PointInShade.X, PointInShade.Y, PointInShade.z, CLng(2), CInt(1), imdata, Pan_mean, Pan_Max, Pan_Min, Pan_SD) InLight = 0 Rem Yellow Call r_transform_ground_to_pixel(i, PointInShade.X, PointInShade.Y, PointInShade.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, -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) End If Form1.Picture1(i).DrawWidth = 1 Rem Assign to public string Cpu Cpu = Cpu & image_info(i).ImageCode & "," Cpu = Cpu & InLight & "," Cpu = Cpu & Format$(SunObjCamAngleXY, "0.00") & "," Cpu = Cpu & Format$(PhaseAngle, "0.00") & "," Cpu = Cpu & Format$(NadirAngle, "0.00") & "," Cpu = Cpu & Format$(imdata.r, "0") & "," Cpu = Cpu & Format$(imdata.G, "0") & "," Cpu = Cpu & Format$(imdata.B, "0") & "," Cpu = Cpu & Format$(imdata.NIR, "0") & "," Cpu = Cpu & Format$(Pan_mean, "0.0") & "," Cpu = Cpu & Format$(Pan_Max, "0.0") & "," Cpu = Cpu & Format$(Pan_Min, "0.0") & "," Cpu = Cpu & Format$(Pan_SD, "0.0") & "," Next i Rem TestOutput 'Open "c:\test.txt" For Append As 1 'Print #1, Cpu 'Close (1) 'Cpu = "" Measurement.cf = Tree.cf Measurement.cont = Tree.cont2 Measurement.pw = Tree.pw Measurement.Crown_RMSE = rmse Measurement.H_LiDAR = HlidMax Measurement.TreeSpecies = LPE_SP(Ntrees) Measurement.z = Z_sol Measurement.X = X_sol Measurement.Y = Y_sol Measurement.Z_dtm = Zgnd Measurement.Foto_h = Z_sol - Zgnd Measurement.num = Measurement.num + 1 Measurement.Lidar_Dcrm = sade(maxA) LidarFitDone = True 'Form1.Label10.Caption = "Plot: " & MPlot & "Dcr: " & Format$(sade(maxA), "0.0 m") & " L-h: " & Format$(HlidMax, "0.0 m") & " F-h: " & Format$(Tree.H, "0.0 m") & " d13: " & Format$(Measurement.dbh_lidar_Dcrm, "0.0 ") Form1.Label10.FontBold = False Form1.Label10.Caption = "Dcr: " & Format$(sade(maxA), "0.0 m") & " H: " & Format$(Measurement.Foto_h, "0.0 m") & " cf: " & Format$(Tree.cf, "0.00") & " pw: " & Format$(Tree.pw, "0.00") & " cont: " & Format$(Tree.cont2 ^ 2, "0.00") & " rms: " & Format$(rmse, "0.00") Call View_Globals NotFound: Close (6) 'Open "c:\mk_lidar.txt" For Append As 6 'Print #6, Format$(Tree.X, "0.00") & "," & Format$(Tree.Y, "0.00") & "," & Format$(Tree.z, "0.00") & "," & Format$((Tree.cf * Tree.h * Sin(4 * (0.3)) ^ Tree.pw + Tree.cont2) * 2, "0.0 ") & "," & Format$(Measurement.RMSE, "0.00") & "," & Format$(Tree.h, "0.00") & "," & Format$(X_sol, "0.00") & "," & Format$(Y_sol, "0.00") & "," & Format$(Z_sol, "0.00") & "," & Format$(getheight(X_sol, Y_sol), "0.00") & "," & Format$(Xm, "0.00") & "," & Format$(Ym, "0.00") & "," & Format$(TreeHeight, "0.00") & "," & Format$(TreeDiam, "0.00") 'Close (6) Rem ******************* Rest of the trees ****************** 'Call Mark_And_Plot(RMSE, Hlimit, Hlowest, CLng(maxA), Ntrees, LPE_SP(), Hftau(), Hfind(), sade(), xydist(), zdist(), Trees(), LPTreeHit()) Exit Sub ErrorInMatchLocalCrown: Close MsgBox ("Error occurred, resuming") End Sub