bounding.glsl

void compute_bounding_points_compact(surface2x3 surf, vec2 u, vec2 v, out bounds2x3 bounds) {

    float surf_c7_c7 = surf.c[7]*surf.c[7];
    float surf_2x_c0 = 2.0*surf.c[0];
    float surf_2x_c1 = 2.0*surf.c[1];
    float surf_2x_c3 = 2.0*surf.c[3];
    float surf_2x_c6 = 2.0*surf.c[6];
    float surf_4x_c1 = 4.0*surf.c[1];
    
    {
	    float surf_4x_c0_c1_sub_c3_c3 = surf_2x_c0*surf_2x_c1 - surf.c[3]*surf.c[3];
        
        // extract contour of quadratic along z plane
        float NXY = surf.c[4]*surf_4x_c1 - surf.c[5]*surf_2x_c3;
        float NY = surf.c[6]*surf_4x_c1 - surf.c[7]*surf_2x_c3;
        float NYY = surf_4x_c0_c1_sub_c3_c3*4.0;

        // compute z bounds of volume
        vec2 vz = surface2x1_bounds(
        	(surf.c[2]*surf_4x_c1 - surf.c[5]*surf.c[5])*NYY - NXY*NXY,
        	(surf.c[8]*surf_4x_c1 - 2.0*surf.c[5]*surf.c[7])*NYY - 2.0*NXY*NY,
        	(surf.c[9]*surf_4x_c1 - surf_c7_c7)*NYY - NY*NY);
        // extract quadratic of plane at volume z bounds
        float slice_z0_c3 = surf.c[4]*vz[0] + surf.c[6];
        float slice_z0_c4 = surf.c[5]*vz[0] + surf.c[7];

        float slice_z1_c3 = surf.c[4]*vz[1] + surf.c[6];
        float slice_z1_c4 = surf.c[5]*vz[1] + surf.c[7];

        // compute position of contact points at volume z bounds
        float f = 1.0/surf_4x_c0_c1_sub_c3_c3;
        bounds.v0 = vec3(
            (surf.c[3]*slice_z0_c4 - surf_2x_c1*slice_z0_c3)*f,
            (surf.c[3]*slice_z0_c3 - surf_2x_c0*slice_z0_c4)*f,
        	vz[0]);
        bounds.v1 = vec3(
            (surf.c[3]*slice_z1_c4 - surf_2x_c1*slice_z1_c3)*f,
            (surf.c[3]*slice_z1_c3 - surf_2x_c0*slice_z1_c4)*f,
        	vz[1]);
    }
    
    // compute z bounds of corner quadratics
    // compute position of contact points at corner z bounds
    {
    	float surf_2x_c7 = 2.0*surf.c[7];
        float surf_c1_v0 = surf.c[1]*v[0];
        float surf_c7_v0 = surf.c[7]*v[0];
        float surf_c1_v1 = surf.c[1]*v[1];
        float surf_c7_v1 = surf.c[7]*v[1];   
        
        // compute z bounds of plane contours    
        // compute position of contact points at plane z bounds
        float surf_2x_c1_neg_rcp = -1.0 / surf_2x_c1;
	    float surf_4x_c9_c1_sub_c7_c7 = surf.c[9]*surf_4x_c1 - surf_c7_c7;
        
        {
            float slice_x0_c0 = (surf.c[0]*u[0] + surf.c[4])*u[0] + surf.c[2];
            float slice_x0_c2 = surf.c[3]*u[0] + surf.c[5];
            float slice_x0_c3 = surf.c[6]*u[0] + surf.c[8];

            {
                vec2 cx0y0 = surface2x1_bounds((surf_c1_v0 + slice_x0_c2)*v[0] + slice_x0_c0, surf_c7_v0 + slice_x0_c3, surf.c[9]);
                bounds.c00_0 = vec3(u[0], v[0], 1.0) * cx0y0[0];
                bounds.c00_1 = vec3(u[0], v[0], 1.0) * cx0y0[1];	
            }
            {
                vec2 cx0y1 = surface2x1_bounds((surf_c1_v1 + slice_x0_c2)*v[1] + slice_x0_c0, surf_c7_v1 + slice_x0_c3, surf.c[9]);
                bounds.c01_1 = vec3(u[0], v[1], 1.0) * cx0y1[1];
                bounds.c01_0 = vec3(u[0], v[1], 1.0) * cx0y1[0];
            }
            {
                vec2 ph0 = surface2x1_bounds(slice_x0_c0*surf_4x_c1 - slice_x0_c2*slice_x0_c2, slice_x0_c3*surf_4x_c1 - slice_x0_c2*surf_2x_c7, surf_4x_c9_c1_sub_c7_c7);
                bounds.eh0_0 = vec3(ph0[0] * u[0], (slice_x0_c2*ph0[0] + surf.c[7]) * surf_2x_c1_neg_rcp, ph0[0]);
                bounds.eh0_1 = vec3(ph0[1] * u[0], (slice_x0_c2*ph0[1] + surf.c[7]) * surf_2x_c1_neg_rcp, ph0[1]);
            }
            
        }
        {
            float slice_x1_c0 = (surf.c[0]*u[1] + surf.c[4])*u[1] + surf.c[2];
            float slice_x1_c2 = surf.c[3]*u[1] + surf.c[5];
            float slice_x1_c3 = surf.c[6]*u[1] + surf.c[8];

            {
                vec2 cx1y0 = surface2x1_bounds((surf_c1_v0 + slice_x1_c2)*v[0] + slice_x1_c0, surf_c7_v0 + slice_x1_c3, surf.c[9]);
                bounds.c10_0 = vec3(u[1], v[0], 1.0) * cx1y0[0];
                bounds.c10_1 = vec3(u[1], v[0], 1.0) * cx1y0[1];
            }
            {
                vec2 cx1y1 = surface2x1_bounds((surf_c1_v1 + slice_x1_c2)*v[1] + slice_x1_c0, surf_c7_v1 + slice_x1_c3, surf.c[9]);
                bounds.c11_1 = vec3(u[1], v[1], 1.0) * cx1y1[1];
                bounds.c11_0 = vec3(u[1], v[1], 1.0) * cx1y1[0];
            }
            {
                vec2 ph1 = surface2x1_bounds(slice_x1_c0*surf_4x_c1 - slice_x1_c2*slice_x1_c2, slice_x1_c3*surf_4x_c1 - slice_x1_c2*surf_2x_c7, surf_4x_c9_c1_sub_c7_c7);        
                bounds.eh1_0 = vec3(ph1[0] * u[1], (slice_x1_c2*ph1[0] + surf.c[7]) * surf_2x_c1_neg_rcp, ph1[0]);
                bounds.eh1_1 = vec3(ph1[1] * u[1], (slice_x1_c2*ph1[1] + surf.c[7]) * surf_2x_c1_neg_rcp, ph1[1]);
            }
            
        }
    }    
    
    {
	    float surf_4x_c0 = 4.0*surf.c[0];        
    	float surf_4x_c9_c0_sub_c6_c6 = surf.c[9]*surf_4x_c0 - surf.c[6]*surf.c[6];
        
        // compute z bounds of plane contours
        // compute position of contact points at plane z bounds
        float surf_2x_c0_neg_rcp = -1.0 / surf_2x_c0;
        {
            float slice_y0_c0 = (surf.c[1]*v[0] + surf.c[5])*v[0] + surf.c[2];
            float slice_y0_c2 = surf.c[3]*v[0] + surf.c[4];
            float slice_y0_c3 = surf.c[7]*v[0] + surf.c[8];
            
            vec2 pv0 = surface2x1_bounds(slice_y0_c0*surf_4x_c0 - slice_y0_c2*slice_y0_c2, slice_y0_c3*surf_4x_c0 - slice_y0_c2*surf_2x_c6, surf_4x_c9_c0_sub_c6_c6);
            bounds.ev0_0 = vec3((slice_y0_c2*pv0[0] + surf.c[6]) * surf_2x_c0_neg_rcp, pv0[0] * v[0], pv0[0]);
            bounds.ev0_1 = vec3((slice_y0_c2*pv0[1] + surf.c[6]) * surf_2x_c0_neg_rcp, pv0[1] * v[0], pv0[1]);
        }
        
        {
            float slice_y1_c0 = (surf.c[1]*v[1] + surf.c[5])*v[1] + surf.c[2];
            float slice_y1_c2 = surf.c[3]*v[1] + surf.c[4];
            float slice_y1_c3 = surf.c[7]*v[1] + surf.c[8];
            
            vec2 pv1 = surface2x1_bounds(slice_y1_c0*surf_4x_c0 - slice_y1_c2*slice_y1_c2, slice_y1_c3*surf_4x_c0 - slice_y1_c2*surf_2x_c6, surf_4x_c9_c0_sub_c6_c6);
            bounds.ev1_0 = vec3((slice_y1_c2*pv1[0] + surf.c[6]) * surf_2x_c0_neg_rcp, pv1[0] * v[1], pv1[0]);
            bounds.ev1_1 = vec3((slice_y1_c2*pv1[1] + surf.c[6]) * surf_2x_c0_neg_rcp, pv1[1] * v[1], pv1[1]);
        }
    }
}