Интерполяция двухмерных данных, кусочно-постоянных на гранях

В MATLAB нет встроенной функции, которая будет делать то, что вы хотите. Вероятно, вы могли бы создать свой собственный алгоритм, используя функцию INPOLYGON как предложено Джонасом, но вы можете создать более быструю реализацию самостоятельно, используя некоторые стандартные алгоритмы для определения того, находится ли точка внутри многоугольника.

Некоторое время назад я написал свой собственный код для поиска точек пересечения между отрезком линии и набором треугольных поверхностей в 3D и нашел эта ссылка на софтсерфера как наиболее полезная для реализации алгоритма. Мой случай был сложнее вашего. Поскольку вы работаете в 2D, вы можете игнорировать первый раздел ссылки о поиске точки, в которой отрезок пересекает плоскость треугольника.

Я включил упрощенную версию моего кода MATLAB ниже для вашего использования. Функция interpolate примет ваши матрицы faces, vertices и values в качестве входных данных и вернет дескриптор функции f, который можно вычислить в заданной точке (x,y), чтобы получить кусочное значение в ограничивающем треугольнике. Вот несколько особенностей этого кода:

• Обработка, которая будет выполнена при оценке f, содержится в вложенная функция evaluate_function. Эта функция имеет доступ к другим переменным в interpolate, поэтому ряд переменных, необходимых для теста внутри треугольника, предварительно вычисляются, чтобы evaluate_function работало как можно быстрее.
• В случаях, когда у вас много треугольников, проверка того, находится ли ваша точка внутри всех из них, может быть дорогостоящей. Таким образом, код сначала находит треугольники, которые находятся в пределах заданного радиуса (т. е. длины самой длинной стороны треугольника) от вашей точки. Только эти соседние треугольники проверяются, чтобы увидеть, находится ли точка внутри них.
• Если точка не попадает ни в одну треугольную область, значение NaN возвращается функцией f.

Есть некоторые вещи, которые не включены в код, которые вы можете добавить в зависимости от того, для чего вы его используете:

• Проверка входных данных. В настоящее время код предполагает, что faces – это матрица размера N на 3, _12  – матрица размером M на 2, а values – вектор длины N. Вы, вероятно, захотите добавить проверку ошибок, чтобы убедиться, что входные данные соответствуют этим требованиям, и выдать ошибку, указывающую, когда один или несколько из них неверны.
• Проверка вырожденного треугольника. Возможно, что один или несколько треугольников, заданных вашими входными данными faces и vertices, могут быть вырожденными (т. е. иметь площадь, равную 0). Это происходит, когда две или более вершин треугольника находятся в одной точке или когда все три вершины треугольника лежат на одной прямой. Вероятно, вы захотите добавить проверку, которая будет игнорировать такие треугольники при оценке f.
• Обработка крайних случаев. Возможно, точка окажется на краю двух или более треугольных областей. Поэтому вам придется решить, какое значение примет точка (т. е. наибольшее из номиналов, среднее номиналов и т. д.). Для подобных крайних случаев приведенный ниже код автоматически выберет значение лица, которое находится ближе к началу списка лиц в вашей переменной faces.

Наконец, вот код:

function f = interpolate(faces,vertices,values)

  %# Precompute some data (helps increase speed):

  triVertex = vertices(faces(:,2),:);              %# Triangles main vertices
  triLegLeft = vertices(faces(:,1),:)-triVertex;   %# Triangles left legs
  triLegRight = vertices(faces(:,3),:)-triVertex;  %# Triangles right legs
  C1 = sum(triLegLeft.*triLegRight,2);  %# Dot product of legs
  C2 = sum(triLegLeft.^2,2);            %# Squared length of left leg
  C3 = sum(triLegRight.^2,2);           %# Squared length of right leg
  triBoundary = max(C2,C3);             %# Squared radius of triangle boundary
  scale = C1.^2-C2.*C3;
  C1 = C1./scale;
  C2 = C2./scale;
  C3 = C3./scale;

  %# Return a function handle to the nested function:

  f = @evaluate_function;

  %# The nested evaluation function:

  function val = evaluate_function(x,y)

    w = [x-triVertex(:,1) y-triVertex(:,2)];
    nearIndex = find(sum(w.^2,2) <= triBoundary);  %# Find nearby triangles
    if isempty(nearIndex)
      val = NaN;           %# Return NaN if no triangles are nearby
      return
    end
    w = w(nearIndex,:);
    wdotu = sum(w.*triLegLeft(nearIndex,:),2);
    wdotv = sum(w.*triLegRight(nearIndex,:),2);
    C = C1(nearIndex);
    s = C.*wdotv-C3(nearIndex).*wdotu;
    t = C.*wdotu-C2(nearIndex).*wdotv;
    inIndex = find((s >= 0) & (t >= 0) & (s+t <= 1),1);
    if isempty(inIndex)
      val = NaN;         %# Return NaN if point is outside all triangles
    else
      val = values(nearIndex(inIndex));
    end

  end

end

Для меня это звучит не столько как интерполяция, сколько как определение того, внутри какой треугольной грани находится точка. Посетите этот сайт, чтобы проверить каждую треугольную грань. Ваша функция просто определит, какая грань находится внутри, и вернет соответствующее значение. Конечно, если у вас много граней или вы делаете это много раз, вам нужно будет найти способы оптимизировать это (сохранить самые дальние точки + и - для каждого треугольника в обоих направлениях x и y и сохранить например, с лицом. Если точка не находится внутри этой ограничивающей рамки, то вы можете не проверять, находится ли она внутри треугольника).

Я очень сомневаюсь, что вы найдете что-то встроенное в Matlab или scipy, чтобы делать именно то, что вы хотите, но я могу ошибаться.

вы можете использовать модуль моста CGAL-python; Если я правильно помню, CGAL предоставляет функциональность для поиска треугольников. Однако это, вероятно, работает наиболее эффективно, если вы строите триангуляцию постепенно, используя их встроенное представление. На скорую руку можно просто найти ближайшую вершину меша к точке запроса либо с помощью диаграммы Вороного (функциональность для этого в Matlab невелика), либо, для одного запроса, вычислив все расстояния и найти минимум, а затем найти все треугольники с этой вершиной.

В Matlab есть встроенная функция inpolygon, которая позволяет проверить, находитесь ли вы внутри треугольника. Я не знаю функции, которая определяла бы, внутри какого лица вы находитесь.

Если бы вы написали такую ​​функцию, я бы сначала проверил, к какой вершине ваша точка ближе всего, а затем оценил inpolygon для всех граней, которые имеют общую вершину, пока вы не найдете совпадение. Это должно быть достаточно быстро.

Взгляните на matplotlib.delaunay.interpolate, хорошо задокументированную оболочку для кода C.
(Однако class LinearInterpolator там написано: "В настоящий момент для интерполяции поддерживаются только обычные прямоугольные сетки".)