Hugintrunk  0.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VertexCoordRemapper.cpp
Go to the documentation of this file.
1 // -*- c-basic-offset: 4 -*-
2 
23 #ifdef __WXMAC__
24 #include "panoinc_WX.h"
25 #include "panoinc.h"
26 #endif
27 
28 #include "VertexCoordRemapper.h"
29 #include <math.h>
30 #include <iostream>
31 #include "ViewState.h"
32 
33 #include "panodata/SrcPanoImage.h"
34 
35 /*******************
36  * Face's bit flags
37  *******************/
38 // split_x is set if the node has been split into two children (subdivided) in x
39 const unsigned short int split_flag_x = 1;
40 // patch_flag_x is set if we are patching a hole / overlap that could be caused
41 // by the subdivision in the previous x location being more dense. It is set
42 // instead of split_flag_x, but we subdivide into faces that do not provide more
43 // detail except to the side where there is higher subdivision.
44 const unsigned short int patch_flag_x = 2;
45 // set if we have subdivided into two halves in the y direction
46 const unsigned short int split_flag_y = 4;
47 // patch misalignment due to higher subdivision a bit lower in the y direction.
48 const unsigned short int patch_flag_y = 8;
49 // We flag vertices in faces that should be two faces over each edge of the
50 // panorama's 180 degree seem, such that by flipping the flagged vertices to
51 // the other side we get a valid face for one side, and by flipping those not
52 // flagged we get the other.
53 const unsigned short int vertex_side_flag_start = 16;
54 // 32, 64 and 128 are for the other vertices.
55 // if the tranformation can't find a result, we set flags to indicate the
56 // vertices we had a problem with.
57 const unsigned short int transform_fail_flag = 256;
58 // 512, 1024, and 2048 are for the other vertices
59 
60 /* where we have a set of flags refering to different vertices, vertex [x][y]
61  * corresponds with first_vertex_flag << (y * 2 + x).
62  */
63 
64 
65 /*************************************
66  * Detail / Acuracy / Speed Tradeoffs
67  *************************************/
68 // Range of acceptable subdivision stopping points:
69 // Depth of forced subdivisions. Increase if there are transformations that show
70 // no signs of increasing detail level until subdivided more times, or cause
71 // problems with the seam detection.
72 const unsigned int min_depth = 4;
73 // Depth at which subdivision stops.
74 // When adjusting also adjust this, also adjust the definition of Tree's nodes
75 // to account for the extra memory required. A high value uses more memory,
76 // a low value limits the detail. The amount of elements in the array should be
77 // the sum from n = 0 to max_depth of 4^n. (so for n=6, 1+4+16+64+256+2048+4096)
78 const unsigned int max_depth = 6;
79 
80 // this is the length in screen pixels under which no subdivision occurs.
81 // higher values are faster, lower ones more accurate.
82 // TODO user preference? Increase during interactive changes?
83 const double min_length = 14.0;
84 // the angle in radians under which no subdivision occurs. Again, higher values
85 // will make it faster, lower ones will give higher accuracy. must be positive.
86 const double min_angle = 0.06;
87 // the distance in absolute screen pixels between twice the length of the
88 // children and the length of the parent nodes, under which no subdivision
89 // occurs. higher values are faster, lower values give higher accuracy. Must be
90 // positive.
91 const double min_length_difference = 3.0;
92 // This is the margin around the edge of the screen in pixels, outside of which
93 // any face is not subdivided. Higher values are less likely to crop small
94 // features from the edge of the display, lower values add speed when there is
95 // faces that are significantly off-screen. It can be any number, but only
96 // positive numbers are recommended.
97 const double offscreen_safety_margin = 10.0;
98 
99 template <class T>
100 inline T sqr(T val)
101 {
102  return val * val;
103 };
104 
107  VisualizationState *visualization_state_in)
108  : MeshRemapper(m_pano_in, image, visualization_state_in)
109 {
110 
111 }
112 
114 {
115  DEBUG_DEBUG("mesh update update reset index");
116  // this sets scale, height, and width.
118  // we need to record the output projection for flipping around 180 degrees.
122  // we want to make a remapped mesh, get the transformation we need:
123 // HuginBase::SrcPanoImage *src = visualization_state->GetSrcImage(image_number);
125  // use the scale to determine edge lengths in pixels for subdivision
126 // DEBUG_INFO("updating mesh for image " << image_number << ", using scale "
127 // << scale << ".\n");
128  // find key points used for +/- 180 degree boundary correction
129  // {x|y}_add_360's are added to a value near the left/top boundary to get
130  // the corresponding point over the right/bottom boundary.
131  // other values are used to check where the boundary is.
133  x_add_360 = info->GetXAdd360();
134  radius = info->GetRadius();
135  y_add_360 = info->GetYAdd360();
136  x_midpoint = info->GetMiddleX();
137  y_midpoint = info->GetMiddleY();
138  lower_bound = info->GetLowerX();
139  upper_bound = info->GetUpperX();
140  lower_bound_h = info->GetLowerY();
141  upper_bound_h = info->GetUpperY();
142  // Find the cropping region of the source image, so we can stick to the
143  // bounding rectangle.
144  SetCrop();
145  // the tree needs to know how to use the croping information for generating
150  // make the tree reflect the transformation
151  RecursiveUpdate(0, 0, 0);
152  // now set up for examining the tree contents.
153  tree.ResetIndex();
154  done_node = true;
155 }
156 
158 {
159 // DEBUG_DEBUG("mesh update get face coords");
160  result->tex_c = tex_coords;
161  result->vertex_c = s_vertex_coords;
162  // if we have some faces left over from a previous clipping operation, give
163  // one of those first:
164  if (GiveClipFaceResult(result)) return true;
165  // when 180 degree bounds correction is working, we'll have two nodes.
166  if (done_node)
167  {
168  do
169  {
170  // this will search the tree for the next leaf node.
172  if (!tree_node_id)
173  {
174  // we've reached last one
175  return false;
176  }
177  }
178  // some of the verticies may have arrived from undefined transformations
179  // if this is one, skip it and try to find another.
180  while ((tree.nodes[tree_node_id].flags & (transform_fail_flag * 15)));
181 
182  // find the coordinates from the tree node
183  result->vertex_c = tree.nodes[tree_node_id].verts;
184  // check that the transformation is properly defined.
185 
187  // if the node has a discontinuity, we want to split it into two
188  // faces and return each one on consecutive calls.
192  {
193  done_node = false;
194  // flip the marked nodes to the other side. copy the coordinates 1st
195  result->vertex_c = s_vertex_coords;
196  for (short unsigned int x = 0; x < 2; x++)
197  {
198  for (short unsigned int y = 0; y < 2; y++)
199  {
200  s_vertex_coords[x][y][0] =
201  tree.nodes[tree_node_id].verts[x][y][0];
202  s_vertex_coords[x][y][1] =
203  tree.nodes[tree_node_id].verts[x][y][1];
204  if (discontinuity_flags & (1 << (x*2 + y)))
205  {
207  }
208  }
209  }
210  }
211  } else {
212  // we flip the other vertices to the ones we did last time.
213  done_node = true;
214  for (short unsigned int x = 0; x < 2; x++)
215  {
216  for (short unsigned int y = 0; y < 2; y++)
217  {
218  s_vertex_coords[x][y][0] =
219  tree.nodes[tree_node_id].verts[x][y][0];
220  s_vertex_coords[x][y][1] =
221  tree.nodes[tree_node_id].verts[x][y][1];
222  if (!(discontinuity_flags & (1 << (x*2 + y))))
223  {
225  }
226  }
227  }
228  }
229  // if we are doing circular cropping, clip the face so it makes the shape
230  if (circle_crop)
231  {
232  // If all points are within the radius, then don't clip
233 // HuginBase::SrcPanoImage *src_img = visualization_state->GetSrcImage(image_number);
234  if ( image->isInside(vigra::Point2D(int(result->tex_c[0][0][0] * width),
235  int(result->tex_c[0][0][1] * height)))
236  && image->isInside(vigra::Point2D(int(result->tex_c[0][1][0] * width),
237  int(result->tex_c[0][1][1] * height)))
238  && image->isInside(vigra::Point2D(int(result->tex_c[1][0][0] * width),
239  int(result->tex_c[1][0][1] * height)))
240  && image->isInside(vigra::Point2D(int(result->tex_c[1][1][0] * width),
241  int(result->tex_c[1][1][1] * height))))
242  {
243  // all inside, doesn't need clipping.
244  return true;
245  }
246  // we do need to clip:
247  ClipFace(result);
248  // if there was anything left, return the first face and leave the rest
249  // for later.
250  if (GiveClipFaceResult(result)) return true;
251  // we clipped to nothing... try and get another face: from the top...
252  return (GetNextFaceCoordinates(result));
253  }
254  return true;
255 }
256 
258 {
259  // we want to flip the given vertex to be the other side of the 180 degree
260  // boundary, whatever the projection format.
261  switch (output_projection)
262  {
264  // There is no 180 degree boundary for rectilinear projections.
265  // Anything containing a vertex beyond 90 degrees (i.e. behind the
267  break;
275  // circular projections. These stretch rather nastily over the
276  // centre, and correcting them doesn't help much, so any image
277  // covering the outer circle is switched to a TexCoordRemapper.
278  break;
285  // flip to the other direction of the other side horizontally.
286  if (vertex_c[0] < x_midpoint) vertex_c[0] += x_add_360;
287  else vertex_c[0] -= x_add_360;
288  break;
291  if (vertex_c[0] < x_midpoint)
292  {
293  vertex_c[0] +=
295  } else {
296  vertex_c[0] -=
298  }
299  break;
301  // flip to the other direction of the other side vertically
302  if (vertex_c[1] < y_midpoint) vertex_c[1] += y_add_360;
303  else vertex_c[1] -= y_add_360;
304  break;
305  default:
306  // all other projection, no special handling
307  break;
308  }
309 }
310 
311 void VertexCoordRemapper::RecursiveUpdate(unsigned int node_num,
312  unsigned int stretch_x, unsigned stretch_y)
313 {
314  // find where we are and what we are mapping
315  // TODO? GetPosition is called by GetInputCoordinates, reuse results?
316  unsigned int x, y, row_size, depth;
317  tree.GetPosition(node_num, x, y, row_size, depth);
319  TreeNode *node = &tree.nodes[node_num],
320  *parent = (depth) ?
321  &tree.nodes[tree.GetParentId(x, y, row_size, depth)] : 0,
322  *left = (x % 2) ?
323  &tree.nodes[tree.GetIndex(x-1, y, row_size, depth)] : 0,
324  *top = (y % 2) ?
325  &tree.nodes[tree.GetIndex(x, y-1, row_size, depth)] : 0;
326  bool valid[2][2];
327  for (unsigned short int i = 0; i < 2; i++)
328  {
329  for (unsigned short int j = 0; j < 2; j++)
330  {
331  if (depth == 0)
332  {
333  // the top level has no parent, so we must calculate all points
334  valid[i][j] =
335  transform.transformImgCoord(node->verts[i][j][0],
336  node->verts[i][j][1],
337  tex_coords[i][j][0] * width,
338  tex_coords[i][j][1] * height);
339  } else
340  // Look up where the point in the tree so far. If this is the first
341  // occurrence of this point, we'll calculate the value.
342  if (i == x %2 && j == y%2 && depth)
343  {
344  // extract a corner from the parent.
345  node->verts[i][j][0] = parent->verts[i][j][0];
346  node->verts[i][j][1] = parent->verts[i][j][1];
347  valid[i][j] = !(parent->flags & (transform_fail_flag << (j*2 +i)));
348  } else if (x % 2 && !i) {
349  // copy from the left
350  node->verts[0][j][0] = left->verts[1][j][0];
351  node->verts[0][j][1] = left->verts[1][j][1];
352  valid[i][j] = !(left->flags & (transform_fail_flag << (j *2 + 1)));
353  } else if (y % 2 && !j) {
354  // copy from the top
355  node->verts[i][0][0] = top->verts[i][1][0];
356  node->verts[i][0][1] = top->verts[i][1][1];
357  valid[i][j] = !(top->flags & (transform_fail_flag << (2 + i)));
358  } else {
359  // We can't find it easily, try a more expensive search.
360  // this will linearly interpolate along the edges where the
361  // subdivision was higher, avoiding gaps when the detail level
362  // was lower above or to the left.
363  if (!tree.GetTransform(x + i * (1 << stretch_x),
364  y + j * (1 << stretch_y),
365  depth, x, y, node->verts[i][j][0],
366  node->verts[i][j][1]))
367  {
368  // We can't find it, so calculate it:
369  valid[i][j] =
370  transform.transformImgCoord(node->verts[i][j][0],
371  node->verts[i][j][1],
372  tex_coords[i][j][0] * width,
373  tex_coords[i][j][1] * height);
374  // If the face results from a patch subdivision (for
375  // aligning a subdivided face, otherwise it need not exist),
376  // use the midpoint of the parent's vertices instead of
377  // calculating a transformed one, so we can use less
378  // subdivision on the other side.
379  // subdivision in x
380  // FIXME there are still gaps. I think there's a logic error
381  if ( depth // not on the top level.
382  // patching in y
383  && (parent->flags & patch_flag_x)
384  // we must be in the middle of the split, the nodes on
385  // the corners of the parent line up anyway.
386  && ((i + x) % 2)
387  // we should be on the side away from the subdivison
388  // (+ve y).
389  && j
390  // If we are alao split in y we can use the middle
391  // node to provide more detail.
392  && (!((parent->flags & split_flag_y) && !(y % 2)))
393  // don't do this if we cross the 180 degree seam
394  && (!(parent->flags & (vertex_side_flag_start * 15))))
395  {
396  node->verts[i][1][0] = (parent->verts[0][1][0]
397  + parent->verts[1][1][0]) / 2.0;
398  node->verts[i][1][1] = (parent->verts[0][1][1]
399  + parent->verts[1][1][1]) / 2.0;
400  }
401  // subdivision in y
402  if ( depth
403  && (parent->flags & patch_flag_y)
404  && ((j + y) % 2)
405  && i
406  && (!((parent->flags & split_flag_x) && !(x % 2)))
407  && (!(parent->flags & (vertex_side_flag_start * 15))))
408  {
409  node->verts[1][j][0] = (parent->verts[1][0][0]
410  + parent->verts[1][1][0]) / 2.0;
411  node->verts[1][j][1] = (parent->verts[1][0][1]
412  + parent->verts[1][1][1]) / 2.0;
413  }
414  } else {
415  // we managed to find it from data already known.
416  valid[i][j] = true;
417  }
418  }
419  }
420  }
421 
422  // now for the recursion
423  // which directions should we divide in?
424  TestSubdivide(node_num);
425  // add the flags for invlaid transormations
426  for (unsigned int i = 0; i < 2; i++)
427  {
428  for (unsigned int j = 0; j < 2; j++)
429  {
430  if (!valid[i][j])
431  {
432  node->flags |= transform_fail_flag << (j * 2 + i);
433  }
434  }
435  }
436  // if the face should be split, now recurse to the child nodes.
437  if (node->flags & (split_flag_x | split_flag_y))
438  {
439  // we will split at least one way.
440  if (!(node->flags & split_flag_x))
441  {
442  // we are not splitting across x, but will will across y.
443  // so the quad will be twice as wide
444  stretch_x++;
445  }
446  else if (!(node->flags & split_flag_y))
447  {
448  stretch_y++;
449  }
450  // find the top left sub-quad
451  x *= 2;
452  y *= 2;
453  row_size *= 2;
454  depth++;
455  // the top left is always generated
456  RecursiveUpdate(tree.GetIndex(x, y, row_size, depth),
457  stretch_x, stretch_y);
458  // split in x
459  if (node->flags & split_flag_x)
460  {
461  RecursiveUpdate(tree.GetIndex(x + 1, y, row_size, depth),
462  stretch_x, stretch_y);
463  }
464  // split in y
465  if (node->flags & split_flag_y)
466  {
467  RecursiveUpdate(tree.GetIndex(x, y + 1, row_size, depth),
468  stretch_x, stretch_y);
469  // if we are splitting in both directions, do the lower right corner
470  if (node->flags & split_flag_x)
471  {
472  RecursiveUpdate(tree.GetIndex(x + 1, y + 1, row_size, depth),
473  stretch_x, stretch_y);
474  }
475  }
476  }
477 }
478 
479 void VertexCoordRemapper::TestSubdivide(unsigned int node_id)
480 {
481  TreeNode *node = &tree.nodes[node_id];
482  unsigned int x, y, row_size, depth;
483  tree.GetPosition(node_id, x, y, row_size, depth);
484  unsigned short int flags = 0;
485  if (depth < min_depth)
486  {
487  // subdivide at least min_depth times
488  // we will need more information for non-trivial children
489  SetLengthAndAngle(node);
490  flags |= split_flag_x | split_flag_y;
491  } else {
492  unsigned int parent_id = tree.GetParentId(node_id);
493  TreeNode *parent = &tree.nodes[parent_id];
494  // minimum length check. We use the length of the top edge to test for
495  // subdivision in x, and the length of the left edge for subdivision in
496  // y.
497  SetLengthAndAngle(node);
498  bool do_not_split_x = node->length_x * scale < min_length,
499  do_not_split_y = node->length_y * scale < min_length;
500  if (depth == max_depth)
501  {
502  // don't subdivide more than max_depth times
503  do_not_split_x = true;
504  do_not_split_y = true;
505  }
506  // if we have stopped splitting in some direction, don't consider
507  // splitting in that direction again.
508  if (!(tree.nodes[parent_id].flags & split_flag_x))
509  {
510  do_not_split_x = true;
511  }
512  else if (!(tree.nodes[parent_id].flags & split_flag_y))
513  {
514  do_not_split_y = true;
515  }
516  // if it has only subdivided to patch up between two subdivision levels,
517  // don't consider subdividing for adding more detail.
518  if (tree.nodes[parent_id].flags & patch_flag_x)
519  {
520  do_not_split_x = true;
521  }
522  else if (tree.nodes[parent_id].flags & patch_flag_y)
523  {
524  do_not_split_y = true;
525  }
526 
527  // If the angles have deviated too much from the parent then we should
528  // add more detail, however if it is fairly flat then we don't need to.
529  // It is possible for the angles to remain constant but the length
530  // of the lines to change dramatically, so we check for a big difference
531  // between the length of the parent node and twice the length of child.
532  // if the child does not change the length much and the angle is small,
533  // then we have enough detail, and we don't split.
534  float ang_x = node->angle_x - parent->angle_x;
535  if (ang_x < 0) ang_x = -ang_x;
536  if (ang_x > M_PI) ang_x = 2 * M_PI - ang_x;
537  float length_difference_x
538  = (parent->length_x - (2.0 * node->length_x)) * scale;
539  if (length_difference_x < 0.0)
540  {
541  length_difference_x = -length_difference_x;
542  }
543  if (ang_x < min_angle && length_difference_x < min_length_difference)
544  {
545  do_not_split_x = true;
546  }
547  float ang_y = node->angle_y - parent->angle_y;
548  if (ang_y < 0) ang_y = -ang_y;
549  if (ang_y > M_PI) ang_y = 2 * M_PI - ang_y;
550  float length_difference_y
551  = (parent->length_y - (2.0 * node->length_y)) * scale;
552  if (length_difference_y < 0.0)
553  {
554  length_difference_y = -length_difference_y;
555  }
556  if (ang_y < min_angle && length_difference_y < min_length_difference)
557  {
558  do_not_split_y = true;
559  }
560  // if the face is entirely off the screen, we should not subdivide it.
561  // get the screen area
562  vigra::Rect2D viewport = visualization_state->GetVisibleArea();
563  // add a margin for safety, we don't want to clip too much stuff that
564  // curls back on to the screen. We add 2 as we need some space around
565  // very small panoramas that have enlarged to fit the preview window,
566  // and even with a fairly large margin rounding to int may lose the
567  // border completely.
568  viewport.addBorder((int) (2.0 + offscreen_safety_margin * scale));
569  bool all_left = true, all_right = true,
570  all_above = true, all_bellow = true;
571  for (unsigned int ix = 0; ix < 2; ix++)
572  {
573  for (unsigned int iy = 0; iy < 2; iy++)
574  {
575  if (node->verts[ix][iy][0] > viewport.left())
576  all_left = false;
577  if (node->verts[ix][iy][0] < viewport.right())
578  all_right = false;
579  if (node->verts[ix][iy][1] > viewport.top())
580  all_above = false;
581  if (node->verts[ix][iy][1] < viewport.bottom())
582  all_bellow = false;
583  }
584  }
585  if (all_left || all_right || all_bellow || all_above)
586  {
587  // all vertices are off a side of the screen. This catches most
588  // cases where the face is off the screen. Don't allow subdivisions:
589  do_not_split_x = true;
590  do_not_split_y = true;
591  }
592  if (!do_not_split_x) flags |= split_flag_x;
593  if (!do_not_split_y) flags |= split_flag_y;
594  }
595  /* Flag the vertices on a different side of the +/-180 degree seam.
596  * We don't want to flag any vertices if the face covers a continuous
597  * area of the transformation.
598  */
599  // We don't need to mark the first few subdivisions, but this is necessary
600  // when patch subdivisions become possible.
601  if (depth >= min_depth)
602  {
603  // determine if it is likely to be non-continuous.
604  // this needs to be false for leaf-node faces in the centre that span
605  // across the '0 degree' point, and true for faces that span the +/-180
606  // degree split. It doesn't really matter what it is set to otherwise.
607  bool noncontinuous = false;
609  switch (output_projection)
610  {
612  // we don't need faces to cross from one side to another. Faces
613  // all / partially 'behind the viewer' are skipped because the
614  // vertices behind the viewer are marked.
615  break;
620  // A mesh covering the extremities of a disk projection should
621  // be using a TexCoordRemapper instead, otherwise, a point
622  // mapping to the border will be stretched across the disk.
623  break;
631  // Cylinderical-like projections have the seam across the left
632  // and right edge. We'll take any face within the middle third
633  // to be continuous, the rest possibly noncontinuous.
634  if ( node->verts[0][0][0] < lower_bound
635  || node->verts[0][0][0] > upper_bound
636  || node->verts[1][0][0] < lower_bound
637  || node->verts[1][0][0] > upper_bound
638  || node->verts[0][1][0] < lower_bound
639  || node->verts[0][1][0] > upper_bound
640  || node->verts[1][1][0] < lower_bound
641  || node->verts[1][1][0] > upper_bound)
642  {
643  noncontinuous = true;
644  // flag nodes on the right hand side.
645  if (node->verts[0][0][0] > x_midpoint)
646  {
647  flags |= vertex_side_flag_start;
648  }
649  if (node->verts[0][1][0] > x_midpoint)
650  {
651  flags |= vertex_side_flag_start * 2;
652  }
653  if (node->verts[1][0][0] > x_midpoint)
654  {
655  flags |= vertex_side_flag_start * 4;
656  }
657  if (node->verts[1][1][0] > x_midpoint)
658  {
659  flags |= vertex_side_flag_start * 8;
660  }
661  }
662  break;
664  // like above, but the bounds change with height
665  if ( node->verts[0][0][0] < i->GetLowerX(node->verts[0][0][1])
666  || node->verts[0][0][0] > i->GetUpperX(node->verts[0][0][1])
667  || node->verts[1][0][0] < i->GetLowerX(node->verts[1][0][1])
668  || node->verts[1][0][0] > i->GetUpperX(node->verts[1][0][1])
669  || node->verts[0][1][0] < i->GetLowerX(node->verts[0][1][1])
670  || node->verts[0][1][0] > i->GetUpperX(node->verts[0][1][1])
671  || node->verts[1][1][0] < i->GetLowerX(node->verts[1][1][1])
672  || node->verts[1][1][0] > i->GetUpperX(node->verts[1][1][1])
673  )
674  {
675  noncontinuous = true;
676  // flag nodes on the right hand side.
677  if (node->verts[0][0][0] > x_midpoint)
678  {
679  flags |= vertex_side_flag_start;
680  }
681  if (node->verts[0][1][0] > x_midpoint)
682  {
683  flags |= vertex_side_flag_start * 2;
684  }
685  if (node->verts[1][0][0] > x_midpoint)
686  {
687  flags |= vertex_side_flag_start * 4;
688  }
689  if (node->verts[1][1][0] > x_midpoint)
690  {
691  flags |= vertex_side_flag_start * 8;
692  }
693  }
694  break;
696  // like the cylindrical ones, but vertically.
697  if ( node->verts[0][0][1] < lower_bound_h
698  || node->verts[0][0][1] > upper_bound_h
699  || node->verts[1][0][1] < lower_bound_h
700  || node->verts[1][0][1] > upper_bound_h
701  || node->verts[0][1][1] < lower_bound_h
702  || node->verts[0][1][1] > upper_bound_h
703  || node->verts[1][1][1] < lower_bound_h
704  || node->verts[1][1][1] > upper_bound_h)
705  {
706  noncontinuous = true;
707  // flag nodes on the bottom side.
708  if (node->verts[0][0][1] > y_midpoint)
709  {
710  flags |= vertex_side_flag_start;
711  }
712  if (node->verts[0][1][1] > y_midpoint)
713  {
714  flags |= vertex_side_flag_start * 2;
715  }
716  if (node->verts[1][0][1] > y_midpoint)
717  {
718  flags |= vertex_side_flag_start * 4;
719  }
720  if (node->verts[1][1][1] > y_midpoint)
721  {
722  flags |= vertex_side_flag_start * 8;
723  }
724  }
725  break;
727  break;
728  default:
729  // all other projections, no special handling
730  break;
731  }
732  if (noncontinuous)
733  {
734  // If all the side flags are set, we only have one face on one side:
735  // Remove all of those flags and we have the same face, but we can
736  // use the presence of any flags to detect when we have two.
737  if ((flags / vertex_side_flag_start) % 16 == 15)
738  {
739  flags &= ~(vertex_side_flag_start * 15);
740  }
741  }
742  }
743 
744  node->flags = flags;
745  // check if the faces to the left are subdivided at a higher level
746  if (x && !(flags & split_flag_y) && (depth < max_depth))
747  {
748  // get the potentially more subdivided version
749  double dest_x, dest_y;
750  unsigned int subdiv_node;
751  subdiv_node = tree.GetTransform(x * 2, y * 2 + 1,
752  depth * 2,
753  x * 2, y * 2, dest_x, dest_y);
754  if (subdiv_node > node_id)
755  {
756  // we should have a more subdivided node on the left.
757  // mark it for patching up to line up with it.
758  node->flags |= split_flag_y | patch_flag_y;
759  }
760  }
761  if (y && !(flags & split_flag_x) && (depth < max_depth))
762  {
763  // get the potentially more subdivided version
764  double dest_x, dest_y;
765  unsigned int subdiv_node;
766  subdiv_node = tree.GetTransform(x * 2 + 1, y * 2,
767  depth * 2,
768  x * 2, y * 2, dest_x, dest_y);
769  if (subdiv_node > node_id)
770  {
771  // we should have a more subdivided node on the left.
772  // mark it for patching up to line up with it.
773  node->flags |= split_flag_x | patch_flag_x;
774  }
775  }
776 }
777 
779 {
780  float xdx = node->verts[0][0][0] - node->verts[1][0][0],
781  xdy = node->verts[0][0][1] - node->verts[1][0][1],
782  ydx = node->verts[0][0][0] - node->verts[0][1][0],
783  ydy = node->verts[0][0][1] - node->verts[0][1][1];
784  // this is the length of the edge with y = 0
785  node->length_x = sqrt(sqr(xdx) + sqr(xdy));
786  // this is the length of the edge with x = 0
787  node->length_y = sqrt(sqr(ydx) + sqr(ydy));
788  // find the angle of the top and left edge.
789  node->angle_x = atan2(xdy, xdx);
790  node->angle_y = atan2(ydy, ydx);
791 }
792 
793 unsigned int VertexCoordRemapper::Tree::GetParentId(const unsigned int nodenum)
794 {
795  // get the parent id of a node specifed by its index.
796  unsigned int x, y, row_size, depth;
797  GetPosition(nodenum, x, y, row_size, depth);
798  return GetParentId(x, y, row_size, depth);
799 }
800 
801 unsigned int VertexCoordRemapper::Tree::GetParentId(unsigned int x,
802  unsigned int y,
803  unsigned int row_size,
804  unsigned int depth)
805 {
806  // get the parent id of a node specified by an exact location.
807  // the parent is the one that takes the group of four elements around this
808  // one in the above depth. All the dimensions are halved.
809  x /= 2;
810  y /= 2;
811  row_size /= 2;
812  depth--;
813  return GetIndex(x, y, row_size, depth);
814 }
815 
816 unsigned int VertexCoordRemapper::Tree::GetDepth(const unsigned int node_num)
817 {
818  unsigned int depth = 0, count = 0;
819  while (node_num > count)
820  {
821  depth++;
822  count += 1 << (depth * 2);
823  }
824  return depth;
825 }
826 
827 void VertexCoordRemapper::Tree::GetPosition(const unsigned int node_num,
828  unsigned int &x, unsigned int &y,
829  unsigned int &row_size,
830  unsigned int &depth)
831 {
832  depth = GetDepth(node_num);
833  row_size = 1 << depth;
834  unsigned int sub = 0;
835  if (depth)
836  {
837  for (unsigned int d = 0; d < depth; d++)
838  {
839  sub += (1 << (d * 2));
840  }
841  }
842  unsigned int position_id = node_num - sub;
843  x = position_id % row_size;
844  y = position_id / row_size;
845 }
846 
847 unsigned int VertexCoordRemapper::Tree::GetIndex(const unsigned int x,
848  const unsigned int y,
849  const unsigned int row_size,
850  unsigned int depth)
851 {
852  unsigned int add = 0;
853  while (depth)
854  {
855  depth--;
856  add += 1 << (depth * 2);
857  }
858  return add + x + y * row_size;
859 }
860 
861 void VertexCoordRemapper::Tree::GetInputCoordinates(const unsigned int node_num,
862  double coords[2][2][2])
863 {
864  // find the coordinates of each point at the vertices in the original image.
865  // this is used to look up the remapped coordinates and provide texture
866  // coordinates.
867  // it halves in size several times depending on depth.
868  unsigned int x, y, row_size, depth;
869  GetPosition(node_num, x, y, row_size, depth);
870  // now we can locate the upper right corner
871  double row_spacing = 1.0 / (double) row_size;
872  coords[0][0][0] = row_spacing * (double) x;
873  coords[0][0][1] = row_spacing * (double) y;
874  // the extent is dependent on the direction of the subdivisions.
875  // at least one direction has an extent of row_spacing. It can double up
876  // if the parent did not subdivide in a direction, recursively up the tree.
877  bool scale_x = false;
878  double opp = 1.0;
879  if (node_num != 0)
880  {
881  unsigned int parent_id = GetParentId(x, y, row_size, depth);
882  if (!(nodes[parent_id].flags & split_flag_x))
883  {
884  while (!(nodes[parent_id].flags & split_flag_x))
885  {
886  parent_id = GetParentId(parent_id);
887  opp *= 2.0;
888  }
889  scale_x = true;
890  } else {
891  while (!(nodes[parent_id].flags & split_flag_y))
892  {
893  parent_id = GetParentId(parent_id);
894  opp *= 2.0;
895  }
896  }
897  }
898  opp *= row_spacing;
899  coords[1][0][0] = coords[0][0][0] + (scale_x ? opp : row_spacing);
900  coords[1][0][1] = coords[0][0][1];
901  coords[0][1][0] = coords[0][0][0];
902  coords[0][1][1] = coords[0][0][1] + (scale_x ? row_spacing : opp);
903  coords[1][1][0] = coords[1][0][0];
904  coords[1][1][1] = coords[0][1][1];
905  // now we transform the results so that we map to the cropped region instead
906  // of the whole image.
907  for (unsigned int i = 0; i < 2; i++)
908  {
909  for (unsigned int j = 0; j < 2; j++)
910  {
911  coords[i][j][0] = coords[i][j][0] * x_crop_scale + x_crop_offs;
912  coords[i][j][1] = coords[i][j][1] * y_crop_scale + y_crop_offs;
913  }
914  }
915 }
916 
918 {
919  cur_tree_node = 0;
920 }
921 
923 {
924  // depth first search for leaf nodes with cur_tree_node
925  unsigned int x, y, row_size, depth;
926  GetPosition(cur_tree_node, x, y, row_size, depth);
927  // we want to backtrack until we find an unexplored sub branch. At this
928  // point, we trace down the tree until we get to a leaf.
929  if (cur_tree_node != 0)
930  {
931  bool done = false;
932  while (!done && cur_tree_node != 0)
933  {
934  unsigned int xd = x % 2;
935  unsigned int yd = y % 2;
936  x /= 2;
937  y /= 2;
938  row_size /= 2;
939  depth--;
940  cur_tree_node = GetIndex(x, y, row_size, depth);
941  if (cur_tree_node == 0 && xd == 1 && yd == 1)
942  {
943  // we've expanded all the options and got back to the top
944  return 0; // no more faces
945  }
946  // where does this split?
947  bool sx = ((nodes[cur_tree_node].flags & split_flag_x) != 0);
948  bool sy = ((nodes[cur_tree_node].flags & split_flag_y) != 0);
949  // have we used all the split options?
950  if (!(((sx && xd) || !sx) && ((sy && yd) || !sy)))
951  {
952  // we've found an unexpanded child, take the next one:
953  done = true;
954  x *= 2;
955  y *= 2;
956  row_size *= 2;
957  depth++;
958  if (sx && !xd && !yd) {
959  x++;
960  } else if ((sx && xd && sy && !yd) || (!sx && sy && !yd)) {
961  y++;
962  } else if (sx && sy && !xd && yd) {
963  x++; y++;
964  }
965  cur_tree_node = GetIndex(x, y, row_size, depth);
966  // if we've made our way to the root node, we've run out out
967  // of options.
968  if (cur_tree_node == 0)
969  {
970  return 0;
971  }
972  };
973  }
974  }
975  // find the first leaf on this subtree, taking top left every time.
976  while (nodes[cur_tree_node].flags & (split_flag_x | split_flag_y))
977  {
978  x *= 2;
979  y *= 2;
980  row_size *= 2;
981  depth++;
982  cur_tree_node = GetIndex(x, y, row_size, depth);
983  }
984  return cur_tree_node;
985 }
986 
987 unsigned int VertexCoordRemapper::Tree::GetTransform(unsigned int src_x, unsigned int src_y,
988  unsigned int max_depth,
989  unsigned int stop_x, unsigned int stop_y,
990  double &dest_x, double &dest_y)
991 {
992  // we start at the top and take children, prefering the lowest numbered one,
993  // until we either have no children, or we have found a child that knows the
994  // exact place we are looking for.
995  unsigned int no_x = 0, no_y = 0, row_size = 1, depth = 0,
996  rem_x = src_x, rem_y = src_y,
997  step_x = 1 << max_depth, step_y = 1 << max_depth,
998  node_id = GetIndex(no_x, no_y, row_size, depth);
999  while ( !( (rem_x == 0 || rem_x == step_x)
1000  && (rem_y == 0 || rem_y == step_y))
1001  && (nodes[node_id].flags & (split_flag_x | split_flag_y)))
1002  {
1003  // split, try and take earliest child node that leads towards goal
1004  depth++;
1005  no_x *= 2;
1006  no_y *= 2;
1007  row_size *= 2;
1008  if (nodes[node_id].flags & split_flag_x)
1009  {
1010  step_x /= 2;
1011  if (rem_x > step_x)
1012  {
1013  no_x ++;
1014  rem_x -= step_x;
1015  }
1016  }
1017  if (nodes[node_id].flags & split_flag_y)
1018  {
1019  step_y /= 2;
1020  if (rem_y > step_y)
1021  {
1022  no_y ++;
1023  rem_y -= step_y;
1024  }
1025  }
1026  // we want to stop if we have got a node at least as high as the
1027  // requested stopping node. Anything at a higher depth is after it.
1028  if (depth > max_depth) return 0;
1029  node_id = GetIndex(no_x, no_y, row_size, depth);
1030  if (depth == max_depth && no_x >= stop_x && no_y >= stop_y) return 0;
1031  }
1032  // if any of the vertices we are want to use are invalid (e.g. behind the
1033  // viewer in a rectilinear projection) we don't want to risk using them:
1034  if (nodes[node_id].flags & (transform_fail_flag * 15)) return 0;
1035  // if this face crosses a discontinuity, we should be using a point off
1036  // screen instead of in the middle. Refuse to use these faces
1037  if (nodes[node_id].flags & (vertex_side_flag_start * 15)) return 0;
1038  // linearly interpolate the node's corners.
1039  // most of the time we only use factors of 0 and 1, we don't want to make
1040  // points up except when trying to connect a point on a highly subdivided
1041  // face to a point that doesn't exist because it is on a less subdivided
1042  // face, in which case it needs to line up with the linear interpolation of
1043  // the less subdivided face. This is along one edge, so the other direction
1044  // should have a blending factor of 0 or 1.
1045  double xf = (double) rem_x / (double) step_x;
1046  double yf = (double) rem_y / (double) step_y;
1047 
1048  double top_x = (1.0 - xf) * nodes[node_id].verts[0][0][0]
1049  + xf * nodes[node_id].verts[1][0][0],
1050  bottom_x = (1-.0 - xf) * nodes[node_id].verts[0][1][0]
1051  + xf * nodes[node_id].verts[1][1][0],
1052  top_y = (1.0 - xf) * nodes[node_id].verts[0][0][1]
1053  + xf * nodes[node_id].verts[1][0][1],
1054  bottom_y = (1.0 - xf) * nodes[node_id].verts[0][1][1]
1055  + xf * nodes[node_id].verts[1][1][1];
1056  dest_x = top_x * (1.0 - yf) + bottom_x * yf;
1057  dest_y = top_y * (1.0 - yf) + bottom_y* yf;
1058  return node_id;
1059 }
1060 
VertexCoordRemapper(HuginBase::Panorama *m_pano, HuginBase::SrcPanoImage *image, VisualizationState *visualization_state)
const double min_length
PanoramaOptions::ProjectionFormat getProjection() const
bool GiveClipFaceResult(Coords *result)
Get a face that was produced by ClipFace.
unsigned int GetDepth(const unsigned int nodenum)
An abstract base class for objects that calculate an approximate remap specified by quadrilatrials...
Definition: MeshRemapper.h:43
const double GetUpperY() const
double crop_x1
Definition: MeshRemapper.h:100
double s_vertex_coords[2][2][2]
double crop_x2
Definition: MeshRemapper.h:100
unsigned int GetIndex(const unsigned int x, const unsigned int y, const unsigned int row_size, unsigned int depth)
double tex_coords[2][2][2]
void SetLengthAndAngle(TreeNode *node)
A class for exchanging pointers to coordinates.
Definition: MeshRemapper.h:60
const unsigned short int split_flag_y
unsigned int getHeight() const
get panorama height
HuginBase::SrcPanoImage * image
Definition: MeshRemapper.h:85
const double GetRadius() const
virtual void UpdateAndResetIndex()
const double GetUpperX() const
void DiscontinuityFlip(double vertex_c[2])
const double min_angle
include file for the hugin project
double(* tex_c)[2][2]
The coordinate in the source image ranging from 0 to 1.
Definition: MeshRemapper.h:64
void RecursiveUpdate(unsigned int node_id, unsigned int stretch_x, unsigned int stretch_y)
bool isInside(vigra::Point2D p, bool ignoreMasks=false) const
check if a coordinate is inside the source image
HuginBase::PanoramaOptions::ProjectionFormat output_projection
const double GetLowerY() const
VisualizationState * visualization_state
Definition: MeshRemapper.h:83
const unsigned short int split_flag_x
const double min_length_difference
const unsigned short int patch_flag_x
Model for a panorama.
Definition: Panorama.h:152
TreeNode nodes[1+4+16+64+256+1024+4096]
void ClipFace(Coords *face)
Crop a face to the source image, return true if there is anything left.
const double GetMiddleY() const
void SetCrop()
Fill the crop values of the MeshRemapper from the source image.
vigra::Rect2D GetVisibleArea()
Definition: ViewState.h:222
void createInvTransform(const vigra::Diff2D &srcSize, VariableMap srcVars, Lens::LensProjectionFormat srcProj, const vigra::Diff2D &destSize, PanoramaOptions::ProjectionFormat destProj, const std::vector< double > &destProjParam, double destHFOV, const vigra::Diff2D &origSrcSize)
create image-&gt;pano transformation
double height
The sizes of the input images in pixels.
Definition: MeshRemapper.h:94
double crop_y1
Definition: MeshRemapper.h:100
unsigned int GetParentId(const unsigned int nodenum)
void GetInputCoordinates(unsigned int node_num, double coords[2][2][2])
float scale
The number number of units between vertex coorinates that gives a pixel in the display.
Definition: MeshRemapper.h:89
const unsigned short int vertex_side_flag_start
const double GetMiddleX() const
T sqr(T val)
virtual OutputProjectionInfo * GetProjectionInfo()
Definition: ViewState.cpp:473
void GetPosition(const unsigned int nodenum, unsigned int &x, unsigned int &y, unsigned int &row_size, unsigned int &depth)
unsigned int GetTransform(unsigned int src_x, unsigned int src_y, unsigned int depth, unsigned int stop_x, unsigned int stop_y, double &dest_x, double &dest_y)
const double GetXAdd360() const
const double GetLowerX() const
HuginBase::PTools::Transform transform
A transform to use to remap the images.
Definition: MeshRemapper.h:96
unsigned int getWidth() const
include file for the hugin project
bool transformImgCoord(double &x_dest, double &y_dest, double x_src, double y_src) const
like transform, but return image coordinates, not cartesian coordinates
const unsigned int min_depth
const unsigned short int transform_fail_flag
const unsigned short int patch_flag_y
virtual HuginBase::PanoramaOptions * GetOptions()
Definition: ViewState.cpp:468
const unsigned int max_depth
double width
Definition: MeshRemapper.h:94
#define DEBUG_DEBUG(msg)
Definition: utils.h:68
virtual bool GetNextFaceCoordinates(Coords *result)
Get the texture and vertex coordinates for the next face.
double crop_y2
Definition: MeshRemapper.h:100
virtual void UpdateAndResetIndex()
static void info(const char *fmt,...)
Definition: svm.cpp:95
unsigned short int discontinuity_flags
const double offscreen_safety_margin
static uint16_t flags
double(* vertex_c)[2][2]
The coordinate in the panorama, in its pixel space.
Definition: MeshRemapper.h:66
All variables of a source image.
Definition: SrcPanoImage.h:194
#define M_PI
Definition: GaborFilter.cpp:34
void TestSubdivide(unsigned int node_id)
const double GetYAdd360() const