반대 선형 보간법?
https://stackoverflow.com/questions/808441
-
03-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian문제
나는 네 2d 점,p0=(x0,y0),p1=(x1,y1),등등.을 형성하는 사각형.나의 케이스에서,쿼드는 직사각형이 아닌,하지만 적어도목하게 되도록 하는 것입니다.
p2 --- p3
| |
t | p |
| |
p0 --- p1
s
내가 사용하여 선형 보간.S 및 T 내에 있는[0..1]보간된 지점에 의해 제공됩니다.
bilerp(s,t) = t*(s*p3+(1-s)*p2) + (1-t)*(s*p1+(1-s)*p0)
여기에 문제가 발생할 수 있다.I have a2d 점 p 나가 알고 있는가 내부에 쿼드.을 찾고 싶은 s,t 는다는 점 때 사용하여 선형 보간.
거기에 간단한 공식을 반대 선형 보간법?
에 대한 감사 솔루션입니다.게시 내의 구현 Naaff 의 솔루션으로 wiki.
해결책
문제를 교차 문제로 생각하는 것이 가장 쉽다고 생각합니다. 지점 P가 P0, P1, P2 및 P3에 의해 정의 된 임의의 2D 이중선 표면과 교차하는 매개 변수 위치 (S, T)는 무엇입니까?
이 문제를 해결하는 데 필요한 접근법은 Tspauld의 제안과 유사합니다.
x와 y 측면에서 두 방정식으로 시작하십시오.
x = (1-s)*( (1-t)*x0 + t*x2 ) + s*( (1-t)*x1 + t*x3 )
y = (1-s)*( (1-t)*y0 + t*y2 ) + s*( (1-t)*y1 + t*y3 )
T를위한 해결 :
t = ( (1-s)*(x0-x) + s*(x1-x) ) / ( (1-s)*(x0-x2) + s*(x1-x3) )
t = ( (1-s)*(y0-y) + s*(y1-y) ) / ( (1-s)*(y0-y2) + s*(y1-y3) )
이제이 두 방정식을 서로 동일하게 설정하여 T를 제거 할 수 있습니다. 모든 것을 왼쪽으로 옮기고 단순화 우리는 양식의 방정식을 얻습니다.
A*(1-s)^2 + B*2s(1-s) + C*s^2 = 0
어디에:
A = (p0-p) X (p0-p2)
B = ( (p0-p) X (p1-p3) + (p1-p) X (p0-p2) ) / 2
C = (p1-p) X (p1-p3)
연산자 X를 사용하여 2D 크로스 제품 (예 : P0 X P1 = X0*Y1 -Y0*x1). 이 방정식을 2 차로 포맷했습니다 번스타인 다항식 이로 인해 물건이 더 우아하고 수치 적으로 안정적입니다. S에 대한 솔루션은이 방정식의 뿌리입니다. Bernstein 다항식에 대한 2 차 공식을 사용하여 뿌리를 찾을 수 있습니다.
s = ( (A-B) +- sqrt(B^2 - A*C) ) / ( A - 2*B + C )
2 차 공식은 +-로 인한 두 가지 답변을 제공합니다. 이중선 표면 내에 P가있는 솔루션에만 관심이 있다면 S가 0과 1 사이의 대답을 폐기 할 수 있습니다. S의 관점에서.
중요한 특별한 경우 하나를 지적해야합니다. 분모 a -2*b + c = 0이면 2 차 다항식은 실제로 선형입니다. 이 경우 S를 찾으려면 훨씬 간단한 방정식을 사용해야합니다.
s = A / (A-C)
AC = 0이 아니라면 정확히 하나의 솔루션을 제공합니다. a = c이면 두 가지 경우가 있습니다. 모두 s의 값은 p를 포함합니다 아니요 s에 대한 값은 p를 포함합니다.
다른 팁
하나는 사용할 수 있습니다 뉴턴의 방법 다음과 같은 비선형 방정식 시스템을 반복적으로 해결하려면 다음과 같습니다.
p = p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t.
두 가지 방정식이 있습니다 (방정식의 X와 Y 구성 요소 모두에서 평등)과 두 개의 미지의 (S 및 T)가 있습니다.
Newton의 방법을 적용하려면 필요합니다 잔여 r, 즉 :
r = p - (p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t),
그리고 자코비아 매트릭스 잔차를 구별하여 발견됩니다. 우리의 문제에 대해 Jacobian은 다음과 같습니다.
J(:,1) = dr/ds = -p0*(1-t) + p1*(1-t) + p2*t - p3*t (first column of matrix)
J(:,2) = dr/dt = -p0*(1-s) - p1*s + p2*s + p3*(1-s) (second column).
Newton의 방법을 사용하려면 초기 추측 (S, T)으로 시작한 다음 다음 반복을 몇 번 수행합니다.
(s,t) = (s,t) - J^-1 r,
S 및 T의 새로운 값으로 각 반복을 재 계산합니다. 각각의 반복에서 지배적 인 비용은 Jacobian의 역수를 잔차 (J^-1 r)에 적용하는데, J^-1 r)는 j를 계수 행렬로, R을 오른쪽으로 풀어야한다.
방법에 대한 직관 :
직관적으로, 사변형이 a 평행 사변형, 문제를 해결하는 것이 훨씬 쉬울 것입니다. 뉴턴의 방법은 연속 평행 사변형 근사치로 사변형 문제를 해결하는 것입니다. 각 반복마다 우리는
평행 사변형으로 사변형을 근사화하려면 지점 (S, T)에서 로컬 미분 정보를 사용하십시오.
선형 시스템을 해결하여 평행 사변형 근사치에서 올바른 (S, T) 값을 찾으십시오.
이 새로운 지점으로 뛰어 들어 반복하십시오.
방법의 장점 :
Newton 유형의 방법에 대해 예상대로 수렴은 매우 빠릅니다. 반복이 진행됨에 따라 방법이 실제 지점에 가까워지고 가까워 질뿐만 아니라 로컬 평행 사변형 근사치도 더 정확 해져 수렴 속도 자체가 증가합니다 (반복적 인 방법 전문 용어에서는 Newton의 방법이 차량으로 수렴). 실제로 이것은 올바른 숫자의 수가 대략 반복을 대략 두 배로 늘린다는 것을 의미합니다.
다음은 내가 한 임의의 시험에서 대표적인 반복 테이블 대 오류입니다 (코드 참조).
Iteration Error
1 0.0610
2 9.8914e-04
3 2.6872e-07
4 1.9810e-14
5 5.5511e-17 (machine epsilon)
두 번의 반복 후에 오류는 효과적으로 눈에 띄지 않고 대부분의 실용적 목적으로 충분히 적을 정도로 작고 5 번의 반복 후 결과는의의 한계에 정확합니다. 기계 정밀도.
반복 횟수 (예 : 대부분의 실제 애플리케이션에 대한 3 개의 반복 및 매우 높은 정확도가 필요한 경우 8 개의 반복)를 고정하면 알고리즘에는 매우 간단하고 간단한 논리가 있으며, 그에 대한 구조가 높아지는 구조가 있습니다. 성능 계산. 모든 종류의 특수 에지 케이스*를 확인하고 결과에 따라 다른 논리를 사용할 필요가 없습니다. 몇 가지 간단한 공식을 포함하는 루프에 대한 것입니다. 아래는 인터넷과 인터넷 주변의 다른 답변에서 발견 된 전통적인 공식 기반 접근법에 대한이 접근법의 몇 가지 장점을 강조합니다.
쉽게 코딩 할 수 있습니다. 루프를 위해 만들고 몇 가지 공식을 입력하십시오.
조건부 나 분기가 없습니다 (if/then), 일반적으로 훨씬 더 나은 것을 허용합니다. 파이프 라인 효율성.
정사각형 뿌리가없고 단 1 개 분할 만 반복 당 필요합니다 (잘 작성된 경우). 다른 모든 작업은 간단한 추가, 빼기 및 곱셈입니다. 정사각형 뿌리와 부문 은닉처 특정 아키텍처의 효율성 (특히 특정 임베디드 시스템에서). 실제로, 당신이 어떻게 후드를 보면 제곱근 그리고 부서 실제로 현대 프로그래밍 언어에 의해 계산되며, 둘 다 Newton의 방법의 변형, 때로는 하드웨어 및 아키텍처에 따라 소프트웨어에 사용됩니다.
쉽게 벡터화 할 수 있습니다 한 번에 많은 수의 사변형으로 배열에서 작동합니다. 이 작업을 수행하는 방법의 예는 아래의 벡터화 된 코드를 참조하십시오.
임의의 치수로 확장됩니다. 이 알고리즘은 여러 차원 (2D, 3D, 4D, ...)에서 다차원 보간을 역전시키는 간단한 방법으로 확장됩니다. 아래에 3D 버전을 포함 시켰으며 간단한 재귀 버전을 작성하거나 자동 차별화 라이브러리를 사용하여 N 차이로 이동하는 것을 상상할 수 있습니다. Newton의 방법은 일반적으로 치수 독립적 수렴 속도를 나타냅니다. 원칙적 으로이 방법은 현재 하드웨어에서 수천 치수 (!)로 확장 할 수 있어야합니다 (그 후 N-By-N 행렬 j를 구성하고 해결하는 것은 아마도 제한적 일 것입니다. 요인).
물론 이러한 것들 대부분은 하드웨어, 컴파일러 및 기타 여러 요인에 따라 다르므로 마일리지는 다를 수 있습니다.
암호:
어쨌든 여기 내 Matlab 코드가 있습니다.
기본 2D 버전 :
function q = bilinearInverse(p,p1,p2,p3,p4,iter)
%Computes the inverse of the bilinear map from [0,1]^2 to the convex
% quadrilateral defined by the ordered points p1 -> p2 -> p3 -> p4 -> p1.
%Uses Newton's method. Inputs must be column vectors.
q = [0.5; 0.5]; %initial guess
for k=1:iter
s = q(1);
t = q(2);
r = p1*(1-s)*(1-t) + p2*s*(1-t) + p3*s*t + p4*(1-s)*t - p;%residual
Js = -p1*(1-t) + p2*(1-t) + p3*t - p4*t; %dr/ds
Jt = -p1*(1-s) - p2*s + p3*s + p4*(1-s); %dr/dt
J = [Js,Jt];
q = q - J\r;
q = max(min(q,1),0);
end
end
예제 사용 :
% Test_bilinearInverse.m
p1=[0.1;-0.1];
p2=[2.2;-0.9];
p3=[1.75;2.3];
p4=[-1.2;1.1];
q0 = rand(2,1);
s0 = q0(1);
t0 = q0(2);
p = p1*(1-s0)*(1-t0) + p2*s0*(1-t0) + p3*s0*t0 + p4*(1-s0)*t0;
iter=5;
q = bilinearInverse(p,p1,p2,p3,p4,iter);
err = norm(q0-q);
disp(['Error after ',num2str(iter), ' iterations: ', num2str(err)])
예제 출력 :
>> test_bilinearInverse
Error after 5 iterations: 1.5701e-16
빠른 벡터화 된 2D 버전 :
function [ss,tt] = bilinearInverseFast(px,py, p1x,p1y, p2x,p2y, p3x,p3y, p4x,p4y, iter)
%Computes the inverse of the bilinear map from [0,1]^2 to the convex
% quadrilateral defined by the ordered points p1 -> p2 -> p3 -> p4 -> p1,
% where the p1x is the x-coordinate of p1, p1y is the y-coordinate, etc.
% Vectorized: if you have a lot of quadrilaterals and
% points to interpolate, then p1x(k) is the x-coordinate of point p1 on the
% k'th quadrilateral, and so forth.
%Uses Newton's method. Inputs must be column vectors.
ss = 0.5 * ones(length(px),1);
tt = 0.5 * ones(length(py),1);
for k=1:iter
r1 = p1x.*(1-ss).*(1-tt) + p2x.*ss.*(1-tt) + p3x.*ss.*tt + p4x.*(1-ss).*tt - px;%residual
r2 = p1y.*(1-ss).*(1-tt) + p2y.*ss.*(1-tt) + p3y.*ss.*tt + p4y.*(1-ss).*tt - py;%residual
J11 = -p1x.*(1-tt) + p2x.*(1-tt) + p3x.*tt - p4x.*tt; %dr/ds
J21 = -p1y.*(1-tt) + p2y.*(1-tt) + p3y.*tt - p4y.*tt; %dr/ds
J12 = -p1x.*(1-ss) - p2x.*ss + p3x.*ss + p4x.*(1-ss); %dr/dt
J22 = -p1y.*(1-ss) - p2y.*ss + p3y.*ss + p4y.*(1-ss); %dr/dt
inv_detJ = 1./(J11.*J22 - J12.*J21);
ss = ss - inv_detJ.*(J22.*r1 - J12.*r2);
tt = tt - inv_detJ.*(-J21.*r1 + J11.*r2);
ss = min(max(ss, 0),1);
tt = min(max(tt, 0),1);
end
end
속도를 위해이 코드는 암시 적으로 2x2 행렬의 역수에 대한 다음 공식을 사용합니다.
[a,b;c,d]^-1 = (1/(ad-bc))[d, -b; -c, a]
예제 사용 :
% test_bilinearInverseFast.m
n_quads = 1e6; % 1 million quads
iter = 8;
% Make random quadrilaterals, ensuring points are ordered convex-ly
n_randpts = 4;
pp_xx = zeros(n_randpts,n_quads);
pp_yy = zeros(n_randpts,n_quads);
disp('Generating convex point ordering (may take some time).')
for k=1:n_quads
while true
p_xx = randn(4,1);
p_yy = randn(4,1);
conv_inds = convhull(p_xx, p_yy);
if length(conv_inds) == 5
break
end
end
pp_xx(1:4,k) = p_xx(conv_inds(1:end-1));
pp_yy(1:4,k) = p_yy(conv_inds(1:end-1));
end
pp1x = pp_xx(1,:);
pp1y = pp_yy(1,:);
pp2x = pp_xx(2,:);
pp2y = pp_yy(2,:);
pp3x = pp_xx(3,:);
pp3y = pp_yy(3,:);
pp4x = pp_xx(4,:);
pp4y = pp_yy(4,:);
% Make random interior points
ss0 = rand(1,n_quads);
tt0 = rand(1,n_quads);
ppx = pp1x.*(1-ss0).*(1-tt0) + pp2x.*ss0.*(1-tt0) + pp3x.*ss0.*tt0 + pp4x.*(1-ss0).*tt0;
ppy = pp1y.*(1-ss0).*(1-tt0) + pp2y.*ss0.*(1-tt0) + pp3y.*ss0.*tt0 + pp4y.*(1-ss0).*tt0;
pp = [ppx; ppy];
% Run fast inverse bilinear interpolation code:
disp('Running inverse bilinear interpolation.')
tic
[ss,tt] = bilinearInverseFast(ppx,ppy, pp1x,pp1y, pp2x,pp2y, pp3x,pp3y, pp4x,pp4y, 10);
time_elapsed = toc;
disp(['Number of quadrilaterals: ', num2str(n_quads)])
disp(['Inverse bilinear interpolation took: ', num2str(time_elapsed), ' seconds'])
err = norm([ss0;tt0] - [ss;tt],'fro')/norm([ss0;tt0],'fro');
disp(['Error: ', num2str(err)])
예제 출력 :
>> test_bilinearInverseFast
Generating convex point ordering (may take some time).
Running inverse bilinear interpolation.
Number of quadrilaterals: 1000000
Inverse bilinear interpolation took: 0.5274 seconds
Error: 8.6881e-16
3D 버전 :
수렴 진행 상황을 표시하는 일부 코드가 포함되어 있습니다.
function ss = trilinearInverse(p, p1,p2,p3,p4,p5,p6,p7,p8, iter)
%Computes the inverse of the trilinear map from [0,1]^3 to the box defined
% by points p1,...,p8, where the points are ordered consistent with
% p1~(0,0,0), p2~(0,0,1), p3~(0,1,0), p4~(1,0,0), p5~(0,1,1),
% p6~(1,0,1), p7~(1,1,0), p8~(1,1,1)
%Uses Gauss-Newton method. Inputs must be column vectors.
tol = 1e-9;
ss = [0.5; 0.5; 0.5]; %initial guess
for k=1:iter
s = ss(1);
t = ss(2);
w = ss(3);
r = p1*(1-s)*(1-t)*(1-w) + p2*s*(1-t)*(1-w) + ...
p3*(1-s)*t*(1-w) + p4*(1-s)*(1-t)*w + ...
p5*s*t*(1-w) + p6*s*(1-t)*w + ...
p7*(1-s)*t*w + p8*s*t*w - p;
disp(['k= ', num2str(k), ...
', residual norm= ', num2str(norm(r)),...
', [s,t,w]= ',num2str([s,t,w])])
if (norm(r) < tol)
break
end
Js = -p1*(1-t)*(1-w) + p2*(1-t)*(1-w) + ...
-p3*t*(1-w) - p4*(1-t)*w + ...
p5*t*(1-w) + p6*(1-t)*w + ...
-p7*t*w + p8*t*w;
Jt = -p1*(1-s)*(1-w) - p2*s*(1-w) + ...
p3*(1-s)*(1-w) - p4*(1-s)*w + ...
p5*s*(1-w) - p6*s*w + ...
p7*(1-s)*w + p8*s*w;
Jw = -p1*(1-s)*(1-t) - p2*s*(1-t) + ...
-p3*(1-s)*t + p4*(1-s)*(1-t) + ...
-p5*s*t + p6*s*(1-t) + ...
p7*(1-s)*t + p8*s*t;
J = [Js,Jt,Jw];
ss = ss - J\r;
end
end
예제 사용 :
%test_trilinearInverse.m
h = 0.25;
p1 = [0;0;0] + h*randn(3,1);
p2 = [0;0;1] + h*randn(3,1);
p3 = [0;1;0] + h*randn(3,1);
p4 = [1;0;0] + h*randn(3,1);
p5 = [0;1;1] + h*randn(3,1);
p6 = [1;0;1] + h*randn(3,1);
p7 = [1;1;0] + h*randn(3,1);
p8 = [1;1;1] + h*randn(3,1);
s0 = rand;
t0 = rand;
w0 = rand;
p = p1*(1-s0)*(1-t0)*(1-w0) + p2*s0*(1-t0)*(1-w0) + ...
p3*(1-s0)*t0*(1-w0) + p4*(1-s0)*(1-t0)*w0 + ...
p5*s0*t0*(1-w0) + p6*s0*(1-t0)*w0 + ...
p7*(1-s0)*t0*w0 + p8*s0*t0*w0;
ss = trilinearInverse(p, p1,p2,p3,p4,p5,p6,p7,p8);
disp(['error= ', num2str(norm(ss - [s0;t0;w0]))])
예제 출력 :
test_trilinearInverse
k= 1, residual norm= 0.38102, [s,t,w]= 0.5 0.5 0.5
k= 2, residual norm= 0.025324, [s,t,w]= 0.37896 0.59901 0.17658
k= 3, residual norm= 0.00037108, [s,t,w]= 0.40228 0.62124 0.15398
k= 4, residual norm= 9.1441e-08, [s,t,w]= 0.40218 0.62067 0.15437
k= 5, residual norm= 3.3548e-15, [s,t,w]= 0.40218 0.62067 0.15437
error= 4.8759e-15
역수가 양의 볼륨을 갖는 경우 역 다중 리어 보간이 잘 정의되기 때문에 입력 지점의 순서에주의를 기울여야합니다. 3D에서는 모양을 그 자체로 내면 내부를 바꾸는 점을 선택하는 것이 훨씬 쉽습니다.
다음은 커뮤니티 위키로 Naaff의 솔루션을 구현했습니다. 다시 한 번 감사드립니다.
이것은 C 구현이지만 C ++에서 작동해야합니다. 퍼즈 테스트 기능이 포함됩니다.
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
int equals( double a, double b, double tolerance )
{
return ( a == b ) ||
( ( a <= ( b + tolerance ) ) &&
( a >= ( b - tolerance ) ) );
}
double cross2( double x0, double y0, double x1, double y1 )
{
return x0*y1 - y0*x1;
}
int in_range( double val, double range_min, double range_max, double tol )
{
return ((val+tol) >= range_min) && ((val-tol) <= range_max);
}
/* Returns number of solutions found. If there is one valid solution, it will be put in s and t */
int inverseBilerp( double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double x, double y, double* sout, double* tout, double* s2out, double* t2out )
{
int t_valid, t2_valid;
double a = cross2( x0-x, y0-y, x0-x2, y0-y2 );
double b1 = cross2( x0-x, y0-y, x1-x3, y1-y3 );
double b2 = cross2( x1-x, y1-y, x0-x2, y0-y2 );
double c = cross2( x1-x, y1-y, x1-x3, y1-y3 );
double b = 0.5 * (b1 + b2);
double s, s2, t, t2;
double am2bpc = a-2*b+c;
/* this is how many valid s values we have */
int num_valid_s = 0;
if ( equals( am2bpc, 0, 1e-10 ) )
{
if ( equals( a-c, 0, 1e-10 ) )
{
/* Looks like the input is a line */
/* You could set s=0.5 and solve for t if you wanted to */
return 0;
}
s = a / (a-c);
if ( in_range( s, 0, 1, 1e-10 ) )
num_valid_s = 1;
}
else
{
double sqrtbsqmac = sqrt( b*b - a*c );
s = ((a-b) - sqrtbsqmac) / am2bpc;
s2 = ((a-b) + sqrtbsqmac) / am2bpc;
num_valid_s = 0;
if ( in_range( s, 0, 1, 1e-10 ) )
{
num_valid_s++;
if ( in_range( s2, 0, 1, 1e-10 ) )
num_valid_s++;
}
else
{
if ( in_range( s2, 0, 1, 1e-10 ) )
{
num_valid_s++;
s = s2;
}
}
}
if ( num_valid_s == 0 )
return 0;
t_valid = 0;
if ( num_valid_s >= 1 )
{
double tdenom_x = (1-s)*(x0-x2) + s*(x1-x3);
double tdenom_y = (1-s)*(y0-y2) + s*(y1-y3);
t_valid = 1;
if ( equals( tdenom_x, 0, 1e-10 ) && equals( tdenom_y, 0, 1e-10 ) )
{
t_valid = 0;
}
else
{
/* Choose the more robust denominator */
if ( fabs( tdenom_x ) > fabs( tdenom_y ) )
{
t = ( (1-s)*(x0-x) + s*(x1-x) ) / ( tdenom_x );
}
else
{
t = ( (1-s)*(y0-y) + s*(y1-y) ) / ( tdenom_y );
}
if ( !in_range( t, 0, 1, 1e-10 ) )
t_valid = 0;
}
}
/* Same thing for s2 and t2 */
t2_valid = 0;
if ( num_valid_s == 2 )
{
double tdenom_x = (1-s2)*(x0-x2) + s2*(x1-x3);
double tdenom_y = (1-s2)*(y0-y2) + s2*(y1-y3);
t2_valid = 1;
if ( equals( tdenom_x, 0, 1e-10 ) && equals( tdenom_y, 0, 1e-10 ) )
{
t2_valid = 0;
}
else
{
/* Choose the more robust denominator */
if ( fabs( tdenom_x ) > fabs( tdenom_y ) )
{
t2 = ( (1-s2)*(x0-x) + s2*(x1-x) ) / ( tdenom_x );
}
else
{
t2 = ( (1-s2)*(y0-y) + s2*(y1-y) ) / ( tdenom_y );
}
if ( !in_range( t2, 0, 1, 1e-10 ) )
t2_valid = 0;
}
}
/* Final cleanup */
if ( t2_valid && !t_valid )
{
s = s2;
t = t2;
t_valid = t2_valid;
t2_valid = 0;
}
/* Output */
if ( t_valid )
{
*sout = s;
*tout = t;
}
if ( t2_valid )
{
*s2out = s2;
*t2out = t2;
}
return t_valid + t2_valid;
}
void bilerp( double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double s, double t, double* x, double* y )
{
*x = t*(s*x3+(1-s)*x2) + (1-t)*(s*x1+(1-s)*x0);
*y = t*(s*y3+(1-s)*y2) + (1-t)*(s*y1+(1-s)*y0);
}
double randrange( double range_min, double range_max )
{
double range_width = range_max - range_min;
double rand01 = (rand() / (double)RAND_MAX);
return (rand01 * range_width) + range_min;
}
/* Returns number of failed trials */
int fuzzTestInvBilerp( int num_trials )
{
int num_failed = 0;
double x0, y0, x1, y1, x2, y2, x3, y3, x, y, s, t, s2, t2, orig_s, orig_t;
int num_st;
int itrial;
for ( itrial = 0; itrial < num_trials; itrial++ )
{
int failed = 0;
/* Get random positions for the corners of the quad */
x0 = randrange( -10, 10 );
y0 = randrange( -10, 10 );
x1 = randrange( -10, 10 );
y1 = randrange( -10, 10 );
x2 = randrange( -10, 10 );
y2 = randrange( -10, 10 );
x3 = randrange( -10, 10 );
y3 = randrange( -10, 10 );
/*x0 = 0, y0 = 0, x1 = 1, y1 = 0, x2 = 0, y2 = 1, x3 = 1, y3 = 1;*/
/* Get random s and t */
s = randrange( 0, 1 );
t = randrange( 0, 1 );
orig_s = s;
orig_t = t;
/* bilerp to get x and y */
bilerp( x0, y0, x1, y1, x2, y2, x3, y3, s, t, &x, &y );
/* invert */
num_st = inverseBilerp( x0, y0, x1, y1, x2, y2, x3, y3, x, y, &s, &t, &s2, &t2 );
if ( num_st == 0 )
{
failed = 1;
}
else if ( num_st == 1 )
{
if ( !(equals( orig_s, s, 1e-5 ) && equals( orig_t, t, 1e-5 )) )
failed = 1;
}
else if ( num_st == 2 )
{
if ( !((equals( orig_s, s , 1e-5 ) && equals( orig_t, t , 1e-5 )) ||
(equals( orig_s, s2, 1e-5 ) && equals( orig_t, t2, 1e-5 )) ) )
failed = 1;
}
if ( failed )
{
num_failed++;
printf("Failed trial %d\n", itrial);
}
}
return num_failed;
}
int main( int argc, char** argv )
{
int num_failed;
srand( 0 );
num_failed = fuzzTestInvBilerp( 100000000 );
printf("%d of the tests failed\n", num_failed);
getc(stdin);
return 0;
}
당신은 2D에서 일하기 때문에 당신의 bilerp 함수는 실제로 2 방정식, x의 경우 1, y의 경우 1입니다. 양식으로 다시 작성할 수 있습니다.
x = t * s * A.x + t * B.x + s * C.x + D.x
y = t * s * A.y + t * B.y + s * C.y + D.y
어디에:
A = p3 - p2 - p1 + p0
B = p2 - p0
C = p1 - p0
D = p0
첫 번째 방정식을 다시 작성하십시오 t 측면에서 s, 두 번째로 대체하고 해결하십시오 s.
이것은 내 구현입니다 ... Tfiniga보다 빠릅니다.
void invbilerp( float x, float y, float x00, float x01, float x10, float x11, float y00, float y01, float y10, float y11, float [] uv ){
// substition 1 ( see. derivation )
float dx0 = x01 - x00;
float dx1 = x11 - x10;
float dy0 = y01 - y00;
float dy1 = y11 - y10;
// substitution 2 ( see. derivation )
float x00x = x00 - x;
float xd = x10 - x00;
float dxd = dx1 - dx0;
float y00y = y00 - y;
float yd = y10 - y00;
float dyd = dy1 - dy0;
// solution of quadratic equations
float c = x00x*yd - y00y*xd;
float b = dx0*yd + dyd*x00x - dy0*xd - dxd*y00y;
float a = dx0*dyd - dy0*dxd;
float D2 = b*b - 4*a*c;
float D = sqrt( D2 );
float u = (-b - D)/(2*a);
// backsubstitution of "u" to obtain "v"
float v;
float denom_x = xd + u*dxd;
float denom_y = yd + u*dyd;
if( abs(denom_x)>abs(denom_y) ){ v = -( x00x + u*dx0 )/denom_x; }else{ v = -( y00y + u*dy0 )/denom_y; }
uv[0]=u;
uv[1]=v;
/*
// do you really need second solution ?
u = (-b + D)/(2*a);
denom_x = xd + u*dxd;
denom_y = yd + u*dyd;
if( abs(denom_x)>abs(denom_y) ){ v = -( x00x + u*dx0 )/denom_x; }else{ v2 = -( y00y + u*dy0 )/denom_y; }
uv[2]=u;
uv[3]=v;
*/
}
그리고 파생
// starting from bilinear interpolation
(1-v)*( (1-u)*x00 + u*x01 ) + v*( (1-u)*x10 + u*x11 ) - x
(1-v)*( (1-u)*y00 + u*y01 ) + v*( (1-u)*y10 + u*y11 ) - y
substition 1:
dx0 = x01 - x00
dx1 = x11 - x10
dy0 = y01 - y00
dy1 = y11 - y10
we get:
(1-v) * ( x00 + u*dx0 ) + v * ( x10 + u*dx1 ) - x = 0
(1-v) * ( y00 + u*dy0 ) + v * ( y10 + u*dy1 ) - y = 0
we are trying to extract "v" out
x00 + u*dx0 + v*( x10 - x00 + u*( dx1 - dx0 ) ) - x = 0
y00 + u*dy0 + v*( y10 - y00 + u*( dy1 - dy0 ) ) - y = 0
substition 2:
x00x = x00 - x
xd = x10 - x00
dxd = dx1 - dx0
y00y = y00 - y
yd = y10 - y00
dyd = dy1 - dy0
// much nicer
x00x + u*dx0 + v*( xd + u*dxd ) = 0
y00x + u*dy0 + v*( yd + u*dyd ) = 0
// this equations for "v" are used for back substition
v = -( x00x + u*dx0 ) / ( xd + u*dxd )
v = -( y00x + u*dy0 ) / ( yd + u*dyd )
// but for now, we eliminate "v" to get one eqution for "u"
( x00x + u*dx0 ) / ( xd + u*dxd ) = ( y00y + u*dy0 ) / ( yd + u*dyd )
put denominators to other side
( x00x + u*dx0 ) * ( yd + u*dyd ) = ( y00y + u*dy0 ) * ( xd + u*dxd )
x00x*yd + u*( dx0*yd + dyd*x00x ) + u^2* dx0*dyd = y00y*xd + u*( dy0*xd + dxd*y00y ) + u^2* dy0*dxd
// which is quadratic equation with these coefficients
c = x00x*yd - y00y*xd
b = dx0*yd + dyd*x00x - dy0*xd - dxd*y00y
a = dx0*dyd - dy0*dxd
P의 단일 값이 P의 단일 값이라면 P가 정사각형의 네 모서리의 최소값과 최대 값 사이에 있으면 일반적으로 단일 솔루션 (S, T)을 찾는 것은 불가능합니다. 이중선 대기업은 당신에게 그 가치를 줄 것입니다.
일반적으로 정사각형 내부에는 무한한 수의 솔루션 (S, T)이 있습니다. 그들은 정사각형을 통해 구부러진 (쌍곡선) 경로를 따라 누워있을 것입니다.
당신의 함수가 벡터 가치가있는 경우, 당신은 정사각형의 알려지지 않은 지점에 두 개의 알려진 값이 있습니까? 정사각형의 각 구석에 알려진 두 매개 변수의 알려진 값이 주어지면 솔루션이 존재할 수 있지만, 이에 대한 보장은 없습니다. 우리는 이것을 두 가지 별도의 독립적 인 문제로 볼 수 있습니다. 그들 각각에 대한 해결책은 정사각형을 통해 쌍곡선 윤곽선을 따라 놓을 것입니다. 윤곽 쌍이 사각형 내부에 교차하면 용액이 존재합니다. 그들이 교차하지 않으면 해결책이 없습니다.
또한 문제를 해결하기 위해 간단한 공식이 존재하는지 묻습니다. 죄송하지만 실제로는 보지 못했습니다. 내가 말했듯이, 곡선은 쌍곡선입니다.
한 가지 해결책은 다른 보간 방법으로 전환하는 것입니다. 따라서 이중선 대신 정사각형을 삼각형으로 나눕니다. 각 삼각형 내에서 우리는 이제 진정한 선형 보간을 사용할 수 있습니다. 이제 각 삼각형 내에서 2 개의 미지의 2 방정식의 선형 시스템을 해결할 수 있습니다. 상응하는 부분 선형 윤곽선이 공동으로 발생하는 희귀 한 퇴화 케이스를 제외하고 각 삼각형에는 하나의 솔루션이있을 수 있습니다.
의 일부 응답에 약간의 해석의 질문입니다.ie.그들은 당신이어 값 의 알 수 없는 삽입 기능보다는 보정 위치로 p(x,y)내부에 쿼드는 당신을 찾으려면(s,t)의 좌표.이것은 단순한 문제가 되도록 보장하는 솔루션의 교차로 두 줄을 통해 쿼드.
라인 중 하나를 통해 세그먼트 p0p1 및 p2p3,다른 것이다 p0p2 및 p1p3 와 유사한 축 정렬된 경우입니다.이러한 라인은 유일하게 정의의 위치로 p(x,y),그리고 물론에서 교차합니다.
고려의 라인을 통해 p0p1 및 p2p3,우리가 상상할 수 있는 가족의 이러한선은 각기 다른 s-값 우리가 선택,각 절단 쿼드에서 다른 폭입니다.면 우리는 해결 s-값을,우리를 찾을 수 있습 두 끝점을 설정하여 t=0,t=1.
그래서 가정한다 라인의 형식은:y=a0*x+b0
그런 다음 우리가 알고 있는 두 끝점 이 라인의 경우,우리는 우리 주 s value.그들은:
(1-s)p0+(s)p1
(1-s)p2+(s)p3
이러한 두 끝점,우리는 결정할 수 있는 가족의 라인을 연결하여으로 그들을 방정식에 대한 라인과 해결을 위한 a0 및 b0 로서의 기능을 s.설정 s-가치를 제공에 대한 공식을 특정한 라인입니다.모든 지금 우리에게 필요한 그 밖으로 하는 라인에서 가족을 우리의점 p(x,y).단순히 연결하의 좌표를 p(x,y)으로 우리 선을 수식,그리고 우리는 그에 대한 해결의 대상 값 s.
해당 접근 할 수 있을 찾 t 니다.
글쎄, p가 2d 포인트라면, 쉽게 얻을 수 있습니다. 이 경우, s는 t에서 사변형의 총 폭의 분수 성분이다.
그러나 P가 스칼라 인 경우, 이중선 보간 함수가 반드시 모 놀리 식 일 필요는 없기 때문에 반드시 가능하지는 않습니다.
