<html><head></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; ">Hi,<div><br></div><div>from here :</div><div><br></div><div><a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.cholesky.html">http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.cholesky.html</a></div><div><br></div><div>seems it is applicable only on a square matrix.</div><div>is that true ?</div><div><br></div><div>reading this pdf :</div><div><br></div><div><a href="http://geomatica.como.polimi.it/corsi/misure_geodetiche/MG_03SistRif2.pdf">http://geomatica.como.polimi.it/corsi/misure_geodetiche/MG_03SistRif2.pdf</a></div><div><br></div><div>scroll to the page : 15</div><div><br></div><div>i can see it isn't a common system : &nbsp;Y = AX &nbsp;</div><div>but it is : &nbsp;Y = AX + B + V</div><div><br></div><div><br></div><div>the points are :</div><div><br></div><div><br></div><div><br></div><div><div>p1:</div><div><div>4402553.569 727053.737 4542823.332&nbsp;</div><div>4399375.518 703845.639 4549214.880&nbsp;</div><div>4412911.336 701094.214 4536517.955</div><div><br></div><div><br></div></div></div><div>p2 :</div><div><div>4402553.334 727053.937 4542823.474&nbsp;</div><div>4399375.347 703845.876 4549215.105&nbsp;</div><div>4412911.150 701094.435 4536518.139</div><div><br></div><div><br></div></div><div><br></div><div><br></div><div>tring to rewrite the code i did :</div><div><br></div><div>L = np.loadtxt(p1)</div><div><br></div><div><div>A = np.zeros((3*L.shape[0],7),float)</div><div>A[ ::3, 0] = 1.0</div><div>A[1::3, 1] = 1.0</div><div>A[2::3, 2] = 1.0</div><div>A[1::3, 3] = L[:,2]</div><div>A[2::3, 3] = -L[:,1]</div><div>A[ ::3, 4] = -L[:,2]</div><div>A[2::3, 4] = L[:,0]</div><div>A[ ::3, 5] = L[:,1]</div><div>A[1::3, 5] = -L[:,0]</div><div>A[ ::3, 6] = L[:,0]</div><div>A[1::3, 6] = L[:,1]</div><div>A[2::3, 6] = L[:,2]</div><div><br></div><div>G = np.loadtxt(p2)</div><div><br></div><div>V = np.zeros((3*G.shape[0],1),float)</div><div>V[ ::3, 0] = G[:,0] - L[:,0]</div><div>V[1::3, 0] = G[:,1] - L[:,1]</div><div>V[2::3, 0] = G[:,2] - L[:,2]</div><div><br></div><div><br></div><div>i used&nbsp;the code in the previouse mail :</div><div><br></div><div><div>N = np.dot(A.T.conj(), A)</div><div>T = np.dot(A.T.conj(), Y)</div><div>C = np.dot(linalg.inv(N), T)</div></div><div><br></div><div>to solve a system like :&nbsp;</div><div><br></div><div>&nbsp;Y = AX &nbsp;</div><div><br></div><div>how can i change the code to solve :&nbsp;</div><div><br></div><div>&nbsp;Y = AX + B + V</div><div><br></div><div>instead ?</div><div><br></div><div>the results in the dpf are :</div><div><br></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 10.5px/normal 'Times New Roman'; "><span style="font: 16.0px 'Times New Roman'"><i>x</i></span><i>0 &nbsp; &nbsp;&nbsp;<span class="Apple-style-span" style="font-size: 16px; font-style: normal; ">-9.256 m&nbsp;</span></i></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 10.5px/normal 'Times New Roman'; "><i></i><span style="font: 16.0px 'Times New Roman'"><i>y</i></span><i>0 &nbsp; &nbsp;&nbsp;<span class="Apple-style-span" style="font-size: 16px; font-style: normal; ">-23.701 m&nbsp;</span></i></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 10.5px/normal 'Times New Roman'; "><i></i><span style="font: 16.0px 'Times New Roman'"><i>z</i></span><i>0 &nbsp; &nbsp;&nbsp;<span class="Apple-style-span" style="font-size: 16px; font-style: normal; ">16.792 m</span></i></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 10.5px/normal 'Times New Roman'; "><i></i><span style="font: 16.0px 'Times New Roman'"><i>R</i></span><i>x &nbsp; &nbsp;&nbsp;<span class="Apple-style-span" style="font-size: 16px; font-style: normal; ">-0.0001990982 rad&nbsp;</span></i></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 10.5px/normal 'Times New Roman'; "><i></i><span style="font: 16.0px 'Times New Roman'"><i>R</i></span><i>y &nbsp; &nbsp;&nbsp;<span class="Apple-style-span" style="font-size: 16px; font-style: normal; ">0.0001778762 rad&nbsp;</span></i></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 10.5px/normal 'Times New Roman'; "><i></i><span style="font: 16.0px 'Times New Roman'"><i>R</i></span><i>z &nbsp; &nbsp;&nbsp;<span class="Apple-style-span" style="font-size: 16px; font-style: normal; ">0.00015&nbsp;rad&nbsp;</span></i></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 10.5px/normal 'Times New Roman'; "><i></i><span style="font: 16.0px Helvetica">μ &nbsp; &nbsp;<span class="Apple-style-span" style="font-family: 'Times New Roman'; ">0.00000046</span>&nbsp;</span></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal 16px/normal 'Times New Roman'; "><br></div><div><font class="Apple-style-span" face="'Times New Roman'" size="4"><span class="Apple-style-span" style="font-size: 16px;"><span class="Apple-style-span" style="font-family: Helvetica; font-size: medium; ">any suggestion, also on how to apply a decomposition can give me a great help!</span></span></font></div><div><br></div><div><div><br></div></div></div><div><br><div><div>Il giorno 02/mar/2010, alle ore 19.00, <a href="mailto:scipy-user-request@scipy.org">scipy-user-request@scipy.org</a> ha scritto:</div><br class="Apple-interchange-newline"><blockquote type="cite">Hi,<br><br>I see your coordinates are very big numbers. These can lead to an<br>ill-conditioned (AtA) matrix which can lead to a wrong inverse. Try<br>using cholesky decomposition instead of using linalg.inv. The method<br>is outlined here<span class="Apple-converted-space">&nbsp;</span><a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">http://en.wikipedia.org/wiki/Cholesky_decomposition</a><br>and numpy has routines for doing just this. The other option is to<br>subtract the centre-of-mass from your coordinates i.e. P1 -&gt; P1x -<br>P1xave, P1y - P1yave etc. Do the same for P2. Calculate the<br>transformation again. The shifts should now be zero but your rotations<br>and scale should be statistically independent and correct. The<br>difference between the two coordinate systems' centres of mass gives<br>the shifts. Transformations far away from the origin rarely behaves<br>well.<br><br>Good luck.<br></blockquote></div><br></div></body></html>