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

SCM Repository

[diderot] Diff of /branches/pure-cfg/test/ridge3d.diderot
ViewVC logotype

Diff of /branches/pure-cfg/test/ridge3d.diderot

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1663, Tue Nov 29 09:48:33 2011 UTC revision 1664, Tue Nov 29 10:43:54 2011 UTC
# Line 1  Line 1 
1    /*! \file ridge3d.diderot
2     *
3     * \author Gordon Kindlmann
4     */
5    
6    /*
7     * COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
8     * All rights reserved.
9     */
10    
11  int gridSize = 120;  int gridSize = 120;
 // HEY this dataset isn't committed yet; GLK's internet connection was bad  
12  field#2(3)[] F = bspln3 ⊛ load("../data/lungcrop.nrrd");  field#2(3)[] F = bspln3 ⊛ load("../data/lungcrop.nrrd");
13  input int stepsMax = 20;  input int stepsMax = 30;
14  real epsilon = 0.005;  input real epsilon = 0.001;
15  real travelMax = 5.0;  input real travelMax = 5.0;
16    input real strnMin = 170.0;
17    
18  strand sample (int ui, int vi, int wi) {  strand sample (int ui, int vi, int wi) {
19      output vec3 pos = [lerp(21.6401, 21.6401 + 140.0*0.5,      output vec3 pos = [lerp(21.6401, 21.6401 + 140.0*0.5,
# Line 16  Line 25 
25      int steps = 0;      int steps = 0;
26      real travel = 0.0;      real travel = 0.0;
27      update {      update {
28          if (!inside(pos, F) || steps > stepsMax || travel > travelMax) {          if (!inside(pos, F) || steps > stepsMax || travel > travelMax)
29              die;              die;
         }  
30          real gmag = |∇F(pos)|;          real gmag = |∇F(pos)|;
31          if (gmag == 0.0) {    // can't compute step if |∇F|, so have to bail          if (gmag == 0.0)    // can't compute step if |∇F|, so have to bail
32              die;              die;
         }  
33          vec3 grad = ∇F(pos);          vec3 grad = ∇F(pos);
34          tensor[3,3] hess = ∇⊗∇F(pos);          tensor[3,3] hess = ∇⊗∇F(pos);
35    
36          real{3} eval = evals(hess);          real{3} eval = evals(hess);
37            real strn = -eval{1};
38          vec3{3} evec = evecs(hess);          vec3{3} evec = evecs(hess);
39          vec3 dir = normalize((evec{2}⊗evec{2} + evec{1}⊗evec{1})•grad);          vec3 dir = normalize((evec{2}⊗evec{2} + evec{1}⊗evec{1})•grad);
40    
41          real fdd = grad•dir;          real fdd = grad•dir;
42          real sdd = dir•hess•dir;          real sdd = dir•hess•dir;
43          vec3 delta = (0.1 if sdd >= 0.0 else -fdd/sdd)*dir;          vec3 delta = (0.1 if sdd >= 0.0 else -fdd/sdd)*dir;
44          if (|delta| < epsilon) {          if (|delta| < epsilon) {
45             if (eval{2} <= -300.0) {             if (strn >= strnMin)
46               stabilize;               stabilize;
            }  
47             die;             die;
48          }          }
         // BUG: can't have:   else {  
49          pos += delta;          pos += delta;
50          steps += 1;          steps += 1;
51          travel += |delta|;          travel += |delta|;
         // BUG: can't have:   }  
52      }      }
53  }  }
54    

Legend:
Removed from v.1663  
changed lines
  Added in v.1664

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