Большая интерполяция карты высот

Question

У меня есть vector<vector<double>> heightmap, который динамически загружается из CSV-файла данных GPS, чтобы быть около 4000x4000. Однако обеспечивает только 140 799 точек.

Создает карту в оттенках серого, как показано ниже:

I wi sh для интерполяции высот между всеми точками создать карту высоты области.

В приведенном ниже коде все известные точки будут смотреть в радиусе 10 м от точки, чтобы найти другие известные точки. Если найдена другая точка, она будет линейно интерполироваться между двумя точками. Интерполированные точки определяются как - высота и неустановленные значения определяются как -1337.

Этот подход невероятно медленный. Я уверен, что есть лучшие способы добиться этого.

bool run_interp = true;

bool interp_interp = false;

int counter = 0;

while (run_interp)
{
    for (auto x = 0; x < map.size(); x++)
    {
        for (auto y = 0; y < map.at(x).size(); y++)
        {
            const auto height = map.at(x).at(y);

            if (height == -1337) continue;

            if (!interp_interp && height < 0) continue;

            //Look in a 10m radius of a known value to see if there
            //Is another known value to linearly interp between
            //Set height to a negative if it has been interped

            const int radius = (1 / resolution) * 10;

            for (auto rxi = 0; rxi < radius * 2; rxi++)
            {
                //since we want to expand outwards
                const int rx = x + ((rxi % 2 == 0) ? rxi / 2 : -(rxi - 1) / 2);

                if (rx < 0 || rx >= map.size()) continue;

                for (auto ryi = 0; ryi < radius * 2; ryi++)
                {
                    const int ry = y + ((rxi % 2 == 0) ? rxi / 2 : -(rxi - 1) / 2);

                    if (ry < 0 || ry >= map.at(x).size()) continue;

                    const auto new_height = map.at(rx).at(ry);

                    if (new_height == -1337) continue;

                    //First go around we don't want to interp
                    //Interps
                    if (!interp_interp && new_height < 0) continue;

                    //We have found a known point within 10m
                    const auto delta = new_height - height;
                    const auto distance = sqrt((rx- x) * (rx - x)
                        + (ry - y) * (ry - y));
                    const auto angle = atan2(ry - y, rx - x);
                    const auto ratio = delta / distance;

                    //Backtrack from found point until we get to know point
                    for (auto radi = 0; radi < distance; radi++)
                    {
                        const auto new_x = static_cast<int>(x + radi * cos(angle));
                        const auto new_y = static_cast<int>(y + radi * sin(angle));

                        if (new_x < 0 || new_x >= map.size()) continue;
                        if (new_y < 0 || new_y >= map.at(new_x).size()) continue;

                        const auto interp_height = map.at(new_x).at(new_y);

                        //If it is a known height don't interp it
                        if (interp_height > 0)
                            continue;

                        counter++;

                        set_height(new_x, new_y, -interp_height);
                    }
                }
            }
        }
        std::cout << x << " " << counter << std::endl;;

    }

    if (interp_interp)
        run_interp = false;

    interp_interp = true;

}



set_height(const int x, const int y, const double height)
{
//First time data being set
if (map.at(x).at(y) == -1337)
{
    map.at(x).at(y) = height;
}
else // Data set already so average it
{
    //While this isn't technically correct and weights
    //Later data significantly more favourablily
    //It should be fine
    //TODO: fix it.
    map.at(x).at(y) += height;
    map.at(x).at(y) /= 2;
}

}

Answer 1

Если вы поместите точки в дерево kd, будет гораздо быстрее найти ближайшую точку (O (nlogn)).

Я не уверен, что это решит все ваши проблемы, но это начало.