1 /**********
2 *
3 *
4 * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
5 * University of Manchester, UK. All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without modification,
8 * are permitted provided that the following conditions are met:
9 *
10 * . Redistributions of source code must retain the above copyright notice,
11 * this list of conditions and the following disclaimer.
12 *
13 * . Redistributions in binary form must reproduce the above copyright notice,
14 * this list of conditions and the following disclaimer in the documentation
15 * and/or other materials provided with the distribution.
16 *
17 * . Neither the name of the University of Manchester nor the names of its
18 * contributors may be used to endorse or promote products derived from this
19 * software without specific prior written permission.
20 *
21 *
22 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
23 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
24 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
25 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
26 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
27 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
28 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
29 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
30 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
31 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
32 * POSSIBILITY OF SUCH DAMAGE.
33 **********
34 *
35 * Program : TINA
36 * File : $Source: /home/tina/cvs/tina-tools/tinatool/tlvision/tlvisTrn_procs.c,v $
37 * Date : $Date: 2009/03/23 16:08:53 $
38 * Version : $Revision: 1.6 $
39 * CVS Id : $Id: tlvisTrn_procs.c,v 1.6 2009/03/23 16:08:53 paul Exp $
40 *
41 * Notes :
42 *
43 *
44 *
45 *********
46 */
47
48
49 #if HAVE_CONFIG_H
50 #include <config.h>
51 #endif
52
53 #include <math.h>
54 #include <stdlib.h>
55
56 #include <tina/sys/sysDef.h>
57 #include <tina/sys/sysPro.h>
58 #include <tina/math/mathDef.h>
59 #include <tina/math/mathPro.h>
60 #include <tina/geometry/geomDef.h>
61 #include <tina/geometry/geomPro.h>
62 #include <tina/image/imgDef.h>
63 #include <tina/image/imgPro.h>
64 #include <tinatool/draw/drawPro.h>
65
66 /* needed for triangle */
67 #define REAL double
68
69 #include <tinatool/tlvision/tlvisTrn_triangle.h>
70 #include "tlvisTrn_procs.h"
71
72
73 #define IMAGE_FILL 0
74 #define SHADE_FILL 1
75 #define MEAN_FILL 2
76
77
78 static Terrain_data *terrain = NULL;
79 static Imrect *im = NULL;
80 static Imrect *shade_im = NULL;
81 static Imrect *mean_im = NULL;
82
83 static float scalez = 0.25;
84 static int samplex = 16, sampley = 16;
85
86 static Bool show_top = true;
87 static Bool show_bot = false;
88 static Bool show_lines = true;
89 static Bool show_fill = false;
90 static Bool show_all = false;
91
92 static int fill_type = IMAGE_FILL;
93 /* static int line_col; */
94
95 void terr_fill_set(int new_fill_type)
96 {
97 fill_type = new_fill_type;
98 }
99
100 /*
101 void terr_line_col_set(int col)
102 {
103 line_col = col;
104 }
105 */
106
107 void terr_show_lines_set(Bool show)
108 {
109 show_lines = show;
110 }
111
112 void terr_show_all_set(Bool show)
113 {
114 show_all = show;
115 }
116
117 void terr_show_fill_set(Bool show)
118 {
119 show_fill = show;
120 }
121
122 void terr_show(Bool top, Bool bot)
123 {
124 show_top = top;
125 show_bot = bot;
126 }
127
128 void terr_im_set(Imrect * new_im)
129 {
130 im = new_im;
131 }
132
133 Imrect *terr_im_get(void)
134 {
135 return (im);
136 }
137
138 void terr_samplex_set(int sx)
139 {
140 samplex = sx;
141 }
142
143 int terr_samplex_get(void)
144 {
145 return (samplex);
146 }
147
148 void terr_sampley_set(int sy)
149 {
150 sampley = sy;
151 }
152
153 int terr_sampley_get(void)
154 {
155 return (sampley);
156 }
157
158 void terr_scalez_set(double sz)
159 {
160 scalez = sz;
161 }
162
163 double terr_scalez_get(void)
164 {
165 return (scalez);
166 }
167
168 void terrain_skeldraw(Tv * tv)
169 {
170 if (terrain == NULL)
171 return;
172
173 tv_screen_double_buffer_start(tv->tv_screen);
174 tv_terrain_skel(tv, terrain->data, terrain->mask, terrain->m, terrain->n);
175 tv_screen_double_buffer_end(tv->tv_screen);
176 }
177
178 static int patch_shade(int i, int j, Tv * tv)
179 {
180 Vec3 p00 = {Vec3_id};
181 Vec3 p01 = {Vec3_id};
182 Vec3 p10 = {Vec3_id};
183 Vec3 n = {Vec3_id};
184
185 if (tv == NULL)
186 return (0);
187 p00 = terrain->data[i][j];
188 p01 = terrain->data[i][j + 1];
189 p10 = terrain->data[i + 1][j];
190 /* BUGFIX Vec3 p11; p11 = terrain->data[i + 1][j + 1]; */
191
192 n = vec3_unit(vec3_cross(vec3_diff(p10, p00), vec3_diff(p01, p00)));
193 return ((int) (64.0 + 191.0 * fabs(vec3_dot(n, tv->ez3))));
194 }
195
196 static int patch_mean(int i, int j)
197 {
198 double sum = 0;
199
200 sum += vec3_z(terrain->data[i][j]);
201 sum += vec3_z(terrain->data[i][j + 1]);
202 sum += vec3_z(terrain->data[i + 1][j]);
203 sum += vec3_z(terrain->data[i + 1][j + 1]);
204 sum /= 4.0 * scalez;
205 return ((sum > 255) ? 255 : sum);
206 }
207
208 static Imrect *get_image_for_surf(int (*intensity_func) ( /* ??? */ ), void *data)
209 {
210 Imregion *region;
211 Imrect *im;
212 int n, m;
213 int i, j;
214 int x, y;
215 int lx, ux, ly, uy;
216 double dx, dy;
217 int g, *row;
218
219 if (terrain == NULL)
220 return (NULL);
221
222 m = terrain->m;
223 n = terrain->n;
224
225 lx = vec3_x(terrain->data[0][0]);
226 ly = vec3_y(terrain->data[0][0]);
227 ux = vec3_x(terrain->data[m - 1][n - 1]) + 1;
228 uy = vec3_y(terrain->data[m - 1][n - 1]) + 1;
229 dx = (double) (ux - lx) / (n - 1);
230 dy = (double) (uy - ly) / (m - 1);
231
232 region = roi_alloc(lx, ly, ux, uy);
233 im = im_alloc(uy, ux, region, uchar_v);
234 rfree((void *) region);
235
236 row = ivector_alloc(lx, ux);
237 for (i = 0; i < m - 1; i++)
238 {
239 for (j = 0; j < n - 1; j++)
240 {
241 g = intensity_func(i, j, data);
242 for (x = j * dx; x < (j + 1) * dx && x + lx < ux; x++)
243 row[lx + x] = g;
244 }
245 for (y = i * dy; y < (i + 1) * dy && y + ly < uy; y++)
246 im_put_row(row, im, y + ly, lx, ux);
247 }
248 ivector_free(row, lx);
249 return (im);
250 }
251
252 void terrain_fulldraw(Tv * tv)
253 {
254 int n, m;
255
256 if (terrain == NULL)
257 return;
258
259 n = terrain->n;
260 m = terrain->m;
261
262 /* tv_set_color(tv, line_col); */
263 (void) tv_set_op(tv, (show_lines == true) ? rop_copy : rop_null);
264 tv_terrain_hidden_init(tv, show_bot, show_top);
265 if (show_fill == true)
266 {
267 switch (fill_type)
268 {
269 case IMAGE_FILL:
270 tv_terrain_fill3(tv, terrain->data, m, n, im);
271 break;
272 case SHADE_FILL:
273 tv_terrain_fill3(tv, terrain->data, m, n, shade_im);
274 break;
275 case MEAN_FILL:
276 tv_terrain_fill3(tv, terrain->data, m, n, mean_im);
277 break;
278 }
279 } else if (show_all)
280 tv_terrain3(tv, terrain->data, m, n);
281 else
282 tv_terrain_hidden3(tv, terrain->data, m, n);
283 (void) tv_set_op(tv, rop_copy);
284 }
285
286 void terrain_backdraw(Tv * tv)
287 {
288 }
289
290 void terr_terrain_set(Terrain_data * new_terrain)
291 {
292 terrain_data_free(terrain);
293 terrain = new_terrain;
294 im_free(shade_im);
295 shade_im = get_image_for_surf(patch_shade, (void *) terrain_tv());
296 im_free(mean_im);
297 mean_im = get_image_for_surf(patch_mean, NULL);
298 }
299
300 void terrain_from_image(void)
301 {
302 if (im == NULL)
303 return;
304
305 terrain_data_free(terrain);
306 terrain = im_surface(im, (Imrect *) NULL, samplex, sampley, scalez);
307 im_free(shade_im);
308 shade_im = get_image_for_surf(patch_shade, (void *) terrain_tv());
309 im_free(mean_im);
310 mean_im = get_image_for_surf(patch_mean, NULL);
311 }
312
313 void terrain_init(Tv * tv)
314 {
315 Vec3 centre = {Vec3_id};
316 Vec3 aim = {Vec3_id};
317 Vec3 down = {Vec3_id};
318 float radius, pscale;
319 int m, n;
320
321 if (tv == NULL || terrain == NULL)
322 return;
323
324 m = terrain->m;
325 n = terrain->n;
326 centre = terrain->data[m / 2][n / 2];
327 radius = 0.5 * vec3_mod(vec3_diff(terrain->data[0][0], terrain->data[m - 1][n - 1]));
328 pscale = 3;
329 /**
330 aim = vec3_unit(vec3(1.0, 1.0, -1.0));
331 down = vec3_minus(vec3_ey());
332 tv_set_axis(tv, vec3_minus(vec3_ez()));
333 **/
334 aim = vec3_unit(vec3(0.0, -0.5, 1.0));
335 down = vec3_ey();
336 tv_camera3(tv, centre, radius, pscale, aim, down);
337 tv_set_axis(tv, vec3_ez());
338 }
339
340 Tv *terrain_tv(void)
341 {
342 static Tv *tv = NULL;
343
344 if (tv != NULL)
345 return (tv);
346
347 tv = tv_create("Terrain");
348 (void) tv_set_backdraw(tv, terrain_backdraw);
349 (void) tv_set_fulldraw(tv, terrain_fulldraw);
350 (void) tv_set_skeldraw(tv, terrain_skeldraw);
351 tv_set_init(tv, terrain_init);
352 (void) tv_set_zoomlevel(tv, ZOOMGR);
353 return (tv);
354 }
355
356 static void triangle_init(struct triangulateio *in, struct triangulateio *out)
357 {
358 /* initialise the input triangle structure */
359
360 in->pointlist = NULL;
361 in->pointattributelist = NULL;
362 in->pointmarkerlist = NULL;
363 in->numberofpoints = 0;
364 in->numberofpointattributes = 0;
365
366 in->trianglelist = NULL;
367 in->triangleattributelist = NULL;
368 in->trianglearealist = NULL;
369 in->neighborlist = NULL;
370 in->numberoftriangles = 0;
371 in->numberofcorners = 0;
372 in->numberoftriangleattributes = 0;
373
374 in->segmentlist = NULL;
375 in->segmentmarkerlist = NULL;
376 in->numberofsegments = 0;
377
378 in->holelist = NULL;
379 in->numberofholes = 0;
380
381 in->regionlist = NULL;
382 in->numberofregions = 0;
383
384 in->edgelist = NULL;
385 in->edgemarkerlist = NULL;
386 in->normlist = NULL;
387 in->numberofedges = 0;
388
389 /* initialise the output triangle structure */
390
391 out->pointlist = NULL;
392 out->pointattributelist = NULL;
393 out->pointmarkerlist = NULL;
394 out->numberofpoints = 0;
395 out->numberofpointattributes = 0;
396
397 out->trianglelist = NULL;
398 out->triangleattributelist = NULL;
399 out->trianglearealist = NULL;
400 out->neighborlist = NULL;
401 out->numberoftriangles = 0;
402 out->numberofcorners = 0;
403 out->numberoftriangleattributes = 0;
404
405 out->segmentlist = NULL;
406 out->segmentmarkerlist = NULL;
407 out->numberofsegments = 0;
408
409 out->holelist = NULL;
410 out->numberofholes = 0;
411
412 out->regionlist = NULL;
413 out->numberofregions = 0;
414
415 out->edgelist = NULL;
416 out->edgemarkerlist = NULL;
417 out->normlist = NULL;
418 out->numberofedges = 0;
419 }
420
421 /* same_side: returns true if pixel in question is on the
422 same side of the line (defined by p1, p2) as p3 */
423
424 static Bool same_side(Vec2 p1, Vec2 p2, Vec2 p3, Vec2 pix)
425 {
426 float m, p3test, pixtest;
427
428 /* vertical line case */
429 if (vec2_x(p1) == vec2_x(p2))
430 {
431 /* both p3 and pix to right of line */
432 if ((vec2_x(p3) > vec2_x(p1)) && (vec2_x(pix) >= vec2_x(p1)))
433 return(true);
434 /* both p3 and pix to left of line */
435 if ((vec2_x(p3) < vec2_x(p1)) && (vec2_x(pix) <= vec2_x(p1)))
436 return(true);
437 }
438
439 /* horizontal line case */
440 else if (vec2_y(p1) == vec2_y(p2))
441 {
442 /* both p3 and pix above line */
443 if ((vec2_y(p3) > vec2_y(p1)) && (vec2_y(pix) >= vec2_y(p1)))
444 return(true);
445 /* both p3 and pix below line */
446 if ((vec2_y(p3) < vec2_y(p1)) && (vec2_y(pix) <= vec2_y(p1)))
447 return(true);
448 }
449
450 /* sloping line case */
451 else
452 {
453 /* find the line gradient */
454 m = (vec2_y(p2) - vec2_y(p1)) / (vec2_x(p2) - vec2_x(p1));
455
456 if (abs(m) <= 1)
457 {
458 /* use x knowns, look at y vals */
459 p3test = m * (vec2_x(p3) - vec2_x(p1)) + vec2_y(p1);
460 pixtest = m * (vec2_x(pix) - vec2_x(p1)) + vec2_y(p1);
461
462 /* both points below line */
463 if ((p3test > vec2_y(p3)) && (pixtest >= vec2_y(pix)))
464 return(true);
465 /* both points above line */
466 if ((p3test < vec2_y(p3)) && (pixtest <= vec2_y(pix)))
467 return(true);
468 }
469 else
470 {
471 /* use y knowns, look at x vals */
472 p3test = vec2_x(p1) + (vec2_y(p3) - vec2_y(p1)) / m;
473 pixtest = vec2_x(p1) + (vec2_y(pix) - vec2_y(p1)) / m;
474
475 /* both points to the left of line */
476 if ((p3test > vec2_x(p3)) && (pixtest >= vec2_x(pix)))
477 return(true);
478 /* both points to the right of line */
479 if ((p3test < vec2_x(p3)) && (pixtest <= vec2_x(pix)))
480 return(true);
481 }
482 }
483
484 return(false);
485 }
486
487 /* interp_val: given a triangle in 3d space (p1, p2, p3)
488 what is the z value of a point on that plane with x y
489 coordinates given by pix */
490
491 static float interp_val(Vec3 p1, Vec3 p2, Vec3 p3, Vec2 pix)
492 {
493 Vec3 normal, interp;
494
495 normal = vec3_cross(vec3_diff(p1, p3), vec3_diff(p2, p3));
496 interp = vec3_inter_line_plane(vec3(vec2_x(pix), vec2_y(pix), 0.0),
497 vec3(0.0, 0.0, 1.0),
498 p1, normal);
499 return(vec3_z(interp));
500 }
501
502 /* triangulate_IM: uses delaunay triangulation routine to
503 divide up a set of points into triangles and then interpolates
504 those triangles to fill in the unknown points in the triangles */
505 /*
506 void triangulate_IM(Imrect *im)
507 {
508 struct triangulateio in, out;
509 Vec2 *lut = NULL, p1, p2, p3;
510 Vec3 t1, t2, t3;
511 int num, x, y, xmin, xmax, ymin, ymax;
512
513 if (im == NULL || im->vtype!=float_v)
514 {
515 error("No IM or IM not of type float", non_fatal);
516 return;
517 }
518
519 triangle_init(&in, &out);
520
521 FOR_IM(im->region, y, x)
522 if (IM_FLOAT(im, y, x))
523 in.numberofpoints++;
524
525 lut = (Vec2 *)calloc(in.numberofpoints, sizeof(Vec2));
526
527
528 in.pointlist = (REAL *)malloc(in.numberofpoints * 2 * sizeof(REAL));
529
530 num = 0;
531
532 FOR_IM(im->region, y, x)
533 {
534 if (IM_FLOAT(im, y, x) == 0.0)
535 continue;
536
537 vec2_x(lut[num]) = x;
538 vec2_y(lut[num]) = y;
539
540 in.pointlist[2*num] = x;
541 in.pointlist[2*num+1] = y;
542 num++;
543 }
544
545 triangulate("zPNQ", &in, &out, NULL);
546
547 for (num=0; num<out.numberoftriangles; num++)
548 {
549 p1 = lut[out.trianglelist[3*num]];
550 p2 = lut[out.trianglelist[3*num+1]];
551 p3 = lut[out.trianglelist[3*num+2]];
552
553 xmin = tina_int(MIN(MIN(vec2_x(p1), vec2_x(p2)), vec2_x(p3)));
554 ymin = tina_int(MIN(MIN(vec2_y(p1), vec2_y(p2)), vec2_y(p3)));
555 xmax = tina_int(MAX(MAX(vec2_x(p1), vec2_x(p2)), vec2_x(p3)));
556 ymax = tina_int(MAX(MAX(vec2_y(p1), vec2_y(p2)), vec2_y(p3)));
557
558 for (y=ymin; y<=ymax; y++)
559 for (x=xmin; x<=xmax; x++)
560 if ((IM_FLOAT(im, y, x) == 0.0) &&
561 same_side(p1, p2, p3, vec2(x, y)) &&
562 same_side(p1, p3, p2, vec2(x, y)) &&
563 same_side(p2, p3, p1, vec2(x, y)))
564 {
565
566 t1 = vec3(vec2_x(p1), vec2_y(p1),
567 IM_FLOAT(im, tina_int(vec2_y(p1)), tina_int(vec2_x(p1))));
568 t2 = vec3(vec2_x(p2), vec2_y(p2),
569 IM_FLOAT(im, tina_int(vec2_y(p2)), tina_int(vec2_x(p2))));
570 t3 = vec3(vec2_x(p3), vec2_y(p3),
571 IM_FLOAT(im, tina_int(vec2_y(p3)), tina_int(vec2_x(p3))));
572
573 IM_FLOAT(im, y, x) = interp_val(t1, t2, t3, vec2(x, y));
574 }
575 }
576
577 cfree(lut);
578 free(in.pointlist);
579 free(out.trianglelist);
580 }
581 */
582
583 /* triangulate_geom: uses delaunay triangulation routine to triangulate
584 a set of 3D points. the main triangulation routine is driven by
585 the x,y coordinates of the points. the results can then be displayed
586 in the threed tv for example */
587
588 /* note: this proc adds the triangulation lines to the existing
589 geomlist, it does NOT free the given geomlist and return a
590 new one */
591 /*
592 List *triangulate_geom(List *geomlist)
593 {
594 struct triangulateio in, out;
595 List *ptr;
596 Point3 *point, **lut = NULL;
597 Vec3 p1, p2;
598 int num;
599
600 if (geomlist == NULL)
601 return(geomlist);
602
603 triangle_init(&in, &out);
604
605 for (ptr=geomlist; ptr!=NULL; ptr=ptr->next)
606 if (ptr->type == POINT3)
607 in.numberofpoints++;
608
609 lut = (Point3 **)calloc(in.numberofpoints, sizeof(Point3 *));
610
611 in.pointlist = (REAL *)malloc(in.numberofpoints * 2 * sizeof(REAL));
612
613 for (num=0, ptr=geomlist; ptr!=NULL; ptr=ptr->next)
614 {
615 if (ptr->type != POINT3)
616 continue;
617
618 point = (Point3 *)ptr->to;
619
620 lut[num] = point;
621
622 in.pointlist[2*num] = vec3_x(point->p);
623 in.pointlist[2*num+1] = vec3_y(point->p);
624 num++;
625 }
626
627 triangulate("zPNQ", &in, &out, NULL);
628
629
630 for (num=0; num<out.numberoftriangles; num++)
631 {
632 p1 = lut[out.trianglelist[3*num]]->p;
633 p2 = lut[out.trianglelist[3*num+1]]->p;
634 geomlist = list_addtostart(geomlist,
635 link_alloc((void *)line3_make(p1, p2, LINE3),
636 LINE3));
637
638 p1 = lut[out.trianglelist[3*num+1]]->p;
639 p2 = lut[out.trianglelist[3*num+2]]->p;
640 geomlist = list_addtostart(geomlist,
641 link_alloc((void *)line3_make(p1, p2, LINE3),
642 LINE3));
643
644 p1 = lut[out.trianglelist[3*num+2]]->p;
645 p2 = lut[out.trianglelist[3*num]]->p;
646 geomlist = list_addtostart(geomlist,
647 link_alloc((void *)line3_make(p1, p2, LINE3),
648 LINE3));
649 }
650
651 cfree(lut);
652 free(in.pointlist);
653 free(out.trianglelist);
654
655 return(geomlist);
656 }
657 */
658
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.