In [1]:
# imports
import math as m
import numpy as np
import pandas as pd
import pickle
import os
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]:
# Example: computing distance between points
A = dfCoordinates[['x_rel','y_rel','z_rel']].values[0]
B = dfCoordinates[['x_rel','y_rel','z_rel']].values[50]
O = np.array([0.0,0.0,0.0]) # origin
In [5]:
distanceAB = np.linalg.norm(B-A)
distanceOA = np.linalg.norm(A-O)
distanceOB = np.linalg.norm(B-O)
print('Distance AB: %s' % distanceAB)
print('Distance OA: %s' % distanceOA)
print('Distance OB: %s' % distanceOB)
In [6]:
# Adding new info to dataframe
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 [7]:
dfCoordinates.head()
Out[7]:
In [8]:
# sorts accordint to distance from the origin.
dfCoordinates = dfCoordinates.sort_values(by=['oDist'],
ascending=True)
In [9]:
pd.__version__
Out[9]:
In [10]:
dfCoordinates.head()
Out[10]:
In [11]:
m.sqrt(m.pow(19.848751,2)+m.pow(7.76,2)) # check!
Out[11]:
In [15]:
# Useful computing. Compute line equation given 2 points A, B.
A = np.array([1,2])
B = np.array([4,5])
# line equation: (y-yo) = m*(x-xo)
m = (A[1]-B[1])/(A[0]/B[0]) # slope
xo = A[0]
yo = A[1]
c = -m*xo+yo
print("Line Solution is y = {m}x + {c}".format(m=m,c=c))
In [9]:
# Defining a generic function to comput the above
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 [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 [6]:
slope, constant = computeLineEquation(np.array([4,5]), np.array([10,1]))
In [7]:
slope, constant = computeLineEquation(np.array([1,2]), np.array([10,1]))
In [46]:
# Example
x = np.linspace(0.0,100.0,num=500) # generate linear space for x values
y = np.linspace(0.0,100.0,num=500)
xGrid, yGrid = np.meshgrid(x,y) # generate mesh grid for plotting sample data
point1 = np.array([[1.0,
2.0,
1.0]])
point2 = np.array([[4.0,
5.0,
0.0]])
point3 = np.array([[10.0,
1.0,
3.0]])
planeVectorOne = point2 - point1 # a vector within the plane
planeVectorTwo = point3 - point2 # a vector within the plane
perpendicularVector = np.cross(planeVectorOne, planeVectorTwo) # plane equation coeficients obtained <a,b,c>
# print(perpendicularVector)
# print(type(perpendicularVector))
# print(perpendicularVector.shape)
# print(point1.shape)
# print(perpendicularVector[0][2])
# saving constants to literals
a = perpendicularVector[0][0]
b = perpendicularVector[0][1]
c = perpendicularVector[0][2]
xo = point1[0][0]
yo = point1[0][1]
zo = point1[0][2]
# Plane equation a(x-xo)+b(y-yo)+c(z-zo)=0 with isolated z:
z = ((-a*(xGrid-xo)-b*(yGrid-yo))/c)+zo
In [47]:
# Steps:
# 1. Compute plane equation from 3 points. 1 plane.
# 2. Compute line equations. 3 equations.
# 3. Sort them according to distance. A >> B >> C
# 4. Volume1 = Double integral xA-xB lineEqAC lineEqAB
# 5. Volume2 = Double integral xB-xC lineEqAC lineEqBC
# 6. VolumeTotal = Volume1 + Volume2
In [49]:
# Setting up sample space
x = np.linspace(0.0,11.0,num=50) # generate linear space for x values
y = np.linspace(0.0,11.0,num=50)
xGrid, yGrid = np.meshgrid(x,y) # generate mesh grid for plotting sample data
In [62]:
# Setting up points in 3D space
A = np.array([1.0,
2.0,
1.0])
B = np.array([4.0,
5.0,
0.0])
C = np.array([10.0,
1.0,
3.0])
In [63]:
# 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
In [71]:
# 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
In [68]:
# computing plane equations
slopeAB, constantAB = computeLineEquation(A,B)
slopeBC, constantBC = computeLineEquation(B,C)
slopeAC, constantAC = computeLineEquation(A,C)
In [78]:
# Compute volumes
from scipy import integrate
volume1, error1 = integrate.dblquad(lambda y, x: ((-a*(x-xo)-b*(y-yo))/c)+zo,
A[0],
B[0],
lambda x: slopeAC*x+constantAC,
lambda x: slopeAB*x+constantAB)
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)
volumeTotal = volume1+volume2
print("Total computed volume is {v}".format(v=volumeTotal))
In [12]:
# Setting up points in 3D space
A = np.array([1.0,
2.0,
10.0])
B = np.array([3.0,
5.0,
10.0])
C = np.array([5.0,
2.0,
10.0])
# A = np.array([1.0,
# 2.0,
# 10.0])
# B = np.array([5.0,
# 1.0,
# 10.0])
# C = np.array([5.0,
# 5.0,
# 10.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
volume1, error1 = integrate.dblquad(lambda y, x: ((-a*(x-xo)-b*(y-yo))/c)+zo,
A[0],
B[0],
lambda x: slopeAC*x+constantAC,
lambda x: slopeAB*x+constantAB)
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)
volumeTotal = volume1+volume2
print("Total computed volume is {v1}+{v2}={v}".format(v1=volume1,v2=volume2,v=volumeTotal))
# Volume Computed Correctly! :)
# Check if it works with line equations like x=3 -- parallel to y axis; Nope! Divide by 0 error. Analysing...
# -- when divide by zero occurs, corresponding volume is zero! :) find a way around this on code.
In [46]:
# Setting up points in 3D space
A = np.array([1.0,
3.0,
10.0])
B = np.array([5.0,
5.0,
10.0])
C = np.array([11.0,
3.0,
10.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))
# Volume Computed Correctly! :)
# Check if it works with line equations like x=3 -- parallel to y axis; Nope! Divide by 0 error. Analysing...
# -- when divide by zero occurs, corresponding volume is zero! :) find a way around this on code.
In [47]:
# Volume pyramid: 3x2x5/3 = 10
In [49]:
# Setting up points in 3D space
A = np.array([0.0,
0.0,
2.0])
B = np.array([0.0,
3.0,
2.0])
C = np.array([5.0,
2.5,
0.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))
In [50]:
# IT WORKED!
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.