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 // panormama'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         
00306     }
00307 }
00308 
00309 void VertexCoordRemapper::RecursiveUpdate(unsigned int node_num,
00310                                    unsigned int stretch_x, unsigned stretch_y)
00311 {
00312     // find where we are and what we are mapping
00313     // TODO? GetPosition is called by GetInputCoordinates, reuse results?
00314     unsigned int x, y, row_size, depth;
00315     tree.GetPosition(node_num, x, y, row_size, depth);
00316     tree.GetInputCoordinates(node_num, tex_coords);
00317     TreeNode *node = &tree.nodes[node_num],
00318              *parent = (depth) ?
00319                        &tree.nodes[tree.GetParentId(x, y, row_size, depth)] : 0,
00320              *left = (x % 2) ?
00321                         &tree.nodes[tree.GetIndex(x-1, y, row_size, depth)] : 0,
00322              *top = (y % 2) ?
00323                         &tree.nodes[tree.GetIndex(x, y-1, row_size, depth)] : 0;
00324     bool valid[2][2];
00325     for (unsigned short int i = 0; i < 2; i++)
00326     {
00327         for (unsigned short int j = 0; j < 2; j++)
00328         {
00329             if (depth == 0)
00330             {
00331                 // the top level has no parent, so we must calculate all points
00332                 valid[i][j] =
00333                       transform.transformImgCoord(node->verts[i][j][0],
00334                                                   node->verts[i][j][1],
00335                                                   tex_coords[i][j][0] * width,
00336                                                   tex_coords[i][j][1] * height);
00337             } else
00338             // Look up where the point in the tree so far. If this is the first
00339             // occurrence of this point, we'll calculate the value.
00340             if (i == x %2 && j == y%2 && depth)
00341             {
00342                 // extract a corner from the parent.
00343                 node->verts[i][j][0] = parent->verts[i][j][0];
00344                 node->verts[i][j][1] = parent->verts[i][j][1];
00345                 valid[i][j] = !(parent->flags & (transform_fail_flag << (j*2 +i)));
00346             } else if (x % 2 && !i) {
00347                 // copy from the left
00348                 node->verts[0][j][0] = left->verts[1][j][0];
00349                 node->verts[0][j][1] = left->verts[1][j][1];
00350                 valid[i][j] = !(left->flags & (transform_fail_flag << (j *2 + 1)));
00351             } else if (y % 2 && !j) {
00352                 // copy from the top
00353                 node->verts[i][0][0] = top->verts[i][1][0];
00354                 node->verts[i][0][1] = top->verts[i][1][1];
00355                 valid[i][j] = !(top->flags & (transform_fail_flag << (2 + i)));
00356             } else {
00357                 // We can't find it easily, try a more expensive search.
00358                 // this will linearly interpolate along the edges where the
00359                 // subdivision was higher, avoiding gaps when the detail level
00360                 // was lower above or to the left.
00361                 if (!tree.GetTransform(x + i * (1 << stretch_x),
00362                                        y + j * (1 << stretch_y),
00363                                        depth, x, y, node->verts[i][j][0],
00364                                        node->verts[i][j][1]))
00365                 {
00366                     // We can't find it, so calculate it:
00367                     valid[i][j] =
00368                       transform.transformImgCoord(node->verts[i][j][0],
00369                                                   node->verts[i][j][1],
00370                                                   tex_coords[i][j][0] * width,
00371                                                   tex_coords[i][j][1] * height);
00372                     // If the face results from a patch subdivision (for
00373                     // aligning a subdivided face, otherwise it need not exist),
00374                     // use the midpoint of the parent's vertices instead of
00375                     // calculating a transformed one, so we can use less
00376                     // subdivision on the other side.
00377                     // subdivision in x
00378                     // FIXME there are still gaps. I think there's a logic error
00379                     if (   depth // not on the top level.
00380                         // patching in y
00381                         && (parent->flags & patch_flag_x)
00382                         // we must be in the middle of the split, the nodes on
00383                         // the corners of the parent line up anyway.
00384                         && ((i + x) % 2)
00385                         // we should be on the side away from the subdivison
00386                         // (+ve y).
00387                         && j
00388                         // If we are alao split in y we can use the middle
00389                         // node to provide more detail.
00390                         && (!((parent->flags & split_flag_y) && !(y % 2)))
00391                         // don't do this if we cross the 180 degree seam
00392                         && (!(parent->flags & (vertex_side_flag_start * 15))))
00393                     {
00394                         node->verts[i][1][0] = (parent->verts[0][1][0]
00395                                                + parent->verts[1][1][0]) / 2.0;
00396                         node->verts[i][1][1] = (parent->verts[0][1][1]
00397                                                + parent->verts[1][1][1]) / 2.0;
00398                     }
00399                     // subdivision in y
00400                     if (   depth
00401                         && (parent->flags & patch_flag_y)
00402                         && ((j + y) % 2)
00403                         && i
00404                         && (!((parent->flags & split_flag_x) && !(x % 2)))
00405                         && (!(parent->flags & (vertex_side_flag_start * 15))))
00406                     {
00407                         node->verts[1][j][0] = (parent->verts[1][0][0]
00408                                                + parent->verts[1][1][0]) / 2.0;
00409                         node->verts[1][j][1] = (parent->verts[1][0][1]
00410                                                + parent->verts[1][1][1]) / 2.0;
00411                     }
00412                 } else {
00413                     // we managed to find it from data already known.
00414                     valid[i][j] = true;
00415                 }
00416             }
00417         }
00418     }
00419     
00420     // now for the recursion
00421     // which directions should we divide in?
00422     TestSubdivide(node_num);
00423     // add the flags for invlaid transormations
00424     for (unsigned int i = 0; i < 2; i++)
00425     {
00426         for (unsigned int j = 0; j < 2; j++)
00427         {
00428             if (!valid[i][j])
00429             {
00430                 node->flags |= transform_fail_flag << (j * 2 + i);
00431             }
00432         }
00433     }
00434     // if the face should be split, now recurse to the child nodes.
00435     if (node->flags & (split_flag_x | split_flag_y))
00436     {
00437         // we will split at least one way.
00438         if (!(node->flags & split_flag_x))
00439         {
00440             // we are not splitting across x, but will will across y.
00441             // so the quad will be twice as wide
00442             stretch_x++;
00443         }
00444         else if (!(node->flags & split_flag_y))
00445         {
00446             stretch_y++;
00447         }
00448         // find the top left sub-quad
00449         x *= 2;
00450         y *= 2;
00451         row_size *= 2;
00452         depth++;
00453         // the top left is always generated
00454         RecursiveUpdate(tree.GetIndex(x, y, row_size, depth),
00455                         stretch_x, stretch_y);
00456         // split in x
00457         if (node->flags & split_flag_x)
00458         {
00459             RecursiveUpdate(tree.GetIndex(x + 1, y, row_size, depth),
00460                             stretch_x, stretch_y);
00461         }
00462         // split in y
00463         if (node->flags & split_flag_y)
00464         {
00465             RecursiveUpdate(tree.GetIndex(x, y + 1, row_size, depth),
00466                             stretch_x, stretch_y);
00467             // if we are splitting in both directions, do the lower right corner
00468             if (node->flags & split_flag_x)
00469             {
00470                 RecursiveUpdate(tree.GetIndex(x + 1, y + 1, row_size, depth),
00471                     stretch_x, stretch_y);
00472             }
00473         }
00474     }
00475 }
00476 
00477 void VertexCoordRemapper::TestSubdivide(unsigned int node_id)
00478 {
00479     TreeNode *node = &tree.nodes[node_id];
00480     unsigned int x, y, row_size, depth;
00481     tree.GetPosition(node_id, x, y, row_size, depth);
00482     unsigned short int flags = 0;
00483     if (depth < min_depth)
00484     {
00485         // subdivide at least min_depth times
00486         // we will need more information for non-trivial children
00487         SetLengthAndAngle(node);
00488         flags |= split_flag_x | split_flag_y;
00489     } else {
00490         unsigned int parent_id = tree.GetParentId(node_id);        
00491         TreeNode *parent = &tree.nodes[parent_id];
00492         // minimum length check. We use the length of the top edge to test for
00493         // subdivision in x, and the length of the left edge for subdivision in
00494         // y.
00495         SetLengthAndAngle(node);
00496         bool do_not_split_x = node->length_x * scale < min_length,
00497              do_not_split_y = node->length_y * scale < min_length;
00498         if (depth == max_depth)
00499         {
00500             // don't subdivide more than max_depth times
00501             do_not_split_x = true;
00502             do_not_split_y = true;
00503         }    
00504         // if we have stopped splitting in some direction, don't consider
00505         // splitting in that direction again.
00506         if (!(tree.nodes[parent_id].flags & split_flag_x))
00507         {
00508             do_not_split_x = true;
00509         }
00510         else if (!(tree.nodes[parent_id].flags & split_flag_y))
00511         {
00512             do_not_split_y = true;
00513         }
00514         // if it has only subdivided to patch up between two subdivision levels,
00515         // don't consider subdividing for adding more detail.
00516         if (tree.nodes[parent_id].flags & patch_flag_x)
00517         {
00518             do_not_split_x = true;
00519         }
00520         else if (tree.nodes[parent_id].flags & patch_flag_y)
00521         {
00522             do_not_split_y = true;
00523         }
00524 
00525         // If the angles have deviated too much from the parent then we should
00526         // add more detail, however if it is fairly flat then we don't need to.
00527         // It is possible for the angles to remain constant but the length
00528         // of the lines to change dramatically, so we check for a big difference
00529         // between the length of the parent node and twice the length of child.
00530         // if the child does not change the length much and the angle is small,
00531         // then we have enough detail, and we don't split.
00532         float ang_x = node->angle_x - parent->angle_x;
00533         if (ang_x < 0) ang_x = -ang_x;
00534         if (ang_x > M_PI) ang_x = 2 * M_PI - ang_x;
00535         float length_difference_x
00536                           = (parent->length_x - (2.0 * node->length_x)) * scale;
00537         if (length_difference_x < 0.0)
00538         {
00539             length_difference_x = -length_difference_x;
00540         }
00541         if (ang_x < min_angle && length_difference_x < min_length_difference)
00542         {
00543             do_not_split_x = true;
00544         }
00545         float ang_y = node->angle_y - parent->angle_y;
00546         if (ang_y < 0) ang_y = -ang_y;
00547         if (ang_y > M_PI) ang_y = 2 * M_PI - ang_y;
00548         float length_difference_y
00549                           = (parent->length_y - (2.0 * node->length_y)) * scale;
00550         if (length_difference_y < 0.0)
00551         {
00552             length_difference_y = -length_difference_y;
00553         }
00554         if (ang_y < min_angle  && length_difference_y < min_length_difference)
00555         {
00556             do_not_split_y = true;
00557         }
00558         // if the face is entirely off the screen, we should not subdivide it.
00559         // get the screen area
00560         vigra::Rect2D viewport = visualization_state->GetVisibleArea();
00561         // add a margin for safety, we don't want to clip too much stuff that
00562         // curls back on to the screen. We add 2 as we need some space around 
00563         // very small panoramas that have enlarged to fit the preview window, 
00564         // and even with a fairly large margin rounding to int may lose the
00565         // border completely.
00566         viewport.addBorder((int) (2.0 + offscreen_safety_margin * scale));
00567         bool all_left = true, all_right = true,
00568              all_above = true, all_bellow = true;
00569         for (unsigned int ix = 0; ix < 2; ix++)
00570         {
00571             for (unsigned int iy = 0; iy < 2; iy++)
00572             {
00573                 if (node->verts[ix][iy][0] > viewport.left())
00574                                                     all_left = false;
00575                 if (node->verts[ix][iy][0] < viewport.right())
00576                                                     all_right = false;
00577                 if (node->verts[ix][iy][1] > viewport.top())
00578                                                     all_above = false;
00579                 if (node->verts[ix][iy][1] < viewport.bottom())
00580                                                     all_bellow = false;
00581             }
00582         }
00583         if (all_left || all_right || all_bellow || all_above)
00584         {
00585             // all vertices are off a side of the screen. This catches most
00586             // cases where the face is off the screen. Don't allow subdivisions:
00587             do_not_split_x = true;
00588             do_not_split_y = true;
00589         }
00590         if (!do_not_split_x) flags |= split_flag_x;
00591         if (!do_not_split_y) flags |= split_flag_y;
00592     }
00593     /* Flag the vertices on a different side of the +/-180 degree seam.
00594      * We don't want to flag any vertices if the face covers a continuous
00595      * area of the transformation.
00596      */
00597     // We don't need to mark the first few subdivisions, but this is necessary
00598     // when patch subdivisions become possible.
00599     if (depth >= min_depth)
00600     {
00601         // determine if it is likely to be non-continuous.
00602         // this needs to be false for leaf-node faces in the centre that span
00603         // across the '0 degree' point, and true for faces that span the +/-180
00604         // degree split. It doesn't really matter what it is set to otherwise.
00605         bool noncontinuous = false;
00606         OutputProjectionInfo *i = visualization_state->GetProjectionInfo();
00607         switch (output_projection)
00608         {
00609             case HuginBase::PanoramaOptions::RECTILINEAR:
00610                 // we don't need faces to cross from one side to another. Faces
00611                 // all / partially 'behind the viewer' are skipped because the
00612                 // vertices behind the viewer are marked.
00613                 break;
00614             case HuginBase::PanoramaOptions::FULL_FRAME_FISHEYE:
00615             case HuginBase::PanoramaOptions::STEREOGRAPHIC:
00616             case HuginBase::PanoramaOptions::LAMBERT_AZIMUTHAL:
00617             case HuginBase::PanoramaOptions::HAMMER_AITOFF:
00618                 // A mesh covering the extremities of a disk projection should
00619                 // be using a TexCoordRemapper instead, otherwise, a point
00620                 // mapping to the border will be stretched across the disk.
00621                 break;
00622             case HuginBase::PanoramaOptions::CYLINDRICAL:
00623             case HuginBase::PanoramaOptions::EQUIRECTANGULAR:
00624             case HuginBase::PanoramaOptions::MERCATOR:
00625             case HuginBase::PanoramaOptions::LAMBERT:
00626             case HuginBase::PanoramaOptions::MILLER_CYLINDRICAL:
00627             case HuginBase::PanoramaOptions::PANINI:
00628                         case HuginBase::PanoramaOptions::BIPLANE:
00629                         case HuginBase::PanoramaOptions::TRIPLANE:
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         }
00729         if (noncontinuous)
00730         {
00731             // If all the side flags are set, we only have one face on one side:
00732             // Remove all of those flags and we have the same face, but we can
00733             // use the presence of any flags to detect when we have two.
00734             if ((flags / vertex_side_flag_start) % 16 == 15)
00735             {
00736                 flags &= ~(vertex_side_flag_start * 15);
00737             }
00738         }
00739     }
00740 
00741     node->flags = flags;
00742     // check if the faces to the left are subdivided at a higher level
00743     if (x && !(flags & split_flag_y) && (depth < max_depth))
00744     {
00745         // get the potentially more subdivided version
00746         double dest_x, dest_y;
00747         unsigned int subdiv_node;
00748         subdiv_node = tree.GetTransform(x * 2, y * 2 + 1,
00749                                         depth * 2,
00750                                         x * 2, y * 2, dest_x, dest_y);
00751         if (subdiv_node > node_id)
00752         {
00753             // we should have a more subdivided node on the left.
00754             // mark it for patching up to line up with it.
00755             node->flags |= split_flag_y | patch_flag_y;
00756         }
00757     }
00758     if (y && !(flags & split_flag_x) && (depth < max_depth))
00759     {
00760         // get the potentially more subdivided version
00761         double dest_x, dest_y;
00762         unsigned int subdiv_node;
00763         subdiv_node = tree.GetTransform(x * 2 +  1, y * 2,
00764                                         depth * 2,
00765                                         x * 2, y * 2, dest_x, dest_y);
00766         if (subdiv_node > node_id)
00767         {
00768             // we should have a more subdivided node on the left.
00769             // mark it for patching up to line up with it.
00770             node->flags |= split_flag_x | patch_flag_x;
00771         }
00772     }
00773 }
00774 
00775 void VertexCoordRemapper::SetLengthAndAngle(TreeNode *node)
00776 {
00777     float xdx = node->verts[0][0][0] - node->verts[1][0][0],
00778           xdy = node->verts[0][0][1] - node->verts[1][0][1],
00779           ydx = node->verts[0][0][0] - node->verts[0][1][0],
00780           ydy = node->verts[0][0][1] - node->verts[0][1][1];
00781     // this is the length of the edge with y = 0
00782     node->length_x = sqrt(sqr(xdx) + sqr(xdy));
00783     // this is the length of the edge with x = 0
00784     node->length_y = sqrt(sqr(ydx) + sqr(ydy));
00785     // find the angle of the top and left edge.
00786     node->angle_x = atan2(xdy, xdx);
00787     node->angle_y = atan2(ydy, ydx);
00788 }
00789 
00790 unsigned int VertexCoordRemapper::Tree::GetParentId(const unsigned int nodenum)
00791 {
00792     // get the parent id of a node specifed by its index.
00793     unsigned int x, y, row_size, depth;
00794     GetPosition(nodenum, x, y, row_size, depth);
00795     return GetParentId(x, y, row_size, depth);
00796 }
00797 
00798 unsigned int VertexCoordRemapper::Tree::GetParentId(unsigned int x,
00799                                              unsigned int y,
00800                                              unsigned int row_size,
00801                                              unsigned int depth)
00802 {
00803     // get the parent id of a node specified by an exact location.
00804     // the parent is the one that takes the group of four elements around this
00805     // one in the above depth. All the dimensions are halved.
00806     x /= 2;
00807     y /= 2;
00808     row_size /= 2;
00809     depth--;    
00810     return GetIndex(x, y, row_size, depth);
00811 }
00812 
00813 unsigned int VertexCoordRemapper::Tree::GetDepth(const unsigned int node_num)
00814 {
00815     unsigned int depth = 0, count = 0;
00816     while (node_num > count)
00817     {
00818         depth++;
00819         count += 1 << (depth * 2);
00820     }
00821     return depth;
00822 }
00823 
00824 void VertexCoordRemapper::Tree::GetPosition(const unsigned int node_num,
00825                                      unsigned int &x, unsigned int &y,
00826                                      unsigned int &row_size,
00827                                      unsigned int &depth)
00828 {
00829     depth = GetDepth(node_num);
00830     row_size = 1 << depth;
00831     unsigned int sub = 0;
00832     if (depth)
00833     {
00834       for (unsigned int d = 0; d < depth; d++)
00835       {
00836           sub += (1 << (d * 2));
00837       }
00838     }
00839     unsigned int position_id = node_num - sub;
00840     x = position_id % row_size;
00841     y = position_id / row_size;
00842 }
00843 
00844 unsigned int VertexCoordRemapper::Tree::GetIndex(const unsigned int x,
00845                                           const unsigned int y,
00846                                           const unsigned int row_size,
00847                                           unsigned int depth)
00848 {
00849     unsigned int add = 0;
00850     while (depth)
00851     {
00852         depth--;
00853         add += 1 << (depth * 2);
00854     }
00855     return add + x + y * row_size;
00856 }
00857 
00858 void VertexCoordRemapper::Tree::GetInputCoordinates(const unsigned int node_num,
00859                                              double coords[2][2][2])
00860 {
00861     // find the coordinates of each point at the vertices in the original image.
00862     // this is used to look up the remapped coordinates and provide texture
00863     // coordinates.
00864     // it halves in size several times depending on depth.
00865     unsigned int x, y, row_size, depth;
00866     GetPosition(node_num, x, y, row_size, depth);
00867     // now we can locate the upper right corner
00868     double row_spacing = 1.0 / (double) row_size;
00869     coords[0][0][0] = row_spacing * (double) x;
00870     coords[0][0][1] = row_spacing * (double) y;
00871     // the extent is dependent on the direction of the subdivisions.
00872     // at least one direction has an extent of row_spacing. It can double up
00873     // if the parent did not subdivide in a direction, recursively up the tree.
00874     bool scale_x = false;
00875     double opp = 1.0;
00876     if (node_num != 0)
00877     {    
00878       unsigned int parent_id = GetParentId(x, y, row_size, depth);
00879       if (!(nodes[parent_id].flags & split_flag_x))
00880       {
00881           while (!(nodes[parent_id].flags & split_flag_x))
00882           {
00883               parent_id = GetParentId(parent_id);
00884               opp *= 2.0;
00885           }
00886           scale_x = true;
00887       } else {
00888           while (!(nodes[parent_id].flags & split_flag_y))
00889           {
00890               parent_id = GetParentId(parent_id);
00891               opp *= 2.0;
00892           }
00893       }
00894     }
00895     opp *= row_spacing;
00896     coords[1][0][0] = coords[0][0][0] + (scale_x ? opp : row_spacing);
00897     coords[1][0][1] = coords[0][0][1];
00898     coords[0][1][0] = coords[0][0][0];
00899     coords[0][1][1] = coords[0][0][1] + (scale_x ? row_spacing : opp);
00900     coords[1][1][0] = coords[1][0][0];
00901     coords[1][1][1] = coords[0][1][1];
00902     // now we transform the results so that we map to the cropped region instead
00903     // of the whole image.
00904     for (unsigned int i = 0; i < 2; i++)
00905     {
00906         for (unsigned int j = 0; j < 2; j++)
00907         {
00908             coords[i][j][0] = coords[i][j][0] * x_crop_scale + x_crop_offs;
00909             coords[i][j][1] = coords[i][j][1] * y_crop_scale + y_crop_offs;
00910         }
00911     }
00912 }
00913 
00914 void VertexCoordRemapper::Tree::ResetIndex()
00915 {
00916     cur_tree_node = 0;
00917 }
00918 
00919 unsigned int VertexCoordRemapper::Tree::GetNext()
00920 {
00921     // depth first search for leaf nodes with cur_tree_node
00922     unsigned int x, y, row_size, depth;
00923     GetPosition(cur_tree_node, x, y, row_size, depth);
00924     // we want to backtrack until we find an unexplored sub branch. At this 
00925     //     point, we trace down the tree until we get to a leaf.
00926     if (cur_tree_node != 0)
00927     {
00928         unsigned int xd = 0, yd = 0;
00929         bool done = false;
00930         while (!done && cur_tree_node != 0)
00931         {
00932             xd = x % 2;
00933             yd = y % 2;
00934             x /= 2;
00935             y /= 2;
00936             row_size /= 2;
00937             depth--;
00938             cur_tree_node = GetIndex(x, y, row_size, depth);
00939             if (cur_tree_node == 0 && xd == 1 && yd == 1)
00940             {
00941                 // we've expanded all the options and got back to the top
00942                 return 0; // no more faces
00943             }
00944             // where does this split?
00945             bool sx = ((nodes[cur_tree_node].flags & split_flag_x) != 0);
00946             bool sy = ((nodes[cur_tree_node].flags & split_flag_y) != 0);
00947             // have we used all the split options?
00948             if (!(((sx && xd) || !sx) && ((sy && yd) || !sy)))
00949             {
00950                 // we've found an unexpanded child, take the next one:
00951                 done = true;
00952                 x *= 2;
00953                 y *= 2;
00954                 row_size *= 2;
00955                 depth++;
00956                 if (sx && !xd && !yd) {
00957                     x++;
00958                 } else if ((sx && xd && sy && !yd) || (!sx && sy && !yd)) {
00959                     y++;
00960                 } else if (sx && sy && !xd && yd) {
00961                     x++; y++;
00962                 }
00963                 cur_tree_node = GetIndex(x, y, row_size, depth);
00964                 // if we've made our way to the root node, we've run out out
00965                 // of options.
00966                 if (cur_tree_node == 0)
00967                 {
00968                     return 0;
00969                 }
00970             };
00971         }
00972     }
00973     // find the first leaf on this subtree, taking top left every time.
00974     while (nodes[cur_tree_node].flags & (split_flag_x | split_flag_y))
00975     {
00976         x *= 2;
00977         y *= 2;
00978         row_size *= 2;
00979         depth++;
00980         cur_tree_node = GetIndex(x, y, row_size, depth);
00981     }
00982     return cur_tree_node;
00983 }
00984 
00985 unsigned int VertexCoordRemapper::Tree::GetTransform(unsigned int src_x, unsigned int src_y,
00986                                       unsigned int max_depth,
00987                                       unsigned int stop_x, unsigned int stop_y,
00988                                       double &dest_x, double &dest_y)
00989 {
00990     // we start at the top and take children, prefering the lowest numbered one,
00991     // until we either have no children, or we have found a child that knows the
00992     // exact place we are looking for.
00993     unsigned int no_x = 0, no_y = 0, row_size = 1, depth = 0,
00994                  rem_x = src_x, rem_y = src_y,
00995                  step_x = 1 << max_depth, step_y = 1 << max_depth,
00996                  node_id = GetIndex(no_x, no_y, row_size, depth);
00997     while  (   !(    (rem_x == 0 || rem_x == step_x)
00998                   && (rem_y == 0 || rem_y == step_y))
00999             &&  (nodes[node_id].flags & (split_flag_x | split_flag_y)))
01000     {
01001         // split, try and take earliest child node that leads towards goal
01002         depth++;
01003         no_x *= 2;
01004         no_y *= 2;
01005         row_size *= 2;
01006         if (nodes[node_id].flags & split_flag_x)
01007         {
01008             step_x /= 2;
01009             if (rem_x > step_x)
01010             {
01011                 no_x ++;
01012                 rem_x -= step_x;
01013             }
01014         }
01015         if (nodes[node_id].flags & split_flag_y)
01016         {
01017             step_y /= 2;
01018             if (rem_y > step_y)
01019             {
01020                 no_y ++;
01021                 rem_y -= step_y;
01022             }
01023         }
01024         // we want to stop if we have got a node at least as high as the
01025         // requested stopping node. Anything at a higher depth is after it.
01026         if (depth > max_depth) return 0;
01027         node_id = GetIndex(no_x, no_y, row_size, depth);
01028         if (depth == max_depth && no_x >= stop_x && no_y >= stop_y) return 0;
01029     }
01030     // if any of the vertices we are want to use are invalid (e.g. behind the
01031     // viewer in a rectilinear projection) we don't want to risk using them:
01032     if (nodes[node_id].flags & (transform_fail_flag * 15)) return 0;
01033     // if this face crosses a discontinuity, we should be using a point off
01034     // screen instead of in the middle. Refuse to use these faces
01035     if (nodes[node_id].flags & (vertex_side_flag_start * 15)) return 0;
01036     // linearly interpolate the node's corners.
01037     // most of the time we only use factors of 0 and 1, we don't want to make
01038     // points up except when trying to connect a point on a highly subdivided
01039     // face to a point that doesn't exist because it is on a less subdivided
01040     // face, in which case it needs to line up with the linear interpolation of
01041     // the less subdivided face. This is along one edge, so the other direction
01042     // should have a blending factor of 0 or 1.
01043     double xf = (double) rem_x / (double) step_x;
01044     double yf = (double) rem_y / (double) step_y;
01045     
01046     double top_x = (1.0 - xf) * nodes[node_id].verts[0][0][0]
01047               + xf * nodes[node_id].verts[1][0][0],
01048            bottom_x = (1-.0 - xf) * nodes[node_id].verts[0][1][0]
01049               + xf * nodes[node_id].verts[1][1][0],
01050            top_y = (1.0 - xf) * nodes[node_id].verts[0][0][1]
01051               + xf * nodes[node_id].verts[1][0][1],
01052            bottom_y = (1.0 - xf) * nodes[node_id].verts[0][1][1]
01053               + xf * nodes[node_id].verts[1][1][1];
01054     dest_x = top_x * (1.0 - yf) + bottom_x * yf;
01055     dest_y = top_y * (1.0 - yf) + bottom_y* yf;
01056     return node_id;
01057 }
01058 

Generated on Sat Aug 2 01:25:38 2014 for Hugintrunk by  doxygen 1.3.9.1