Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] View of /branches/ein16/test/lic-turb2d.diderot
ViewVC logotype

View of /branches/ein16/test/lic-turb2d.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3682 - (download) (annotate)
Thu Feb 18 20:13:18 2016 UTC (3 years, 3 months ago) by cchiw
File size: 2476 byte(s)
creating stable branch that represents ein ir
// lic-turb.diderot
// 
// demo of line integral convolution (LIC) on 2D turbulent flow
//
// First, need to create padded vector dataset v.nrrd with:
// unu pad -i ../data/turb2d.nrrd -min 0 -5 -5 -max M M+5 M+5 -b pad -v 0.0001 -o v.nrrd
// We pad with value 0.0001; we padded with zero than normalize() would fail
// But if you try take too many steps (increasing stepNum) then eventually
// you will leave the vector field domain, which will crash the program,
// so really the padding should be done with a vector field that 
// pointwards inwards, towards the center of the domain, so that anything
// leaving the domain will get pushed back in.  OR, really, the code
// below should be smart and use inside() or some other test to stop
// the integration when we get too close to the domain boundary.
//
// Then, need to create noise texture with:
// unu slice -i ../data/turb2d.nrrd -a 0 -p 0 | unu resample -s 1020 561 | unu 1op nrand -o R.nrrd
// where "1020 561" is copied from imgSizeX and imgSizeY below; we start
// with unu slice in order to inherit the image orientation info, to
// create a noise texture at the same world-space position as the data
//
// process output with:
// unu reshape -i mip.txt -s 1020 561 | unu quantize -b 8  -o lic.png


int imgSizeX = 1020;
int imgSizeY = 561;
real h = 0.005;    // step size of integration
int stepNum = 10; // take this many steps both upstream and downstream

field#1(2)[2] V = image("v.nrrd") ⊛ ctmr;
field#0(2)[] R = image("R.nrrd") ⊛ tent;

strand LIC (int xi, int yi) {
    real xx = lerp(0.0, 6.78, 0.0, real(xi), real(imgSizeX)-1.0);
    real yy = lerp(0.0, 3.74, 0.0, real(yi), real(imgSizeY)-1.0);
    vec2 pos0 = [xx,yy];
    vec2 forw = pos0;
    vec2 back = pos0;
    output real sum = R(pos0);
    int step = 0;

    update {
        // Euler integration step
        // forw = forw + h*V(forw);
        // back = back - h*V(back);
        // Midpoint method step
        forw += h*normalize(V(forw + 0.5*h*normalize(V(forw))));
        back -= h*normalize(V(back - 0.5*h*normalize(V(back))));
        sum += R(forw) + R(back);
        step += 1;
        if (step == stepNum) {
          // the pow() is a way to adjust how much velocity modulates
          // contrast; lower values increase contrast at small velocity
          sum = pow(|V(pos0)|,0.6)*sum/real(1 + 2*stepNum);
          stabilize;
        }
    }
}

initially [ LIC(xi, yi) | yi in 0..(imgSizeY-1), xi in 0..(imgSizeX-1) ];

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0