Index out of range in lattice classes overload operator ()

[lukadt]

Dear @ddemidov, I am using your library in my current project and found an indexing issue in the following code:

https://github.com/ddemidov/mba/blob/e2003d8381f05633d68172d319a054bfd591181f/mba/mba.hpp#L373-L380

I changed to the code below by adding a guard condition since floor(u) - 1can take negative value assigned to an unsigned type.

double operator()(const point<NDim> &p) const {
            index<NDim> i;
            point<NDim> s;

            for(unsigned d = 0; d < NDim; ++d) {
                double u = (p[d] - cmin[d]) * hinv[d];
                i[d] = floor(u) - 1;                       
                s[d] = u - floor(u);
                if (floor(u) - 1 < 0)
                {
                    i[d] = 0;
                }
            }

Let me know if this fix is correct, in case I can provide you the dataset to reproduce the issue. Thanks in advance, Luca

[ddemidov]

I believe the method expects the point p to be within the bounding box specified in the constructor. If you want the method to work with points outside of the bounding box, you should also check for out of bounds error at the right hand side of the interval (something like i[d] + 4 < grid[d]).

Frankly, I would just leave the responsibility with the library user, so that no compute cycles are wasted unnecessarily. Or possibly hide the check behind #ifndef NDEBUG preprocessor guard, similar to how C++ deals with the out-of-bounds accesses to std::vector<>.

[lukadt]

Dear @ddemidov, thanks a lot for your explanation. I would like to have more explanations on the role of data points bounding box. The things I don't like and understand is that using my function interpolateData below the error I was referring previously exhibit only on particular values of oversampling param ( equal to 3.0 I got the error with oversampling equal to 5.0 the interpolation works fine).

So to sum up, how the bounding box reflects on data interpolation? It should fit closely or can also provide a bigger one?

    bool interpolateData(cv::Mat& xCoordinates, cv::Mat& yCoordinates, cv::Mat& values, double oversampling, uint maxlevels, double tool)
    {

        std::vector<mba::point<2>> referenceCoordinates;
        std::vector<double> xCoords, yCoords;
        std::vector<double> dataValues;

        // Bounding box containing the data points.
        mba::point<2> lo = { 0, 0 };
        mba::point<2> hi = { xCoordinates.cols * oversampling ,  xCoordinates.rows * oversampling  };

        int gridSize = xCoordinates.rows * oversampling;
        mba::index<2> grid = { 2, 2 };

        for (size_t j = 0; j < gridSize / oversampling; ++j)
        {
            for (size_t i = 0; i < gridSize / oversampling; ++i)
            {

                referenceCoordinates.push_back({ j * oversampling + (int)oversampling / 2,i * oversampling + (int)oversampling / 2 });

                xCoords.push_back(xCoordinates.at<double>(i, j));
                yCoords.push_back(yCoordinates.at<double>(i, j));
                dataValues.push_back(values.at<double>(i, j));

            }
        }

        mba::MBA<2> phiXCoord(lo, hi, grid, std::begin(referenceCoordinates), std::end(referenceCoordinates), std::begin(xCoords), maxlevels, tool);
        mba::MBA<2> phiYCoord(lo, hi, grid, std::begin(referenceCoordinates), std::end(referenceCoordinates), std::begin(yCoords), maxlevels, tool);
        mba::MBA<2> phiValue(lo, hi, grid, std::begin(referenceCoordinates), std::end(referenceCoordinates), std::begin(dataValues), maxlevels, tool);

        cv::Mat interpXCoord = cv::Mat::zeros(gridSize, gridSize, CV_64FC1);
        cv::Mat interpYCoord = cv::Mat::zeros(gridSize, gridSize, CV_64FC1);
        cv::Mat interpValues = cv::Mat::zeros(gridSize, gridSize, CV_64FC1);

        for (int i = 0; i < gridSize; i++)
        {
            for (int j = 0; j < gridSize; j++)
            {
                interpXCoord.at<double>(i, j) = phiXCoord({ (double)j,(double)i });
                interpYCoord.at<double>(i, j) = phiYCoord({ (double)j,(double)i });
                interpValues.at<double>(i, j) = phiValue({ (double)j,(double)i });
            }
        }

        xCoordinates = interpXCoord;
        yCoordinates = interpYCoord;
        values = interpValues;

        return true;

    }

[ddemidov]

You should provide a box that is slightly larger than the bounding box for your point cloud. Something like

span = max(p) - min(p);
lo = min(p) - 1e-3 * span;
hi = max(p) + 1e-3 * span;

[lukadt]

Dear @ddemidov , thanks a lot for the explanation. Maybe as a further improvement you could automatically calculate for input data points the correct bounding box to get a reliable interpolation.