समेकित समीकरणों के लिए कम स्क्वायर समाधान

Question

मैं निर्देशांक के एक सेट से दूसरे में एक रूपांतरण को फिट करने की कोशिश कर रहा हूं।

x' = R + Px + Qy 
y' = S - Qx + Py 
Where P,Q,R,S are constants, P = scale*cos(rotation). Q=scale*sin(rotation) 

पी, क्यू, आर, एस को संबंधित बिंदुओं के एक सेट में फ़िट करने के लिए एक प्रसिद्ध 'हाथ से' सूत्र है। लेकिन मुझे फिट पर एक त्रुटि अनुमान होना चाहिए - इसलिए मुझे कम से कम वर्ग समाधान की आवश्यकता है।

'संख्यात्मक व्यंजनों' पढ़ें लेकिन मुझे यह पता लगाने में परेशानी हो रही है कि एक्स और वाई के साथ डेटा सेट के लिए इसे कैसे किया जाए।

क्या कोई यह कैसे कर सकता है इसके उदाहरण/ट्यूटोरियल/कोड नमूना को इंगित कर सकता है?
भाषा के बारे में भी परेशान नहीं है।
लेकिन - बस मैटलैब/लैपैक/numpy/R की सुविधा में निर्मित उपयोग शायद उपयोगी नहीं है!

संपादित करें: मेरे पास फिट करने के लिए पुराना (एक्स, वाई) नया (एक्स, वाई) का एक बड़ा सेट है। समस्या अति निर्धारित है (अज्ञात से अधिक डेटा अंक) इतना आसान मैट्रिक्स उलटा पर्याप्त नहीं है - और जैसा कि मैंने कहा है कि मुझे वास्तव में फिट पर त्रुटि की आवश्यकता है।

Answer 1

निम्नलिखित कोड को चाल करना चाहिए। मैं बच के लिए निम्न सूत्र का प्रयोग किया:

residual[i] = (computed_x[i] - actual_x[i])^2 
       + (computed_y[i] - actual_y[i])^2 

और फिर व्युत्पन्न कम से कम वर्गों general procedure Wolfram के मैथवर्ल्ड पर दिए गए वर्णन के आधार पर सूत्रों।

मैंने एक्सेल में इस एल्गोरिदम का परीक्षण किया और यह अपेक्षित प्रदर्शन करता है। मैंने दस यादृच्छिक बिंदुओं का संग्रह किया जो बाद में घूर्णन, अनुवादित और यादृच्छिक रूप से जेनरेट किए गए रूपांतरण मैट्रिक्स द्वारा स्केल किए गए थे।

कोई यादृच्छिक शोर उत्पादन डेटा पर लागू के साथ

, इस कार्यक्रम चार मानकों (P, Q, R, और S) जो इनपुट पैरामीटर के समान हैं, और शून्य का एक rSquared मूल्य पैदा करता है।

आउटपुट बिंदुओं पर अधिक से अधिक यादृच्छिक शोर लागू होता है, इसलिए स्थिरांक सही मानों से दूर जाने लगते हैं, और rSquared मूल्य तदनुसार बढ़ता है। ,

// test data 
const int N = 1000; 
float oldPoints_x[N] = { ... }; 
float oldPoints_y[N] = { ... }; 
float newPoints_x[N] = { ... }; 
float newPoints_y[N] = { ... }; 

// compute various sums and sums of products 
// across the entire set of test data 
float Ex = Sum(oldPoints_x, N); 
float Ey = Sum(oldPoints_y, N); 
float Exn = Sum(newPoints_x, N); 
float Eyn = Sum(newPoints_y, N); 
float Ex2 = SumProduct(oldPoints_x, oldPoints_x, N); 
float Ey2 = SumProduct(oldPoints_y, oldPoints_y, N); 
float Exxn = SumProduct(oldPoints_x, newPoints_x, N); 
float Exyn = SumProduct(oldPoints_x, newPoints_y, N); 
float Eyxn = SumProduct(oldPoints_y, newPoints_x, N); 
float Eyyn = SumProduct(oldPoints_y, newPoints_y, N); 

// compute the transformation constants 
// using least-squares regression 
float divisor = Ex*Ex + Ey*Ey - N*(Ex2 + Ey2); 
float P = (Exn*Ex + Eyn*Ey - N*(Exxn + Eyyn))/divisor; 
float Q = (Exn*Ey + Eyn*Ex + N*(Exyn - Eyxn))/divisor; 
float R = (Exn - P*Ex - Q*Ey)/N; 
float S = (Eyn - P*Ey + Q*Ex)/N; 

// compute the rSquared error value 
// low values represent a good fit 
float rSquared = 0; 
float x; 
float y; 
for (int i = 0; i < N; i++) 
{ 
    x = R + P*oldPoints_x[i] + Q*oldPoints_y[i]; 
    y = S - Q*oldPoints_x[i] + P*oldPoints_y[i]; 
    rSquared += (x - newPoints_x[i])^2; 
    rSquared += (y - newPoints_y[i])^2; 
} 

Answer 2

एक मुद्दा यह है कि इस तरह की संख्यात्मक सामग्री अक्सर मुश्किल होती है। यहां तक ​​कि जब एल्गोरिदम सरल होते हैं, तब भी अक्सर ऐसी गणना होती है जो वास्तविक गणना में दिखाई देती हैं।

इसी कारण से, यदि कोई सिस्टम है तो आप आसानी से प्राप्त कर सकते हैं जिसमें एक अंतर्निहित सुविधा है, इसका उपयोग करना सबसे अच्छा हो सकता है।

Answer 3

पी, क्यू, आर, और एस खोजने के लिए, तो आप कम से कम वर्गों का उपयोग कर सकते हैं। मुझे लगता है कि उलझन में बात यह है कि कम से कम वर्गों का सामान्य विवरण एक्स और वाई का उपयोग करता है, लेकिन वे आपकी समस्या में एक्स और वाई से मेल नहीं खाते हैं। आपको बस अपनी समस्या को कम से कम वर्ग ढांचे में सावधानी से अनुवाद करने की आवश्यकता है। आपके मामले में स्वतंत्र चर अविभाज्य समन्वय x और y हैं, निर्भर चर परिवर्तनीय समन्वय x 'और y' हैं, और समायोज्य पैरामीटर पी, क्यू, आर, और एस हैं (यदि यह पर्याप्त स्पष्ट नहीं है, मुझे बताएं और मैं अधिक जानकारी पोस्ट करूंगा।)

एक बार जब आपको पी, क्यू, आर, और एस मिल जाए, तो स्केल = sqrt (पी^2 + क्यू^2) और फिर आप रोटेशन पा सकते हैं पाप रोटेशन = क्यू/स्केल और कॉस रोटेशन = पी/स्केल से।

Answer 4

3x3 मैट्रिक्स टी (पी, क्यू, आर, एस) को परिभाषित करें जैसे (x',y',1) = T (x,y,1)। फिर

A = \sum_i |(T (x_i,y_i,1)) - (x'_i,y'_i,1)|^2 

और ए (पी, क्यू, आर, एस) के खिलाफ ए को कम करें।

यह कोडिंग स्वयं बड़े आकार के प्रोजेक्ट का माध्यम है जब तक कि आप यह गारंटी नहीं दे सकते कि डेटा अच्छी तरह से वातानुकूलित है, खासकर जब आप प्रक्रिया से अच्छा त्रुटि अनुमान चाहते हैं। आप शायद एक मौजूदा minimizer कि त्रुटि अनुमानों का समर्थन करता का उपयोग कर बंद सबसे अच्छा कर रहे हैं ..

कण भौतिकी प्रकार minuit का प्रयोग करेंगे तो सीधे CERNLIB से (कोडिंग सबसे आसानी से फ़ोरट्रॉन 77 में किया साथ), या ROOT से (ग में कोडिंग के साथ ++ , या यह पाइथन बाइंडिंग के बावजूद सुलभ होना चाहिए)। लेकिन अगर आपके पास पहले से ही इनमें से कोई भी उपकरण नहीं है तो यह एक बड़ी स्थापना है।

मुझे यकीन है कि अन्य अन्य मिनीमाइज़र सुझा सकते हैं।

Answer 5

आप इसकी गणना करने के लिए levmar प्रोग्राम का उपयोग कर सकते हैं। यह मेरा सहित कई उत्पादों में परीक्षण और एकीकृत है।इसका जीपीएल के तहत लाइसेंस प्राप्त है, लेकिन यदि यह एक गैर-ओपनसोर्स प्रोजेक्ट है, तो वह आपके लिए लाइसेंस (शुल्क के लिए)

Answer 6

धन्यवाद eJames लगभग exaclty मैं क्या है thats:

यहाँ कोड है। मैंने इसे एक पुराने सेना सर्वेक्षण मैनुअल से कोडित किया जो पहले "सर्वेक्षकों के लिए निर्देश" नोट पर आधारित था जो 100 साल पुराना होना चाहिए! (यह एक्स/वाई के बजाए उत्तर और पूर्व के लिए एन और ई का उपयोग करता है)

फिट पैरामीटर की भलाई बहुत उपयोगी होगी - अगर वे फिट खराब कर देते हैं तो मैं चुने गए बिंदुओं को अंतःस्थापित कर सकता हूं।

FindTransformation(vector<Point2D> known,vector<Point2D> unknown) { 
{ 
    // sums 
    for (unsigned int ii=0;ii<known.size();ii++) { 
     sum_e += unknown[ii].x; 
     sum_n += unknown[ii].y; 
     sum_E += known[ii].x; 
     sum_N += known[ii].y;        
     ++n;   
    } 

    // mean position 
    me = sum_e/(double)n; 
    mn = sum_n/(double)n; 
    mE = sum_E/(double)n; 
    mN = sum_N/(double)n; 

    // differences 
    for (unsigned int ii=0;ii<known.size();ii++) { 

     de = unknown[ii].x - me; 
     dn = unknown[ii].y - mn; 

     // for P 
     sum_deE += (de*known[ii].x); 
     sum_dnN += (dn*known[ii].y); 
     sum_dee += (de*unknown[ii].x); 
     sum_dnn += (dn*unknown[ii].y); 

     // for Q 
     sum_dnE += (dn*known[ii].x); 
     sum_deN += (de*known[ii].y);      
    } 

double P = (sum_deE + sum_dnN)/(sum_dee + sum_dnn); 
double Q = (sum_dnE - sum_deN)/(sum_dee + sum_dnn); 

double R = mE - (P*me) - (Q*mn); 
double S = mN + (Q*me) - (P*mn); 
}