657

edits

Jump to navigation
Jump to search

→Distance between points

Line 98:
Line 98:

=== Distance between points ===

=== Distance between points ===

==== Problem ====

==== Problem ====

+Given a set of points in space we want to calculate the distance between all of them. Each point is described by its components <math> (x_i,y_i,\ldots)</math>. Asusme that the points are saved in a matrix '''<tt>P</tt>''' with '''<tt>N</tt>''' rows (one for each point) and '''<tt>D</tt>''' columns, one for each component.

==== Solution ====

==== Solution ====

+One way of proceeding is to use the broadcast properties of operators in GNU Octave. The square distance between the points can be calculated with the code

+<!-- {{SyntaxHighlight| -->

+{{Code|Calculate square distance between points|<syntaxhighlight lang="octave" style="font-size:13px">

+[N, dim] = size (P);

+Dsq = zeros (N);

+for i = 1:dim

+ Dsq += (P(:,i) - P(:,i)').^2;

+endfor

+</syntaxhighlight>}}

+This matrix is symmetric with zero diagonal.

+
+Similarly the vectors pointing from one point to the another is

+<!-- {{SyntaxHighlight| -->

+{{Code|Calculate radius vector between points|<syntaxhighlight lang="octave" style="font-size:13px">

+R = zeros (N,N,dim);

+for i = 1:dim

+ R(:,:,i) = P(:,i) - P(:,i)';

+endfor

+</syntaxhighlight>}}

+
+The relation between <tt>Dsq</tt> and <tt>R</tt> is

+<!-- {{SyntaxHighlight| -->

+{{Code|Calculate square distance from radius vector between points|<syntaxhighlight lang="octave" style="font-size:13px">

+Dsq = sumsq (R,3);

+</syntaxhighlight>}}

+
+
==== Discussion ====

==== Discussion ====

+The calculation can be implemented using functions like <tt>cellfun</tt> and avoid the loop over components of the points. However in most cases we will have more points than components and the improvement, if any, will be minimal.

+
+Another observation is that the matrix Dsq is symmetric and we could store only the lower or upper triangular part. To use this optimization in a practical way check the help of the functions <tt>vech</tt> and <tt>unvech</tt> (this one is in the Forge package ''general''). Two functions that haven't seen the light yet are <tt>sub2ind_tril</tt> and <tt>ind2sub_tril</tt> (currently private functions in the [[Mechanics_package | Forge package mechanics]]) that are useful to index the elements of a vector constructed with the function <tt>vech</tt>. Each page (the third index) of the multidimensional array <tt>R</tt> is an anti-symmetric matrix and we could also save some memory by keeping only one of the triangular submatrices.

== Plotting ==

== Plotting ==

== User input ==

== User input ==

657

edits

Retrieved from "https://wiki.octave.org/Special:MobileDiff/1852"