# Lazarus

## Programming => General => Topic started by: ranny on May 24, 2019, 10:28:02 am

Title: Matrix unit not giving expected results
Post by: ranny on May 24, 2019, 10:28:02 am
Hello,

I am new to programming and trying out a few maths related programs that will help in my work.  One such idea is solving and displaying 3D data.  Working with matrices is best so I found the 'matrix'unit that should have been a help.

However, it seems that when multiplying a matrix and vector dose not work properly, the results are not correct when compared to simple hand calculations to prove what the result should be.  I thought it might be the convention for RxC or CxR in the matrix when defining the values so I tried the transpose function and it did not transpose the matrix.

One point is that my vectors are always part of an array for example:-

vo[1]:=m1*v[1];

This does not give the correct results for a simple translation matrix so I have had to write my own which seems unnecessary.

Anyone else had the same issue?
Title: Re: Matrix unit not giving expected results
Post by: wp on May 24, 2019, 12:27:31 pm
vo[1]:=m1*v[1];
Seeing a "1" here makes me think that you might not have caught the point that the indexes used by the matrix and vector objects begin with "0".

I've never used this unit, but in my tests for this post it seems that it is working correctly, at least as long as matrix*matrix and matrix*vector multiplication are concerned. What is strange is that the routine for the Inverse of a matrix requires the determinant as a parameter - I've never seen this.

As some help to get started see the attached demo:
Code: Pascal  [Select][+][-]
1. program Project1;
2.
3. {\$mode objfpc}{\$H+}
4.
5. uses
6.   Matrix;
7.
8. var
9.   M: TMatrix3_double;
10.   v: TVector3_double;
11.
12.   procedure PrintMatrix(Caption: String; M: TMatrix3_double);
13.   var
14.     i, j: Integer;
15.   begin
16.     WriteLn(Caption, ' =');
17.     for j := 0 to 2 do begin
18.       for i := 0 to 2 do
19.         Write(M.Data[i, j]:10:4);
20.       WriteLn;
21.     end;
22.     WriteLn;
23.   end;
24.
25.   procedure PrintVector(Caption: String; V: TVector3_double);
26.   var
27.     i: Integer;
28.   begin
29.     WriteLn(Caption, ' =');
30.     for i := 0 to 2 do
31.       WriteLn(V.Data[i]:10:4);
32.     WriteLn;
33.   end;
34.
35. begin
36.   M.Init(
37.     11, 12, 13,
38.     21, 22, 23,
39.     31, 32, 33);
40.
41.   v.Init(
42.     2,
43.     3,
44.     1
45.   );
46.
47.   PrintMatrix('M', M);
48.   PrintVector('M*v', M*v);
49.
51. end.
Title: Re: Matrix unit not giving expected results
Post by: devEric69 on May 24, 2019, 01:45:07 pm
You didn't post the results of your calculations (PrintMatrix, M)  :( ?
Don't you find [line 1, col. 1] = 84?
Title: Re: Matrix unit not giving expected results
Post by: ranny on May 24, 2019, 02:39:16 pm
WP - Thanks for the test run.

The [1] is not the data, that is my array index.  For data it would be vector.data[0] for the first entry and in my case being an array the first entry would be vector[1].data[0].  That part seems to work okay.  I made this earlier....

Code: Pascal  [Select][+][-]
1. procedure TForm1.Button1Click(Sender: TObject);
2.     var
3.     m1:TMatrix4_single;
4.     vi,vo:TVector4_single;
5.     i,j:integer;
6.
7.     begin
8.
9.
10.     stringgrid1.Cells[1,0]:='c0';
11.     stringgrid1.Cells[2,0]:='c1';
12.     stringgrid1.Cells[3,0]:='c2';
13.     stringgrid1.Cells[4,0]:='c3';
14.     stringgrid1.Cells[5,0]:='vi';
15.     stringgrid1.Cells[6,0]:='vo';
16.
17.     stringgrid1.Cells[0,1]:='r0';
18.     stringgrid1.Cells[0,2]:='r1';
19.     stringgrid1.Cells[0,3]:='r2';
20.     stringgrid1.Cells[0,4]:='r3';
21.
22.     stringgrid1.Cells[0,6]:='r0';
23.     stringgrid1.Cells[0,7]:='r1';
24.     stringgrid1.Cells[0,8]:='r2';
25.     stringgrid1.Cells[0,9]:='r3';
26.
27.     m1.init_identity;
28.
29.     m1.data[3,0]:=1;
30.     m1.data[3,1]:=1;
31.
32.     vo.init_zero;
33.
34.     vi.data[0]:=2;
35.     vi.data[1]:=3;
36.     vi.data[2]:=4;
37.     vi.data[3]:=1;
38.
39.     vo:=matvec(m1,vi);
40.     for i:=0 to 3 do
41.         begin
42.         for j:=0 to 3 do
43.             begin
44.             stringgrid1.cells[i+1,j+1]:=floattostr(m1.data[i,j]);
45.             end;
46.         stringgrid1.cells[5,i+1]:=floattostr(vi.data[i]);
47.         stringgrid1.cells[6,i+1]:=floattostr(vo.data[i]);
48.         end;
49.
50.     m1.transpose;
51.     vo:=m1*vi;
52.     for i:=0 to 3 do
53.         begin
54.         for j:=0 to 3 do
55.             begin
56.             stringgrid1.cells[i+1,j+6]:=floattostr(m1.data[i,j]);
57.             end;
58.         stringgrid1.cells[5,i+6]:=floattostr(vi.data[i]);
59.         stringgrid1.cells[6,i+6]:=floattostr(vo.data[i]);
60.         end;
61.
62.
63.     end;

The results are shown below with the first table displaying the results of a translation matrix as expected from my simple matrix x vector function.  The table below is using the matrix unit function of matrix x vector and the results are wrong, with the m1 matrix transposed or not.

I'm new to this forum so apologies if I am not passing on the info in the right format.....

Title: Re: Matrix unit not giving expected results
Post by: tetrastes on May 24, 2019, 07:08:04 pm

Code: Pascal  [Select][+][-]
1.      m1.transpose;
2.     vo:=m1*vi;
3.     ....
4.

The table below is using the matrix unit function of matrix x vector and the results are wrong, with the m1 matrix transposed or not.

Try
Code: Pascal  [Select][+][-]
1.      m1 := m1.transpose;
2.     vo:=m1*vi;
3.     ....
4.

;)
Title: Re: Matrix unit not giving expected results
Post by: wp on May 24, 2019, 07:38:38 pm
I assume from your screenshot that you want to multiply the matrix
Code: [Select]
`  1  0  0  1  0  1  0  1  0  0  1  0  0  0  0  1`with the vector [2, 3, 4, 1].

Using SciLab it can be shown that this product is [3, 4, 4, 1]:
Code: [Select]
`M =    1.0000    0.0000    0.0000    1.0000    0.0000    1.0000    0.0000    1.0000    0.0000    0.0000    1.0000    0.0000    0.0000    0.0000    0.0000    1.0000v =    2.0000    3.0000    4.0000    1.0000M*v =    3.0000    4.0000    4.0000    1.0000`
And here matrix.pas:
Code: Pascal  [Select][+][-]
1. program Project1;
2.
3. {\$mode objfpc}{\$H+}
4.
5. uses
6.   Matrix;
7.
8. var
9.   M: TMatrix4_double;
10.   v: TVector4_double;
11.
12.   procedure PrintMatrix(Caption: String; M: TMatrix4_double);
13.   var
14.     i, j: Integer;
15.   begin
16.     WriteLn(Caption, ' =');
17.     for i := 0 to 3 do begin
18.       for j := 0 to 3 do
19.         Write(M.Data[i, j]:10:4);
20.       WriteLn;
21.     end;
22.     WriteLn;
23.   end;
24.
25.   procedure PrintVector(Caption: String; V: TVector4_double);
26.   var
27.     i: Integer;
28.   begin
29.     WriteLn(Caption, ' =');
30.     for i := 0 to 3 do
31.       WriteLn(V.Data[i]:10:4);
32.     WriteLn;
33.   end;
34.
35. begin
36.   M.Init(
37.     1, 0, 0, 1,
38.     0, 1, 0, 1,
39.     0, 0, 1, 0,
40.     0, 0, 0, 1
41.   );
42.
43.   v.Init(
44.     2,
45.     3,
46.     4,
47.     1
48.   );
49.
50.   PrintMatrix('M', M);
51.   PrintVector('v', v);
52.   PrintVector('M*v', M*v);
53.
55. end.

Output:
Code: [Select]
`M =    1.0000    0.0000    0.0000    1.0000    0.0000    1.0000    0.0000    1.0000    0.0000    0.0000    1.0000    0.0000    0.0000    0.0000    0.0000    1.0000v =    2.0000    3.0000    4.0000    1.0000M*v =    3.0000    4.0000    4.0000    1.0000`The same result! This confirms that the calculation of the matrix unit is correct.

So, what is wrong with your example?

I see the discrepancy between the screenshot (which shows the matrix above) and the instructions how you create it:
Code: Pascal  [Select][+][-]
1.     m1.init_identity;
2.
3.     m1.data[3,0]:=1;
4.     m1.data[3,1]:=1;
You start with an identity matrix and put "1" into the elements [3,0] and [3,1]. However, according to math conventions the first index of a matrix is the index of the row and the second index is the index of the column. Therefore you do not add the 1's to the upper right corner of the matrix, but to the lower left corner! And this, effectively, is the transpose of the matrix that you want. According to SciLab, the product of the transposed matrix with the vector is [2, 3, 4, 6]. This is the result that you show in the lower part of the screenshot.

What's wrong with the screenshot? You did the same error that I did in my previous post: confused the order of the index. As already said math identifies the rows with the first index and the columns with the second index, but the stringgrid uses the first index for the column and the second index for the row; therefore, in your code the matrix indexes i, j show appear in reverse order for the grid cells. Since they are in the same order you effectively draw the transposed matrix!

Another issue: You do have a "m1.Transpose" in your code (I don't know why), but this is useless, because Transpose is a function: it does not in-place transpose the matrix data but transfers the transposed data to another matrix - but you do not assign the result of the function to a variable; therefore you discard the result of the operation.
Title: Re: Matrix unit not giving expected results
Post by: ranny on May 24, 2019, 08:09:47 pm
WP....

Thank you very much for your reply.  Its very interesting and enlightening because if you look at the original post I thought it might be a mix up with row and column but when transpose didn't work I thought there was something else afoot.  if m1.transpose doesn't work then I am surprised I didn't get a compile error, if I did I may have stumbled across the answer myself!

My other code in my program is (unwittingly I guess) structured to read in the way that is the convention shown in many maths and wikipedia pages and I guess i got it to work through my own coding the way I expected it to.  Now as you point out  if the structure is of R versus C is opposite to my thinking then I can deal with it.

I could leave my original function for m*v that I have written but was hoping that the matrix unit functions would be faster.  I will muck about with the code and do some tests to see how it works both ways now you have pointed me in a good direction.

Thank you very much for your time on this, very much appreciated.  For your information, this is what I am doing... solving 3D polynomials and showing the result graphically so that one can see how the degree of fit is with respect to the original data.  Its getting into shape slowly....

Thanks again.

UPDATE:-  Actually, now that I look at it, I realise how daft I have been, should have seen it!!
Title: Re: Matrix unit not giving expected results
Post by: wp on May 24, 2019, 10:20:10 pm
Just to make sure: look at the image at the top of https://en.wikipedia.org/wiki/Matrix_(mathematics). It shows the 1st index of the matrix elements running down along the columns - this means: the first index is the index of the ROWs. And the second index runs along the rows, thus it identifies the COLUMNs. A matrix element therefore is addressed as M[r, c] where r and c are the row and column indexes, respectivley.

In a stringgrid, all cell-related functions and properties use the opposite order: Grid.Cells[c, r].

And I also would like to draw your attention to TAChart which contains a self-contained, linear least-squares library (unit TAFitLib, http://wiki.lazarus.freepascal.org/TAChart_documentation#Fit_series); it is 2D only, however. (Run the demo in (lazarus)/components/tachart/demo/fit).
Title: Re: Matrix unit not giving expected results
Post by: devEric69 on May 31, 2019, 09:42:51 pm
wp solved everything well.

@ranny: beautiful graphic.
For your information, the convention that I find obsolete inherited from the early days of Borland, calling everything "ABlob", does not bring much semantics into a program.
You can use the Hungarian notation to compensate for this and no longer search, nor where are the rows and columns, and what type is their indexes...
Code: Pascal  [Select][+][-]
1. procedure TfrmMain.PrintMatrix(sCaption: String; tdM: TMatrix4_double);
2.  var
3.    iLine, jCol: Integer;
4.    sTemp: string;
5.  begin
6.    memForResults.Append(sCaption+ ' =');
7.    for iLine := 0 to 3 do begin
8.      for jCol := 0 to 3 do begin
9.        sTemp:= sTemp + FloatToStr(tdM.Data[iLine, jCol]) + #9;
10.      end;
11.     sTemp:= sTemp + #13#10;
12.    end;
13.    memForResults.Append(sTemp + #13#10);
14.  end;
15.
16.  procedure TfrmMain.PrintVector(sCaption: String; tdV: TVector4_double);
17.  var
18.    i: Integer;
19.    sTemp: string;
20.  begin
21.    memForResults.Append(sCaption + ' =');
22.    for i := 0 to 3 do
23.      sTemp:= sTemp + FloatToStr(tdV.Data[i]) +#13#10;
24.    memForResults.Append(sTemp + #13#10);
25.  end;
26.
27. procedure TfrmMain.btnCalculateClick(Sender: TObject);
28. var
29.   tdM: TMatrix4_double;
30.   tdV: TVector4_double;
31. begin
32.   tdM.Init(
33.           1, 0, 0, 1,
34.           0, 1, 0, 1,
35.           0, 0, 1, 0,
36.           0, 0, 0, 1
37.         );
38.
39.   tdV.Init(
40.           2,
41.           3,
42.           4,
43.           1
44.         );
45.
46.     PrintMatrix('M', tdM);
47.     PrintVector('v', tdV);
48.     PrintVector('M*v', tdM*tdV);
49. end;
50.

Title: Re: Matrix unit not giving expected results
Post by: ranny on June 10, 2019, 08:32:41 am
Thanks DevEric69, that's worth looking into.

I added the ability to shade and apply colour gradient now, its working well now but I will look into the matrix stuff in a more detailed way.
Title: Re: Matrix unit not giving expected results
Post by: devEric69 on June 10, 2019, 01:05:07 pm
Nice rendering of your polynomial interpolation (I imagine these are spline in ^3, or ^⁴)  8) .

Your 3D rendering graphic engine is also written in Pascal (if so, would you tell the solution used if it's not indiscreet)?
Title: Re: Matrix unit not giving expected results
Post by: ranny on June 10, 2019, 01:31:53 pm
@devEric69

Thanks.  It is a straight polynomial in x and z so y:=x^n.z^n+x^n.z^n-1 etc. and the order of fit can be changed from 1 to a maximum of 7 (just to stop it going a bit crazy) so you get a single equation solution.  Must have enough data points to solve so npoints >(degree of fit +1)^2.  The purpose of the program is to make solutions to 3D maps but show the result so that you can see if the equation is good for your purpose (another program perhaps).  Often values outside the original data bounds will be dodgy so with this you can see whats happening.

Regarding the display of the image, the best thing to do is look for OneLoneCoder videos on YouTube.  This guy shows code development for doing 3D graphics (in C++ but it can be re-written).  I took his basic code and removed the viewer viewing position and also disabled the selection of which faces are drawn since you want to see under the image.  Its very good stuff.  The gradient to go from blue to red based on value is just a colour map that is picked from depending on the y value.

Have a look at the video https://www.youtube.com/watch?v=ih20l3pJoeU and its later videos.

My code is a mess but if after you have looked at the videos you want some more info then we can think about how to do that in a way it does not look a mess!
Title: Re: Matrix unit not giving expected results
Post by: devEric69 on June 12, 2019, 06:18:10 pm
Thank you for this link: I looked. This person is really very, very didactic (3D is a theme in which I have been interested in dotted lines, but I have never found anything so complete and applied. All application and transformation matrices are very well explained).
But... because there is a but: this person calls a point vector (x, y, z), which is in fact a single-line matrix and that he uses as such, i. e. a transposed vector :) .
==> hence a quiporoquo, from the beginning of his videos. But it's a small detail considering all his explanations. And given his mastery of the subject, he must have a reason to leave in this way - perhaps, a facility with the grammar or the functioning of the operators, his development tool?  - that he forgot to mention.
Title: Re: Matrix unit not giving expected results
Post by: ranny on June 13, 2019, 09:00:12 am
I thought his videos were excellent mainly because of his explanation of what he was doing.  The (x,y,z) is a vertex and a vector depending on the context of what is being applied but as a single row (or column, read on!) it allows matrix multiplications. Watch the videos and get an understanding of the maths and have a look at fundamental 3D .

However, this whole thread I started was because of the matrix unit in Lazarus not giving expected results.  The reasons for this are clearly down to his videos showing that his matrices are defined as column then row whereas wikipedia and the like show row then column.  You can see the difference when he explains the transformation matrix then compare it to others online.  So if you use his code throughout it works, but when I tried to use the matrix.vector routines from the matrix unit it was all over the place.

If you are interested in generated 3D curves and simple images then its worth trying to write your own based on his tutorials, and again looking online shows different methods but consistency is key.  I will see if I can smarten mine up a bit and pass it on.