- Main Page
- Namespaces
- Classes
- Files

# VCRA-SZCZERBACKI-PointBasedRendering/analytic.inl

    00001 void analyticPipeline( const std::vector<Splat> &splats, ABuffer *abuffer, const Matrix4x4 &objectToCamera, const Matrix4x4 &projection, const Vector4 &viewport, float front, float back )
    00002 {
    00003         SS_ASSERT(abuffer);
    00004         SS_ASSERT(subpixels == 1);      //the analytic pipeline supports only 1 sample per pixel
    00005
    00006         const int                       exptablesize = 256;
    00007         static float            exptable[exptablesize];
    00008         static bool                     exptableinitialized = false;
    00009         const float                     exptablescale = float(exptablesize)/F;
    00010         const float                     oosqrtF = 1.0f / sqrt(F);
    00011
    00012         if( !exptableinitialized )
    00013         {
    00014                 for(int i=0;i<exptablesize;++i)
    00015                 {
    00016                         exptable[i] = exp(-0.5f * F * float(i)/float(exptablesize));
    00017                 }
    00018         }
    00019
    00020         int width = abuffer->getWidth();
    00021         int     height = abuffer->getHeight();
    00022         int     subpixels = abuffer->getSubpixels();
    00023
    00024         int scissorxmin = intFloor(viewport.x);
    00025         int scissorymin = intFloor(viewport.y);
    00026         int scissorxmax = intFloor(viewport.z);
    00027         int scissorymax = intFloor(viewport.w);
    00028         if( scissorxmin < 0 ) scissorxmin = 0;
    00029         if( scissorymin < 0 ) scissorymin = 0;
    00030         if( scissorxmax > width ) scissorxmax = width;
    00031         if( scissorymax > height ) scissorymax = height;
    00032
    00033         Vector2 vpScale( 0.5f * (viewport.z - viewport.x), 0.5f * (viewport.w - viewport.y) );
    00034         Vector2 vpMid( 0.5f * (viewport.z + viewport.x), 0.5f * (viewport.w + viewport.y) );
    00035
    00036         Vector2 oovpScale(1.0f/vpScale.x, 1.0f/vpScale.y);
    00037
    00038         //prefilter
    00039         Matrix2x2 Vp( prefiltersize*oovpScale.x*oovpScale.x, 0,
    00040                                   0, prefiltersize*oovpScale.y*oovpScale.y );
    00041
    00042         Matrix4x4 objectToScreen;
    00043         objectToScreen = objectToCamera;
    00044         objectToScreen *= projection;
    00045
    00046         Matrix4x4 screenToObject;
    00047         screenToObject = objectToScreen;
    00048         bool invertible = screenToObject.invert();
    00049         SS_ASSERT(invertible);
    00050         if( screenToObject[2].w != 0.0f ) screenToObject[2] *= 1.0f/screenToObject[2].w;
    00051         else screenToObject[2] = -screenToObject[2];
    00052         Vector4 c = screenToObject[2];
    00053
    00054         Vector4 light(10,10,-20,0);
    00055         Matrix4x4 cameraToObject = objectToCamera;
    00056         cameraToObject.invert();
    00057         light *= cameraToObject;
    00058         Vector3 lightInObject(light.x, light.y, light.z);
    00059         lightInObject.normalize();
    00060
    00061         //clip plane setup
    00062         Vector4 clipOrigin[6],clipNormal[6];
    00063         int clipPlanes = 6;
    00064         {
    00065                 float s = (float)sqrt(2.0)/2.0f;
    00066                 clipOrigin[0].set(0,0,0,0);
    00067                 clipNormal[0].set(0,0,1,0);
    00068
    00069                 clipOrigin[1].set(0,0,back,0);
    00070                 clipNormal[1].set(0,0,-1,0);
    00071
    00072                 clipOrigin[2].set(0,0,0,0);
    00073                 clipNormal[2].set(-s,0,0,s);
    00074
    00075                 clipOrigin[3].set(0,0,0,0);
    00076                 clipNormal[3].set(s,0,0,s);
    00077
    00078                 clipOrigin[4].set(0,0,0,0);
    00079                 clipNormal[4].set(0,-s,0,s);
    00080
    00081                 clipOrigin[5].set(0,0,0,0);
    00082                 clipNormal[5].set(0,s,0,s);
    00083         }
    00084
    00085         int subpixelsTested = 0;
    00086         int subpixelsPassed = 0;
    00087
    00088         for(int i=0;i<splats.size();++i)
    00089         {
    00090                 Vector3 p0 = splats[i].m_position;
    00091                 Vector3 n = splats[i].m_normal;
    00092
    00093                 //backface cull (in object-space)
    00094                 Vector4 pleq(n.x, n.y, n.z, -dot(p0,n));
    00095                 int dotSign = extractSign( dot(c, pleq) );
    00096                 if( dotSign ) continue;
    00097
    00098                 //set up the conic matrix
    00099                 Matrix2x2 Vr( splats[i].m_size, 0,
    00100                                           0,                splats[i].m_size );
    00101                 Matrix2x2 VrT = Vr;
    00102                 VrT.transpose();
    00103                 Matrix2x2 Q = Vr * VrT;
    00104                 Q.invert();
    00105
    00106                 //make the tangent plane (in object-space)
    00107                 Vector3 tu = perpendicular( n );
    00108                 tu.normalize();
    00109                 Vector3 tv = cross(n, tu);
    00110
    00111                 //transform the tangents and
    00112                 // the splat centre to the screen-space
    00113                 Vector4 tuh(tu.x, tu.y, tu.z, 0);
    00114                 Vector4 tvh(tv.x, tv.y, tv.z, 0);
    00115                 Vector4 p0h(p0.x, p0.y, p0.z, 1);
    00116                 tuh *= objectToScreen;
    00117                 tvh *= objectToScreen;
    00118                 p0h *= objectToScreen;
    00119
    00120                 //viewport culling
    00121                 int cp;
    00122                 for(cp=0;cp<clipPlanes;cp++) {
    00123                         float a = dot(tuh, clipNormal[cp]);
    00124                         float b = dot(tvh, clipNormal[cp]);
    00125                         float c = dot(p0h, clipNormal[cp]) - dot(clipOrigin[cp], clipNormal[cp]);
    00126                         if( a == 0.0f && b == 0.0f ) { //the planes are parallel
    00127                                 if( c < 0.0f ) break; //behind the clip plane
    00128                         } else {
    00129                                 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 )
    00130                                 { //the clip plane intersects the splat
    00131                                         if( cp == 0 ) break; //front plane discards the splat
    00132                                 } else { //the clip plane doesn't intersect the splat
    00133                                         if( c < 0.0f ) break; //behind the clip plane
    00134                                 }
    00135                         }
    00136                 }
    00137                 if( cp < clipPlanes ) continue; //the splat was behind some of the planes
    00138
    00139                 //mapping from the splat to the screen-space
    00140                 Matrix3x3 M(tuh.x, tuh.y, tuh.w,
    00141                                     tvh.x, tvh.y, tvh.w,
    00142                                     p0h.x, p0h.y, p0h.w);
    00143
    00144                 //mapping from the screen-space to the splat
    00145                 if( M.det() == 0.0f ) continue; //don't draw if the matrix is singular
    00146                 Matrix3x3 Mi = M;
    00147                 Mi.invert();
    00148
    00149                 //discard the splat if the mapping matrix's condition number is too high
    00150                 float conditionNumber = norm_inf(M) * norm_inf(Mi);
    00151                 if( conditionNumber > conditionThreshold ) continue;
    00152
    00153                 //affine approximation
    00154                 Matrix3x3 Qh( Q[0][0], Q[0][1], 0,
    00155                                           Q[1][0], Q[1][1], 0,
    00156                                           0,       0,       -F );
    00157                 Matrix3x3 MiT = Mi;
    00158                 MiT.transpose();
    00159                 Matrix3x3 Qhs = Mi * Qh * MiT;
    00160
    00161                 float delta = Qhs[0][0] * Qhs[1][1] - Qhs[0][1] * Qhs[1][0];
    00162                 if( delta <= 0.0f ) continue; //draw only ellipses
    00163                 float oodelta = 1.0f/delta;
    00164
    00165                 float xt = (Qhs[0][1] * Qhs[1][2] - Qhs[1][1] * Qhs[0][2]) * oodelta;
    00166                 float yt = (Qhs[0][1] * Qhs[0][2] - Qhs[0][0] * Qhs[1][2]) * oodelta;
    00167
    00168                 Matrix2x2 Qs( Qhs[0][0], Qhs[0][1],
    00169                                           Qhs[1][0], Qhs[1][1] );
    00170
    00171                 float f = -Qhs[2][2] - Qhs[0][2] * xt - Qhs[1][2] * yt;
    00172                 if( f <= 0.0f ) continue;
    00173                 Qs *= F / f;
    00174
    00175                 //make the resampling filter, the code below is equivalent to
    00176                 float detQs = Qs.det();
    00177                 if( detQs <= 0.0f ) continue;
    00178                 Qs[0][0] = (Qs[0][0] + detQs*Vp[1][1]);
    00179                 Qs[0][1] = (Qs[0][1] - detQs*Vp[0][1]);
    00180                 Qs[1][0] = (Qs[1][0] - detQs*Vp[1][0]);
    00181                 Qs[1][1] = (Qs[1][1] + detQs*Vp[0][0]);
    00182                 Qs *= detQs/Qs.det();
    00183
    00184                 //the bounding rectangle
    00185                 delta = Qs.det();
    00186                 if( delta <= 0.0f ) continue;
    00187                 oodelta = 1.0f / delta;
    00188                 float bx = sqrt(Qs[1][1] * F * oodelta);
    00189                 float by = sqrt(Qs[0][0] * F * oodelta);
    00190                 float xmin = xt - bx, xmax = xt + bx;
    00191                 float ymin = yt - by, ymax = yt + by;
    00192
    00193                 //determine the integer bounding rectangle and do scissoring
    00194                 xmin = xmin * vpScale.x + vpMid.x;
    00195                 xmax = xmax * vpScale.x + vpMid.x;
    00196                 ymin = ymin * vpScale.y + vpMid.y;
    00197                 ymax = ymax * vpScale.y + vpMid.y;
    00198
    00199                 int sx = intFloor( xmin ), ex = intCeil( xmax );
    00200                 if( sx >= width || ex <= 0 ) continue;
    00201                 sx = max2(sx,scissorxmin);
    00202                 ex = min2(ex,scissorxmax);
    00203                 if( ex - sx == 0 ) continue;
    00204
    00205                 int sy = intFloor( ymin ), ey = intCeil( ymax );
    00206                 if( sy >= height || ey <= 0 ) continue;
    00207                 sy = max2(sy,scissorymin);
    00208                 ey = min2(ey,scissorymax);
    00209                 if( ey - sy == 0 ) continue;
    00210
    00211                 //depth interpolation setup
    00212                 float zmin = p0h.z,zmax = p0h.z;
    00213                 float zx = 0.0f,zy = 0.0f;
    00214                 float zt = p0h.z;
    00215                 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;
    00216                 if( l != 0.0f ) {
    00217                         //calculate zmin and zmax
    00218                         l = sqrt(max2(0.0f,F/l));
    00219                         zmin = p0h.z - l * (tuh.z*tuh.z + tvh.z*tvh.z);
    00220                         zmax = p0h.z + l * (tuh.z*tuh.z + tvh.z*tvh.z);
    00221                         zt = (zmin + zmax)*0.5f;
    00222
    00223                         //calculate the depth gradient, not needed for our surface reconstruction
    00224                         zx = -tvh.z * M[0][1] + tuh.z * M[1][1];
    00225                         zy =  tvh.z * M[0][0] - tuh.z * M[1][0];
    00226                         float dz = F*(zy*zy*Qs[0][0] - 2*zy*zx*Qs[0][1] + zx*zx*Qs[1][1]);
    00227                         if( dz != 0.0f ) dz = (zmax-zt) * sqrt(max2(0.0f,delta / dz));
    00228                         zx *= dz;
    00229                         zy *= dz;
    00230                 }
    00231
    00232                 //call the splat shader here
    00233                 Vector4 color(0,0,0,1);
    00234                 {
    00235                         //------------------- template ---------------------------------
    00236                         //diffuse
    00237                         float d = dot(splats[i].m_normal, lightInObject);
    00238                         d = max2(0.5f,d);
    00239                         color = splats[i].m_color;
    00240                         color *= d;
    00241
    00242                         //specular
    00243                         Vector3 halfway(cameraToObject[3].x, cameraToObject[3].y, cameraToObject[3].z);
    00244                         halfway -= splats[i].m_position;
    00245                         halfway.normalize();
    00246                         halfway += lightInObject;
    00247                         halfway.normalize();
    00248                         d = dot(halfway, splats[i].m_normal);
    00249                         if( d > 0.0f )
    00250                         {
    00251                                 d *= d;
    00252                                 d *= d;
    00253                                 d *= d;
    00254                                 d *= d;
    00255                                 d *= d;
    00256                                 color.x += d;
    00257                                 color.y += d;
    00258                                 color.z += d;
    00259                         }
    00260                         color.w = 1.0f;
    00261                         //------------------- template ---------------------------------
    00262                 }
    00263
    00264                 if( rasterization )
    00265                 {
    00266                         //rasterization loop
    00267                         ABuffer::Fragment fragment;
    00268                         fragment.zmin = zmin;
    00269                         fragment.zmax = zmax;
    00270                         fragment.dotSign = dotSign;
    00271                         fragment.color = color;
    00272                         for(int y=sy;y<ey;++y)
    00273                         {
    00274                                 for(int x=sx;x<ex;++x)
    00275                                 {
    00276                                         Vector2 xp(((float)x - vpMid.x)*oovpScale.x - xt,
    00277                                                                  ((float)y - vpMid.y)*oovpScale.y - yt);
    00278                                         subpixelsTested++;
    00279                                         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;
    00280                                         if( radius < F )
    00281                                         {
    00282                                                 subpixelsPassed++;
    00283                                                 float z = zt + xp.x * zx + xp.y * zy; //needed, if the surface reconstruction uses z
    00284                                                 fragment.weight = exptable[max2(0,(int)intFloor(radius*exptablescale))];
    00285                                                 fragment.z = z;
    00286                                                 abuffer->appendFragment(x,y,0,fragment);
    00287                                         }
    00288                                 }
    00289                         }
    00290                 }
    00291
    00292                 if( pointvisualization )
    00293                 {
    00294                         float w0 = 1.0f / p0h.w;
    00295                         float midx = p0h.x * w0;
    00296                         float midy = p0h.y * w0;
    00297                         glPointSize(3.0f);
    00298                         glBegin(GL_POINTS);
    00299                         glVertex3f(midx*vpScale.x*figsx, midy*vpScale.y*figsy, -1);
    00300                         glEnd();
    00301                 }
    00302
    00303                 if( reconstructionvisualization )
    00304                 {
    00305 /*                      //draw bounding box
    00306                         glBegin(GL_LINE_LOOP);
    00307                         glVertex3f((xmin-vpMid.x)*figsx, (ymax-vpMid.y)*figsy, -1);
    00308                         glVertex3f((xmax-vpMid.x)*figsx, (ymax-vpMid.y)*figsy, -1);
    00309                         glVertex3f((xmax-vpMid.x)*figsx, (ymin-vpMid.y)*figsy, -1);
    00310                         glVertex3f((xmin-vpMid.x)*figsx, (ymin-vpMid.y)*figsy, -1);
    00311                         glEnd();
    00312 */
    00313                         //draw splat boundary
    00314                         float w0 = 1.0f / p0h.w;
    00315                         float midx = p0h.x * w0;
    00316                         float midy = p0h.y * w0;
    00317
    00318                         int segments = 25;
    00319                         float a = 0.0f;
    00320                         float aa = 2.0f * 3.1415f / float(segments);
    00321                         Matrix2x2 J( Mi[0][0], Mi[0][1],
    00322                                                  Mi[1][0], Mi[1][1] );
    00323                         Matrix2x2 Jt = J;
    00324                         Jt.transpose();
    00325                         J *= Q;
    00326                         J *= Jt;
    00327                         float sk = sqrt(F);
    00328                         glBegin(GL_LINE_LOOP);
    00329                         for(int i=0;i<segments;++i)
    00330                         {
    00331                                 Vector2 uv;
    00332                                 float f,g,t;
    00333
    00334                                 //perspective correct position
    00335                                 // a*x^2+2*b*x*y+c*y^2 = k * (oow0 + dwx*x + dwy*y)^2
    00336                                 // sqrt(a*x^2+2*b*x*y+c*y^2) = sqrt(k) * (oow0 + dwx*x + dwy*y)
    00337                                 // y = f/g*x =>
    00338                                 // 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)
    00339                                 // 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)
    00340                                 // 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
    00341                                 // (sqrt(a*g^2+2*b*f*g+c*f^2) - sqrt(k)*dwx*g - sqrt(k)*dwy*f) * x = sqrt(k)*oow0*g
    00342                                 // x = sqrt(k)*oow0*g / (sqrt(a*g^2+2*b*f*g+c*f^2) - sqrt(k)*dwx*g - sqrt(k)*dwy*f)
    00343                                 // y = sqrt(k)*oow0*f / (sqrt(a*g^2+2*b*f*g+c*f^2) - sqrt(k)*dwx*g - sqrt(k)*dwy*f)
    00344
    00345                                 f = sin(a);
    00346                                 g = cos(a);
    00347                                 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]));
    00348                                 glVertex3f(((midx + g * t)*vpScale.x)*figsx, ((midy + f * t)*vpScale.y)*figsy, -1);
    00349
    00350                                 a += aa;
    00351                         }
    00352                         glEnd();
    00353                 }
    00354         }
    00355 }