# MatrixOP help

sjr51

Mon, 01/05/2015 - 07:52 am

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

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

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

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

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.

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.

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.

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

January 5, 2015 at 11:31 pm - Permalink