tadded optional lower cutoff for particles rendered with a colorbar - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
 (HTM) git clone git://src.adamsgaard.dk/sphere
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
 (DIR) commit fa965228a4037d63d6dba70c89b922a1d6b054d1
 (DIR) parent 45188dc7dbc00a30c9eeab280bbd087f055c5360
 (HTM) Author: Anders Damsgaard <adc@geo.au.dk>
       Date:   Fri,  4 Jan 2013 13:14:46 +0100
       
       added optional lower cutoff for particles rendered with a colorbar
       
       Diffstat:
         M src/main.cpp                        |      19 +++++++++++++++----
         M src/raytracer.cuh                   |      65 ++++++++++++++++++++++++-------
         M src/sphere.h                        |       5 +++--
       
       3 files changed, 69 insertions(+), 20 deletions(-)
       ---
 (DIR) diff --git a/src/main.cpp b/src/main.cpp
       t@@ -33,7 +33,8 @@ int main(const int argc, const char *argv[])
            int render = 0; // whether to render an image
            int method = 0; // visualization method
            int nfiles = 0; // number of input files
       -    float max_val = 0.0f;
       +    float max_val = 0.0f;       // max value of colorbar
       +    float lower_cutoff = 0.0f;  // lower cutoff, particles below will not be rendered
        
            // Process input parameters
            int i;
       t@@ -52,9 +53,12 @@ int main(const int argc, const char *argv[])
                        << "-r, --render\t\trender input files instead of simulating temporal evolution\n"
                        << "-dc, --dont-check\tdon't check values before running\n" 
                        << "\nRaytracer (-r) specific options:\n"
       -                << "-m <method> <maxval>, --method <method> <maxval>\n\tcolor visualization method, possible values:\n"
       +                << "-m <method> <maxval> [-l <lower cutoff val>], or\n"
       +                << "--method <method> <maxval> [-l <lower cutoff val>]\n"
       +                << "\tcolor visualization method, possible values:\n"
                        << "\tnormal, pres, vel, angvel, xdisp, angpos\n"
                        << "\t'normal' is the default mode\n"
       +                << "\tif -l is appended, don't render particles with value below\n"
                        << std::endl;
                    return 0; // Exit with success
                }
       t@@ -114,7 +118,14 @@ int main(const int argc, const char *argv[])
                    // Read max. value of colorbar as next argument
                    if (method != 0) {
                        max_val = atof(argv[i+2]);
       -                i += 2; // skip ahead
       +
       +                // Check if a lower cutoff value was specified
       +                if (std::string(argv[i+3]) == "-l") {
       +                    lower_cutoff = atof(argv[i+4]);
       +                    i += 4; // skip ahead
       +                } else {
       +                    i += 2; // skip ahead
       +                }
                    } else {
                        i += 1;
                    }
       t@@ -133,7 +144,7 @@ int main(const int argc, const char *argv[])
        
                    // Render image if requested
                    if (render == 1)
       -                dem.render(method, max_val);
       +                dem.render(method, max_val, lower_cutoff);
        
                    // Otherwise, start iterating through time
                    else
 (DIR) diff --git a/src/raytracer.cuh b/src/raytracer.cuh
       t@@ -143,6 +143,7 @@ __global__ void rayIntersectSpheresColormap(float4* dev_ray_origo,
                Float4* dev_vel,
                Float*  dev_linarr,
                float max_val,
       +        float lower_cutoff,
                unsigned char* dev_img)
        {
            // Compute pixel position from threadIdx/blockIdx
       t@@ -168,6 +169,8 @@ __global__ void rayIntersectSpheresColormap(float4* dev_ray_origo,
            // Iterate through all particles
            for (unsigned int i=0; i<devC_np; ++i) {
        
       +        __syncthreads();
       +
                // Read sphere coordinate and radius
                Float4 x = dev_x[i];
                float3 c = make_float3(x.x, x.y, x.z);
       t@@ -178,32 +181,65 @@ __global__ void rayIntersectSpheresColormap(float4* dev_ray_origo,
                    - 4.0f*dot(d,d)        // -4*A
                    * (dot((e-c),(e-c)) - R*R);  // C
        
       +
       +
                // If the determinant is positive, there are two solutions
                // One where the line enters the sphere, and one where it exits
       -        if (Delta > 0.0f) {
       +        if (lower_cutoff > 0.0) {
        
       -            // Calculate roots, Shirley 2009 p. 77
       -            float t_minus = ((dot(-d,(e-c)) - sqrt( dot(d,(e-c))*dot(d,(e-c)) - dot(d,d)
       -                            * (dot((e-c),(e-c)) - R*R) ) ) / dot(d,d));
       +            // value on colorbar
       +            float val = dev_linarr[i];
        
       -            // Check wether intersection is closer than previous values
       -            if (fabs(t_minus) < tdist) {
       -                p = e + t_minus*d;
       -                tdist = fabs(t_minus);
       -                n = normalize(2.0f * (p - c));   // Surface normal
       -                hitidx = i;
       -            }
       +            // particle is fixed if value > 0
       +            float fixvel = dev_vel[i].w;
        
       -        } // End of solution branch
       +            // only render particles which are above the lower cutoff
       +            // and which are not fixed at a velocity
       +            if (Delta > 0.0f && val > lower_cutoff && fixvel == 0.f) {
       +
       +                // Calculate roots, Shirley 2009 p. 77
       +                float t_minus = ((dot(-d,(e-c)) - sqrt( dot(d,(e-c))*dot(d,(e-c)) - dot(d,d)
       +                                * (dot((e-c),(e-c)) - R*R) ) ) / dot(d,d));
       +
       +                // Check wether intersection is closer than previous values
       +                if (fabs(t_minus) < tdist) {
       +                    p = e + t_minus*d;
       +                    tdist = fabs(t_minus);
       +                    n = normalize(2.0f * (p - c));   // Surface normal
       +                    hitidx = i;
       +                }
       +
       +            } // End of solution branch
       +
       +        } else {
       +
       +            // render particle if it intersects with the ray
       +            if (Delta > 0.0f) {
       +
       +                // Calculate roots, Shirley 2009 p. 77
       +                float t_minus = ((dot(-d,(e-c)) - sqrt( dot(d,(e-c))*dot(d,(e-c)) - dot(d,d)
       +                                * (dot((e-c),(e-c)) - R*R) ) ) / dot(d,d));
       +
       +                // Check wether intersection is closer than previous values
       +                if (fabs(t_minus) < tdist) {
       +                    p = e + t_minus*d;
       +                    tdist = fabs(t_minus);
       +                    n = normalize(2.0f * (p - c));   // Surface normal
       +                    hitidx = i;
       +                }
       +
       +            } // End of solution branch
       +        }
        
            } // End of particle loop
        
            // Write pixel color
            if (tdist < 1e10) {
        
       +        __syncthreads();
                // Fetch particle data used for color
       -        float ratio = dev_linarr[hitidx] / max_val;
                float fixvel = dev_vel[hitidx].w;
       +        float ratio = dev_linarr[hitidx] / max_val;
        
                // Make sure the ratio doesn't exceed the 0.0-1.0 interval
                if (ratio < 0.01f)
       t@@ -362,6 +398,7 @@ __host__ float3 DEM::maxPos()
        __host__ void DEM::render(
                const int method,
                const float maxval,
       +        const float lower_cutoff,
                const float focalLength,
                const unsigned int img_width,
                const unsigned int img_height)
       t@@ -504,7 +541,7 @@ __host__ void DEM::render(
                rayIntersectSpheresColormap<<< blocksPerGrid, threadsPerBlock >>>(
                        dev_ray_origo, dev_ray_direction,
                        dev_x, dev_vel,
       -                dev_linarr, maxval,
       +                dev_linarr, maxval, lower_cutoff,
                        dev_img);
        
            }
 (DIR) diff --git a/src/sphere.h b/src/sphere.h
       t@@ -165,8 +165,9 @@ class DEM {
                // Render particles using raytracing
                void render(
                        const int method = 1,
       -                const float maxval = 1.0e3,
       -                const float focalLength = 1.0,
       +                const float maxval = 1.0e3f,
       +                const float lower_cutoff = 0.0f,
       +                const float focalLength = 1.0f,
                        const unsigned int img_width = 800,
                        const unsigned int img_height = 800);