source: GTP/trunk/App/Demos/Vis/FriendlyCulling/src/ShadowMapping.cpp @ 2942

#include "ShadowMapping.h"
#include "FrameBufferObject.h"
#include "RenderState.h"
#include "RenderTraverser.h"
#include "Light.h"
#include "Polygon3.h"
#include "Polyhedron.h"

#include <IL/il.h>
#include <assert.h>


using namespace std;


namespace CHCDemoEngine
{

static CGprogram sCgShadowProgram;
static CGparameter sShadowParam;


static Polyhedron *polyhedron = NULL;
static Polyhedron *lightPoly = NULL;
static Vector3 dummyPt;
static Matrix4x4 dummyMat;


static void PrintGLerror(char *msg)
{
        GLenum errCode;
        const GLubyte *errStr;
        
        if ((errCode = glGetError()) != GL_NO_ERROR) 
        {
                errStr = gluErrorString(errCode);
                fprintf(stderr,"OpenGL ERROR: %s: %s\n", errStr, msg);
        }
}


static Polyhedron *CreatePolyhedron(const Matrix4x4 &lightMatrix, 
                                                                        const AxisAlignedBox3 &sceneBox)
{
        Frustum frustum(lightMatrix);

        vector<Plane3> clipPlanes;

        for (int i = 0; i < 6; ++ i)
        {
                ////////////
                //-- normalize the coefficients

                // the clipping planes look outward the frustum,
                // so distances > 0 mean that a point is outside
                const float invLength = -1.0f / Magnitude(frustum.mClipPlanes[i].mNormal);

                frustum.mClipPlanes[i].mD *= invLength;
                frustum.mClipPlanes[i].mNormal *= invLength;
        }

        // first create near plane because of precision issues
        clipPlanes.push_back(frustum.mClipPlanes[4]);

        clipPlanes.push_back(frustum.mClipPlanes[0]);
        clipPlanes.push_back(frustum.mClipPlanes[1]);
        clipPlanes.push_back(frustum.mClipPlanes[2]);
        clipPlanes.push_back(frustum.mClipPlanes[3]);
        clipPlanes.push_back(frustum.mClipPlanes[5]);

        return Polyhedron::CreatePolyhedron(clipPlanes, sceneBox);
}


static void GrabDepthBuffer(float *data, GLuint depthTexture)
{
        glEnable(GL_TEXTURE_2D);
        glBindTexture(GL_TEXTURE_2D, depthTexture);

        glGetTexImage(GL_TEXTURE_2D, 0, GL_DEPTH_COMPONENT, GL_FLOAT, data);

        glBindTexture(GL_TEXTURE_2D, 0);
        glDisable(GL_TEXTURE_2D);
}


static void ExportDepthBuffer(float *data, int size)
{
        ilInit();
        assert(ilGetError() == IL_NO_ERROR);

        ILstring filename = ILstring("shadow.tga");
        ilRegisterType(IL_FLOAT);

        const int depth = 1;
        const int bpp = 1;

        if (!ilTexImage(size, size, depth, bpp, IL_LUMINANCE, IL_FLOAT, data))
        {
                cerr << "IL error " << ilGetError() << endl;
        
                ilShutDown();
                assert(ilGetError() == IL_NO_ERROR);

                return;
        }

        if (!ilSaveImage(filename))
        {
                cerr << "TGA write error " << ilGetError() << endl;
        }

        ilShutDown();
        assert(ilGetError() == IL_NO_ERROR);

        cout << "exported depth buffer" << endl;
}



static AxisAlignedBox3 GetExtremalPoints(const Matrix4x4 &m,
                                                                                 const VertexArray &pts)
{
        AxisAlignedBox3 extremalPoints;
        extremalPoints.Initialize();

        VertexArray::const_iterator it, it_end = pts.end();
                
        for (it = pts.begin(); it != it_end; ++ it)
        {
                Vector3 pt = *it;
                pt = m * pt;

                extremalPoints.Include(pt);
        }

        return extremalPoints;
}


ShadowMap::ShadowMap(Light *light, int size, const AxisAlignedBox3 &sceneBox, Camera *cam):
mSceneBox(sceneBox), mSize(size), mCamera(cam), mLight(light)
{
        mFbo = new FrameBufferObject(size, size, FrameBufferObject::DEPTH_32, true);
        // the diffuse color buffer
        mFbo->AddColorBuffer(ColorBufferObject::BUFFER_UBYTE, ColorBufferObject::WRAP_CLAMP_TO_EDGE, ColorBufferObject::FILTER_LINEAR, false);
        //mFbo->AddColorBuffer(ColorBufferObject::BUFFER_FLOAT_32, ColorBufferObject::WRAP_CLAMP_TO_EDGE, ColorBufferObject::FILTER_LINEAR, false);

        mShadowCam = new Camera(mSize, mSize);//mSceneBox.Size().x * 0.5f, mSceneBox.Size().y * 0.5f);
        mShadowCam->SetOrtho(true);
}


ShadowMap::~ShadowMap()
{
        DEL_PTR(mFbo);
        DEL_PTR(mShadowCam);
}


void ShadowMap::DrawPoly(Polyhedron *poly, const Vector3 &color, const Vector3 &color2)
{
        if (!poly) return;

        for (size_t i = 0; i < poly->NumPolygons(); ++ i)
        {
                glColor3f(color.x, color.y, color.z);

                if (i == poly->NumPolygons() - 1) // hack: different color for near
                {
                        //glLineWidth(2);
                        glColor3f(color2.x, color2.y, color2.z);
                
                }
                glBegin(GL_LINE_LOOP);

                Polygon3 *p = poly->GetPolygons()[i];

                for (size_t j = 0; j < p->mVertices.size(); ++ j)
                {
                        Vector3 v = p->mVertices[j];
                        glVertex3d(v.x, v.y, v.z);
                }

                glEnd();
                
                //if (i == poly->NumPolygons() - 1) glLineWidth(1);
        }
}


void ShadowMap::DrawPolys()
{
        DrawPoly(lightPoly, Vector3(1, 0, 1), Vector3(1, 1, 1));
        DrawPoly(polyhedron, Vector3(0, 1, 0), Vector3(0, 1, 1));

        glPointSize(10.0f);

        Vector3 pt = Vector3::ZERO();
        //Vector3 pt = dummyPt;

        Matrix4x4 myMat = Invert(dummyMat);
        pt = myMat * pt;

        glBegin(GL_POINTS);
        glVertex3f(pt.x, pt.y, pt.z);
        glEnd();
}


float ShadowMap::ComputeN(const AxisAlignedBox3 &extremalPoints) const 
{
        const float n = mCamera->GetNear();
        
        const float d = fabs(extremalPoints.Max()[2] - extremalPoints.Min()[2]); 
        
        const float dotProd = DotProd(mCamera->GetDirection(), mShadowCam->GetDirection());
        const float sinGamma = sin(fabs(acos(dotProd)));

        return (n + sqrt(n * (n + d * sinGamma))) /  sinGamma;
}


Matrix4x4 ShadowMap::CalcLispSMTransform(const Matrix4x4 &lightSpace, 
                                                                                 const AxisAlignedBox3 &extremalPoints,
                                                                                 const VertexArray &body
                                                                                 )
{
        //return IdentityMatrix();

        AxisAlignedBox3 bounds_ls = GetExtremalPoints(lightSpace, body);

        ///////////////
        //-- We apply the lispsm algorithm in order to calculate an optimal light projection matrix
        //-- first find the free parameter values n, and P (the projection center), and the projection depth

        //const float n = 1e2f;
        //const float n = 1e6f;
        const float n = ComputeN(bounds_ls);

        cout << "n: " << n << endl;

        const Vector3 nearPt = GetNearCameraPointE(body);
        
        //get the coordinates of the near camera point in light space
        const Vector3 lsNear = lightSpace * nearPt;

        // the start point has the x and y coordinate of e, the z coord of the near plane of the light volume
        const Vector3 startPt = Vector3(lsNear.x, lsNear.y, bounds_ls.Max().z);
        //const Vector3 schas = lightSpace * mCamera->GetPosition();
        //const Vector3 startPt = Vector3(schas.x, schas.y, bounds_ls.Max().z);

        cout << "mx: " <<  bounds_ls.Max() << endl;
        cout << "mn: " << bounds_ls.Min() << endl;

        // the new projection center
        const Vector3 projCenter = startPt + Vector3::UNIT_Z() * n;

        cout <<"start: " << startPt << " " << projCenter << " " << Distance(lightSpace * mCamera->GetPosition(), startPt) << endl;

        dummyPt = startPt;

        //construct a translation that moves to the projection center
        const Matrix4x4 projectionCenter = TranslationMatrix(-projCenter);

        // light space y size
        const float d = fabs(bounds_ls.Max()[2] - bounds_ls.Min()[2]);

        const float dy = fabs(bounds_ls.Max()[1] - bounds_ls.Min()[1]);
        const float dx = fabs(bounds_ls.Max()[0] - bounds_ls.Min()[0]);

        cout << "d: " << d << " dy: " << dy << " dx: " << dx << endl;

        

        //////////
        //-- now apply these values to construct the perspective lispsm matrix

        Matrix4x4 matLispSM;
        
        matLispSM = GetFrustum(-1.0, 1.0, -1.0, 1.0, n, n + d);

        //cout << "lispsm\n" << matLispSM << endl;

        // translate to the projection center
        matLispSM = projectionCenter * matLispSM;

        //cout << "new\n" << matLispSM << endl;

        // transform into OpenGL right handed system
        Matrix4x4 refl = ScaleMatrix(1.0f, 1.0f, -1.0f);
        matLispSM *= refl;
        
        return matLispSM;
}

#if 0
Vector3 ShadowMap::GetNearCameraPointE(const VertexArray &pts) const
{
        float maxDist = -1e25f;
        Vector3 nearest = Vector3::ZERO();

        Matrix4x4 eyeView;
        mCamera->GetModelViewMatrix(eyeView);

        VertexArray newPts;
        polyhedron->CollectVertices(newPts);
        
        //the LVS volume is always in front of the camera
        VertexArray::const_iterator it, it_end = pts.end();     

        for (it = pts.begin(); it != it_end; ++ it)
        {
                Vector3 pt = *it;
                Vector3 ptE = eyeView * pt;
//cout<<"i"<< pt.z;
                if (ptE.z > 0) cerr <<"should not happen " << ptE.z << endl;
                else
                if (ptE.z > maxDist)
                {
                        cout << " d " << ptE.z;
        
                        maxDist = ptE.z;
                        nearest = pt;
                }
        }

        //      return Invert(eyeView) * nearest;
        return nearest;
}

#else

Vector3 ShadowMap::GetNearCameraPointE(const VertexArray &pts) const
{
        VertexArray newPts;
        polyhedron->CollectVertices(newPts);

        Vector3 nearest = Vector3::ZERO();
        float minDist = 1e25f;

        const Vector3 camPos = mCamera->GetPosition();

        VertexArray::const_iterator it, it_end = newPts.end();

        for (it = newPts.begin(); it != it_end; ++ it)
        {
                Vector3 pt = *it;

                const float dist = SqrDistance(pt, camPos);

                if (dist < minDist)
                {
                        minDist = dist;
                        nearest = pt;
                }
        }

        return nearest;
}

#endif

Vector3 ShadowMap::GetProjViewDir(const Matrix4x4 &lightSpace, 
                                                                  const VertexArray &pts) const 
{
        //get the point in the LVS volume that is nearest to the camera
        const Vector3 e = GetNearCameraPointE(pts);

        //construct edge to transform into light-space
        const Vector3 b = e + mCamera->GetDirection();
        //transform to light-space
        const Vector3 e_lp = lightSpace * e;
        const Vector3 b_lp = lightSpace * b;

        Vector3 projDir(b_lp - e_lp);

        Matrix4x4 dummy = Invert(lightSpace);
        Invert(dummy);
        Vector3 dummyVec = dummy * e_lp;
        Vector3 dummyVec2 = dummy * b_lp;

        //projDir.z = -projDir.z;

        cout << "dummy: " << Normalize(dummyVec2 - dummyVec) << endl;
        //project the view direction into the shadow map plane
        projDir.y = .0f;

        return Normalize(projDir);
        //return projDir;
}


bool ShadowMap::CalcLightProjection(Matrix4x4 &lightProj)
{
        ///////////////////
        //-- First step: calc frustum clipped by scene box

        DEL_PTR(polyhedron);
        polyhedron = CalcClippedFrustum(mSceneBox);

        if (!polyhedron) return false; // something is wrong

        // include the part of the light volume that "sees" the frustum
        // we only require frustum vertices

        VertexArray frustumPoints;
        //Vector3 lightDir = Vector3(0, 0, -1);
        IncludeLightVolume(*polyhedron, frustumPoints, mLight->GetDirection(), mSceneBox);


        ///////////////
        //-- transform points from world view to light view and calculate extremal points

        Matrix4x4 lightView;
        mShadowCam->GetModelViewMatrix(lightView);

        const AxisAlignedBox3 extremalPoints = GetExtremalPoints(lightView, frustumPoints);

        // we use directional lights, so the projection can be set to identity
        lightProj = IdentityMatrix(); 

        // switch coordinate system to that used in the lispsm algorithm for calculations
        Matrix4x4 transform2LispSM = ZeroMatrix();

        transform2LispSM.x[0][0] =  1.0f;
        transform2LispSM.x[1][2] =  -1.0f; // y => -z
        transform2LispSM.x[2][1] =  1.0f; // z => y
        transform2LispSM.x[3][3] =  1.0f;


        //switch to the lightspace used in the article
        lightProj *= transform2LispSM;

        const Vector3 projViewDir = GetProjViewDir(lightView * lightProj, frustumPoints);

        cout << "projViewDir: " << projViewDir << " orig " << mCamera->GetDirection() << endl;
                                                                
        //do Light Space Perspective shadow mapping
        //rotate the lightspace so that the projected light view always points upwards
        //calculate a frame matrix that uses the projViewDir[lightspace] as up vector
        //look(from position, into the direction of the projected direction, with unchanged up-vector)
        //const Matrix4x4 frame = MyLookAt2(Vector3::ZERO(), projViewDir, Vector3::UNIT_Y());
        const Matrix4x4 frame = MyLookAt(Vector3::ZERO(), projViewDir, Vector3::UNIT_Y());

        cout << "frame\n " << frame << endl;
        lightProj *= frame;

        const Matrix4x4 matLispSM = 
                CalcLispSMTransform(lightView * lightProj, extremalPoints, frustumPoints);

        lightProj *= matLispSM;

        // change back to GL coordinate system
        Matrix4x4 transformToGL = ZeroMatrix();
        
        transformToGL.x[0][0] =   1.0f;
        transformToGL.x[1][2] =   1.0f; // z => y
        transformToGL.x[2][1] =  -1.0f; // y => -z
        transformToGL.x[3][3] =   1.0f;

        lightProj *= transformToGL;

        AxisAlignedBox3 lightPts = GetExtremalPoints(lightView * lightProj, frustumPoints);

        cout << "max: " << lightPts.Max() << " min: " << lightPts.Min() << endl;
        // focus projection matrix on the extremal points => scale to unit cube
        Matrix4x4 scaleTranslate = GetFittingProjectionMatrix(lightPts);

        cout << "scaleTranslate:\n" << scaleTranslate << endl;

        lightProj = lightProj * scaleTranslate;

        Matrix4x4 mymat = lightView * lightProj;

        AxisAlignedBox3 lightPtsNew = GetExtremalPoints(mymat, frustumPoints);

        cout << "newmax: " << lightPtsNew.Max() << endl;
        cout << "newmin: " << lightPtsNew.Min() << endl;

        // we have to flip the signs in order to tranform to opengl right handed system
        Matrix4x4 refl = ScaleMatrix(1, 1, -1); 
        lightProj *= refl;
        
        return true;
}


Polyhedron *ShadowMap::CalcClippedFrustum(const AxisAlignedBox3 &box) const
{
        Polyhedron *p = mCamera->ComputeFrustum();
        
        Polyhedron *clippedPolyhedron = box.CalcIntersection(*p);

        DEL_PTR(p);
        
        return clippedPolyhedron;
}


//calculates the up vector for the light coordinate frame
static Vector3 CalcUpVec(const Vector3 viewDir, const Vector3 lightDir) 
{
        //we do what gluLookAt does...
        //left is the normalized vector perpendicular to lightDir and viewDir
        //this means left is the normalvector of the yz-plane from the paper
        Vector3 left = CrossProd(lightDir, viewDir);
        
        //we now can calculate the rotated(in the yz-plane) viewDir vector
        //and use it as up vector in further transformations
        Vector3 up = CrossProd(left, lightDir);

        return Normalize(up);
}


void ShadowMap::GetTextureMatrix(Matrix4x4 &m) const
{
        m = mTextureMatrix;
}

 
unsigned int ShadowMap::GetDepthTexture() const
{
        return mFbo->GetDepthTex();
}

unsigned int ShadowMap::GetShadowColorTexture() const
{
        return mFbo->GetColorBuffer(0)->GetTexture();
        
}


void ShadowMap::IncludeLightVolume(const Polyhedron &polyhedron, 
                                                                   VertexArray &frustumPoints, 
                                                                   const Vector3 lightDir,
                                                                   const AxisAlignedBox3 &sceneBox
                                                                   ) 
{
        // we don't need closed form anymore => just store vertices
        VertexArray vertices;
        polyhedron.CollectVertices(vertices);

        // we 'look' at each point and calculate intersections of rays with scene bounding box
        VertexArray::const_iterator it, it_end = vertices.end();

        for (it = vertices.begin(); it != it_end; ++ it)
        {
                Vector3 v  = *it;

                frustumPoints.push_back(v);
                
                // hack: start at point which is guaranteed to be outside of box
                v -= Magnitude(mSceneBox.Diagonal()) * lightDir;

                SimpleRay ray(v, lightDir);

                float tNear, tFar;

                if (sceneBox.Intersects(ray, tNear, tFar))
                {
                        Vector3 newpt = ray.Extrap(tNear);
                        frustumPoints.push_back(newpt);                 
                }
        }
}


void ShadowMap::ComputeShadowMap(RenderTraverser *renderer, const Matrix4x4 &projView)
{
        mFbo->Bind();
        
        glDrawBuffers(1, mrt);

        glPushAttrib(GL_VIEWPORT_BIT);
        glViewport(0, 0, mSize, mSize);

        glDisable(GL_LIGHTING);
        glDisable(GL_TEXTURE_2D);
        glColorMask(GL_FALSE, GL_FALSE, GL_FALSE, GL_FALSE);

        glShadeModel(GL_FLAT);


        /////////////
        //-- render scene into shadow map

        _Render(renderer);


        //////////////
        //-- compute texture matrix

        static Matrix4x4 biasMatrix(0.5f, 0.0f, 0.0f, 0.5f,
                                                                0.0f, 0.5f, 0.0f, 0.5f,
                                                                0.0f, 0.0f, 0.5f, 0.5f,
                                                                0.0f, 0.0f, 0.0f, 1.0f);

        mTextureMatrix = mLightProjView * biasMatrix; 

        glPopAttrib();
        
        glEnable(GL_LIGHTING);
        glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);

#if 0
        float *data = new float[mSize * mSize];

        GrabDepthBuffer(data, mFbo->GetDepthTex());
        ExportDepthBuffer(data, mSize);

        delete [] data;
        
        PrintGLerror("shadow map");
#endif
        
        FrameBufferObject::Release();
}


void ShadowMap::RenderShadowView(RenderTraverser *renderer, const Matrix4x4 &projView)
{
        glEnable(GL_LIGHTING);
        glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);

        _Render(renderer);
        
        glDisable(GL_POLYGON_OFFSET_FILL);
        glDisable(GL_LIGHTING);
}


void ShadowMap::_Render(RenderTraverser *renderer)
{
        const Vector3 dir = mLight->GetDirection();
        //const Vector3 dir(0, 0, -1);

        cout << "dir: " << dir << endl;
        mShadowCam->SetDirection(dir);

        // set position so that we can see the whole scene
        Vector3 pos = mSceneBox.Center();
        pos -= dir * Magnitude(mSceneBox.Diagonal() * 0.5f);

        mShadowCam->SetPosition(mCamera->GetPosition());

        Vector3 upVec = CalcUpVec(mCamera->GetDirection(), dir);
        Matrix4x4 lightView = MyLookAt(mShadowCam->GetPosition(), dir, upVec);

        mShadowCam->mViewOrientation = lightView;

        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

        glPolygonOffset(1.0f, 2000.0f);
        glEnable(GL_POLYGON_OFFSET_FILL);
        
        Matrix4x4 lightProj;
        CalcLightProjection(lightProj);

        glMatrixMode(GL_PROJECTION);
        glPushMatrix();
        glLoadMatrixf((float *)lightProj.x);

        mLightProjView = lightView * lightProj;
        dummyMat = mLightProjView;

        DEL_PTR(lightPoly);
        lightPoly = CreatePolyhedron(mLightProjView, mSceneBox);

        glMatrixMode(GL_MODELVIEW);
        glPushMatrix();
        glLoadIdentity();

        //glDisable(GL_CULL_FACE);
        mShadowCam->SetupCameraView();

        
        /////////////
        //-- render scene into shadow map

        renderer->RenderScene();

        glMatrixMode(GL_MODELVIEW);
        glPopMatrix();

        glMatrixMode(GL_PROJECTION);
        glPopMatrix();

        
        glDisable(GL_POLYGON_OFFSET_FILL);
}
} // namespace