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

SCM Repository

[diderot] Annotation of /tests/examples/vimg/vimg.diderot
ViewVC logotype

Annotation of /tests/examples/vimg/vimg.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4640 - (view) (download)

1 : glk 4640 /* ==========================================
2 :     ## vimg.diderot: 2D image sampler/viewer
3 :    
4 :     This program needs a dataset to render, and a colormap, for example:
5 :    
6 :     ln -s ../data/sscand.nrrd img.nrrd
7 :     ln -s ../cmap/spiral.nrrd cmap.nrrd
8 :    
9 :     Then compile this program; assuming the directions at
10 :     https://github.com/Diderot-Language/examples you can:
11 :    
12 :     ../../vis12/bin/diderotc --exec vimg.diderot
13 :    
14 :     If the needed `img.nrrd` file is missing, the error message looks something like:
15 :    
16 :     uncaught exception Fail [Fail: Nrrd file "img.nrrd" does not exist]
17 :     raised at common/phase-timer.sml:78.50-78.52
18 :     raised at common/phase-timer.sml:78.50-78.52
19 :     raised at nrrd/nrrd-info.sml:146.15-146.74
20 :    
21 :     in which case you should run the `ln -s` command above, or link `img.nrrd`
22 :     to some other 2D scalar nrrd file to view. The same applies to the need
23 :     for `cmap.nrrd` to link to a colormap.
24 :    
25 :     The `-which` option will determine which function is sampled; look
26 :     for `(0 == which)` below to see the start of the function definitions.
27 :     Assuming the symbolic links given above (of `img.nrrd`
28 :     and `cmap.nrrd` to `../data/sscand.nrrd` and `../cmap/spiral.nrrd`
29 :     respectively), some examples usages are:
30 :     * Grayscale image of reconstructed field:
31 :     `./vimg -cent 290 414 -fov 45 -which 0`
32 :     * Grayscale image of gradient magnitude:
33 :     `./vimg -cent 290 414 -fov 45 -which 1`
34 :     * Color image of gradient vector (blue is always zero):
35 :     `./vimg -cent 290 414 -fov 45 -which 2`
36 :     * Colormapped field, with naive pseudo-isocontour:
37 :     `./vimg -cent 290 414 -fov 45 -which 3 -cmin -500 -cmax 1900 -iso 1210 -th 20`
38 :     * Colormapped field, with smarter pseudo-isocontour:
39 :     `./vimg -cent 290 414 -fov 45 -which 4 -cmin -500 -cmax 1900 -iso 1210 -th 0.2`
40 :     * Colormapped field, with ridge lines:
41 :     `./vimg -cent 290 414 -fov 45 -which 5 -cmin -500 -cmax 1900 -th 0.2 -sthr 2`
42 :     * Colormapped field, with valley lines:
43 :     `./vimg -cent 290 414 -fov 45 -which 6 -cmin -500 -cmax 1900 -th 0.2 -sthr 2`
44 :     * Colormapped field, with blue maxima:
45 :     `./vimg -cent 290 414 -fov 45 -which 7 -cmin -500 -cmax 1900 -th 0.25 -sthr 25 -fcol 0 0 0.8`
46 :     * Colormapped field, with green saddle points:
47 :     `./vimg -cent 290 414 -fov 45 -which 8 -cmin -500 -cmax 1900 -th 0.25 -fcol 0 1 0`
48 :     * Colormapped field, with cyan minima:
49 :     `./vimg -cent 290 414 -fov 45 -which 9 -cmin -500 -cmax 1900 -th 0.25 -sthr 25 -fcol 0 1 1`
50 :    
51 :     Each command can be followed by `unu quantize -b 8 -i rgb.nrrd -o rgb.png` to create
52 :     an 8-bit image version of the output. The `-which 3` and `-which 4` commands
53 :     show an important comparison, demonstrating how knowing the gradient permits
54 :     drawing of equal-thickness isocontours (according to the first-order Taylor
55 :     expansion).
56 :    
57 :     Viewing `../data/sscand.nrrd` with the parameters above gives a roughly
58 :     100km view of the area around Geilo, Norway, site of the
59 :     [Winter School](http://www.sintef.no/projectweb/geilowinterschool/2016-scientific-visualization/)
60 :     for which this program was originally created.
61 :     ========================================== */
62 :    
63 :     input image(2)[] img ("2D image dataset to view") = image("../data/sscand.nrrd");
64 :     input image(1)[3] cmimg ("colormap to use") = image("../cmap/spiral.nrrd");
65 :     input vec2 cent ("center of viewing window") = [0,0];
66 :     input real fov ("height of viewing window") = 2;
67 :     input int sx ("# of horizontal samples in viewing window") = 640;
68 :     input int sy ("# of vertical samples in viewing window") = 480;
69 :     input real angle ("orientation (in counter-clockwise degrees) of viewing window") = 0;
70 :     input vec3 bkgd ("RGB color to show outside image domain") = [0.2,0.1,0.1];
71 :     input vec3 fcol ("RGB color to show inside features") = [0.0,0.0,0.0];
72 :     input int which ("what to show about the image, currently from 0 to 9") = 0;
73 :     input real iso ("value at which to show approximate isocontour") = 0;
74 :     input real th ("apparent thickness of feature (meaning varies depending on which)") = 1;
75 :     input real cmin ("value to use at min end of colormap") = 0;
76 :     input real cmax ("value to use at max end of colormap") = 1;
77 :     input real seps ("epsilon >= 0 value to use in divisor when computing strength") = 0.001;
78 :     input real sthr ("strength threshold for showing ridges and valleys") = 0;
79 :    
80 :     /*
81 :     The field to sample and view; needs to be more or less differentiable
82 :     depending on what we want to do with it.
83 :     C1 kernels include: c1tent, ctmr, bspln3, c4hexic
84 :     C2 kernels include: c2ctmr, bspln3, c4hexic
85 :     */
86 :     field#2(2)[] F = c4hexic ⊛ img;
87 :    
88 :     /*
89 :     The colormap can always use tent; it is not differentiated. Note the
90 :     clamp() here: it means that before evaluating cmap(x), x will be clamped to
91 :     inside the domain of cmap, as defined by the location of the samples in the
92 :     cmimg, and the support of the reconstruction filter tent. From a
93 :     functional sense, it might make more sense for the notation to use
94 :     composition: "tent⊛cming◦clamp": clamp first, then lookup into cmimg and
95 :     reconstruct with tent, but this looks a bit weird.
96 :     */
97 :     field#0(1)[3] cmap = tent ⊛ clamp(cmimg);
98 :    
99 :     // compute the image-to-world transform based on specification
100 :     // of the viewing window
101 :     real wdth = fov*sx/sy;
102 :     real phi = angle*π/180;
103 :     tensor[2,2] rot = [[cos(phi),sin(phi)],[-sin(phi),cos(phi)]];
104 :     vec2 spc = [wdth/(sx-1), fov/(sy-1)];
105 :     vec2 dir0 = rot•[spc[0], 0];
106 :     vec2 dir1 = rot•[0, spc[1]];
107 :     vec2 orig = cent - (dir0*(sx-1) + dir1*(sy-1))/2;
108 :    
109 :     // A function that is 1 at 0, and goes down linearly to 0 according to
110 :     // thickness parameter th. Note that if a function is a simple expression like
111 :     // this it can be defined with a single `=`, with no need for a `return`
112 :     // inside a sequence of statements. You might want to call this function
113 :     // "tent" but that name is already taken by a kernel.
114 :     function real bump(real z) = max(0, 1 - |z|/th);
115 :    
116 :     // 2x2 matrix inversion should be part of the language; until
117 :     // then it is pretty easy to create a function for it.
118 :     function tensor[2,2] inv(tensor[2,2] m) {
119 :     real d = m[0,0]*m[1,1] - m[0,1]*m[1,0];
120 :     return [[m[1,1],-m[0,1]],[-m[1,0],m[0,0]]]/d;
121 :     }
122 :    
123 :     strand sample(int ui, int vi) {
124 :     // The output type for a Diderot program is fixed at compile-time,
125 :     // regardless of input variable values or execution outcomes. So
126 :     // this program sometimes saves what is really just grayscale
127 :     // information in an RGB color image.
128 :     output vec3 rgb = bkgd;
129 :     real ll = 0;
130 :     update {
131 :     vec2 x = orig + ui*dir0 + vi*dir1;
132 :     if (inside(x, F)) {
133 :     if (0 == which) {
134 :     // the image value itself
135 :     real y = F(x);
136 :     rgb = [y,y,y];
137 :     } else if (1 == which) {
138 :     // the gradient magnitude
139 :     real y = |∇F(x)|;
140 :     rgb = [y,y,y];
141 :     } else if (2 == which) {
142 :     // the gradient vector, saved in
143 :     // first two components of output
144 :     vec2 gg = ∇F(x);
145 :     rgb = [gg[0],gg[1],0];
146 :     } else if (3 == which) {
147 :     /* Colormapped, with simple value-based isocontour. Note that
148 :     there isn't (currently) a way to learn within the language what
149 :     the domain of colormap (or of any field) is, so rescaling the
150 :     function range to [0,1] here reflects either prior knowledge
151 :     or an assumption about the orientation meta-data in whatever
152 :     cmap.nrrd links to. There is no need to clamp before evaluating
153 :     cmap(); that is handled by the clamp in cmap's declaration */
154 :     rgb = cmap(lerp(0,1,cmin,F(x),cmax));
155 :     ll = bump(F(x) - iso);
156 :     } else if (4 == which) {
157 :     // colormapped, with Newton-based isocontour
158 :     rgb = cmap(lerp(0,1,cmin,F(x),cmax));
159 :     ll = bump((F(x) - iso)/|∇F(x)|);
160 :     } else if (5 == which) {
161 :     // colormapped, with Newton-based ridge line
162 :     rgb = cmap(lerp(0,1,cmin,F(x),cmax));
163 :     vec2{2} evec = evecs(∇⊗∇F(x));
164 :     real{2} eval = evals(∇⊗∇F(x));
165 :     real strn = -eval{1}/(|∇F(x)| + seps);
166 :     if (eval{1} < 0 && strn > sthr) {
167 :     vec2 delta = (evec{1}⊗evec{1}/eval{1})•∇F(x);
168 :     ll = bump(|delta|);
169 :     }
170 :     } else if (6 == which) {
171 :     // colormapped, with Newton-based valley line
172 :     rgb = cmap(lerp(0,1,cmin,F(x),cmax));
173 :     vec2{2} evec = evecs(∇⊗∇F(x));
174 :     real{2} eval = evals(∇⊗∇F(x));
175 :     real strn = eval{0}/(|∇F(x)| + seps);
176 :     if (eval{0} > 0 && strn > sthr) {
177 :     vec2 delta = (evec{0}⊗evec{0}/eval{0})•∇F(x);
178 :     ll = bump(|delta|);
179 :     }
180 :     } else if (7 == which || 8 == which || 9 == which) {
181 :     // colormapped, with critical points
182 :     rgb = cmap(lerp(0,1,cmin,F(x),cmax));
183 :     vec2 delta = inv(∇⊗∇F(x))•∇F(x);
184 :     real{2} eval = evals(∇⊗∇F(x));
185 :     if (7 == which) { // maxima
186 :     if (eval{0} < 0 && eval{1} < 0 && eval{0} < -sthr) {
187 :     ll = bump(|delta|);
188 :     }
189 :     } else if (8 == which) { // saddle
190 :     if (eval{0} > 0 && eval{1} < 0) {
191 :     ll = bump(|delta|);
192 :     }
193 :     } else { // 9 == which: minima
194 :     if (eval{0} > 0 && eval{1} > 0 && eval{1} > sthr) {
195 :     ll = bump(|delta|);
196 :     }
197 :     }
198 :     } else {
199 :     // update "input int which" annotation above as cases are added
200 :     print("Sorry, no function defined for which = ", which, "\n");
201 :     }
202 :     rgb = lerp(rgb, fcol, ll);
203 :     }
204 :     stabilize;
205 :     }
206 :     }
207 :    
208 :     initially [ sample(ui, vi) | vi in 0..(sy-1),
209 :     ui in 0..(sx-1) ];

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