معكوس خطيه الاستيفاء?
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)
هنا المشكلة..لدي 2d النقطة p التي أعلم أنها داخل رباعية.أريد أن أجد s,t من شأنها أن تعطي لي هذه النقطة عند استخدام الصورة الاستيفاء.
هل هناك صيغة بسيطة عكس الصورة الاستيفاء?
شكرا على الحلول.نشرت لي تنفيذ Naaff الحل كما ويكي.
المحلول
وأعتقد أنه من الأسهل أن نفكر في المشكلة على النحو مشكلة تقاطع: ما هو الموقع المعلمة (ق، ر) حيث يتقاطع النقطة ص التعسفي سطح 2D المترابط يحددها P0، P1، P2 و P3
وهذا النهج سوف تأخذ في حل هذه المشكلة مشابه للاقتراح 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 = ( (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) )
ويمكننا الآن تعيين هذه المعادلات اثنين متساوية مع بعضها البعض لقضاء ر. نقل كل شيء إلى الجانب الأيسر وتبسيط نحصل على معادلة النموذج:
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). لقد تنسيق هذه المعادلة باعتباره متعدد الحدود بيرنشتاين متعدد الحدود لأن ذلك يجعل أشياء أكثر أناقة وأكثر استقرارا من الناحية العددية. الحلول إلى s هي جذور هذه المعادلة. يمكن أن نجد جذور باستخدام الصيغة التربيعية لمتعددو الحدود بيرنشتاين:
s = ( (A-B) +- sqrt(B^2 - A*C) ) / ( A - 2*B + C )
والصيغة التربيعية يعطي إجابتين نتيجة ل+ -. إذا كنت ترغب فقط في حلول أين يكمن ص داخل سطح المترابط ثم يمكنك تجاهل أي إجابة حيث الصورة ليست بين 0 و 1. للعثور ر، ببساطة استبدال الصورة مرة أخرى إلى واحد من اثنين من المعادلات أعلاه حيث حللنا لر من حيث الصورة.
وأود أن أشير إلى حالة واحدة استثنائية هامة. إذا كان القاسم A - 2 * B + C = 0 ثم متعدد الحدود من الدرجة الثانية هو في الواقع الخطي. في هذه الحالة، يجب استخدام المعادلة أبسط من ذلك بكثير للعثور على الصورة:
s = A / (A-C)
وهذا سيعطي لك بالضبط حل واحد، إلا إذا AC = 0. إذا A = C ثم لديك حالتين: A = C = 0 وسائل <م> جميع قيم الصورة احتواء ص، وإلا <م> لا تحتوي على قيم الصورة ع.
نصائح أخرى
يمكن للمرء أن استخدام طريقة نيوتن إلى تكرارا حل التالية الخطية نظام من المعادلات:
p = p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t.
علما أن هناك اثنين من المعادلات (المساواة في كل من x و y مكونات المعادلة) و اثنين من المجاهيل (s و t).
تطبيق طريقة نيوتن ، نحن في حاجة إلى المتبقية r ، الذي هو:
r = p - (p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t),
و Jacobian مصفوفة ي ، التي وجدت قبل التفريق المتبقية.لدينا مشكلة 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).
استخدام طريقة نيوتن واحد يبدأ مع تخمين الأولي (s,t), ثم ينفذ التالية التكرار حفنة من الأوقات:
(s,t) = (s,t) - J^-1 r,
إعادة حساب J و r كل التكرار مع قيم جديدة من s و t.في كل التكرار السائد التكلفة في تطبيق معكوس Jacobian المتبقية (J^-1 ص), من خلال حل 2x2 نظام خطي مع ي كما معامل مصفوفة r كما الجانب الأيمن.
الحدس بالنسبة لطريقة:
حدسي ، إذا كانت الرباعي متوازي الاضلاع, سيكون أسهل بكثير من حل المشكلة.طريقة نيوتن هو حل الرباعي مشكلة مع المتعاقبة متوازي الاضلاع تقريبية.في كل تكرار نحن
استخدام المحلية مشتق المعلومات عند النقطة (s,t) لتقريب الرباعي مع متوازي الاضلاع.
العثور على الصحيح (s,t) القيم تحت متوازي الاضلاع التقريب ، من خلال حل نظام الخطي.
الانتقال إلى هذه المرحلة الجديدة وتكرار.
مزايا هذه الطريقة:
كما يتوقع نيوتن-نوع أساليب التقارب هو سريع للغاية.كما التكرارات المضي قدما ، ليس فقط طريقة الحصول على أوثق وأقرب إلى نقطة الحقيقية ، ولكن المحلية متوازي الاضلاع تقريبية تصبح أكثر دقة أيضا ، وبالتالي فإن معدل التقارب في حد ذاته يزيد (في أساليب متكررة المصطلحات ، ونحن نقول أن نيوتن الطريقة quadratically متقاربة).في الواقع هذا يعني أن عدد تصحيح أرقام حوالي يتضاعف كل التكرار.
وهنا ممثل جدول التكرار مقابلخطأ عشوائي المحاكمة فعلت (انظر أدناه للحصول على التعليمات البرمجية):
Iteration Error
1 0.0610
2 9.8914e-04
3 2.6872e-07
4 1.9810e-14
5 5.5511e-17 (machine epsilon)
بعد اثنين من تكرار الخطأ صغيرة بما يكفي لتكون فعالة unnoticable جيدة بما فيه الكفاية بالنسبة لمعظم الأغراض العملية, و بعد 5 التكرار النتيجة دقيقة إلى حدود آلة الدقة.
إذا كان يمكنك إصلاح عدد التكرارات (3 التكرارات بالنسبة لمعظم التطبيقات العملية و 8 التكرار إذا كنت بحاجة إلى دقة عالية جدا) ، ثم خوارزمية بسيطة جدا ومباشرة المنطق ، مع الهيكل الذي يفسح المجال بشكل جيد عالية الأداء حساب.ليست هناك حاجة لفحص كل أنواع خاصة الحالات حافة* باستخدام منطق مختلف اعتمادا على النتائج.بل هو مجرد حلقة تحتوي على عدد قليل من صيغ بسيطة.فيما يلي تسليط الضوء على عدة مزايا هذا النهج على الصيغة التقليدية القائمة على نهج وجدت في إجابات أخرى هنا وفي جميع أنحاء شبكة الإنترنت:
من السهل رمز.جعل مجرد حلقة واكتب في بعض الصيغ.
لا الشرطية أو المتفرعة (إذا/إذن) ، والتي عادة ما يسمح أفضل بكثير pipelining الكفاءة.
لا التربيعي 1 شعبة مطلوب في التكرار (إذا كانت مكتوبة بشكل جيد).جميع العمليات الأخرى هي إضافات بسيطة, الطرح, الضرب.مربع جذور الانقسامات عادة عدة مرات أبطأ من الإضافات/الضرب/الضرب ، و يمكن أن المسمار ذاكرة التخزين المؤقت الكفاءة على بعض أبنية (وأبرزها على بعض الأنظمة المضمنة).في الواقع, إذا كنت تبحث تحت غطاء محرك السيارة في كيفية الجذور التربيعية و الشعب هي في الواقع تحسب لغات البرمجة الحديثة, كلاهما استخدام المتغيرات من طريقة نيوتن ، وأحيانا في الأجهزة في بعض الأحيان في البرامج اعتمادا على العمارة.
يمكن بسهولة vectorized للعمل على المصفوفات مع عدد كبير من رباعيات في آن واحد.أرى vectorized البرمجية أدناه للحصول على مثال عن كيفية القيام بذلك.
يمتد إلى التعسفي الأبعاد.الخوارزمية يمتد في طريقة مباشرة معكوس multilinear الاستيفاء في أي عدد من الأبعاد (2d, 3d, 4d, ...).لقد تضمنت نسخة 3D أدناه ، يمكن للمرء أن يتصور كتابة بسيطة متكررة الإصدار ، أو باستخدام التلقائي التمايز المكتبات للذهاب إلى n-الأبعاد.طريقة نيوتن عادة المعارض البعد مستقلة تقارب أسعار الفائدة ، وذلك في مبدأ الأسلوب يجب أن تكون قابلة للتطوير إلى بضعة آلاف الأبعاد (!) على الأجهزة الحالية (بعد مرحلة بناء وحل ن-ب-ن مصفوفة ي من المحتمل أن يكون العامل المحدد).
بالطبع معظم هذه الأمور يعتمد أيضا على الأجهزة, مترجم, ومجموعة من العوامل الأخرى ، لذلك الأميال الخاص بك قد تختلف.
كود:
على أي حال, هنا هو بلدي Matlab code:(أنا الإفراج عن كل شيء هنا إلى المجال العام)
الأساسية 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
سريع vectorized 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
سرعة هذا القانون ضمنا تستخدم الصيغة التالية معكوس من مصفوفة 2 × 2:
[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
على المرء أن يكون حذرا حول ترتيب المدخلات نقطة منذ معكوس multilinear الاستيفاء فقط واضحة المعالم إذا كان شكل إيجابية حجم و في 3D هو أسهل بكثير من اختيار النقاط التي تجعل شكل بدوره داخل-خارج نفسه.
وهنا تنفيذ نظري حل Naaff، كما ويكي المجتمع. شكرا مرة أخرى.
وهذا هو تطبيق 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 المعادلات، 1 ل x و 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 غير بين الحد الأدنى والحد الأقصى للقيم في الأركان الأربعة للمربع، ثم لا، ليس من الممكن بشكل عام لإيجاد حل واحد (ق، ر ) بحيث interpolant المترابط سوف اعطيكم تلك القيمة.
في عام، وسوف يكون هناك عدد لانهائي من الحلول (ق، ر) داخل الساحة. وسوف تقع على طول (القطعي) مسار منحني من خلال مربع.
إذا وظيفة الخاص بك هو متجه قيمتها واحد، بحيث يكون لديك اثنين من القيم المعروفة في بعض نقطة غير معروفة في الساحة؟ نظرا القيم المعروفة للمعلمتين في كل زاوية من الساحة، ثم حل قد تكون موجودة، ولكن لا يوجد ضمان لذلك. نتذكر أننا يمكن أن نرى في ذلك مشكلتين منفصلتين ومستقلة. فإن الحل لكل واحد منهم يكذب على طول خط كفاف القطعي من خلال مربع. إذا كان الزوج من معالم عبور داخل الساحة، ثم وجود حل. إذا لم يكن لعبور، ثم لا يوجد حل.
وتسأل أيضا في حالة وجود صيغة بسيطة لحل المشكلة. آسف، ولكن ليس حقا أن أرى. كما قلت، منحنيات هي القطعي.
وحل واحد هو للتبديل إلى طريقة مختلفة من الاستيفاء. وذلك بدلا من المترابط، وكسر مربع في زوج من مثلثات. داخل كل مثلث، يمكننا الآن استخدام الخطية حقا. وحتى الآن لا يمكننا حل النظام الخطي من 2 المعادلات في 2 المجهولة داخل كل مثلث. قد يكون هناك حل واحد في كل مثلث، باستثناء حالة منحطة النادرة التي يحدث الخطوط الكنتورية دالة متعددة التعريف الخطية المقابلة أن يكون شارك في الحادث.
ولقد يساء تفسيرها بعض الاستجابات قليلا سؤالك. بمعنى آخر. فهي افتراض يتم منحك <م> قيمة من وظيفة غير معروف محرف، بدلا من وضع ص محرف (س، ص) داخل رباعية أنك تريد العثور على (ق، ر) إحداثيات. هذه هي مشكلة أبسط وهناك ضمان أن يكون الحل الذي هو تقاطع مستقيمين من خلال رباعية.
واحد من خطوط ستخفض من خلال p0p1 قطاعات وp2p3، والآخر ستخفض من خلال p0p2 وp1p3، على غرار حالة الانحياز المحور. وتعرف هذه السطور فريد من موقف ص (س، ص)، وسوف تتقاطع الواضح في هذه المرحلة.
وبالنظر فقط الخط الذي يمر عبر p0p1 وp2p3، يمكننا أن نتصور عائلة من هذه الخطوط، ولكل مختلفة ق القيمة نختار، كل قطع رباعية في عرض مختلفة. إذا كنا إصلاح وذات القيمة الرمزية، يمكن أن نجد نقاط النهاية اثنين عن طريق تحديد ر = 0 و t = 1.
وهكذا لأول مرة تحمل خط ديه النموذج: ذ = A0 * س + B0
وبعد ذلك نحن نعرف اثنين من النهاية من هذا الخط، إذا نصلح قيمة معينة الصورة. وهم:
و(1-ق) P0 + (ق) P1
و(1-ق) P2 + (ق) P3
ونظرا لهذه نقطة نهاية اثنين، يمكننا تحديد عائلة الخطوط عن طريق توصيل لهم في المعادلة للخط، وحل لA0 وB0 <م> وظائف من الصورة . وضع على القيمة الرمزية يعطي صيغة للحصول على خط معين. كل ما نحتاج إليه الآن هو معرفة أي سطر في هذه العائلة يضرب وجهة نظرنا ص (س، ص). ببساطة توصيل إحداثيات ص (س، ص) في المعادلة الخطية لدينا، ويمكن بعد ذلك حل لقيمة الهدف من الصورة.
ويمكن أن يتم هذا النهج المقابل لإيجاد ر كذلك.
حسنا، إذا ص هو نقطة 2D، نعم، يمكنك الحصول عليه بسهولة. في هذه الحالة، S هو العنصر كسري من إجمالي عرض الرباعي في T، T هو أيضا عنصر الكسري للارتفاع إجمالي الرباعي في S.
إذا، ومع ذلك، ف هو العددية، انها لا يمكن بالضرورة، لأن وظيفة المترابط الاستيفاء ليست متجانسة بالضرورة.
