VertexCoordRemapper.cpp

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
00002 
00023 #ifdef __WXMAC__
00024 #include "panoinc_WX.h"
00025 #include "panoinc.h"
00026 #endif
00027 
00028 #include "VertexCoordRemapper.h"
00029 #include <math.h>
00030 #include <iostream>
00031 #include "ViewState.h"
00032 
00033 #include "panodata/SrcPanoImage.h"
00034 
00035 /*******************
00036  * Face's bit flags 
00037  *******************/
00038 // split_x is set if the node has been split into two children (subdivided) in x
00039 const unsigned short int split_flag_x = 1;
00040 // patch_flag_x is set if we are patching a hole / overlap that could be caused
00041 // by the subdivision in the previous x location being more dense. It is set
00042 // instead of split_flag_x, but we subdivide into faces that do not provide more
00043 // detail except to the side where there is higher subdivision. 
00044 const unsigned short int patch_flag_x = 2;
00045 // set if we have subdivided into two halves in the y direction
00046 const unsigned short int split_flag_y = 4;
00047 // patch misalignment due to higher subdivision a bit lower in the y direction.
00048 const unsigned short int patch_flag_y = 8;
00049 // We flag vertices in faces that should be two faces over each edge of the
00050 // panorama's 180 degree seem, such that by flipping the flagged vertices to
00051 // the other side we get a valid face for one side, and by flipping those not
00052 // flagged we get the other.
00053 const unsigned short int vertex_side_flag_start = 16;
00054 // 32, 64 and 128 are for the other vertices.
00055 // if the tranformation can't find a result, we set flags to indicate the
00056 // vertices we had a problem with.
00057 const unsigned short int transform_fail_flag = 256;
00058 // 512, 1024, and 2048 are for the other vertices
00059 
00060 /* where we have a set of flags refering to different vertices, vertex [x][y]
00061  * corresponds with first_vertex_flag << (y * 2 + x).
00062  */
00063 
00064 
00065 /*************************************
00066  * Detail / Acuracy / Speed Tradeoffs
00067  *************************************/
00068 // Range of acceptable subdivision stopping points:
00069 // Depth of forced subdivisions. Increase if there are transformations that show
00070 // no signs of increasing detail level until subdivided more times, or cause
00071 // problems with the seam detection.
00072 const unsigned int min_depth = 4;
00073 // Depth at which subdivision stops.
00074 // When adjusting also adjust this, also adjust the definition of Tree's nodes
00075 // to account for the extra memory required. A high value uses more memory,
00076 // a low value limits the detail. The amount of elements in the array should be
00077 // the sum from n = 0 to max_depth of 4^n. (so for n=6, 1+4+16+64+256+2048+4096)
00078 const unsigned int max_depth = 6;
00079 
00080 // this is the length in screen pixels under which no subdivision occurs.
00081 // higher values are faster, lower ones more accurate.
00082 // TODO user preference? Increase during interactive changes?
00083 const double min_length = 14.0;
00084 // the angle in radians under which no subdivision occurs. Again, higher values
00085 // will make it faster, lower ones will give higher accuracy. must be positive.
00086 const double min_angle = 0.06;
00087 // the distance in absolute screen pixels between twice the length of the
00088 // children and the length of the parent nodes, under which no subdivision
00089 // occurs. higher values are faster, lower values give higher accuracy. Must be
00090 // positive.
00091 const double min_length_difference = 3.0;
00092 // This is  the margin around the edge of the screen in pixels, outside of which
00093 // any face is not subdivided. Higher values are less likely to crop small
00094 // features from the edge of the display, lower values add speed when there is
00095 // faces that are significantly off-screen. It can be any number, but only
00096 // positive numbers are recommended.
00097 const double offscreen_safety_margin = 10.0;
00098 
00099 template <class T>
00100 inline T sqr(T val)
00101 {
00102     return val * val;
00103 };
00104 
00105 VertexCoordRemapper::VertexCoordRemapper(HuginBase::Panorama *m_pano_in,
00106                                          HuginBase::SrcPanoImage *image,
00107                                          VisualizationState *visualization_state_in)
00108     : MeshRemapper(m_pano_in, image, visualization_state_in)
00109 {
00110     
00111 }
00112 
00113 void VertexCoordRemapper::UpdateAndResetIndex()
00114 {
00115     DEBUG_DEBUG("mesh update update reset index");
00116     // this sets scale, height, and width.
00117     MeshRemapper::UpdateAndResetIndex();
00118     // we need to record the output projection for flipping around 180 degrees.
00119     output_projection = visualization_state->GetOptions()->getProjection();
00120     o_width = visualization_state->GetOptions()->getWidth();
00121     o_height = visualization_state->GetOptions()->getHeight();
00122     // we want to make a remapped mesh, get the transformation we need:
00123 //    HuginBase::SrcPanoImage *src = visualization_state->GetSrcImage(image_number);
00124     transform.createInvTransform(*image, *(visualization_state->GetOptions()));
00125     // use the scale to determine edge lengths in pixels for subdivision
00126 //    DEBUG_INFO("updating mesh for image " << image_number << ", using scale "
00127 //              << scale << ".\n");
00128     // find key points used for +/- 180 degree boundary correction
00129     // {x|y}_add_360's are added to a value near the left/top boundary to get
00130     // the corresponding point over the right/bottom boundary.
00131     // other values are used to check where the boundary is.
00132     OutputProjectionInfo *info = visualization_state->GetProjectionInfo();
00133     x_add_360 = info->GetXAdd360();
00134     radius = info->GetRadius();
00135     y_add_360 = info->GetYAdd360();
00136     x_midpoint = info->GetMiddleX();
00137     y_midpoint = info->GetMiddleY();
00138     lower_bound = info->GetLowerX();
00139     upper_bound = info->GetUpperX();
00140     lower_bound_h = info->GetLowerY();
00141     upper_bound_h = info->GetUpperY();
00142     // Find the cropping region of the source image, so we can stick to the
00143     // bounding rectangle.
00144     SetCrop();
00145     // the tree needs to know how to use the croping information for generating
00146     tree.x_crop_scale = crop_x2 - crop_x1;
00147     tree.x_crop_offs = crop_x1;
00148     tree.y_crop_scale = crop_y2 - crop_y1;
00149     tree.y_crop_offs = crop_y1;
00150     // make the tree reflect the transformation
00151     RecursiveUpdate(0, 0, 0);
00152     // now set up for examining the tree contents.
00153     tree.ResetIndex();
00154     done_node = true;
00155 }
00156 
00157 bool VertexCoordRemapper::GetNextFaceCoordinates(Coords *result)
00158 {
00159 //    DEBUG_DEBUG("mesh update get face coords");
00160     result->tex_c = tex_coords;
00161     result->vertex_c = s_vertex_coords;
00162     // if we have some faces left over from a previous clipping operation, give
00163     // one of those first:
00164     if (GiveClipFaceResult(result)) return true;
00165     // when 180 degree bounds correction is working, we'll have two nodes.
00166     if (done_node)
00167     {
00168         do
00169         {
00170             // this will search the tree for the next leaf node.
00171             tree_node_id = tree.GetNext();
00172             if (!tree_node_id)
00173             {
00174                 // we've reached last one
00175                 return false;
00176             }
00177         }
00178         // some of the verticies may have arrived from undefined transformations
00179         // if this is one, skip it and try to find another.
00180         while ((tree.nodes[tree_node_id].flags & (transform_fail_flag * 15)));
00181             
00182         // find the coordinates from the tree node
00183         result->vertex_c = tree.nodes[tree_node_id].verts;
00184         // check that the transformation is properly defined.
00185 
00186         tree.GetInputCoordinates(tree_node_id, tex_coords);
00187         // if the node has a discontinuity, we want to split it into two
00188         // faces and return each one on consecutive calls.
00189         discontinuity_flags =
00190                  (tree.nodes[tree_node_id].flags / vertex_side_flag_start) % 16;
00191         if (discontinuity_flags)
00192         {
00193             done_node = false;
00194             // flip the marked nodes to the other side. copy the coordinates 1st
00195             result->vertex_c = s_vertex_coords;            
00196             for (short unsigned int x = 0; x < 2; x++)
00197             {
00198                 for (short unsigned int y = 0; y < 2; y++)
00199                 {
00200                       s_vertex_coords[x][y][0] =
00201                                         tree.nodes[tree_node_id].verts[x][y][0];
00202                       s_vertex_coords[x][y][1] =
00203                                         tree.nodes[tree_node_id].verts[x][y][1];
00204                       if (discontinuity_flags & (1 << (x*2 + y)))
00205                       {
00206                           DiscontinuityFlip(s_vertex_coords[x][y]);
00207                       }
00208                 }
00209             }
00210         }
00211     } else {
00212         // we flip the other vertices to the ones we did last time.
00213         done_node = true;
00214         for (short unsigned int x = 0; x < 2; x++)
00215         {
00216             for (short unsigned int y = 0; y < 2; y++)
00217             {
00218                   s_vertex_coords[x][y][0] =
00219                                         tree.nodes[tree_node_id].verts[x][y][0];
00220                   s_vertex_coords[x][y][1] =
00221                                         tree.nodes[tree_node_id].verts[x][y][1];
00222                   if (!(discontinuity_flags & (1 << (x*2 + y))))
00223                   {
00224                       DiscontinuityFlip(s_vertex_coords[x][y]);
00225                   }
00226             }
00227         }
00228     }
00229     // if we are doing circular cropping, clip the face so it makes the shape
00230     if (circle_crop)
00231     {
00232         // If all points are within the radius, then don't clip
00233 //        HuginBase::SrcPanoImage *src_img = visualization_state->GetSrcImage(image_number);
00234         if (   image->isInside(vigra::Point2D(int(result->tex_c[0][0][0] * width),
00235                                                 int(result->tex_c[0][0][1] * height)))
00236             && image->isInside(vigra::Point2D(int(result->tex_c[0][1][0] * width),
00237                                                 int(result->tex_c[0][1][1] * height)))
00238             && image->isInside(vigra::Point2D(int(result->tex_c[1][0][0] * width),
00239                                                 int(result->tex_c[1][0][1] * height)))
00240             && image->isInside(vigra::Point2D(int(result->tex_c[1][1][0] * width),
00241                                                 int(result->tex_c[1][1][1] * height))))
00242         {
00243             // all inside, doesn't need clipping.
00244             return true;
00245         }
00246         // we do need to clip:     
00247         ClipFace(result);
00248         // if there was anything left, return the first face and leave the rest
00249         // for later.
00250         if (GiveClipFaceResult(result)) return true;
00251         // we clipped to nothing... try and get another face: from the top...
00252         return (GetNextFaceCoordinates(result));
00253     }
00254     return true;
00255 }
00256 
00257 void VertexCoordRemapper::DiscontinuityFlip(double vertex_c[2])
00258 {
00259     // we want to flip the given vertex to be the other side of the 180 degree
00260     // boundary, whatever the projection format.
00261     switch (output_projection)
00262     {
00263         case HuginBase::PanoramaOptions::RECTILINEAR:
00264             // There is no 180 degree boundary for rectilinear projections.
00265             // Anything containing a vertex beyond 90 degrees (i.e. behind the
00267             break;
00268         case HuginBase::PanoramaOptions::FULL_FRAME_FISHEYE:
00269         case HuginBase::PanoramaOptions::STEREOGRAPHIC:
00270         case HuginBase::PanoramaOptions::LAMBERT_AZIMUTHAL:
00271         case HuginBase::PanoramaOptions::HAMMER_AITOFF:
00272         case HuginBase::PanoramaOptions::ORTHOGRAPHIC:
00273         case HuginBase::PanoramaOptions::EQUISOLID:
00274         case HuginBase::PanoramaOptions::THOBY_PROJECTION:
00275             // circular projections. These stretch rather nastily over the
00276             // centre, and correcting them doesn't help much, so any image
00277             // covering the outer circle is switched to a TexCoordRemapper.
00278             break;
00279         case HuginBase::PanoramaOptions::CYLINDRICAL:
00280         case HuginBase::PanoramaOptions::EQUIRECTANGULAR:
00281         case HuginBase::PanoramaOptions::MERCATOR:
00282         case HuginBase::PanoramaOptions::LAMBERT:
00283         case HuginBase::PanoramaOptions::MILLER_CYLINDRICAL:
00284         case HuginBase::PanoramaOptions::ARCHITECTURAL:
00285             // flip to the other direction of the other side horizontally.
00286             if (vertex_c[0] < x_midpoint) vertex_c[0] += x_add_360;
00287             else vertex_c[0] -= x_add_360;
00288             break;
00289         case HuginBase::PanoramaOptions::SINUSOIDAL:
00290         case HuginBase::PanoramaOptions::ALBERS_EQUAL_AREA_CONIC:
00291             if (vertex_c[0] < x_midpoint)
00292             {
00293                 vertex_c[0] +=
00294                       visualization_state->GetProjectionInfo()->GetXAdd360(vertex_c[1]);
00295             } else {
00296                 vertex_c[0] -=
00297                       visualization_state->GetProjectionInfo()->GetXAdd360(vertex_c[1]);
00298             }
00299             break;
00300         case HuginBase::PanoramaOptions::TRANSVERSE_MERCATOR:
00301             // flip to the other direction of the other side vertically
00302             if (vertex_c[1] < y_midpoint) vertex_c[1] += y_add_360;
00303             else vertex_c[1] -= y_add_360;
00304             break;
00305         default:
00306             // all other projection, no special handling
00307             break;
00308     }
00309 }
00310 
00311 void VertexCoordRemapper::RecursiveUpdate(unsigned int node_num,
00312                                    unsigned int stretch_x, unsigned stretch_y)
00313 {
00314     // find where we are and what we are mapping
00315     // TODO? GetPosition is called by GetInputCoordinates, reuse results?
00316     unsigned int x, y, row_size, depth;
00317     tree.GetPosition(node_num, x, y, row_size, depth);
00318     tree.GetInputCoordinates(node_num, tex_coords);
00319     TreeNode *node = &tree.nodes[node_num],
00320              *parent = (depth) ?
00321                        &tree.nodes[tree.GetParentId(x, y, row_size, depth)] : 0,
00322              *left = (x % 2) ?
00323                         &tree.nodes[tree.GetIndex(x-1, y, row_size, depth)] : 0,
00324              *top = (y % 2) ?
00325                         &tree.nodes[tree.GetIndex(x, y-1, row_size, depth)] : 0;
00326     bool valid[2][2];
00327     for (unsigned short int i = 0; i < 2; i++)
00328     {
00329         for (unsigned short int j = 0; j < 2; j++)
00330         {
00331             if (depth == 0)
00332             {
00333                 // the top level has no parent, so we must calculate all points
00334                 valid[i][j] =
00335                       transform.transformImgCoord(node->verts[i][j][0],
00336                                                   node->verts[i][j][1],
00337                                                   tex_coords[i][j][0] * width,
00338                                                   tex_coords[i][j][1] * height);
00339             } else
00340             // Look up where the point in the tree so far. If this is the first
00341             // occurrence of this point, we'll calculate the value.
00342             if (i == x %2 && j == y%2 && depth)
00343             {
00344                 // extract a corner from the parent.
00345                 node->verts[i][j][0] = parent->verts[i][j][0];
00346                 node->verts[i][j][1] = parent->verts[i][j][1];
00347                 valid[i][j] = !(parent->flags & (transform_fail_flag << (j*2 +i)));
00348             } else if (x % 2 && !i) {
00349                 // copy from the left
00350                 node->verts[0][j][0] = left->verts[1][j][0];
00351                 node->verts[0][j][1] = left->verts[1][j][1];
00352                 valid[i][j] = !(left->flags & (transform_fail_flag << (j *2 + 1)));
00353             } else if (y % 2 && !j) {
00354                 // copy from the top
00355                 node->verts[i][0][0] = top->verts[i][1][0];
00356                 node->verts[i][0][1] = top->verts[i][1][1];
00357                 valid[i][j] = !(top->flags & (transform_fail_flag << (2 + i)));
00358             } else {
00359                 // We can't find it easily, try a more expensive search.
00360                 // this will linearly interpolate along the edges where the
00361                 // subdivision was higher, avoiding gaps when the detail level
00362                 // was lower above or to the left.
00363                 if (!tree.GetTransform(x + i * (1 << stretch_x),
00364                                        y + j * (1 << stretch_y),
00365                                        depth, x, y, node->verts[i][j][0],
00366                                        node->verts[i][j][1]))
00367                 {
00368                     // We can't find it, so calculate it:
00369                     valid[i][j] =
00370                       transform.transformImgCoord(node->verts[i][j][0],
00371                                                   node->verts[i][j][1],
00372                                                   tex_coords[i][j][0] * width,
00373                                                   tex_coords[i][j][1] * height);
00374                     // If the face results from a patch subdivision (for
00375                     // aligning a subdivided face, otherwise it need not exist),
00376                     // use the midpoint of the parent's vertices instead of
00377                     // calculating a transformed one, so we can use less
00378                     // subdivision on the other side.
00379                     // subdivision in x
00380                     // FIXME there are still gaps. I think there's a logic error
00381                     if (   depth // not on the top level.
00382                         // patching in y
00383                         && (parent->flags & patch_flag_x)
00384                         // we must be in the middle of the split, the nodes on
00385                         // the corners of the parent line up anyway.
00386                         && ((i + x) % 2)
00387                         // we should be on the side away from the subdivison
00388                         // (+ve y).
00389                         && j
00390                         // If we are alao split in y we can use the middle
00391                         // node to provide more detail.
00392                         && (!((parent->flags & split_flag_y) && !(y % 2)))
00393                         // don't do this if we cross the 180 degree seam
00394                         && (!(parent->flags & (vertex_side_flag_start * 15))))
00395                     {
00396                         node->verts[i][1][0] = (parent->verts[0][1][0]
00397                                                + parent->verts[1][1][0]) / 2.0;
00398                         node->verts[i][1][1] = (parent->verts[0][1][1]
00399                                                + parent->verts[1][1][1]) / 2.0;
00400                     }
00401                     // subdivision in y
00402                     if (   depth
00403                         && (parent->flags & patch_flag_y)
00404                         && ((j + y) % 2)
00405                         && i
00406                         && (!((parent->flags & split_flag_x) && !(x % 2)))
00407                         && (!(parent->flags & (vertex_side_flag_start * 15))))
00408                     {
00409                         node->verts[1][j][0] = (parent->verts[1][0][0]
00410                                                + parent->verts[1][1][0]) / 2.0;
00411                         node->verts[1][j][1] = (parent->verts[1][0][1]
00412                                                + parent->verts[1][1][1]) / 2.0;
00413                     }
00414                 } else {
00415                     // we managed to find it from data already known.
00416                     valid[i][j] = true;
00417                 }
00418             }
00419         }
00420     }
00421     
00422     // now for the recursion
00423     // which directions should we divide in?
00424     TestSubdivide(node_num);
00425     // add the flags for invlaid transormations
00426     for (unsigned int i = 0; i < 2; i++)
00427     {
00428         for (unsigned int j = 0; j < 2; j++)
00429         {
00430             if (!valid[i][j])
00431             {
00432                 node->flags |= transform_fail_flag << (j * 2 + i);
00433             }
00434         }
00435     }
00436     // if the face should be split, now recurse to the child nodes.
00437     if (node->flags & (split_flag_x | split_flag_y))
00438     {
00439         // we will split at least one way.
00440         if (!(node->flags & split_flag_x))
00441         {
00442             // we are not splitting across x, but will will across y.
00443             // so the quad will be twice as wide
00444             stretch_x++;
00445         }
00446         else if (!(node->flags & split_flag_y))
00447         {
00448             stretch_y++;
00449         }
00450         // find the top left sub-quad
00451         x *= 2;
00452         y *= 2;
00453         row_size *= 2;
00454         depth++;
00455         // the top left is always generated
00456         RecursiveUpdate(tree.GetIndex(x, y, row_size, depth),
00457                         stretch_x, stretch_y);
00458         // split in x
00459         if (node->flags & split_flag_x)
00460         {
00461             RecursiveUpdate(tree.GetIndex(x + 1, y, row_size, depth),
00462                             stretch_x, stretch_y);
00463         }
00464         // split in y
00465         if (node->flags & split_flag_y)
00466         {
00467             RecursiveUpdate(tree.GetIndex(x, y + 1, row_size, depth),
00468                             stretch_x, stretch_y);
00469             // if we are splitting in both directions, do the lower right corner
00470             if (node->flags & split_flag_x)
00471             {
00472                 RecursiveUpdate(tree.GetIndex(x + 1, y + 1, row_size, depth),
00473                     stretch_x, stretch_y);
00474             }
00475         }
00476     }
00477 }
00478 
00479 void VertexCoordRemapper::TestSubdivide(unsigned int node_id)
00480 {
00481     TreeNode *node = &tree.nodes[node_id];
00482     unsigned int x, y, row_size, depth;
00483     tree.GetPosition(node_id, x, y, row_size, depth);
00484     unsigned short int flags = 0;
00485     if (depth < min_depth)
00486     {
00487         // subdivide at least min_depth times
00488         // we will need more information for non-trivial children
00489         SetLengthAndAngle(node);
00490         flags |= split_flag_x | split_flag_y;
00491     } else {
00492         unsigned int parent_id = tree.GetParentId(node_id);        
00493         TreeNode *parent = &tree.nodes[parent_id];
00494         // minimum length check. We use the length of the top edge to test for
00495         // subdivision in x, and the length of the left edge for subdivision in
00496         // y.
00497         SetLengthAndAngle(node);
00498         bool do_not_split_x = node->length_x * scale < min_length,
00499              do_not_split_y = node->length_y * scale < min_length;
00500         if (depth == max_depth)
00501         {
00502             // don't subdivide more than max_depth times
00503             do_not_split_x = true;
00504             do_not_split_y = true;
00505         }    
00506         // if we have stopped splitting in some direction, don't consider
00507         // splitting in that direction again.
00508         if (!(tree.nodes[parent_id].flags & split_flag_x))
00509         {
00510             do_not_split_x = true;
00511         }
00512         else if (!(tree.nodes[parent_id].flags & split_flag_y))
00513         {
00514             do_not_split_y = true;
00515         }
00516         // if it has only subdivided to patch up between two subdivision levels,
00517         // don't consider subdividing for adding more detail.
00518         if (tree.nodes[parent_id].flags & patch_flag_x)
00519         {
00520             do_not_split_x = true;
00521         }
00522         else if (tree.nodes[parent_id].flags & patch_flag_y)
00523         {
00524             do_not_split_y = true;
00525         }
00526 
00527         // If the angles have deviated too much from the parent then we should
00528         // add more detail, however if it is fairly flat then we don't need to.
00529         // It is possible for the angles to remain constant but the length
00530         // of the lines to change dramatically, so we check for a big difference
00531         // between the length of the parent node and twice the length of child.
00532         // if the child does not change the length much and the angle is small,
00533         // then we have enough detail, and we don't split.
00534         float ang_x = node->angle_x - parent->angle_x;
00535         if (ang_x < 0) ang_x = -ang_x;
00536         if (ang_x > M_PI) ang_x = 2 * M_PI - ang_x;
00537         float length_difference_x
00538                           = (parent->length_x - (2.0 * node->length_x)) * scale;
00539         if (length_difference_x < 0.0)
00540         {
00541             length_difference_x = -length_difference_x;
00542         }
00543         if (ang_x < min_angle && length_difference_x < min_length_difference)
00544         {
00545             do_not_split_x = true;
00546         }
00547         float ang_y = node->angle_y - parent->angle_y;
00548         if (ang_y < 0) ang_y = -ang_y;
00549         if (ang_y > M_PI) ang_y = 2 * M_PI - ang_y;
00550         float length_difference_y
00551                           = (parent->length_y - (2.0 * node->length_y)) * scale;
00552         if (length_difference_y < 0.0)
00553         {
00554             length_difference_y = -length_difference_y;
00555         }
00556         if (ang_y < min_angle  && length_difference_y < min_length_difference)
00557         {
00558             do_not_split_y = true;
00559         }
00560         // if the face is entirely off the screen, we should not subdivide it.
00561         // get the screen area
00562         vigra::Rect2D viewport = visualization_state->GetVisibleArea();
00563         // add a margin for safety, we don't want to clip too much stuff that
00564         // curls back on to the screen. We add 2 as we need some space around 
00565         // very small panoramas that have enlarged to fit the preview window, 
00566         // and even with a fairly large margin rounding to int may lose the
00567         // border completely.
00568         viewport.addBorder((int) (2.0 + offscreen_safety_margin * scale));
00569         bool all_left = true, all_right = true,
00570              all_above = true, all_bellow = true;
00571         for (unsigned int ix = 0; ix < 2; ix++)
00572         {
00573             for (unsigned int iy = 0; iy < 2; iy++)
00574             {
00575                 if (node->verts[ix][iy][0] > viewport.left())
00576                                                     all_left = false;
00577                 if (node->verts[ix][iy][0] < viewport.right())
00578                                                     all_right = false;
00579                 if (node->verts[ix][iy][1] > viewport.top())
00580                                                     all_above = false;
00581                 if (node->verts[ix][iy][1] < viewport.bottom())
00582                                                     all_bellow = false;
00583             }
00584         }
00585         if (all_left || all_right || all_bellow || all_above)
00586         {
00587             // all vertices are off a side of the screen. This catches most
00588             // cases where the face is off the screen. Don't allow subdivisions:
00589             do_not_split_x = true;
00590             do_not_split_y = true;
00591         }
00592         if (!do_not_split_x) flags |= split_flag_x;
00593         if (!do_not_split_y) flags |= split_flag_y;
00594     }
00595     /* Flag the vertices on a different side of the +/-180 degree seam.
00596      * We don't want to flag any vertices if the face covers a continuous
00597      * area of the transformation.
00598      */
00599     // We don't need to mark the first few subdivisions, but this is necessary
00600     // when patch subdivisions become possible.
00601     if (depth >= min_depth)
00602     {
00603         // determine if it is likely to be non-continuous.
00604         // this needs to be false for leaf-node faces in the centre that span
00605         // across the '0 degree' point, and true for faces that span the +/-180
00606         // degree split. It doesn't really matter what it is set to otherwise.
00607         bool noncontinuous = false;
00608         OutputProjectionInfo *i = visualization_state->GetProjectionInfo();
00609         switch (output_projection)
00610         {
00611             case HuginBase::PanoramaOptions::RECTILINEAR:
00612                 // we don't need faces to cross from one side to another. Faces
00613                 // all / partially 'behind the viewer' are skipped because the
00614                 // vertices behind the viewer are marked.
00615                 break;
00616             case HuginBase::PanoramaOptions::FULL_FRAME_FISHEYE:
00617             case HuginBase::PanoramaOptions::STEREOGRAPHIC:
00618             case HuginBase::PanoramaOptions::LAMBERT_AZIMUTHAL:
00619             case HuginBase::PanoramaOptions::HAMMER_AITOFF:
00620                 // A mesh covering the extremities of a disk projection should
00621                 // be using a TexCoordRemapper instead, otherwise, a point
00622                 // mapping to the border will be stretched across the disk.
00623                 break;
00624             case HuginBase::PanoramaOptions::CYLINDRICAL:
00625             case HuginBase::PanoramaOptions::EQUIRECTANGULAR:
00626             case HuginBase::PanoramaOptions::MERCATOR:
00627             case HuginBase::PanoramaOptions::LAMBERT:
00628             case HuginBase::PanoramaOptions::MILLER_CYLINDRICAL:
00629             case HuginBase::PanoramaOptions::PANINI:
00630                         case HuginBase::PanoramaOptions::BIPLANE:
00631                         case HuginBase::PanoramaOptions::TRIPLANE:
00632             case HuginBase::PanoramaOptions::GENERAL_PANINI:
00633                 // Cylinderical-like projections have the seam across the left
00634                 // and right edge. We'll take any face within the middle third
00635                 // to be continuous, the rest possibly noncontinuous.
00636                 if (   node->verts[0][0][0] < lower_bound
00637                     || node->verts[0][0][0] > upper_bound
00638                     || node->verts[1][0][0] < lower_bound
00639                     || node->verts[1][0][0] > upper_bound
00640                     || node->verts[0][1][0] < lower_bound
00641                     || node->verts[0][1][0] > upper_bound
00642                     || node->verts[1][1][0] < lower_bound
00643                     || node->verts[1][1][0] > upper_bound)
00644                 {
00645                     noncontinuous = true;
00646                     // flag nodes on the right hand side.
00647                     if (node->verts[0][0][0] > x_midpoint)
00648                     {
00649                         flags |= vertex_side_flag_start;
00650                     }
00651                     if (node->verts[0][1][0] > x_midpoint)
00652                     {
00653                         flags |= vertex_side_flag_start * 2;
00654                     }
00655                     if (node->verts[1][0][0] > x_midpoint)
00656                     {
00657                         flags |= vertex_side_flag_start * 4;
00658                     }
00659                     if (node->verts[1][1][0] > x_midpoint)
00660                     {
00661                         flags |= vertex_side_flag_start * 8;
00662                     }
00663                 }
00664                 break;
00665             case HuginBase::PanoramaOptions::SINUSOIDAL:
00666                 // like above, but the bounds change with height
00667                 if (   node->verts[0][0][0] < i->GetLowerX(node->verts[0][0][1])
00668                     || node->verts[0][0][0] > i->GetUpperX(node->verts[0][0][1])
00669                     || node->verts[1][0][0] < i->GetLowerX(node->verts[1][0][1])
00670                     || node->verts[1][0][0] > i->GetUpperX(node->verts[1][0][1])
00671                     || node->verts[0][1][0] < i->GetLowerX(node->verts[0][1][1])
00672                     || node->verts[0][1][0] > i->GetUpperX(node->verts[0][1][1])
00673                     || node->verts[1][1][0] < i->GetLowerX(node->verts[1][1][1])
00674                     || node->verts[1][1][0] > i->GetUpperX(node->verts[1][1][1])
00675                    )
00676                 {
00677                     noncontinuous = true;
00678                     // flag nodes on the right hand side.
00679                     if (node->verts[0][0][0] > x_midpoint)
00680                     {
00681                         flags |= vertex_side_flag_start;
00682                     }
00683                     if (node->verts[0][1][0] > x_midpoint)
00684                     {
00685                         flags |= vertex_side_flag_start * 2;
00686                     }
00687                     if (node->verts[1][0][0] > x_midpoint)
00688                     {
00689                         flags |= vertex_side_flag_start * 4;
00690                     }
00691                     if (node->verts[1][1][0] > x_midpoint)
00692                     {
00693                         flags |= vertex_side_flag_start * 8;
00694                     }
00695                 }
00696                 break;
00697             case HuginBase::PanoramaOptions::TRANSVERSE_MERCATOR:
00698                 // like the cylindrical ones, but vertically.
00699                 if (   node->verts[0][0][1] < lower_bound_h
00700                     || node->verts[0][0][1] > upper_bound_h
00701                     || node->verts[1][0][1] < lower_bound_h
00702                     || node->verts[1][0][1] > upper_bound_h
00703                     || node->verts[0][1][1] < lower_bound_h
00704                     || node->verts[0][1][1] > upper_bound_h
00705                     || node->verts[1][1][1] < lower_bound_h
00706                     || node->verts[1][1][1] > upper_bound_h)
00707                 {
00708                     noncontinuous = true;
00709                     // flag nodes on the bottom side.
00710                     if (node->verts[0][0][1] > y_midpoint)
00711                     {
00712                         flags |= vertex_side_flag_start;
00713                     }
00714                     if (node->verts[0][1][1] > y_midpoint)
00715                     {
00716                         flags |= vertex_side_flag_start * 2;
00717                     }
00718                     if (node->verts[1][0][1] > y_midpoint)
00719                     {
00720                         flags |= vertex_side_flag_start * 4;
00721                     }
00722                     if (node->verts[1][1][1] > y_midpoint)
00723                     {
00724                         flags |= vertex_side_flag_start * 8;
00725                     }
00726                 }
00727                 break;
00728             case HuginBase::PanoramaOptions::ALBERS_EQUAL_AREA_CONIC:
00729                 break;
00730             default:
00731                 // all other projections, no special handling
00732                 break;
00733         }
00734         if (noncontinuous)
00735         {
00736             // If all the side flags are set, we only have one face on one side:
00737             // Remove all of those flags and we have the same face, but we can
00738             // use the presence of any flags to detect when we have two.
00739             if ((flags / vertex_side_flag_start) % 16 == 15)
00740             {
00741                 flags &= ~(vertex_side_flag_start * 15);
00742             }
00743         }
00744     }
00745 
00746     node->flags = flags;
00747     // check if the faces to the left are subdivided at a higher level
00748     if (x && !(flags & split_flag_y) && (depth < max_depth))
00749     {
00750         // get the potentially more subdivided version
00751         double dest_x, dest_y;
00752         unsigned int subdiv_node;
00753         subdiv_node = tree.GetTransform(x * 2, y * 2 + 1,
00754                                         depth * 2,
00755                                         x * 2, y * 2, dest_x, dest_y);
00756         if (subdiv_node > node_id)
00757         {
00758             // we should have a more subdivided node on the left.
00759             // mark it for patching up to line up with it.
00760             node->flags |= split_flag_y | patch_flag_y;
00761         }
00762     }
00763     if (y && !(flags & split_flag_x) && (depth < max_depth))
00764     {
00765         // get the potentially more subdivided version
00766         double dest_x, dest_y;
00767         unsigned int subdiv_node;
00768         subdiv_node = tree.GetTransform(x * 2 +  1, y * 2,
00769                                         depth * 2,
00770                                         x * 2, y * 2, dest_x, dest_y);
00771         if (subdiv_node > node_id)
00772         {
00773             // we should have a more subdivided node on the left.
00774             // mark it for patching up to line up with it.
00775             node->flags |= split_flag_x | patch_flag_x;
00776         }
00777     }
00778 }
00779 
00780 void VertexCoordRemapper::SetLengthAndAngle(TreeNode *node)
00781 {
00782     float xdx = node->verts[0][0][0] - node->verts[1][0][0],
00783           xdy = node->verts[0][0][1] - node->verts[1][0][1],
00784           ydx = node->verts[0][0][0] - node->verts[0][1][0],
00785           ydy = node->verts[0][0][1] - node->verts[0][1][1];
00786     // this is the length of the edge with y = 0
00787     node->length_x = sqrt(sqr(xdx) + sqr(xdy));
00788     // this is the length of the edge with x = 0
00789     node->length_y = sqrt(sqr(ydx) + sqr(ydy));
00790     // find the angle of the top and left edge.
00791     node->angle_x = atan2(xdy, xdx);
00792     node->angle_y = atan2(ydy, ydx);
00793 }
00794 
00795 unsigned int VertexCoordRemapper::Tree::GetParentId(const unsigned int nodenum)
00796 {
00797     // get the parent id of a node specifed by its index.
00798     unsigned int x, y, row_size, depth;
00799     GetPosition(nodenum, x, y, row_size, depth);
00800     return GetParentId(x, y, row_size, depth);
00801 }
00802 
00803 unsigned int VertexCoordRemapper::Tree::GetParentId(unsigned int x,
00804                                              unsigned int y,
00805                                              unsigned int row_size,
00806                                              unsigned int depth)
00807 {
00808     // get the parent id of a node specified by an exact location.
00809     // the parent is the one that takes the group of four elements around this
00810     // one in the above depth. All the dimensions are halved.
00811     x /= 2;
00812     y /= 2;
00813     row_size /= 2;
00814     depth--;    
00815     return GetIndex(x, y, row_size, depth);
00816 }
00817 
00818 unsigned int VertexCoordRemapper::Tree::GetDepth(const unsigned int node_num)
00819 {
00820     unsigned int depth = 0, count = 0;
00821     while (node_num > count)
00822     {
00823         depth++;
00824         count += 1 << (depth * 2);
00825     }
00826     return depth;
00827 }
00828 
00829 void VertexCoordRemapper::Tree::GetPosition(const unsigned int node_num,
00830                                      unsigned int &x, unsigned int &y,
00831                                      unsigned int &row_size,
00832                                      unsigned int &depth)
00833 {
00834     depth = GetDepth(node_num);
00835     row_size = 1 << depth;
00836     unsigned int sub = 0;
00837     if (depth)
00838     {
00839       for (unsigned int d = 0; d < depth; d++)
00840       {
00841           sub += (1 << (d * 2));
00842       }
00843     }
00844     unsigned int position_id = node_num - sub;
00845     x = position_id % row_size;
00846     y = position_id / row_size;
00847 }
00848 
00849 unsigned int VertexCoordRemapper::Tree::GetIndex(const unsigned int x,
00850                                           const unsigned int y,
00851                                           const unsigned int row_size,
00852                                           unsigned int depth)
00853 {
00854     unsigned int add = 0;
00855     while (depth)
00856     {
00857         depth--;
00858         add += 1 << (depth * 2);
00859     }
00860     return add + x + y * row_size;
00861 }
00862 
00863 void VertexCoordRemapper::Tree::GetInputCoordinates(const unsigned int node_num,
00864                                              double coords[2][2][2])
00865 {
00866     // find the coordinates of each point at the vertices in the original image.
00867     // this is used to look up the remapped coordinates and provide texture
00868     // coordinates.
00869     // it halves in size several times depending on depth.
00870     unsigned int x, y, row_size, depth;
00871     GetPosition(node_num, x, y, row_size, depth);
00872     // now we can locate the upper right corner
00873     double row_spacing = 1.0 / (double) row_size;
00874     coords[0][0][0] = row_spacing * (double) x;
00875     coords[0][0][1] = row_spacing * (double) y;
00876     // the extent is dependent on the direction of the subdivisions.
00877     // at least one direction has an extent of row_spacing. It can double up
00878     // if the parent did not subdivide in a direction, recursively up the tree.
00879     bool scale_x = false;
00880     double opp = 1.0;
00881     if (node_num != 0)
00882     {    
00883       unsigned int parent_id = GetParentId(x, y, row_size, depth);
00884       if (!(nodes[parent_id].flags & split_flag_x))
00885       {
00886           while (!(nodes[parent_id].flags & split_flag_x))
00887           {
00888               parent_id = GetParentId(parent_id);
00889               opp *= 2.0;
00890           }
00891           scale_x = true;
00892       } else {
00893           while (!(nodes[parent_id].flags & split_flag_y))
00894           {
00895               parent_id = GetParentId(parent_id);
00896               opp *= 2.0;
00897           }
00898       }
00899     }
00900     opp *= row_spacing;
00901     coords[1][0][0] = coords[0][0][0] + (scale_x ? opp : row_spacing);
00902     coords[1][0][1] = coords[0][0][1];
00903     coords[0][1][0] = coords[0][0][0];
00904     coords[0][1][1] = coords[0][0][1] + (scale_x ? row_spacing : opp);
00905     coords[1][1][0] = coords[1][0][0];
00906     coords[1][1][1] = coords[0][1][1];
00907     // now we transform the results so that we map to the cropped region instead
00908     // of the whole image.
00909     for (unsigned int i = 0; i < 2; i++)
00910     {
00911         for (unsigned int j = 0; j < 2; j++)
00912         {
00913             coords[i][j][0] = coords[i][j][0] * x_crop_scale + x_crop_offs;
00914             coords[i][j][1] = coords[i][j][1] * y_crop_scale + y_crop_offs;
00915         }
00916     }
00917 }
00918 
00919 void VertexCoordRemapper::Tree::ResetIndex()
00920 {
00921     cur_tree_node = 0;
00922 }
00923 
00924 unsigned int VertexCoordRemapper::Tree::GetNext()
00925 {
00926     // depth first search for leaf nodes with cur_tree_node
00927     unsigned int x, y, row_size, depth;
00928     GetPosition(cur_tree_node, x, y, row_size, depth);
00929     // we want to backtrack until we find an unexplored sub branch. At this 
00930     //     point, we trace down the tree until we get to a leaf.
00931     if (cur_tree_node != 0)
00932     {
00933         bool done = false;
00934         while (!done && cur_tree_node != 0)
00935         {
00936             unsigned int xd = x % 2;
00937             unsigned int yd = y % 2;
00938             x /= 2;
00939             y /= 2;
00940             row_size /= 2;
00941             depth--;
00942             cur_tree_node = GetIndex(x, y, row_size, depth);
00943             if (cur_tree_node == 0 && xd == 1 && yd == 1)
00944             {
00945                 // we've expanded all the options and got back to the top
00946                 return 0; // no more faces
00947             }
00948             // where does this split?
00949             bool sx = ((nodes[cur_tree_node].flags & split_flag_x) != 0);
00950             bool sy = ((nodes[cur_tree_node].flags & split_flag_y) != 0);
00951             // have we used all the split options?
00952             if (!(((sx && xd) || !sx) && ((sy && yd) || !sy)))
00953             {
00954                 // we've found an unexpanded child, take the next one:
00955                 done = true;
00956                 x *= 2;
00957                 y *= 2;
00958                 row_size *= 2;
00959                 depth++;
00960                 if (sx && !xd && !yd) {
00961                     x++;
00962                 } else if ((sx && xd && sy && !yd) || (!sx && sy && !yd)) {
00963                     y++;
00964                 } else if (sx && sy && !xd && yd) {
00965                     x++; y++;
00966                 }
00967                 cur_tree_node = GetIndex(x, y, row_size, depth);
00968                 // if we've made our way to the root node, we've run out out
00969                 // of options.
00970                 if (cur_tree_node == 0)
00971                 {
00972                     return 0;
00973                 }
00974             };
00975         }
00976     }
00977     // find the first leaf on this subtree, taking top left every time.
00978     while (nodes[cur_tree_node].flags & (split_flag_x | split_flag_y))
00979     {
00980         x *= 2;
00981         y *= 2;
00982         row_size *= 2;
00983         depth++;
00984         cur_tree_node = GetIndex(x, y, row_size, depth);
00985     }
00986     return cur_tree_node;
00987 }
00988 
00989 unsigned int VertexCoordRemapper::Tree::GetTransform(unsigned int src_x, unsigned int src_y,
00990                                       unsigned int max_depth,
00991                                       unsigned int stop_x, unsigned int stop_y,
00992                                       double &dest_x, double &dest_y)
00993 {
00994     // we start at the top and take children, prefering the lowest numbered one,
00995     // until we either have no children, or we have found a child that knows the
00996     // exact place we are looking for.
00997     unsigned int no_x = 0, no_y = 0, row_size = 1, depth = 0,
00998                  rem_x = src_x, rem_y = src_y,
00999                  step_x = 1 << max_depth, step_y = 1 << max_depth,
01000                  node_id = GetIndex(no_x, no_y, row_size, depth);
01001     while  (   !(    (rem_x == 0 || rem_x == step_x)
01002                   && (rem_y == 0 || rem_y == step_y))
01003             &&  (nodes[node_id].flags & (split_flag_x | split_flag_y)))
01004     {
01005         // split, try and take earliest child node that leads towards goal
01006         depth++;
01007         no_x *= 2;
01008         no_y *= 2;
01009         row_size *= 2;
01010         if (nodes[node_id].flags & split_flag_x)
01011         {
01012             step_x /= 2;
01013             if (rem_x > step_x)
01014             {
01015                 no_x ++;
01016                 rem_x -= step_x;
01017             }
01018         }
01019         if (nodes[node_id].flags & split_flag_y)
01020         {
01021             step_y /= 2;
01022             if (rem_y > step_y)
01023             {
01024                 no_y ++;
01025                 rem_y -= step_y;
01026             }
01027         }
01028         // we want to stop if we have got a node at least as high as the
01029         // requested stopping node. Anything at a higher depth is after it.
01030         if (depth > max_depth) return 0;
01031         node_id = GetIndex(no_x, no_y, row_size, depth);
01032         if (depth == max_depth && no_x >= stop_x && no_y >= stop_y) return 0;
01033     }
01034     // if any of the vertices we are want to use are invalid (e.g. behind the
01035     // viewer in a rectilinear projection) we don't want to risk using them:
01036     if (nodes[node_id].flags & (transform_fail_flag * 15)) return 0;
01037     // if this face crosses a discontinuity, we should be using a point off
01038     // screen instead of in the middle. Refuse to use these faces
01039     if (nodes[node_id].flags & (vertex_side_flag_start * 15)) return 0;
01040     // linearly interpolate the node's corners.
01041     // most of the time we only use factors of 0 and 1, we don't want to make
01042     // points up except when trying to connect a point on a highly subdivided
01043     // face to a point that doesn't exist because it is on a less subdivided
01044     // face, in which case it needs to line up with the linear interpolation of
01045     // the less subdivided face. This is along one edge, so the other direction
01046     // should have a blending factor of 0 or 1.
01047     double xf = (double) rem_x / (double) step_x;
01048     double yf = (double) rem_y / (double) step_y;
01049     
01050     double top_x = (1.0 - xf) * nodes[node_id].verts[0][0][0]
01051               + xf * nodes[node_id].verts[1][0][0],
01052            bottom_x = (1-.0 - xf) * nodes[node_id].verts[0][1][0]
01053               + xf * nodes[node_id].verts[1][1][0],
01054            top_y = (1.0 - xf) * nodes[node_id].verts[0][0][1]
01055               + xf * nodes[node_id].verts[1][0][1],
01056            bottom_y = (1.0 - xf) * nodes[node_id].verts[0][1][1]
01057               + xf * nodes[node_id].verts[1][1][1];
01058     dest_x = top_x * (1.0 - yf) + bottom_x * yf;
01059     dest_y = top_y * (1.0 - yf) + bottom_y* yf;
01060     return node_id;
01061 }
01062 

Generated on 28 Jul 2016 for Hugintrunk by  doxygen 1.4.7