VCRA-SZCZERBACKI-PointBasedRendering/analytic.inl

void analyticPipeline( const std::vector<Splat> &splats, ABuffer *abuffer, const Matrix4x4 &objectToCamera, const Matrix4x4 &projection, const Vector4 &viewport, float front, float back )
{
        SS_ASSERT(abuffer);
        SS_ASSERT(subpixels == 1);      //the analytic pipeline supports only 1 sample per pixel

        const int                       exptablesize = 256;
        static float            exptable[exptablesize];
        static bool                     exptableinitialized = false;
        const float                     exptablescale = float(exptablesize)/F;
        const float                     oosqrtF = 1.0f / sqrt(F);

        if( !exptableinitialized )
        {
                for(int i=0;i<exptablesize;++i)
                {
                        exptable[i] = exp(-0.5f * F * float(i)/float(exptablesize));
                }
        }

        int width = abuffer->getWidth();
        int     height = abuffer->getHeight();
        int     subpixels = abuffer->getSubpixels();

        int scissorxmin = intFloor(viewport.x);
        int scissorymin = intFloor(viewport.y);
        int scissorxmax = intFloor(viewport.z);
        int scissorymax = intFloor(viewport.w);
        if( scissorxmin < 0 ) scissorxmin = 0;
        if( scissorymin < 0 ) scissorymin = 0;
        if( scissorxmax > width ) scissorxmax = width;
        if( scissorymax > height ) scissorymax = height;

        Vector2 vpScale( 0.5f * (viewport.z - viewport.x), 0.5f * (viewport.w - viewport.y) );
        Vector2 vpMid( 0.5f * (viewport.z + viewport.x), 0.5f * (viewport.w + viewport.y) );

        Vector2 oovpScale(1.0f/vpScale.x, 1.0f/vpScale.y);

        //prefilter
        Matrix2x2 Vp( prefiltersize*oovpScale.x*oovpScale.x, 0,
                                  0, prefiltersize*oovpScale.y*oovpScale.y );

        Matrix4x4 objectToScreen;
        objectToScreen = objectToCamera;
        objectToScreen *= projection;

        Matrix4x4 screenToObject;
        screenToObject = objectToScreen;
        bool invertible = screenToObject.invert();
        SS_ASSERT(invertible);
        if( screenToObject[2].w != 0.0f ) screenToObject[2] *= 1.0f/screenToObject[2].w;
        else screenToObject[2] = -screenToObject[2];
        Vector4 c = screenToObject[2];

        Vector4 light(10,10,-20,0);
        Matrix4x4 cameraToObject = objectToCamera;
        cameraToObject.invert();
        light *= cameraToObject;
        Vector3 lightInObject(light.x, light.y, light.z);
        lightInObject.normalize();

        //clip plane setup
        Vector4 clipOrigin[6],clipNormal[6];
        int clipPlanes = 6;
        {
                float s = (float)sqrt(2.0)/2.0f;
                clipOrigin[0].set(0,0,0,0);
                clipNormal[0].set(0,0,1,0);

                clipOrigin[1].set(0,0,back,0);
                clipNormal[1].set(0,0,-1,0);

                clipOrigin[2].set(0,0,0,0);
                clipNormal[2].set(-s,0,0,s);

                clipOrigin[3].set(0,0,0,0);
                clipNormal[3].set(s,0,0,s);

                clipOrigin[4].set(0,0,0,0);
                clipNormal[4].set(0,-s,0,s);

                clipOrigin[5].set(0,0,0,0);
                clipNormal[5].set(0,s,0,s);
        }

        int subpixelsTested = 0;
        int subpixelsPassed = 0;

        for(int i=0;i<splats.size();++i)
        {
                Vector3 p0 = splats[i].m_position;
                Vector3 n = splats[i].m_normal;

                //backface cull (in object-space)
                Vector4 pleq(n.x, n.y, n.z, -dot(p0,n));
                int dotSign = extractSign( dot(c, pleq) );
                if( dotSign ) continue;

                //set up the conic matrix
                Matrix2x2 Vr( splats[i].m_size, 0,
                                          0,                splats[i].m_size );
                Matrix2x2 VrT = Vr;
                VrT.transpose();
                Matrix2x2 Q = Vr * VrT;
                Q.invert();

                //make the tangent plane (in object-space)
                Vector3 tu = perpendicular( n );
                tu.normalize();
                Vector3 tv = cross(n, tu);

                //transform the tangents and
                // the splat centre to the screen-space
                Vector4 tuh(tu.x, tu.y, tu.z, 0);
                Vector4 tvh(tv.x, tv.y, tv.z, 0);
                Vector4 p0h(p0.x, p0.y, p0.z, 1);
                tuh *= objectToScreen;
                tvh *= objectToScreen;
                p0h *= objectToScreen;

                //viewport culling
                int cp;
                for(cp=0;cp<clipPlanes;cp++) {
                        float a = dot(tuh, clipNormal[cp]);
                        float b = dot(tvh, clipNormal[cp]);
                        float c = dot(p0h, clipNormal[cp]) - dot(clipOrigin[cp], clipNormal[cp]);
                        if( a == 0.0f && b == 0.0f ) { //the planes are parallel
                                if( c < 0.0f ) break; //behind the clip plane
                        } else {
                                if( F * (Q[0][0]*b*b - 2.0f*Q[0][1]*a*b + Q[1][1]*a*a) - c*c*Q.det() >= 0.0f )
                                { //the clip plane intersects the splat
                                        if( cp == 0 ) break; //front plane discards the splat
                                } else { //the clip plane doesn't intersect the splat
                                        if( c < 0.0f ) break; //behind the clip plane
                                }
                        }
                }
                if( cp < clipPlanes ) continue; //the splat was behind some of the planes

                //mapping from the splat to the screen-space
                Matrix3x3 M(tuh.x, tuh.y, tuh.w,
                                    tvh.x, tvh.y, tvh.w,
                                    p0h.x, p0h.y, p0h.w);

                //mapping from the screen-space to the splat
                if( M.det() == 0.0f ) continue; //don't draw if the matrix is singular
                Matrix3x3 Mi = M;
                Mi.invert();

                //discard the splat if the mapping matrix's condition number is too high
                float conditionNumber = norm_inf(M) * norm_inf(Mi);
                if( conditionNumber > conditionThreshold ) continue;

                //affine approximation
                Matrix3x3 Qh( Q[0][0], Q[0][1], 0,
                                          Q[1][0], Q[1][1], 0,
                                          0,       0,       -F );
                Matrix3x3 MiT = Mi;
                MiT.transpose();
                Matrix3x3 Qhs = Mi * Qh * MiT;

                float delta = Qhs[0][0] * Qhs[1][1] - Qhs[0][1] * Qhs[1][0];
                if( delta <= 0.0f ) continue; //draw only ellipses
                float oodelta = 1.0f/delta;

                float xt = (Qhs[0][1] * Qhs[1][2] - Qhs[1][1] * Qhs[0][2]) * oodelta;
                float yt = (Qhs[0][1] * Qhs[0][2] - Qhs[0][0] * Qhs[1][2]) * oodelta;

                Matrix2x2 Qs( Qhs[0][0], Qhs[0][1],
                                          Qhs[1][0], Qhs[1][1] );

                float f = -Qhs[2][2] - Qhs[0][2] * xt - Qhs[1][2] * yt;
                if( f <= 0.0f ) continue;
                Qs *= F / f;

                //make the resampling filter, the code below is equivalent to
                float detQs = Qs.det();
                if( detQs <= 0.0f ) continue;
                Qs[0][0] = (Qs[0][0] + detQs*Vp[1][1]);
                Qs[0][1] = (Qs[0][1] - detQs*Vp[0][1]);
                Qs[1][0] = (Qs[1][0] - detQs*Vp[1][0]);
                Qs[1][1] = (Qs[1][1] + detQs*Vp[0][0]);
                Qs *= detQs/Qs.det();

                //the bounding rectangle
                delta = Qs.det();
                if( delta <= 0.0f ) continue;
                oodelta = 1.0f / delta;
                float bx = sqrt(Qs[1][1] * F * oodelta);
                float by = sqrt(Qs[0][0] * F * oodelta);
                float xmin = xt - bx, xmax = xt + bx;
                float ymin = yt - by, ymax = yt + by;

                //determine the integer bounding rectangle and do scissoring
                xmin = xmin * vpScale.x + vpMid.x;
                xmax = xmax * vpScale.x + vpMid.x;
                ymin = ymin * vpScale.y + vpMid.y;
                ymax = ymax * vpScale.y + vpMid.y;

                int sx = intFloor( xmin ), ex = intCeil( xmax );
                if( sx >= width || ex <= 0 ) continue;
                sx = max2(sx,scissorxmin);
                ex = min2(ex,scissorxmax);
                if( ex - sx == 0 ) continue;

                int sy = intFloor( ymin ), ey = intCeil( ymax );
                if( sy >= height || ey <= 0 ) continue;
                sy = max2(sy,scissorymin);
                ey = min2(ey,scissorymax);
                if( ey - sy == 0 ) continue;

                //depth interpolation setup
                float zmin = p0h.z,zmax = p0h.z;
                float zx = 0.0f,zy = 0.0f;
                float zt = p0h.z;
                float l = Q[0][0]*tuh.z*tuh.z + 2.0f*Q[0][1]*tuh.z*tvh.z + Q[1][1]*tvh.z*tvh.z;
                if( l != 0.0f ) {
                        //calculate zmin and zmax
                        l = sqrt(max2(0.0f,F/l));
                        zmin = p0h.z - l * (tuh.z*tuh.z + tvh.z*tvh.z);
                        zmax = p0h.z + l * (tuh.z*tuh.z + tvh.z*tvh.z);
                        zt = (zmin + zmax)*0.5f;

                        //calculate the depth gradient, not needed for our surface reconstruction
                        zx = -tvh.z * M[0][1] + tuh.z * M[1][1];
                        zy =  tvh.z * M[0][0] - tuh.z * M[1][0];
                        float dz = F*(zy*zy*Qs[0][0] - 2*zy*zx*Qs[0][1] + zx*zx*Qs[1][1]);
                        if( dz != 0.0f ) dz = (zmax-zt) * sqrt(max2(0.0f,delta / dz));
                        zx *= dz;
                        zy *= dz;
                }

                //call the splat shader here
                Vector4 color(0,0,0,1);
                {
                        //------------------- template ---------------------------------
                        //diffuse
                        float d = dot(splats[i].m_normal, lightInObject);
                        d = max2(0.5f,d);
                        color = splats[i].m_color;
                        color *= d;

                        //specular
                        Vector3 halfway(cameraToObject[3].x, cameraToObject[3].y, cameraToObject[3].z);
                        halfway -= splats[i].m_position;
                        halfway.normalize();
                        halfway += lightInObject;
                        halfway.normalize();
                        d = dot(halfway, splats[i].m_normal);
                        if( d > 0.0f )
                        {
                                d *= d;
                                d *= d;
                                d *= d;
                                d *= d;
                                d *= d;
                                color.x += d;
                                color.y += d;
                                color.z += d;
                        }
                        color.w = 1.0f;
                        //------------------- template ---------------------------------
                }

                if( rasterization )
                {
                        //rasterization loop
                        ABuffer::Fragment fragment;
                        fragment.zmin = zmin;
                        fragment.zmax = zmax;
                        fragment.dotSign = dotSign;
                        fragment.color = color;
                        for(int y=sy;y<ey;++y)
                        {
                                for(int x=sx;x<ex;++x)
                                {
                                        Vector2 xp(((float)x - vpMid.x)*oovpScale.x - xt,
                                                                 ((float)y - vpMid.y)*oovpScale.y - yt);
                                        subpixelsTested++;
                                        float radius = Qs[0][0]*xp.x*xp.x + 2.0f*Qs[0][1]*xp.x*xp.y + Qs[1][1]*xp.y*xp.y;
                                        if( radius < F )
                                        {
                                                subpixelsPassed++;
                                                float z = zt + xp.x * zx + xp.y * zy; //needed, if the surface reconstruction uses z
                                                fragment.weight = exptable[max2(0,(int)intFloor(radius*exptablescale))];
                                                fragment.z = z;
                                                abuffer->appendFragment(x,y,0,fragment);
                                        }
                                }
                        }
                }

                if( pointvisualization )
                {
                        float w0 = 1.0f / p0h.w;
                        float midx = p0h.x * w0;
                        float midy = p0h.y * w0;
                        glPointSize(3.0f);
                        glBegin(GL_POINTS);
                        glVertex3f(midx*vpScale.x*figsx, midy*vpScale.y*figsy, -1);
                        glEnd();
                }

                if( reconstructionvisualization )
                {
/*                      //draw bounding box
                        glBegin(GL_LINE_LOOP);
                        glVertex3f((xmin-vpMid.x)*figsx, (ymax-vpMid.y)*figsy, -1);
                        glVertex3f((xmax-vpMid.x)*figsx, (ymax-vpMid.y)*figsy, -1);
                        glVertex3f((xmax-vpMid.x)*figsx, (ymin-vpMid.y)*figsy, -1);
                        glVertex3f((xmin-vpMid.x)*figsx, (ymin-vpMid.y)*figsy, -1);
                        glEnd();
*/
                        //draw splat boundary
                        float w0 = 1.0f / p0h.w;
                        float midx = p0h.x * w0;
                        float midy = p0h.y * w0;

                        int segments = 25;
                        float a = 0.0f;
                        float aa = 2.0f * 3.1415f / float(segments);
                        Matrix2x2 J( Mi[0][0], Mi[0][1],
                                                 Mi[1][0], Mi[1][1] );
                        Matrix2x2 Jt = J;
                        Jt.transpose();
                        J *= Q;
                        J *= Jt;
                        float sk = sqrt(F);
                        glBegin(GL_LINE_LOOP);
                        for(int i=0;i<segments;++i)
                        {
                                Vector2 uv;
                                float f,g,t;

                                //perspective correct position
                                // a*x^2+2*b*x*y+c*y^2 = k * (oow0 + dwx*x + dwy*y)^2
                                // sqrt(a*x^2+2*b*x*y+c*y^2) = sqrt(k) * (oow0 + dwx*x + dwy*y)
                                // y = f/g*x =>
                                // sqrt(a*x^2+2*b*x*f/g*x+c*f/g*x*f/g*x) = sqrt(k) * (oow0 + dwx*x + dwy*f/g*x)
                                // sqrt(a*x^2*g^2+2*b*x^2*f*g+c*f^2*x^2) = sqrt(k) * (oow0*g + dwx*x*g + dwy*f*x)
                                // sqrt(a*g^2+2*b*f*g+c*f^2) * x - sqrt(k)*dwx*g*x - sqrt(k)*dwy*f*x = sqrt(k)*oow0*g
                                // (sqrt(a*g^2+2*b*f*g+c*f^2) - sqrt(k)*dwx*g - sqrt(k)*dwy*f) * x = sqrt(k)*oow0*g
                                // x = sqrt(k)*oow0*g / (sqrt(a*g^2+2*b*f*g+c*f^2) - sqrt(k)*dwx*g - sqrt(k)*dwy*f)
                                // y = sqrt(k)*oow0*f / (sqrt(a*g^2+2*b*f*g+c*f^2) - sqrt(k)*dwx*g - sqrt(k)*dwy*f)

                                f = sin(a);
                                g = cos(a);
                                t = sk * w0 / (sqrt(J[0][0] * g * g + 2.0f * J[0][1] * f * g + J[1][1] * f * f) - sk * (g * Mi[0][2] + f * Mi[1][2]));
                                glVertex3f(((midx + g * t)*vpScale.x)*figsx, ((midy + f * t)*vpScale.y)*figsy, -1);

                                a += aa;
                        }
                        glEnd();
                }
        }
}

---