Generate a Finite Plane of a Given Size and Orientation Using its Normal Vector. OpenGL, Open Scene Graph.

*) Given a Plane with equation Ax+By+Cz+D=0 find two director vectors.
*) Once obtained, vectors VD1 and VD2, find a third vector in the same plane VP that is perpendicular to VD1

plane ecuation: Ax+By+Cz+D=0, being (A,B,C) the director vector of the plane
Getting two points in the plane:
Given x=0, y=0 => Cz+D=0 -> Cz=-D/C ==> PointA=(0,0,-D/C)
Given x=1, y=0 => A+Cz+D=0 -> Cz= -D-A -> z= (-D-A)/C   ==> PointB=(1,0,(-D-A)/C)
Given x=0, y=1 => B+Cz+D=0 -> Cz= -D-B -> z= (-D-B)/C   ==> PointC=(0,1,(-D-B)/C)

This was wrong, check code comments, solved there

directorVectorOne= VectorFromPointAToPointB= (1,0, ((-D-A)/C)-(-D/C) )
= (1,0, (-D-A+D)/C ) = (1,0,-A/C) )
directorVectorTwo= VectorFromPointAToPointC= (0,1, ((-D-B)/C)-(-D/C) )
= (0,1, (-D-B+D)/C ) = (0,1, -B/C) )

TestVector_BC=(-1,1,(-D-B)/C-(-D-A)/C)) = (-1,1,(-B+A)/C)

Example: 2x+y+8z-21=0
A=2 B=1 C=8  : Plane Normal Vector= (2,1,8)
DVector1=(1,0,-A/C) = (1,0,-0.25)
DVector2=(0,1,-B/C) = (0,1,-0.125)

TestVector_BC= (-1,1,(-B+A)/C) = (-1,1,0.125)

Property of determinants: If tree vectors are in the same plane, then their
determinant is 0:
Det(DVector1,DVector2,TestVector_BC)=0 (checked)
*/
double A,B,C,D;
osg::Vec4 pvector= _plane.asVec4();
A=pvector[0];
B=pvector[1];
C=pvector[2];
D=pvector[3];
_directorVectorOne=osg::Vec3(1.0f,0.0f,-A/C);
_directorVectorOne.normalize();
_directorVectorTwo=osg::Vec3(0,1.0f,-B/C);
_directorVectorTwo.normalize();
}

void calculateQuadranglePoints()
{
/*We got two director vectors for the plane
VD1=(1,0,-A/C)
VD2=(0,1,-B/C) already calculated in calculateDirectorVectors()
I want a third vector VDP in the same plane thanVD1 and VD2 but
perpendicular to VD1, that is, VD1 dot VDP is zero.

Ax+By+Cz+D=0 for  x=1,y=0, z=-A/C obtained before

det(VD1,VD2,VDP) =0 as they are coplanar

VD1 dot VDP =0 as they are perpendicular.

Lets find VDP
VD1=(f,g,h) VD2=(u,v,w) VDP=(p,q,r)
p,q,r are our unknowns

Ap+Bq+Cr+D=0;
det(VD1,VD2,VD3)= det( (g h),(v w) )*p - det( (f h),(u w) )*q+ det( (f g),(u v) )*r=0
fp+gq+hr=0;

That is, the system to solve for p,q,r is:

(gw-hv)*p-(fw-hu)*q+(fv-gu)*r=0
Ap+Bq+Cr+D=0
fp+gq+hr=0

Solve by Kramer Rule


| (gw-hv) (fw-hu) (fv-gu) |
|A|= | A         B               C |
| f         g                h |

R=  ( 0, -D , 0)



        | 0 (fw-hu) (fv-gu) |
|Ap|= | -D B                 C |
        | 0 g                h |



        | (qw-hv) 0 (fv-gu) |
|Aq|= | A         -D         C |
        | f         0           h |




| (qw-hv) (fw-hu) 0 |
|Ar|= | A            B         -D |
| f            g           0 |

---->
|A| = -A*( h(fw-hu) - g(fv-gu) ) + B* (h(gw-hv) - f(fv-gu)) - C(g(gw-hv)-f(fw-hu))
|Ap|= D*(h(fw-hu)-g(fv-gu))
|Aq|= (-D)*(h(gw-hv)-f(fv-gu))
|Ar|=D*(g(gw-hv)-f(fw-hu))

so that:
p=|Ap| / |A|
q=|Aq| / |A|
r=|Ar| / |A|

The resulting vector (p,r,q) is coplanar to VD1 and perpendicular to it.

*/

//////////////////////////////////////////////
// A limited plane calculated from its normal
// Used for testing, no real use now
// EXAMPLE USAGE:
//osg::Vec3 pnormal=osg::Vec3(1,0.001,0.001);
// pnormal.normalize();
// osg::Plane plane=osg::Plane(pnormal,20);
// GeoPlane::GeoPlane(plane,50);

#pragma once
#include <iostream>
#include <osg/ref_ptr>
#include <osgViewer/Viewer>
#include <osg/ShapeDrawable>
#include <OSGGeneralUtils.hh>
class GeoPlane
{
public:

GeoPlane(osg::Plane plane,double size)
{
osg::NotifyHandler *notifyhandler=new LogFileHandler("geoplanelog.txt");
osg::setNotifyLevel( osg::INFO );
osg::setNotifyHandler(notifyhandler);
_plane=plane;
calculateDirectorVectors();
calculateSecondPerpendicularDirectorVector();

_u=_directorVectorOne;
_v=_directorVectorPerpendicular;
_size=size;
calculatePoints();
}

virtual ~GeoPlane(void)
{

}

/**
Returns the first director vector u, perpendicular to v
**/
osg::Vec3 getU()
{
return _u;
}

/**
Returns the second director vector v, perpendicular to u
**/
osg::Vec3 getV()
{
return _v;
}

/**
Returns the plane normal.
**/
osg::Vec3 getNormal()
{
return _plane.getNormal();
}

/**
Returns the plane D factor.
**/
double getD()
{
return _plane.asVec4()[3];
}

/**
Returns the plane Drawing size factor.
**/
double getSize()
{
return _size;
}


private:


void calculateDirectorVectors()
{
/* plane ecuation: Ax+By+Cz+D=0, being (A,B,C) the director vector of the plane
Getting two points in the plane:
Given x=0, y=0 => Cz+D=0 -> Cz=-D/C ==> PointA=(0,0,-D/C)
Given x=1, y=0 => A+Cz+D=0 -> Cz= -D-A -> z= (-D-A)/C   ==> PointB=(1,0,(-D-A)/C)
Given x=0, y=1 => B+Cz+D=0 -> Cz= -D-B -> z= (-D-B)/C   ==> PointC=(0,1,(-D-B)/C)

//Becareful, if C is 0 then other points should be chosen, for example
// x=0,y=?, z=1 we will solve that later.

directorVectorOne= VectorFromPointAToPointB=
(1,0,z1-zo)=(1,0,-A/C)
//Wrong:(1,0, ((-D-A)/C)-(-D/C) )= (1,0, (-D-A+D)/C ) = (1,0,-A/C) )
directorVectorTwo= VectorFromPointAToPointC= (0,1.0f,-B/C)

//Wrong:(0,1, ((-D-B)/C)-(-D/C) )= (0,1, (-D-B+D)/C ) = (0,1, -B/C) )

directorVectorThree=(-1,1,(-B+A)/C);
Wrong:(-1,1,(-D-B)/C-(-D-A)/C)) = (-1,1,(-B+A)/C)

Example: 2x+y+8z-21=0
A=2 B=1 C=8  : Plane Normal Vector= (2,1,8)
DVector1=(1,0,-A/C) = (1,0,-0.25)
DVector2=(0,1,-B/C) = (0,1,-0.125)

TestVector_BC= (-1,1,(-B+A)/C) = (-1,1,0.125)

Property of determinants: If tree vectors are in the same plane, then their
determinant is 0:
Det(DVector1,DVector2,TestVector_BC)=0 (checked)
*/
double A,B,C,D;
osg::Vec4 pvector= _plane.asVec4();
A=pvector[0];
B=pvector[1];
C=pvector[2];
D=pvector[3];
/* Example:
p0=(x0,0,0)
p1=(x1,1,0)
p2=(x2,0,1)
p0p1=(x1-x0,1,0)
p0p2=(x2-x0,0,1)
p1p2=(x2-x1,-1,1)
x0: Ax0+D=0 --> x0=-D/A
x1: Ax1+B+D=0 --> x1=(-B-D)/A
x2: Ax2+C+D=0 --> x2=(-C-D)/A
==> p0p1=(-B/A,1,0);
p0p2=(-C/A,0,1)
p1p2=((-C+B)/A)
*/


if (C>0.0001)
{
_origin=osg::Vec3(0,0,-D/C);

_directorVectorOne=osg::Vec3(1.0f,0.0f,-A/C); //-(A+D)/C);
_directorVectorOne.normalize();
_directorVectorTwo=osg::Vec3(0,1.0f,-B/C); //-(B+D)/C);
_directorVectorTwo.normalize();
_directorVectorThree=osg::Vec3(-1,1,(-B+A)/C); //(A-B-D)/C);
_directorVectorThree.normalize();
} else if (B>0.0001)
{
_origin=osg::Vec3(0,-B/D,0);
_directorVectorOne=osg::Vec3(1.0f,-A/B,0); //(1.0f,-(A+D)/B,0);
_directorVectorOne.normalize();
_directorVectorTwo=osg::Vec3(0,-C/B,1); //(0,-(C+D)/B,1.0f);
_directorVectorTwo.normalize();
_directorVectorThree=osg::Vec3(-1,(-C+A)/B,1); //(-1,(A-C-D)/B,1);
_directorVectorThree.normalize();
} else
{
_origin=osg::Vec3(-D/A,0,0);
_directorVectorOne=osg::Vec3(-B/A,1,0); //(-(B+D)/A,1.0f,0);
_directorVectorOne.normalize();
_directorVectorTwo=osg::Vec3(-C/A,0,1);//(-(C+D)/A,0,1.0f);
_directorVectorTwo.normalize();
_directorVectorThree=osg::Vec3((-C+B)/A,-1,1);//((B-C-D)/A,-1,1);
_directorVectorThree.normalize();

}
OSG_INFO<<"Distance origin to plane:"<<pointPlaneDistance(_origin)<<std::endl;


}

void calculationsTestOne()
{
double dot=_directorVectorOne*_directorVectorTwo;
std::cout<<dot<<std::endl; //Normally should not be 0.
double det=OSGGeneralUtils::determinant(_directorVectorOne,_directorVectorTwo,_directorVectorThree);
std::cout<<det<<std::endl; // Should be 0 as they are coplanar
}

void calculateSecondPerpendicularDirectorVector()
{
 /*We got three director vectors for the plane
  VD1=(1,0,-A/C)
  VD2=(0,1,-B/C)
  VD3=(-1,1,(-B+A)/C)
  already calculated in calculateDirectorVectors()
  I want a third vector VDP in the same plane than VD1,VD2,VD3 but
  perpendicular to VD1, that is, VD1 dot VDP is zero. and
  det(VD1,VD3,VDP)=0;

  A*p+B*q+C*r+D=0; for being in the plane

  det(VD1,VD2,VDP) =0 as they are coplanar

  VD1 dot VDP =0 as they are perpendicular.

  Lets find VDP
  VD1=(f,g,h) VD2=(u,v,w) VDP=(p,q,r)
  p,q,r are our unknowns

  Ap+Bq+Cr+D=0;
  det(VD1,VD2,VDP)= det( (g h),(v w) )*p - det( (f h),(u w) )*q+ det( (f g),(u v) )*r=0
  fp+gq+hr=0;

  That is, the system to solve for p,q,r is:

 p*(g*w-h*v)-q*(f*w-h*u)+r*(f*v-g*u)=0;
  A*p+B*q+C*r+D=0
  f*p+g*q+h*r=0

 Solve by Kramer Rule


| (qw-hv)  (fw-hu) (fv-gu) |
|A|= | A          B               C |
| f          g               h |

R=  ( 0, -D , 0)



         | 0 (fw-hu)  (fv-gu) |
|Ap|= | -D B              C |
         | 0 g              h |



| (qw-hv) 0 (fv-gu) |
|Aq|= | A          -D C   |
| f          0 h |




| (qw-hv) (fw-hu) 0 |
|Ar|= | A             B          -D |
| f             g           0 |

---->
|A| = -A*( h(fw-hu) - g(fv-gu) ) + B* (h(gw-hv) - f(fv-gu)) - C(g(gw-hv)-f(fw-hu))
|Ap|= D*(h(fw-hu)-g(fv-gu))
|Aq|= (-D)*(h(gw-hv)-f(fv-gu))
|Ar|=D*(g(gw-hv)-f(fw-hu))

so that:
p=|Ap| / |A|
q=|Aq| / |A|
r=|Ar| / |A|

The resulting vector (p,r,q) is coplanar to VD1 and perpendicular to it.
*/
double A,B,C,D;
double DetA,DetAp,DetAq,DetAr;
double f,g,h,u,v,w,p,q,r;
osg::Vec4 pvector= _plane.asVec4();
A=pvector[0];
B=pvector[1];
C=pvector[2];
D=pvector[3];

f=_directorVectorOne[0];
g=_directorVectorOne[1];
h=_directorVectorOne[2];

u=_directorVectorTwo[0];
v=_directorVectorTwo[1];
w=_directorVectorTwo[2];

/*

DetA= -A*( h*(f*w-h*u) - g*(f*v-g*u) ) +
B* (h*(g*w-h*v) - f*(f*v-g*u)) -
C*(g*(g*w-h*v)-f*(f*w-h*u));
DetAp= D*(h*(f*w-h*u)-g*(f*v-g*u));
DetAq= (-D)*(h*(g*w-h*v)-f*(f*v-g*u));
DetAr= D*(g*(g*w-h*v)-f*(f*w-h*u));

p=DetAp/DetA;
q=DetAq/DetA;
r=DetAr/DetA;*/

DetA=((g*g+f*f)*w-g*h*v-f*h*u)*C+(-g*h*w+(h*h+f*f)*v-f*g*u)*B+(-f*h*w-f*g*v+(h*h+g*g)*u)*A;

p=((f*h*w+f*g*v+(-(h*h)-(g*g))*u)*D)/DetA;
q=((g*h*w+(-(h*h)-f*f)*v+f*g*u)*D)/DetA;
r=-(((g*g+f*f)*w-g*h*v-f*h*u)*D)/DetA;

_directorVectorPerpendicular[0]=p;
_directorVectorPerpendicular[1]=q;
_directorVectorPerpendicular[2]=r;
_directorVectorPerpendicular.normalize();
calculationsTestTwo();
}

void calculatePoints()
{
osg::Vec3 upr=(_u+_v)*(_size/2.0);
osg::Vec3 dwr=(_u-_v)*(_size/2.0);
osg::Vec3 ulf=(-_u+_v)*(_size/2.0);
osg::Vec3 dlf=(-_u-_v)*(_size/2.0);
//OSG IS CCW
_pointOne=_origin+upr;
_pointTwo=_origin+ulf;
_pointThree=_origin+dlf;
_pointFour=_origin+dwr;


/*OSG_INFO<<"Distance to plane p1:"<<pointPlaneDistance(_pointOne)<<std::endl;
OSG_INFO<<"Distance to plane p2:"<<pointPlaneDistance(_pointTwo)<<std::endl;
OSG_INFO<<"Distance to plane p3:"<<pointPlaneDistance(_pointThree)<<std::endl;
OSG_INFO<<"Distance to plane p4:"<<pointPlaneDistance(_pointFour)<<std::endl;*/
 
OSG_INFO<<std::endl;
// Now translate this beast to the origin to make it visible:
_pointOne-=_origin;
_pointTwo-=_origin;
_pointThree-=_origin;
_pointFour-=_origin;

OSG_INFO<<"p1:["<<_pointOne[0]<<","<<_pointOne[1]<<","<<_pointOne[2]<<"];"<<std::endl;
OSG_INFO<<"p2:["<<_pointTwo[0]<<","<<_pointTwo[1]<<","<<_pointTwo[2]<<"];"<<std::endl;
OSG_INFO<<"p3:["<<_pointThree[0]<<","<<_pointThree[1]<<","<<_pointThree[2]<<"];"<<std::endl;
OSG_INFO<<"p4:["<<_pointFour[0]<<","<<_pointFour[1]<<","<<_pointFour[2]<<"];"<<std::endl;

OSG_INFO<<"plane1: "<<_plane.getNormal()[0]<<"*x+"<<_plane.getNormal()[1]<<"*y+"<<
_plane.getNormal()[2]<<"*z+"<<_plane.asVec4()[3]<<";"<<std::endl;

OSG_INFO<<"draw3d(enhanced3d=true,implicit(plane1,x,-100,100,y,-100,100,z,-30,30));"<<
std::endl;

//OSG_INFO<<"draw3d(color = red, point_size = 2, point_type = 3,points([p1,p2,p3,p4]));"<<std::endl;

OSG_INFO<<"thequad : mesh([p1,p2], [p4,p3]);"<<std::endl;

OSG_INFO<<"draw3d(line_width = 3,color = red,thequad);"<<std::endl;
}

void calculationsTestTwo()
{
// check that DirectorVectorOne and Perpendicular are
// perpendicular:
double dot=_directorVectorOne*_directorVectorPerpendicular;
std::cout<<dot<<std::endl;
// check that DirectorVectorOne,DirectorVectorTwo and Perpendicular are
// coplanar:
double det=OSGGeneralUtils::determinant(_directorVectorOne,_directorVectorTwo,_directorVectorPerpendicular);
std::cout<<det<<std::endl;
}

double pointPlaneDistance(osg::Vec3 p)
{
// used for testing
/*The distance from P(x,y,z) to the plane ax+by+cz=d is:

|ax+by+cz-d| / sqrt(a^2 + b^2 + c^2)

Therefore |ax+by+cz-d| = |3(2) + (-3) + (-2) - 15| = |-14| = 14.

sqrt(a^2 + b^2 + c^2) = sqrt (9 + 1 + 4) = sqrt14.

The distance is 14/sqrt14 or sqrt 14.*/
double A,B,C,D;
osg::Vec4 pvector= _plane.asVec4();
A=pvector[0];
B=pvector[1];
C=pvector[2];
D=pvector[3];

double numerator=A*p[0]+B*p[1]+C*p[2]+D;
double denom=sqrt(A*A+B*B+C*C);
double dist=numerator/denom;
return dist;
}

osg::Vec3 _pointOne,_pointTwo,_pointThree,_pointFour;
osg::Vec3 _origin;
osg::Plane _plane;
double _size;
osg::Vec3 _u,_v;
osg::Vec3 _directorVectorOne,_directorVectorTwo,_directorVectorThree,_directorVectorPerpendicular;
};