In [29]:
# imports
import math as m
import numpy as np
import pandas as pd
import pickle
import os
from scipy import integrate
In [2]:
# creating dataframe from csv file
print('Loding survey data')
dfCoordinates = pd.read_csv('gpsDataUTM.csv') # reads data into dataframe
In [3]:
dfCoordinates.head()
Out[3]:
In [4]:
# Adding new info to dataframe: distance from origin (common reference to all points)
def distanceFromOrigin(x,y,z):
# compute distance from origin <0,0,0>
point = np.array([x,y,z]) # place columns into an array <x,y,z>
origin = np.zeros(3) # create origin <0,0,0>
return np.linalg.norm(point-origin) # compute norm
dfCoordinates['oDist'] = np.vectorize(distanceFromOrigin)(dfCoordinates['x_rel'],
dfCoordinates['y_rel'],
dfCoordinates['z_rel'])
In [5]:
dfCoordinates.head()
Out[5]:
In [6]:
# sorts accordint to distance from the origin.
dfCoordinates = dfCoordinates.sort_values(by=['oDist'],
ascending=True)
In [11]:
dfCoordinates.head()
Out[11]:
In [30]:
def computeLineEquation(A,B):
# computes line equation given 2 points: A and B
# Both A and B need to be numpy arrays
if (B[0]-A[0]) == 0: # Check if line is parallel to y axis (zero divide)
#print("Zero division. Line is parallel to y axis. x={x}".format(x=A[0]))
return np.array([None,None])
else: # detect slope and intersection
m = (B[1]-A[1])/(B[0]-A[0]) # slope
xo = A[0]
yo = A[1]
c = -m*xo+yo
#print("The line equation is y = {m}x + {c}".format(m=m,c=c))
return np.array([m,c]) # returns slope and constant in numpy array format
In [150]:
def computeVolume(A,B,C):
# computes volume given 3 points using double integral, line and plane equations
# A, B and C are numpy arrays
# Setting up plane equation
AB = B-A # vector within plane
BC = C-B # vector within plane
NV = np.cross(AB,BC) # vector that is normal to plane
# Plane values; plane equation is a*(x-xo)+b*(y-yo)+c*(z-zo)=0
a = NV[0]
b = NV[1]
c = NV[2]
xo = A[0]
yo = A[1]
zo = A[2]
# Compute line equations:
slopeAB, constantAB = computeLineEquation(A,B)
slopeBC, constantBC = computeLineEquation(B,C)
slopeAC, constantAC = computeLineEquation(A,C)
# Compute double integral
if ((slopeAC is not None) and (slopeAB is not None)): # if slopes are defined
volume1, error1 = integrate.dblquad(lambda y, x: ((-a*(x-xo)-b*(y-yo))/c)+zo,
A[0],
B[0],
lambda x: slopeAB*x+constantAB,
lambda x: slopeAC*x+constantAC)
else:
volume1 = 0.0
if ((slopeAC is not None) and (slopeBC is not None)): # if slopes are defined
volume2, error2 = integrate.dblquad(lambda y, x: ((-a*(x-xo)-b*(y-yo))/c)+zo,
B[0],
C[0],
lambda x: slopeAC*x+constantAC,
lambda x: slopeBC*x+constantBC)
else:
volume2 = 0.0
return abs(volume1)+abs(volume2) # return total volume
In [155]:
# Required number of loops to compute total volume = number of points - 1
totalVolume = 0.0 # initialize volume summation holder
for i in range(dfCoordinates['x'].count()-2): # iterate thru all groups of 3 points in ascending order from origin.
s1 = dfCoordinates.iloc[i] # get series (row entry)
s2 = dfCoordinates.iloc[i+1]
s3 = dfCoordinates.iloc[i+2]
dfAnalysis = pd.DataFrame([]) # build dataframe composed of these 3 entry set
dfAnalysis = dfAnalysis.append(s1,ignore_index=True)
dfAnalysis = dfAnalysis.append(s2,ignore_index=True)
dfAnalysis = dfAnalysis.append(s3,ignore_index=True)
dfAnalysis = dfAnalysis.sort_values(by=['x_rel'], ascending=True)# order by lowest x_rel >> A,B,C
A = np.array([dfAnalysis.iloc[0]['x_rel'],
dfAnalysis.iloc[0]['y_rel'],
dfAnalysis.iloc[0]['z_rel']])
B = np.array([dfAnalysis.iloc[1]['x_rel'],
dfAnalysis.iloc[1]['y_rel'],
dfAnalysis.iloc[1]['z_rel']])
C = np.array([dfAnalysis.iloc[2]['x_rel'],
dfAnalysis.iloc[2]['y_rel'],
dfAnalysis.iloc[2]['z_rel']])
# Compute line equations, plane equations and volume
volume = computeVolume(A,B,C)
totalVolume += volume # add to total volume
print('Total volume is: {v} cubic meters.'.format(v=totalVolume))
In [92]:
A = np.array([1.2462402103701609, -3.6996389972046018, -0.83000000000004093])
B = np.array([17.149493328179236, -10.234135422855616, -0.87000000000000455])
In [93]:
np.cross(A,B)
Out[93]:
In [11]:
m.sqrt(m.pow(19.848751,2)+m.pow(7.76,2)) # check!
Out[11]:
In [11]:
slope, constant = computeLineEquation(np.array([1,2]), np.array([1,15]))
In [5]:
# the same but different...
slope, constant = computeLineEquation(np.array([1,2]), np.array([4,5]))
In [51]:
# Known volume: 25
In [52]:
# Setting up points in 3D space
A = np.array([0.0,
3.0,
4.0])
B = np.array([0.0,
0.0,
4.0])
C = np.array([5.0,
2.5,
2.0])
# Setting up plane equation
AB = B-A # vector within plane
BC = C-B # vector within plane
NV = np.cross(AB,BC) # vector that is normal to plane
# Plane values
a = NV[0]
b = NV[1]
c = NV[2]
xo = A[0]
yo = A[1]
zo = A[2]
print("The plane equation is {a}*(x-{xo})+{b}*(y-{yo})+{c}*(z-{zo})=0".format(a=a,b=b,c=c,xo=xo,yo=yo,zo=zo))
# With these values, you can fully define the plane equation
# computing plane equations
slopeAB, constantAB = computeLineEquation(A,B)
slopeBC, constantBC = computeLineEquation(B,C)
slopeAC, constantAC = computeLineEquation(A,C)
# Compute volumes
from scipy import integrate
if ((slopeAC is not None) and (slopeAB is not None)): # if slopes are defined
volume1, error1 = integrate.dblquad(lambda y, x: ((-a*(x-xo)-b*(y-yo))/c)+zo,
A[0],
B[0],
lambda x: slopeAB*x+constantAB,
lambda x: slopeAC*x+constantAC)
else:
volume1 = 0.0
if ((slopeAC is not None) and (slopeBC is not None)): # if slopes are defined
volume2, error2 = integrate.dblquad(lambda y, x: ((-a*(x-xo)-b*(y-yo))/c)+zo,
B[0],
C[0],
lambda x: slopeAC*x+constantAC,
lambda x: slopeBC*x+constantBC)
else:
volume2 = 0.0
volumeTotal = abs(volume1)+abs(volume2)
print("Volume 1: {v}".format(v=volume1))
print("Volume 2: {v}".format(v=volume2))
print("Total computed volume is {v}".format(v=volumeTotal))
So at the end of this chapter, we know how to calculate the volume of a solid object formed by a set of 3 points in space. Now it is a matter of ordering how this computations are to be made for a point cloud and sum all of them for total volume.