0

In trying to use the midpoint method to create geodesics on a surface, but I am not comfortable with the Matlab syntax needed to iteratively update a matrix containing points on a surface embedded in R3, and, therefore, composed of 3 rows and n columns with the first element given as, say, A[:,1], and the last element, A[:,n].

Initially two points would be chosen on the XY plane along a line. Multiple points along the line segment on the XY plane running through A[:,1] and A[:,n] would be projected onto the surface using a function f(x,y), so that we would end up with the initial matrix with the first row containing a number of coordinates in the X axis; the second row, the coordinates in the Y axis; and the third row, the height of the curve at each point of f(x,y). Each column is one point.

So if the first two columns are A[:,1] and A[:,2], I would like to obtain an updated matrix where the first column would be unchanged, A[:,1], and the second column would be the entrywise average of A[:,1] and A[:,2]. In other words, the first entry of the column (x) will be the average of the first entry of the first and second columns; same for the second (y) and third (z) entries. At the other end the matrix would end with the average of A[:,n - 1] and A[:,n] in the penultimate column, and A[:,n] in the very last column.

Then a minimization step would take place to project each column to the surface f(x,y), and the same process of averaging would start again.

Note that with every step the matrix would grow by one column.

Logically, this process would be within a loop in a function.

I would like to ask how to iteratively achieve the averaging step, which, for example would in the first step go from

A =

   0   2   1   4
   1   3   3   2
   1   2   2   2

to

A =
     0   1     1.5   2.5    4    
     1   2     3     2.5    2    
     1   1.5   2     2      2

After the accepted answer, this is the code:

A = [0, 2, 1, 4; 1, 3, 3, 2; 1, 2, 2, 2]
A = [A(:,1), (A(:,1:end-1) + A(:,2:end))/2, A(:,end)]
Antoni Parellada
  • 4,253
  • 6
  • 49
  • 114

1 Answers1

2

The following line updates A based on your description. However, since A's number of column increases (as you already know), you'd better preallocate it and slightly adjust the following assignment (which would go in the for loop).

A = [A(:,1), (A(:,1:end-1) + A(:,2:end))/2, A(:,end)];

However, you might want to check some function from the Image Processing Toolbox, such as imfilter.

Enlico
  • 23,259
  • 6
  • 48
  • 102