Missing Points in Container RangeSearch

[marekdam]

I am using a simple custom version of your octree that is something like:

// Adaptor
using int3D = std::array<long int, 3>;
using CustomBoundingBox = std::array<int3D ,2>;
struct AdaptorBasicsCustom
{
    static inline long int& point_comp(int3D & pt, OrthoTree::dim_type iDimension)
    {
        return pt[iDimension];
    }

    static constexpr long int point_comp_c(int3D const& pt, OrthoTree::dim_type iDimension)
    {
        return pt[iDimension];
    }

    static inline int3D& box_min(CustomBoundingBox& box) { return box[0]; }
    static inline int3D& box_max(CustomBoundingBox& box) { return box[1]; }
    static constexpr int3D const& box_min_c(CustomBoundingBox const& box) { return box[0]; }
    static constexpr int3D const& box_max_c(CustomBoundingBox const& box) { return box[1]; }
};

using AdaptorCustom = OrthoTree::AdaptorGeneralBase<3, int3D, CustomBoundingBox, AdaptorBasicsCustom, long int>;

// Tailored Quadtree objects
using QuadtreePointCustom = OrthoTree::OrthoTreePoint<3, int3D, CustomBoundingBox, AdaptorCustom, long int>;
using QuadtreeContainerCustom = OrthoTree::OrthoTreeContainerPoint<QuadtreePointCustom, int3D>;

I have tried the RangeSearch function, but when I compare it to a simple slow method I use now it misses certain points. It happens on a octree with 423 points. I have not been able to reproduce the issue on a smaller problem.

However, I think the lines [1368-1374] might be an issue. Debugging the code I found that a nodeChild bounding box containing the missing point should have been giving a true statement for overlap with the range but it wasn't. Are these comparisons supposed to be strictly lesser than and strictly greater than? I think it will work if it is changed to <= or >=. I am a beginner when it comes to octtrees and the morton code indexing, so I am not sure if I am missing something.

        for (dim_type iDim = 0; iDim < nDimension && bOverlap; ++iDim)
        {
          if (IsValidKey(keyChild & (morton_node_id_type{ 1 } << iDim)))
            bOverlap &= AD::point_comp_c(AD::box_min_c(nodeChild.box), iDim) < AD::point_comp_c(AD::box_max_c(range), iDim);
          else
            bOverlap &= AD::point_comp_c(AD::box_max_c(nodeChild.box), iDim) > AD::point_comp_c(AD::box_min_c(range), iDim);
        }

[attcs]

Hi! Thank you for the bug report. The issue is more sophisticated than this, the solution requires more caution. The issue is connected with the point classification: the point on the edge of the grid will be placed in the higher grid number, which could eventually be a different Octree node, but if the range is fit into the smaller grid, it will not be able to reach this node.

I think both your suggestion and the following code replacement together could help you temporarily, while I am working on a proper solution:

• replace in the rangeSearchRoot(): from auto keyNodeSmallest = this->FindSmallestNodeKey(keyRange); to auto keyNodeSmallest = GetRootKey();
  • And the code fix what you mentioned:

    for (dim_type iDim = 0; iDim < nDimension && bOverlap; ++iDim)
    {
      if (IsValidKey(keyChild & (morton_node_id_type{ 1 } << iDim)))
        bOverlap &= AD::point_comp_c(AD::box_min_c(nodeChild.box), iDim) <= AD::point_comp_c(AD::box_max_c(range), iDim);
      else
        bOverlap &= AD::point_comp_c(AD::box_max_c(nodeChild.box), iDim) > AD::point_comp_c(AD::box_min_c(range), iDim);
    }

    It will be slower because it has to go down on the tree from the root.

[marekdam]

Nice thanks! That seems to work for this smaller problem at the moment. I'll leave the issue open if you plan on implementing the rigorous solution.

[attcs]

I pushed up the fix.

[marekdam]

I tried the fix you pushed, but it is still missing some points. I see that you did not modify the greater than line for bOverlap. If you make that >= I do get the correct result again.

[attcs]

I did! I resolved this issue in two commits, one of the relation fix and the other is the smallest node. This is the current state:

        auto bOverlap = true;
        for (dim_type iDim = 0; iDim < nDimension && bOverlap; ++iDim)
        {
          if (IsValidKey(keyChild & (morton_node_id_type{ 1 } << iDim)))
            bOverlap &= AD::point_comp_c(AD::box_min_c(nodeChild.box), iDim) <= AD::point_comp_c(AD::box_max_c(range), iDim);
          else
            bOverlap &= AD::point_comp_c(AD::box_max_c(nodeChild.box), iDim) > AD::point_comp_c(AD::box_min_c(range), iDim);
        }
        if (!bOverlap)
          continue;

[marekdam]

Ok, but I mean the ">" in the else branch. Shouldn't that be ">=", so in total:

if (IsValidKey(keyChild & (morton_node_id_type{ 1 } << iDim)))
  bOverlap &= AD::point_comp_c(AD::box_min_c(nodeChild.box), iDim) <= AD::point_comp_c(AD::box_max_c(range), iDim);
else
  bOverlap &= AD::point_comp_c(AD::box_max_c(nodeChild.box), iDim) >= AD::point_comp_c(AD::box_min_c(range), iDim);

That is the change I made in addition to your recent commits and it works for me.

[attcs]

You are right, thank you! I made the fixup.

[marekdam]

Great, yes it works. Thanks again

[marekdam]

Can you reopen this issue? I have been trying a larger example and missing points again.

[attcs]

Of course!

[marekdam]

I don't think I understand your octree well enough to debug this time but I made some tests and exported the data to see if maybe you could figure it out. The attached text file contains the point data found in octree.GetData(). I created the Tree with nDepthMax=3.

In one test case I used a bounding box defined by: {39,43,72} and {49,53,76}. My brute force algorithm finds 44 points in the box, while the RangeSearch function finds 40 points. The missing points have ids of {164,165,375,1549}

The missing points have coordinates: 164: {46, 43, 76} 165: {46, 43, 75} 375: {41, 43, 76} 1549: {41, 43, 75}

Let me know if that helps!

octree_data.txt

[attcs]

It is great help! I reproduced the bug, I'll try to solve at evening!

[attcs]

I solved. The issue was the child node's bounding box were fastly accumulated numerical errors if the geometry_type is an integral type, therefore highly diverging from the theoretical raster which is used for the point classification. The solution is a new method for the node's bounding box calculation which does not inherit the sizes from the parent node.

I attached a unit test (BruteForceRangeSearch_UsingPredefinedData_IfAvailable), but not your dataset, to make it easier to check.

[marekdam]

Great, that works now!