MatrixOP help

Share icon

Hi, I would like to use a MatrixOP to speed up my code but I'm not sure how to do it. I only partially understand how matrix operations work, even after reading the help.

My problem: I have a bunch of 3D objects that I am rotating about the z-axis and then by the y-axis (to do this I use MatrixMultiply and this is fine). What I'm doing is performing this for all possible angles about z and y (and computing distances to find the lowest total value - this is not the problem). Right now my code makes two 1D waves that correspond to all possible angles, I then step through the combinations using For-EndFor loops (see code). I am sure that I must be able to make a matrix/2D wave of these values and use a MatrixOp to do this and I know that this must be faster, but I can't see how to do it.

Any help would be very much appreciated.

Function LowPoint()
   
    variable theta,phi //in radians
    Make /o /n=180 MatThetaWave
    MatThetaWave =x/(180/PI)
    Make /o /n=360 MatPhiWave
    MatPhiWave =x/(180/PI)
    Make /o /n=(360,180) MatResWave

    String wList=wavelist("MT*",";","")
    String wName
    Variable i,j,k
    Variable px //pixel value - how much black would be seen in gizmo en face

    For (j = 0; j < numpnts(MatThetaWave); j +=1)
        Theta=MatThetaWave(j)
        For (k = 0; k < numpnts(MatPhiWave);  k +=1)
            px=0
            Phi=MatPhiWave(k)
            Make/o zRotationMatrix={{cos(phi),-sin(phi),0},{sin(phi),cos(phi),0},{0,0,1}}
            Make/o yRotationMatrix={{cos(theta),0,sin(theta)},{0,1,0},{-sin(theta),0,cos(theta)}}
                For (i = 0; i < ItemsInList(wList); i += 1)
                    wName = StringFromList(i, wList)
                    Wave/Z w = $wName
                    MatrixMultiply w,zRotationMatrix
                    Wave M_Product
                    MatrixMultiply M_Product,yRotationMatrix
                    px +=sqrt(((M_Product[1][0]-M_Product[0][0])^2)+((M_Product[1][1]-M_Product[0][1])^2))
                EndFor
            MatResWave[k][j]=px
        EndFor
    EndFor
    WaveStats MatResWave
    Print "Rotating z by phi =",V_minRowLoc/(180/PI)," and then rotating y by theta = ",V_minColLoc/(180/PI)
    Straight(V_minRowLoc/(180/PI),V_minColLoc/(180/PI))
End
I'll start by suggesting that this is probably not the best way to find a minimum (if that is your goal).

Looking at your code here are a few quick comments:

1. it would be useful to know the dimensions of your MT* waves. This is important if you want to use MatrixOP in a meaningful way. The general idea would be to try and pack your MT* waves into a single matrix that would be multiplied by the various rotations.
2. two rotations can be concatenated into a single matrix multiplication.
3. it is not clear why you need both angle waves as they are directly generated from the index.

A.G.
WaveMetrics, Inc.

January 5, 2015 at 11:03 am - Permalink

Thank you for the advice. Yes, I'm sure there are better ways to find the minimum but this seemed like a good place to start!

In answer to your comments:
1. The MT* waves are three columns corresponding to X,Y,Z coordinates. They have a start and end i.e. two rows. The idea is to find a rotation condition that will maximally align all waves to the zenith. This will then allow us to compare MT* waves from many experiments.
2. I didn't realise that. Do you mean have the second matrix (yRotationMatrix) as second layer to the first matrix?
3. Ah yes I see. But would you still advocate using this looped index approach rather than a matrix to set these values?

Thank you again for the help. I'd appreciate further input.

January 5, 2015 at 12:45 pm - Permalink

sjr51 wrote:

1. The MT* waves are three columns corresponding to X,Y,Z coordinates. They have a start and end i.e. two rows. The idea is to find a rotation condition that will maximally align all waves to the zenith. This will then allow us to compare MT* waves from many experiments.


I am not sure I understand the details but my hunch is that you might do better by running PCA on the data and then performing a single rotation.

Quote:

2. I didn't realise that.

This is a general rule. You can find the details e.g., http://en.wikipedia.org/wiki/Euler%27s_rotation_theorem. See the third paragraph.

Quote:
Do you mean have the second matrix (yRotationMatrix) as second layer to the first matrix?

No. Just write out the matrix representation of rotating any vector and then applying another rotation to the result. You should then see what I mean.
Quote:

3. Ah yes I see. But would you still advocate using this looped index approach rather than a matrix to set these values?

I'd be tempted to store each MT wave as a layer of a 3D wave say MT3D. Next, for each proposed rotation you find the equivalent matrix RM and then the rotation of all MT's is simply:
MatrixOP/O aa=MT3D x RM
Finally, a tricky way to evaluate the square root expression is to appropriately transpose (see ImageTransform transposeVol) the resulting aa so that you are subtracting layers.

I hope this helps,

A.G.
WaveMetrics, Inc.

January 5, 2015 at 05:36 pm - Permalink

Thank you for your input A.G. I'll try 2, then 3. I'll explore 1 as an alternative solution.

January 5, 2015 at 11:31 pm - Permalink