Python – projecting points onto a plane given by a normal and a point

mathnumpypython

I am using the following function to try to project bunch of vertices onto a plane. Basically mapping a polyhedron to polygon. My plane is defined by a normal vector and a point (called centroid here).

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy import linalg as la
#######################################################################################################
#this part is for drawing the vector arrow
#copy pasted
#http://stackoverflow.com/questions/22867620/putting-arrowheads-on-vectors-in-matplotlibs-3d-plot
#######################################################################################################
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import matplotlib.cm as cm#for color selection
import itertools
class Arrow3D(FancyArrowPatch):
    def __init__(self, xs, ys, zs, *args, **kwargs):
        FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
        self._verts3d = xs, ys, zs

    def draw(self, renderer):
        xs3d, ys3d, zs3d = self._verts3d
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
        self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
        FancyArrowPatch.draw(self, renderer)
#######################################################################################################
#function to project onto the plane
#takes normal vector of the plane, the centroid of the plane and vertices to be projected as argument
#returns the position vector of the projected point on the plane
#######################################################################################################
def planeprojection(normalvector,centroid,vertices):
    shape = vertices.shape#shape of vertex array, can be one vertex or multiple vertices to project
    if len(shape)==1:#meaning there is only one vertex
        vertex = vertices
        #dot product of position vector to the vertex from plane and normal vector
        dotscalar = np.dot(np.subtract(vertex,centroid),normalvector)
        #now returning the position vector of the projection onto the plane
        return np.subtract(vertex,dotscalar*normalvector)
    else:
        #array to store projectedvectors
        projectedvectors = np.zeros((shape[0],shape[1]))
        #now projecting onto plane, one by one
        for counter in range(shape[0]):
            vertex = vertices[counter,:]
            print vertex
            dotscalar = np.dot(np.subtract(centroid,vertex),normalvector)
            print dotscalar
            #now returning the position vector of the projection onto the plane
            projectedvectors[counter,:] = np.subtract(vertex,dotscalar*normalvector)
        #now returning the vectors projected
        return projectedvectors
#######################################################################################################
###normalising function
#######################################################################################################
def normalise(v):
    norm=np.linalg.norm(v)
    if norm==0: 
       return v
    return v/norm
##############################################################################################################
##########points to project####################
##############################################################################################################
xcod = np.array([1,2,1,3,-1,1])
ycod = np.array([2,1,4.5,5.,6,2])
zcod = np.array([1,-2,0,2,3,1])
num = len(xcod)-1
#centroid of the cell
centroid = np.array([np.mean(xcod[0:num]),np.mean(ycod[0:num]),np.mean(zcod[0:num])])
#getting tuples of x,y,z
verts = [zip(xcod,ycod,zcod)]
#numpy array of vertices
vertices =np.array(verts[0])
#normal to the plane
averagenormal = np.array([ 0.91008281, -0.24978471,  0.3306915 ])
#Projecting the vertices now
projectedvertices = planeprojection(averagenormal,centroid,vertices)
#changing the format to tuple for plotting polyhedron surface
projectedverts = [map(tuple,projectedvertices)]
################################################################################
######plotting #################################################################
################################################################################
#also defining the plot in 3d for start plotting
fig = plt.figure()
ax = fig.add_subplot(111,projection = '3d')
#plotting all the points
ax.plot(xcod,ycod,zcod,'x-')
#plotting polyhedron surface
ax.add_collection3d(Poly3DCollection(projectedverts,alpha=0.5,color="k"))
#adding labels for vertice
for i in range(num):
    ax.text(xcod[i],ycod[i],zcod[i],'%d(%.2f,%.2f,%.2f)'%(i,xcod[i],ycod[i],zcod[i]))
#drawing vector arrow for averal normal vector
drawvec=Arrow3D([centroid[0],averagenormal[0]+centroid[0]],[centroid[1],averagenormal[1]+centroid[1]],[centroid[2],averagenormal[2]+centroid[2]],mutation_scale=8,lw=1,color='k')
ax.add_artist(drawvec)
#draw average vector from origin as originally found out
drawvec = Arrow3D([0,averagenormal[0]],[0,averagenormal[1]],[0,averagenormal[2]],mutation_scale=8,lw=1,color='g')
ax.add_artist(drawvec)
#plotting averagenormal vector point
ax.scatter(averagenormal[0],averagenormal[1],averagenormal[2],marker = 'o',color="g")
#plotting centroid 
ax.scatter(centroid[0],centroid[1],centroid[2],marker = 'o',color="g")
#plotting averagenormal vector point from centroid
ax.scatter(averagenormal[0]+centroid[0],averagenormal[1]+centroid[1],averagenormal[2]+centroid[2],marker = 'o',color="g")#plot show
ax.set_xlim([-1,5])
ax.set_ylim([-1,6])
ax.set_zlim([-5,6])
plt.show()

Here, in the figure below, black shaded part is supposed to be the plane. Something is not working right shouldnt this be the right way to go?
I am basically using following formula for projection, where Proj(P) is the projection onto plane, Centriod is the point on the plane and n is the normal vector.

Proj(P) = P-((P-Centroid).n)n

EDIT

Changing the mistake suggested by Nico below, dotscalar = np.dot(np.subtract(vertex,centroid),normalvector) it works !

Best Answer

Your formula

Proj(P) = P-((P-Centroid).n)n

is correct.

But your implementation

dotscalar = np.dot(np.subtract(centroid,vertex),normalvector)

is not. You swapped centroid and vertex. Thus, instead of moving the point onto the plane, your algorithm doubles its distance from the plane.

Related Solutions

Python – How to pretty-print a numpy.array without scientific notation and with given precision

You can use set_printoptions to set the precision of the output:

import numpy as np
x=np.random.random(10)
print(x)
# [ 0.07837821  0.48002108  0.41274116  0.82993414  0.77610352  0.1023732
#   0.51303098  0.4617183   0.33487207  0.71162095]

np.set_printoptions(precision=3)
print(x)
# [ 0.078  0.48   0.413  0.83   0.776  0.102  0.513  0.462  0.335  0.712]

And suppress suppresses the use of scientific notation for small numbers:

y=np.array([1.5e-10,1.5,1500])
print(y)
# [  1.500e-10   1.500e+00   1.500e+03]
np.set_printoptions(suppress=True)
print(y)
# [    0.      1.5  1500. ]

See the docs for set_printoptions for other options.

To apply print options locally, using NumPy 1.15.0 or later, you could use the numpy.printoptions context manager. For example, inside the with-suite precision=3 and suppress=True are set:

x = np.random.random(10)
with np.printoptions(precision=3, suppress=True):
    print(x)
    # [ 0.073  0.461  0.689  0.754  0.624  0.901  0.049  0.582  0.557  0.348]

But outside the with-suite the print options are back to default settings:

print(x)    
# [ 0.07334334  0.46132615  0.68935231  0.75379645  0.62424021  0.90115836
#   0.04879837  0.58207504  0.55694118  0.34768638]

If you are using an earlier version of NumPy, you can create the context manager yourself. For example,

import numpy as np
import contextlib

@contextlib.contextmanager
def printoptions(*args, **kwargs):
    original = np.get_printoptions()
    np.set_printoptions(*args, **kwargs)
    try:
        yield
    finally: 
        np.set_printoptions(**original)

x = np.random.random(10)
with printoptions(precision=3, suppress=True):
    print(x)
    # [ 0.073  0.461  0.689  0.754  0.624  0.901  0.049  0.582  0.557  0.348]

To prevent zeros from being stripped from the end of floats:

np.set_printoptions now has a formatter parameter which allows you to specify a format function for each type.

np.set_printoptions(formatter={'float': '{: 0.3f}'.format})
print(x)

which prints

[ 0.078  0.480  0.413  0.830  0.776  0.102  0.513  0.462  0.335  0.712]

instead of

[ 0.078  0.48   0.413  0.83   0.776  0.102  0.513  0.462  0.335  0.712]

Improved Area Lighting in WebGL & ThreeJS

The good news is there is a solution; but first the bad news.

Your approach of using the point that maximizes the dot product is fundamentally flawed, and not physically plausible.

In your first illustration above, suppose that your area light consisted of only the left half.

The "purple" point -- the one that maximizes the dot-product for the left half -- is the same as the point that maximizes the dot-product for both halves combined.

Therefore, if one were to use your proposed solution, one would conclude that the left half of the area light emits the same radiation as the entire light. Obviously, that is impossible.

The solution for computing the total amount of light that the area light casts on a given point is rather complicated, but for reference, you can find an explanation in the 1994 paper The Irradiance Jacobian for Partially Occluded Polyhedral Sources here.

I suggest you look at Figure 1, and a few paragraphs of Section 1.2 -- and then stop. :-)

To make it easy, I have coded a very simple shader that implements the solution using the three.js WebGLRenderer -- not the deferred one.

EDIT: Here is an updated fiddle: http://jsfiddle.net/hh74z2ft/1/

enter image description here

The core of the fragment shader is quite simple

// direction vectors from point to area light corners

for( int i = 0; i < NVERTS; i ++ ) {

    lPosition[ i ] = viewMatrix * lightMatrixWorld * vec4( lightverts[ i ], 1.0 ); // in camera space

    lVector[ i ] = normalize( lPosition[ i ].xyz + vViewPosition.xyz ); // dir from vertex to areaLight

}

// vector irradiance at point

vec3 lightVec = vec3( 0.0 );

for( int i = 0; i < NVERTS; i ++ ) {

    vec3 v0 = lVector[ i ];
    vec3 v1 = lVector[ int( mod( float( i + 1 ), float( NVERTS ) ) ) ]; // ugh...

    lightVec += acos( dot( v0, v1 ) ) * normalize( cross( v0, v1 ) );

}

// irradiance factor at point

float factor = max( dot( lightVec, normal ), 0.0 ) / ( 2.0 * 3.14159265 );

More Good News:

This approach is physically correct.
Attenuation is handled automatically. ( Be aware that smaller lights will require a larger intensity value. )
In theory, this approach should work with arbitrary polygons, not just rectangular ones.

Caveats:

I have only implemented the diffuse component, because that is what your question addresses.
You will have to implement the specular component using a reasonable heuristic -- similar to what you already have coded, I expect.
This simple example does not handle the case where the area light is "partially below the horizon" -- i.e. not all 4 vertices are above the plane of the face.
Since WebGLRenderer does not support area lights, you can't "add the light to the scene" and expect it to work. This is why I pass all necessary data into the custom shader. ( WebGLDeferredRenderer does support area lights, of course. )
Shadows are not supported.

three.js r.73

EDIT

Best Answer

Related Solutions

Python – How to pretty-print a numpy.array without scientific notation and with given precision

Improved Area Lighting in WebGL & ThreeJS

Related Topic