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::GENERAL_PANINI:
00631                 // Cylinderical-like projections have the seam across the left
00632                 // and right edge. We'll take any face within the middle third
00633                 // to be continuous, the rest possibly noncontinuous.
00634                 if (   node->verts[0][0][0] < lower_bound
00635                     || node->verts[0][0][0] > upper_bound
00636                     || node->verts[1][0][0] < lower_bound
00637                     || node->verts[1][0][0] > upper_bound
00638                     || node->verts[0][1][0] < lower_bound
00639                     || node->verts[0][1][0] > upper_bound
00640                     || node->verts[1][1][0] < lower_bound
00641                     || node->verts[1][1][0] > upper_bound)
00642                 {
00643                     noncontinuous = true;
00644                     // flag nodes on the right hand side.
00645                     if (node->verts[0][0][0] > x_midpoint)
00646                     {
00647                         flags |= vertex_side_flag_start;
00648                     }
00649                     if (node->verts[0][1][0] > x_midpoint)
00650                     {
00651                         flags |= vertex_side_flag_start * 2;
00652                     }
00653                     if (node->verts[1][0][0] > x_midpoint)
00654                     {
00655                         flags |= vertex_side_flag_start * 4;
00656                     }
00657                     if (node->verts[1][1][0] > x_midpoint)
00658                     {
00659                         flags |= vertex_side_flag_start * 8;
00660                     }
00661                 }
00662                 break;
00663             case HuginBase::PanoramaOptions::SINUSOIDAL:
00664                 // like above, but the bounds change with height
00665                 if (   node->verts[0][0][0] < i->GetLowerX(node->verts[0][0][1])
00666                     || node->verts[0][0][0] > i->GetUpperX(node->verts[0][0][1])
00667                     || node->verts[1][0][0] < i->GetLowerX(node->verts[1][0][1])
00668                     || node->verts[1][0][0] > i->GetUpperX(node->verts[1][0][1])
00669                     || node->verts[0][1][0] < i->GetLowerX(node->verts[0][1][1])
00670                     || node->verts[0][1][0] > i->GetUpperX(node->verts[0][1][1])
00671                     || node->verts[1][1][0] < i->GetLowerX(node->verts[1][1][1])
00672                     || node->verts[1][1][0] > i->GetUpperX(node->verts[1][1][1])
00673                    )
00674                 {
00675                     noncontinuous = true;
00676                     // flag nodes on the right hand side.
00677                     if (node->verts[0][0][0] > x_midpoint)
00678                     {
00679                         flags |= vertex_side_flag_start;
00680                     }
00681                     if (node->verts[0][1][0] > x_midpoint)
00682                     {
00683                         flags |= vertex_side_flag_start * 2;
00684                     }
00685                     if (node->verts[1][0][0] > x_midpoint)
00686                     {
00687                         flags |= vertex_side_flag_start * 4;
00688                     }
00689                     if (node->verts[1][1][0] > x_midpoint)
00690                     {
00691                         flags |= vertex_side_flag_start * 8;
00692                     }
00693                 }
00694                 break;
00695             case HuginBase::PanoramaOptions::TRANSVERSE_MERCATOR:
00696                 // like the cylindrical ones, but vertically.
00697                 if (   node->verts[0][0][1] < lower_bound_h
00698                     || node->verts[0][0][1] > upper_bound_h
00699                     || node->verts[1][0][1] < lower_bound_h
00700                     || node->verts[1][0][1] > upper_bound_h
00701                     || node->verts[0][1][1] < lower_bound_h
00702                     || node->verts[0][1][1] > upper_bound_h
00703                     || node->verts[1][1][1] < lower_bound_h
00704                     || node->verts[1][1][1] > upper_bound_h)
00705                 {
00706                     noncontinuous = true;
00707                     // flag nodes on the bottom side.
00708                     if (node->verts[0][0][1] > y_midpoint)
00709                     {
00710                         flags |= vertex_side_flag_start;
00711                     }
00712                     if (node->verts[0][1][1] > y_midpoint)
00713                     {
00714                         flags |= vertex_side_flag_start * 2;
00715                     }
00716                     if (node->verts[1][0][1] > y_midpoint)
00717                     {
00718                         flags |= vertex_side_flag_start * 4;
00719                     }
00720                     if (node->verts[1][1][1] > y_midpoint)
00721                     {
00722                         flags |= vertex_side_flag_start * 8;
00723                     }
00724                 }
00725                 break;
00726             case HuginBase::PanoramaOptions::ALBERS_EQUAL_AREA_CONIC:
00727                 break;
00728             default:
00729                 // all other projections, no special handling
00730                 break;
00731         }
00732         if (noncontinuous)
00733         {
00734             // If all the side flags are set, we only have one face on one side:
00735             // Remove all of those flags and we have the same face, but we can
00736             // use the presence of any flags to detect when we have two.
00737             if ((flags / vertex_side_flag_start) % 16 == 15)
00738             {
00739                 flags &= ~(vertex_side_flag_start * 15);
00740             }
00741         }
00742     }
00743 
00744     node->flags = flags;
00745     // check if the faces to the left are subdivided at a higher level
00746     if (x && !(flags & split_flag_y) && (depth < max_depth))
00747     {
00748         // get the potentially more subdivided version
00749         double dest_x, dest_y;
00750         unsigned int subdiv_node;
00751         subdiv_node = tree.GetTransform(x * 2, y * 2 + 1,
00752                                         depth * 2,
00753                                         x * 2, y * 2, dest_x, dest_y);
00754         if (subdiv_node > node_id)
00755         {
00756             // we should have a more subdivided node on the left.
00757             // mark it for patching up to line up with it.
00758             node->flags |= split_flag_y | patch_flag_y;
00759         }
00760     }
00761     if (y && !(flags & split_flag_x) && (depth < max_depth))
00762     {
00763         // get the potentially more subdivided version
00764         double dest_x, dest_y;
00765         unsigned int subdiv_node;
00766         subdiv_node = tree.GetTransform(x * 2 +  1, y * 2,
00767                                         depth * 2,
00768                                         x * 2, y * 2, dest_x, dest_y);
00769         if (subdiv_node > node_id)
00770         {
00771             // we should have a more subdivided node on the left.
00772             // mark it for patching up to line up with it.
00773             node->flags |= split_flag_x | patch_flag_x;
00774         }
00775     }
00776 }
00777 
00778 void VertexCoordRemapper::SetLengthAndAngle(TreeNode *node)
00779 {
00780     float xdx = node->verts[0][0][0] - node->verts[1][0][0],
00781           xdy = node->verts[0][0][1] - node->verts[1][0][1],
00782           ydx = node->verts[0][0][0] - node->verts[0][1][0],
00783           ydy = node->verts[0][0][1] - node->verts[0][1][1];
00784     // this is the length of the edge with y = 0
00785     node->length_x = sqrt(sqr(xdx) + sqr(xdy));
00786     // this is the length of the edge with x = 0
00787     node->length_y = sqrt(sqr(ydx) + sqr(ydy));
00788     // find the angle of the top and left edge.
00789     node->angle_x = atan2(xdy, xdx);
00790     node->angle_y = atan2(ydy, ydx);
00791 }
00792 
00793 unsigned int VertexCoordRemapper::Tree::GetParentId(const unsigned int nodenum)
00794 {
00795     // get the parent id of a node specifed by its index.
00796     unsigned int x, y, row_size, depth;
00797     GetPosition(nodenum, x, y, row_size, depth);
00798     return GetParentId(x, y, row_size, depth);
00799 }
00800 
00801 unsigned int VertexCoordRemapper::Tree::GetParentId(unsigned int x,
00802                                              unsigned int y,
00803                                              unsigned int row_size,
00804                                              unsigned int depth)
00805 {
00806     // get the parent id of a node specified by an exact location.
00807     // the parent is the one that takes the group of four elements around this
00808     // one in the above depth. All the dimensions are halved.
00809     x /= 2;
00810     y /= 2;
00811     row_size /= 2;
00812     depth--;    
00813     return GetIndex(x, y, row_size, depth);
00814 }
00815 
00816 unsigned int VertexCoordRemapper::Tree::GetDepth(const unsigned int node_num)
00817 {
00818     unsigned int depth = 0, count = 0;
00819     while (node_num > count)
00820     {
00821         depth++;
00822         count += 1 << (depth * 2);
00823     }
00824     return depth;
00825 }
00826 
00827 void VertexCoordRemapper::Tree::GetPosition(const unsigned int node_num,
00828                                      unsigned int &x, unsigned int &y,
00829                                      unsigned int &row_size,
00830                                      unsigned int &depth)
00831 {
00832     depth = GetDepth(node_num);
00833     row_size = 1 << depth;
00834     unsigned int sub = 0;
00835     if (depth)
00836     {
00837       for (unsigned int d = 0; d < depth; d++)
00838       {
00839           sub += (1 << (d * 2));
00840       }
00841     }
00842     unsigned int position_id = node_num - sub;
00843     x = position_id % row_size;
00844     y = position_id / row_size;
00845 }
00846 
00847 unsigned int VertexCoordRemapper::Tree::GetIndex(const unsigned int x,
00848                                           const unsigned int y,
00849                                           const unsigned int row_size,
00850                                           unsigned int depth)
00851 {
00852     unsigned int add = 0;
00853     while (depth)
00854     {
00855         depth--;
00856         add += 1 << (depth * 2);
00857     }
00858     return add + x + y * row_size;
00859 }
00860 
00861 void VertexCoordRemapper::Tree::GetInputCoordinates(const unsigned int node_num,
00862                                              double coords[2][2][2])
00863 {
00864     // find the coordinates of each point at the vertices in the original image.
00865     // this is used to look up the remapped coordinates and provide texture
00866     // coordinates.
00867     // it halves in size several times depending on depth.
00868     unsigned int x, y, row_size, depth;
00869     GetPosition(node_num, x, y, row_size, depth);
00870     // now we can locate the upper right corner
00871     double row_spacing = 1.0 / (double) row_size;
00872     coords[0][0][0] = row_spacing * (double) x;
00873     coords[0][0][1] = row_spacing * (double) y;
00874     // the extent is dependent on the direction of the subdivisions.
00875     // at least one direction has an extent of row_spacing. It can double up
00876     // if the parent did not subdivide in a direction, recursively up the tree.
00877     bool scale_x = false;
00878     double opp = 1.0;
00879     if (node_num != 0)
00880     {    
00881       unsigned int parent_id = GetParentId(x, y, row_size, depth);
00882       if (!(nodes[parent_id].flags & split_flag_x))
00883       {
00884           while (!(nodes[parent_id].flags & split_flag_x))
00885           {
00886               parent_id = GetParentId(parent_id);
00887               opp *= 2.0;
00888           }
00889           scale_x = true;
00890       } else {
00891           while (!(nodes[parent_id].flags & split_flag_y))
00892           {
00893               parent_id = GetParentId(parent_id);
00894               opp *= 2.0;
00895           }
00896       }
00897     }
00898     opp *= row_spacing;
00899     coords[1][0][0] = coords[0][0][0] + (scale_x ? opp : row_spacing);
00900     coords[1][0][1] = coords[0][0][1];
00901     coords[0][1][0] = coords[0][0][0];
00902     coords[0][1][1] = coords[0][0][1] + (scale_x ? row_spacing : opp);
00903     coords[1][1][0] = coords[1][0][0];
00904     coords[1][1][1] = coords[0][1][1];
00905     // now we transform the results so that we map to the cropped region instead
00906     // of the whole image.
00907     for (unsigned int i = 0; i < 2; i++)
00908     {
00909         for (unsigned int j = 0; j < 2; j++)
00910         {
00911             coords[i][j][0] = coords[i][j][0] * x_crop_scale + x_crop_offs;
00912             coords[i][j][1] = coords[i][j][1] * y_crop_scale + y_crop_offs;
00913         }
00914     }
00915 }
00916 
00917 void VertexCoordRemapper::Tree::ResetIndex()
00918 {
00919     cur_tree_node = 0;
00920 }
00921 
00922 unsigned int VertexCoordRemapper::Tree::GetNext()
00923 {
00924     // depth first search for leaf nodes with cur_tree_node
00925     unsigned int x, y, row_size, depth;
00926     GetPosition(cur_tree_node, x, y, row_size, depth);
00927     // we want to backtrack until we find an unexplored sub branch. At this 
00928     //     point, we trace down the tree until we get to a leaf.
00929     if (cur_tree_node != 0)
00930     {
00931         bool done = false;
00932         while (!done && cur_tree_node != 0)
00933         {
00934             unsigned int xd = x % 2;
00935             unsigned int yd = y % 2;
00936             x /= 2;
00937             y /= 2;
00938             row_size /= 2;
00939             depth--;
00940             cur_tree_node = GetIndex(x, y, row_size, depth);
00941             if (cur_tree_node == 0 && xd == 1 && yd == 1)
00942             {
00943                 // we've expanded all the options and got back to the top
00944                 return 0; // no more faces
00945             }
00946             // where does this split?
00947             bool sx = ((nodes[cur_tree_node].flags & split_flag_x) != 0);
00948             bool sy = ((nodes[cur_tree_node].flags & split_flag_y) != 0);
00949             // have we used all the split options?
00950             if (!(((sx && xd) || !sx) && ((sy && yd) || !sy)))
00951             {
00952                 // we've found an unexpanded child, take the next one:
00953                 done = true;
00954                 x *= 2;
00955                 y *= 2;
00956                 row_size *= 2;
00957                 depth++;
00958                 if (sx && !xd && !yd) {
00959                     x++;
00960                 } else if ((sx && xd && sy && !yd) || (!sx && sy && !yd)) {
00961                     y++;
00962                 } else if (sx && sy && !xd && yd) {
00963                     x++; y++;
00964                 }
00965                 cur_tree_node = GetIndex(x, y, row_size, depth);
00966                 // if we've made our way to the root node, we've run out out
00967                 // of options.
00968                 if (cur_tree_node == 0)
00969                 {
00970                     return 0;
00971                 }
00972             };
00973         }
00974     }
00975     // find the first leaf on this subtree, taking top left every time.
00976     while (nodes[cur_tree_node].flags & (split_flag_x | split_flag_y))
00977     {
00978         x *= 2;
00979         y *= 2;
00980         row_size *= 2;
00981         depth++;
00982         cur_tree_node = GetIndex(x, y, row_size, depth);
00983     }
00984     return cur_tree_node;
00985 }
00986 
00987 unsigned int VertexCoordRemapper::Tree::GetTransform(unsigned int src_x, unsigned int src_y,
00988                                       unsigned int max_depth,
00989                                       unsigned int stop_x, unsigned int stop_y,
00990                                       double &dest_x, double &dest_y)
00991 {
00992     // we start at the top and take children, prefering the lowest numbered one,
00993     // until we either have no children, or we have found a child that knows the
00994     // exact place we are looking for.
00995     unsigned int no_x = 0, no_y = 0, row_size = 1, depth = 0,
00996                  rem_x = src_x, rem_y = src_y,
00997                  step_x = 1 << max_depth, step_y = 1 << max_depth,
00998                  node_id = GetIndex(no_x, no_y, row_size, depth);
00999     while  (   !(    (rem_x == 0 || rem_x == step_x)
01000                   && (rem_y == 0 || rem_y == step_y))
01001             &&  (nodes[node_id].flags & (split_flag_x | split_flag_y)))
01002     {
01003         // split, try and take earliest child node that leads towards goal
01004         depth++;
01005         no_x *= 2;
01006         no_y *= 2;
01007         row_size *= 2;
01008         if (nodes[node_id].flags & split_flag_x)
01009         {
01010             step_x /= 2;
01011             if (rem_x > step_x)
01012             {
01013                 no_x ++;
01014                 rem_x -= step_x;
01015             }
01016         }
01017         if (nodes[node_id].flags & split_flag_y)
01018         {
01019             step_y /= 2;
01020             if (rem_y > step_y)
01021             {
01022                 no_y ++;
01023                 rem_y -= step_y;
01024             }
01025         }
01026         // we want to stop if we have got a node at least as high as the
01027         // requested stopping node. Anything at a higher depth is after it.
01028         if (depth > max_depth) return 0;
01029         node_id = GetIndex(no_x, no_y, row_size, depth);
01030         if (depth == max_depth && no_x >= stop_x && no_y >= stop_y) return 0;
01031     }
01032     // if any of the vertices we are want to use are invalid (e.g. behind the
01033     // viewer in a rectilinear projection) we don't want to risk using them:
01034     if (nodes[node_id].flags & (transform_fail_flag * 15)) return 0;
01035     // if this face crosses a discontinuity, we should be using a point off
01036     // screen instead of in the middle. Refuse to use these faces
01037     if (nodes[node_id].flags & (vertex_side_flag_start * 15)) return 0;
01038     // linearly interpolate the node's corners.
01039     // most of the time we only use factors of 0 and 1, we don't want to make
01040     // points up except when trying to connect a point on a highly subdivided
01041     // face to a point that doesn't exist because it is on a less subdivided
01042     // face, in which case it needs to line up with the linear interpolation of
01043     // the less subdivided face. This is along one edge, so the other direction
01044     // should have a blending factor of 0 or 1.
01045     double xf = (double) rem_x / (double) step_x;
01046     double yf = (double) rem_y / (double) step_y;
01047     
01048     double top_x = (1.0 - xf) * nodes[node_id].verts[0][0][0]
01049               + xf * nodes[node_id].verts[1][0][0],
01050            bottom_x = (1-.0 - xf) * nodes[node_id].verts[0][1][0]
01051               + xf * nodes[node_id].verts[1][1][0],
01052            top_y = (1.0 - xf) * nodes[node_id].verts[0][0][1]
01053               + xf * nodes[node_id].verts[1][0][1],
01054            bottom_y = (1.0 - xf) * nodes[node_id].verts[0][1][1]
01055               + xf * nodes[node_id].verts[1][1][1];
01056     dest_x = top_x * (1.0 - yf) + bottom_x * yf;
01057     dest_y = top_y * (1.0 - yf) + bottom_y* yf;
01058     return node_id;
01059 }
01060 

Generated on 20 Apr 2018 for Hugintrunk by  doxygen 1.4.7