Interpolation des données 2d qui sont constantes par morceaux sur les faces

Question

J'ai un maillage irrégulier qui est décrit par deux variables - un tableau de faces qui stocke les indices des sommets constituant chaque face, et un tableau vert qui stocke les coordonnées de chaque sommet . J'ai aussi une fonction qui est supposée être constante par morceaux sur chaque face, et elle est stockée sous la forme d'un tableau de valeurs par face.

Je cherche un moyen de construire une fonction f à partir de ces données. Quelque chose le long des lignes suivantes:

faces = [[0,1,2], [1,2,3], [2,3,4] ...] 
verts = [[0,0], [0,1], [1,0], [1,1],....] 
vals = [0.0, 1.0, 0.5, 3.0,....] 

f = interpolate(faces, verts, vals) 

f(0.2, 0.2) = 0.0 # point inside face [0,1,2] 
f(0.6, 0.6) = 1.0 # point inside face [1,2,3] 

La méthode manuelle d'évaluation f(x,y) serait de trouver la face correspondante que le point x,y se trouve dans, et retourner la valeur qui est stockée dans ce visage. Y at-il une fonction qui l'implémente déjà dans scipy (ou dans matlab)?

Answer 1

Il n'y a pas fonction intégrée dans MATLAB qui fera ce que vous voulez. Vous pourriez probablement construire votre propre algorithme en utilisant la fonction INPOLYGON en tant que suggested by Jonas, mais vous pourriez être en mesure de créer vous-même une implémentation plus rapide en utilisant des algorithmes standards pour trouver si un point est dans un polygone. Il y a quelque temps, j'ai écrit mon propre code pour trouver les points d'intersection entre un segment de ligne et un ensemble de surfaces triangulaires en 3D, et j'ai trouvé que this softsurfer link était le plus utile pour implémenter l'algorithme. Mon cas était plus compliqué que le tien.Puisque vous travaillez en 2-D, vous pouvez ignorer la première section du lien sur la recherche du point où le segment croise le plan du triangle.

J'ai inclus une version simplifiée de mon code MATLAB ci-dessous pour que vous puissiez l'utiliser. La fonction interpolate prendra en entrée les matrices faces, vertices et values et renverra une poignée de fonction f pouvant être évaluée à un point donné (x, y) pour obtenir la valeur par morceaux dans le triangle de délimitation. Voici quelques caractéristiques de ce code:

• Le traitement qui sera fait lorsque vous évaluez f est contenu dans le nested functionevaluate_function. Cette fonction a accès aux autres variables dans interpolate, donc un certain nombre de variables nécessaires pour le test en triangle sont précalculés afin que evaluate_function s'exécute aussi vite que possible.
• Dans les cas où vous avez beaucoup de triangles, tester si votre point est à l'intérieur de tous peut être coûteux. Par conséquent, le code trouve d'abord des triangles qui se trouvent dans un rayon donné (c'est-à-dire la longueur de la jambe la plus longue des triangles) de votre point. Seuls ces triangles proches sont testés pour voir si le point est en eux.
• Si un point ne se trouve pas dans une région triangulaire, la valeur NaN est renvoyée par f.

Il y a des choses qui ne sont pas inclus dans le code, que vous pouvez ajouter en fonction de ce que vous utilisez pour:

• Contrôle des entrées: Le code suppose actuellement que faces est une matrice N-by-3, vertices est une matrice M-par-2, et values est une longueur N vecteur. Vous voudrez probablement ajouter une vérification d'erreur pour vous assurer que les entrées sont conformes à ces exigences et émettre une erreur indiquant quand un ou plusieurs d'entre eux sont incorrects.
• contrôle triangle dégénérée: Il est possible que l'un ou plusieurs des triangles définis par vos peut-être dégénéré (à savoir qu'ils peuvent avoir une zone de 0) entrées faces et vertices. Cela se produit lorsque deux ou plusieurs des sommets des triangles sont le même point exact, ou lorsque les trois sommets du triangle se trouvent sur une ligne droite. Vous voudrez probablement ajouter un contrôle qui ignorera ces triangles quand il s'agit d'évaluer f.
• Manipulation des bords: Il est possible qu'un point puisse se retrouver sur le bord de deux ou plusieurs régions triangulaires. Vous devrez donc décider quelle valeur prendra le point (c'est-à-dire la plus grande des valeurs faciales, la moyenne des valeurs faciales, etc.). Pour les cas de bords comme celui-ci, le code ci-dessous choisira automatiquement la valeur du visage qui est plus proche du début de la liste des visages dans votre variable faces.

Enfin, voici le code:

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 

Answer 2

Cela ne me semble pas tellement comme interpolation que de trouver quelle face triangulaire le point est à l'intérieur. Vérifiez this site sur un moyen de tester chaque visage triangulaire. Votre fonction déterminera simplement la face à l'intérieur et retournera la valeur correspondante. Bien sûr, si vous avez beaucoup de visages ou si vous faites beaucoup de choses, alors vous voudrez trouver des moyens d'optimiser cela (stockez les points + et - les plus éloignés pour chaque triangle dans les directions x et y et stockez par exemple, si le point n'est pas dans cette zone de délimitation, alors vous pouvez aussi ne pas vérifier s'il s'agit de l'intérieur du triangle). Je doute vraiment que vous trouverez quelque chose qui est intégré à Matlab ou Scipy pour faire exactement ce que vous voulez, mais je peux me tromper.

Answer 3

vous pouvez utiliser le module de pont CGAL-python; Si je me souviens bien, CGAL fournit des fonctionnalités pour la recherche de triangle. Cependant, cela fonctionne probablement plus efficacement si vous construisez la triangulation de façon incrémentielle en utilisant leur représentation intégrée. Pour un problème rapide, vous pouvez simplement trouver le sommet de maillage le plus proche du point de la requête, soit par un diagramme de Voronoi (la fonctionnalité dans Matlab n'est pas géniale), soit, pour une seule requête, calculer tous les distances, et trouver le minimum, puis rechercher tous les triangles avec ce sommet.

Answer 4

Matlab a la fonction intégrée inpolygon, qui vous permet de tester si vous êtes dans un triangle. Je ne connais pas de fonction qui identifierait à l'intérieur de quel visage vous êtes. Si vous deviez écrire une telle fonction, je testerais d'abord le sommet le plus proche de votre point, puis j'évaluerais inpolygon pour toutes les faces qui partagent le sommet jusqu'à ce que vous trouviez la correspondance. Cela devrait être raisonnablement rapide.

Answer 5

Jetez un oeil à matplotlib.delaunay.interpolate, qui est un wrapper bien documenté pour le code C.
(Cependant class LinearInterpolator il dit « À l'heure actuelle, seuls sont pris en charge régulièrement des grilles rectangulaires pour l'interpolation. »)