d2::image::defined_scale_by_half(): Replace incorrect operators for geometric mean.
[Ale.git] / d3 / pt.h
blobd9e646cf688cfbb36c5e49c7df1c4f167db1b739
1 // Copyright 2005 David Hilvert <dhilvert@auricle.dyndns.org>,
2 // <dhilvert@ugcs.caltech.edu>
4 /* This file is part of the Anti-Lamenessing Engine.
6 The Anti-Lamenessing Engine is free software; you can redistribute it
7 and/or modify it under the terms of the GNU General Public License as
8 published by the Free Software Foundation; either version 2 of the License,
9 or (at your option) any later version.
11 The Anti-Lamenessing Engine is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with the Anti-Lamenessing Engine; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 * d3/et.h: Represent 3D->2D projective transformations.
25 #ifndef __pt_h__
26 #define __pt_h__
28 #include "space.h"
29 #include "et.h"
32 * Structure to describe a 3D->2D projective transformation. 3D information is
33 * preserved by adding a depth element to the result.
35 * The following coordinate systems are considered:
37 * P: projective
38 * C: local cartesian
39 * W: world
42 struct pt {
43 private:
44 d2::transformation t;
45 et euclidean;
46 ale_real _view_angle; /* XXX: should be ale_pos */
47 ale_pos scale_factor;
48 mutable ale_pos diag_per_depth;
50 public:
53 * Constructor
56 pt() {
57 _view_angle = M_PI / 4;
58 scale_factor = 1;
59 diag_per_depth = 0;
62 pt(d2::transformation t, et e, ale_real va, ale_pos sf = 1) {
63 this->t = t;
64 euclidean = e;
65 _view_angle = va;
66 scale_factor = sf;
67 diag_per_depth = 0;
71 * Output function
74 void debug_output() {
75 t.debug_output();
76 euclidean.debug_output();
77 fprintf(stderr, "[pt.do va=%f sf=%f dpd=%f\n]",
78 _view_angle, scale_factor, diag_per_depth);
82 * Get euclidean transformation reference.
85 et &e() {
86 return euclidean;
90 * Modify scale factor
92 void scale(ale_pos sf) {
93 scale_factor = sf;
97 * Modify or get view angle
99 void view_angle(ale_pos va) {
100 diag_per_depth = 0;
101 _view_angle = va;
104 ale_pos view_angle() {
105 return _view_angle;
109 * Get the 2D scale factor
111 ale_pos scale_2d() const {
112 return t.scale();
116 * Transform W to C.
118 point wc(point p) const {
119 return euclidean(p);
123 * Transform C to P for given width and height.
125 point cp_generic(point p, ale_pos w, ale_pos h) const {
127 * Divide x and y by negative z
130 p[0] /= -p[2];
131 p[1] /= -p[2];
134 * Scale x and y
137 ale_pos scaling_factor = sqrt(w*w + h*h) / (2 * tan(_view_angle / 2));
138 p[0] *= scaling_factor;
139 p[1] *= scaling_factor;
142 * Add an offset so that the upper-left corner is the origin.
145 p[0] += h / 2;
146 p[1] += w / 2;
148 return p;
152 * Transform point p.
154 struct point wp_generic(struct point p, ale_pos w, ale_pos h) const {
155 return cp_generic(wc(p), w, h);
159 * Width and height
162 ale_pos scaled_width() const {
163 return t.scaled_width() * scale_factor;
166 ale_pos scaled_height() const {
167 return t.scaled_height() * scale_factor;
170 int scaled_in_bounds(point p) const {
171 return (p[0] >= 0 && p[0] <= scaled_height() - 1
172 && p[1] >= 0 && p[1] <= scaled_width() - 1);
175 int unscaled_in_bounds(point p) const {
176 return (p[0] >= 0 && p[0] <= unscaled_height() - 1
177 && p[1] >= 0 && p[1] <= unscaled_width() - 1);
180 ale_pos unscaled_width() const {
181 return t.unscaled_width() * scale_factor;
184 ale_pos unscaled_height() const {
185 return t.unscaled_height() * scale_factor;
189 * Scaled transforms
192 point cp_scaled(point p) const {
193 return cp_generic(p, scaled_width(), scaled_height());
196 point wp_scaled(point p) const {
197 return wp_generic(p, scaled_width(), scaled_height());
201 * Unscaled transforms
204 point cp_unscaled(point p) const {
205 return cp_generic(p, unscaled_width(), unscaled_height());
208 point wp_unscaled(point p) const {
209 return wp_generic(p, unscaled_width(), unscaled_height());
213 * Transform P to C.
215 point pc_generic(point p, ale_pos w, ale_pos h) const {
217 * Subtract offset
220 p[0] -= h / 2;
221 p[1] -= w / 2;
224 * Scale x and y
227 ale_pos scaling_factor = sqrt(w*w + h*h) / (2 * tan(_view_angle / 2));
228 p[0] /= scaling_factor;
229 p[1] /= scaling_factor;
232 * Multiply x and y by negative z
235 p[0] *= -p[2];
236 p[1] *= -p[2];
238 return p;
242 * Transform C to W
244 point cw(point p) const {
245 return euclidean.inverse_transform(p);
249 * Transform P to W
251 point pw_generic(point p, ale_pos w, ale_pos h) const {
252 return cw(pc_generic(p, w, h));
256 * Inverse transforms for scaled points.
259 point pc_scaled(point p) const {
260 return pc_generic(p, scaled_width(), scaled_height());
263 point pw_scaled(point p) const {
264 return pw_generic(p, scaled_width(), scaled_height());
268 * Inverse transforms for unscaled points.
271 point pc_unscaled(point p) const {
272 return pc_generic(p, unscaled_width(), unscaled_height());
275 point pw_unscaled(point p) const {
276 return pw_generic(p, unscaled_width(), unscaled_height());
280 * Density calculation
283 ale_pos c_density_scaled(point p) const {
284 ale_pos one_density = 1 / (pc_scaled(point(0, 0, -1)).lengthto(pc_scaled(point(0, 1, -1)))
285 * pc_scaled(point(0, 0, -1)).lengthto(pc_scaled(point(1, 0, -1))));
287 return one_density / (p[2] * p[2]);
290 ale_pos w_density_scaled(point p) const {
291 return c_density_scaled(wc(p));
294 ale_pos w_density_scaled_max(point w0, point w1, point w2) {
295 point c0 = wc(w0);
296 point c1 = wc(w1);
297 point c2 = wc(w2);
300 * Select the point closest to the camera.
303 if (c0[2] > c1[2] && c0[2] > c2[2])
304 return c_density_scaled(c0);
305 else if (c1[2] > c2[2])
306 return c_density_scaled(c1);
307 else
308 return c_density_scaled(c2);
311 private:
312 void calculate_diag_per_depth() const {
313 if (diag_per_depth != 0)
314 return;
315 ale_pos w = unscaled_width();
316 ale_pos h = unscaled_height();
318 diag_per_depth = sqrt(2) * (2 * tan(_view_angle / 2)) / sqrt(w*w + h*h);
321 public:
324 * Get a trilinear coordinate for a given depth.
326 ale_pos trilinear_coordinate(ale_pos depth, ale_pos diagonal) {
327 calculate_diag_per_depth();
329 return log(diagonal / (fabs(depth) * diag_per_depth)) / log(2);
334 * Get a trilinear coordinate for a given position in the world and
335 * a given 2D diagonal distance.
337 ale_pos trilinear_coordinate(point w, ale_pos diagonal) {
338 return trilinear_coordinate(wc(w)[2], diagonal);
342 * Get a trilinear coordinate for a given subspace.
344 ale_pos trilinear_coordinate(const space::traverse &st) {
345 point min = st.get_min();
346 point max = st.get_max();
347 point avg = (min + max) / (ale_pos) 2;
349 ale_pos diagonal = min.lengthto(max) * sqrt(2) / sqrt(3);
351 return trilinear_coordinate(avg, diagonal);
355 * Get a diagonal distance for a given position in the world
356 * and a given trilinear coordinate.
358 ale_pos diagonal_distance(point w, ale_pos coordinate) const {
359 calculate_diag_per_depth();
361 ale_pos depth = fabs(wc(w)[2]);
362 ale_pos diagonal = pow(2, coordinate) * depth * diag_per_depth;
364 return diagonal;
368 * Get the 3D diagonal for a given depth and trilinear coordinate.
370 ale_pos diagonal_distance_3d(ale_pos depth, ale_pos coordinate) const {
371 calculate_diag_per_depth();
372 return pow(2, coordinate) * fabs(depth) * diag_per_depth * sqrt(3) / sqrt(2);
376 * Get the 1D distance for a given depth and trilinear coordinate.
378 ale_pos distance_1d(ale_pos depth, ale_pos coordinate) const {
379 calculate_diag_per_depth();
380 return pow(2, coordinate) * fabs(depth) * diag_per_depth / sqrt(2);
383 ale_pos distance_1d(point iw, ale_pos coordinate) const {
384 if (wc(iw)[2] >= 0)
385 return point::undefined()[0];
386 return distance_1d(-wc(iw)[2], coordinate);
390 * Check for inclusion of a point in the bounding box of projected
391 * vertices. This function returns non-zero when a point is included,
392 * when one of the vertices is infinite or undefined, or when a vertex
393 * is behind the point of projection.
395 * WBB is assumed to contain {volume_min, volume_max}.
398 int check_inclusion(const point *wbb, const d2::point &pc_min, const d2::point &pc_max, int scaled) const {
400 assert(pc_min[0] <= pc_max[0]);
401 assert(pc_min[1] <= pc_max[1]);
403 int test[2][2] = {{0, 0}, {0, 0}};
405 for (int x = 0; x < 2; x++)
406 for (int y = 0; y < 2; y++)
407 for (int z = 0; z < 2; z++) {
409 point p = scaled ? wp_scaled(point(wbb[x][0], wbb[y][1], wbb[z][2]))
410 : wp_unscaled(point(wbb[x][0], wbb[y][1], wbb[z][2]));
412 if (!p.finite())
413 return 1;
415 if (p[2] > 0)
416 return 1;
418 for (int d = 0; d < 2; d++) {
419 if (p[d] <= pc_max[d])
420 test[d][0] = 1;
421 if (p[d] >= pc_min[d])
422 test[d][1] = 1;
426 for (int d = 0; d < 2; d++)
427 for (int c = 0; c < 2; c++)
428 if (test[d][c] == 0)
429 return 0;
431 return 1;
434 int check_inclusion_scaled(const point *wbb, d2::point pc_min, d2::point pc_max) const {
435 return check_inclusion(wbb, pc_min, pc_max, 1);
438 int check_inclusion_scaled(const space::traverse &st, d2::point pc_min, d2::point pc_max) const {
439 return check_inclusion_scaled(st.get_bounds(), pc_min, pc_max);
442 int check_inclusion_scaled(const space::traverse &st, d2::point pc) {
443 return check_inclusion_scaled(st, pc, pc);
447 * Get bounding box for projection of a subspace.
450 void get_view_local_bb(point volume_min, point volume_max, point bb[2], int scaled) const {
452 point min = point::posinf();
453 point max = point::neginf();
455 point wbb[2] = { volume_min, volume_max };
457 for (int x = 0; x < 2; x++)
458 for (int y = 0; y < 2; y++)
459 for (int z = 0; z < 2; z++) {
460 point p = scaled ? wp_scaled(point(wbb[x][0], wbb[y][1], wbb[z][2]))
461 : wp_unscaled(point(wbb[x][0], wbb[y][1], wbb[z][2]));
463 for (int d = 0; d < 3; d++) {
464 if (p[d] < min[d])
465 min[d] = p[d];
466 if (p[d] > max[d])
467 max[d] = p[d];
472 * Clip bounding box to image extents.
475 if (min[0] < 0)
476 min[0] = 0;
477 if (min[1] < 0)
478 min[1] = 0;
479 if (max[0] > scaled_height() - 1)
480 max[0] = scaled_height() - 1;
481 if (max[1] > scaled_width() - 1)
482 max[1] = scaled_width() - 1;
484 bb[0] = min;
485 bb[1] = max;
488 void get_view_local_bb_unscaled(point volume_min, point volume_max, point bb[2]) const {
489 get_view_local_bb(volume_min, volume_max, bb, 0);
492 void get_view_local_bb_scaled(point volume_min, point volume_max, point bb[2]) const {
493 get_view_local_bb(volume_min, volume_max, bb, 1);
496 void get_view_local_bb_scaled(const space::traverse &st, point bb[2]) const {
497 get_view_local_bb_scaled(st.get_min(), st.get_max(), bb);
500 void get_view_local_bb_unscaled(const space::traverse &t, point bb[2]) const {
501 get_view_local_bb_unscaled(t.get_min(), t.get_max(), bb);
505 * Get the in-bounds centroid for a subspace, if one exists.
508 point centroid(point volume_min, point volume_max) const {
510 point bb[2];
512 get_view_local_bb_unscaled(volume_min, volume_max, bb);
514 point min = bb[0];
515 point max = bb[1];
517 for (int d = 0; d < 2; d++)
518 if (min[d] > max[d])
519 return point::undefined();
521 if (min[2] >= 0)
522 return point::undefined();
524 if (max[2] > 0)
525 max[2] = 0;
527 return (max + min) / 2;
530 point centroid(const space::traverse &t) {
531 return centroid(t.get_min(), t.get_max());
535 * Get the local space origin in world space.
538 point origin() {
539 return cw(point(0, 0, 0));
543 #endif