3D matrix rearrangement
L.K.Chris
Fri, 06/21/2013 - 02:03 am
I'm having some difficulties writing a procedure to automatize analysis of fluorescence spectra. I initially start out with a 3D dataset, which I then rearrange to 3D>2D>1D>3D, so that I end up with a matrix where my frames are aligned. I now have the problem that every layer represents an episode of 10 frames but with alternating spectra: 1st layer=CFP fluorescence, 2nd layer = YFP fluorescence, 3rd layer= CFP fluorescence and so on. The question is: is there a clever way to shuffle the layers to be able to average all the YFP/CFP spectra? F.ex if i have a matrix of 512(fluorescence intensity)x10(frames)x30(episodes) how do I end up with a 512x10x2(CFP and YFP average) matrix?
The initial wave is coming from the the Speloader function which I found here (// Originally from Sjors Wurpel 2005
// This procedure file contains the LoadPrincetonSPE() function to import binary winspec (.spe) files). The code I have so far is:
Function GetRow(wName,colToTake)
//Takes the Row containing the fluorescence intensity and creates a 2D wave.
wave wName
variable colToTake
variable Nrows
variable Ncols
variable Nlayers
Nrows = DimSize(wName, 0)
Ncols = DimSize(wName, 1)
Nlayers = DimSize(wName, 2)
string NewWname = nameofwave(wName)+"_R"+num2str(coltotake)
make/O/N =(Nrows,Nlayers) wave2d
wave2d= wName[p][colToTake][q]
rename wave2d $NewWname
end
Function MakeMatrix(w2d, FramePerEpisode)
//The function will take a 2D wave containing all the frames as columns and redimension it into a 3d matrix
//with each frame of each episode aligned in the z-plane.
wave w2d
variable FramePerEpisode
variable Ncols
Variable NEps
variable Nrows
variable sizeOf1D
duplicate/O w2d, tempWave
Nrows = DimSize(w2d, 0)
Ncols = DimSize(w2d, 1)
NEps = Ncols / FramePerEpisode
sizeof1D = Nrows*Ncols
redimension/N=(sizeof1d) tempwave
redimension/N=(512, frameperepisode, Neps) tempwave
string newWave = nameofwave(w2d)+"Matrix"
rename tempwave $newWave
setdimlabel 1,0, Frame1, $newWave
setdimlabel 1,1, Frame2, $newWave
setdimlabel 1,2, Frame3, $newWave
setdimlabel 1,3, Frame4, $newWave
setdimlabel 1,4, Frame5, $newWave
setdimlabel 1,5, Frame6, $newWave
setdimlabel 1,6, Frame7, $newWave
setdimlabel 1,7, Frame8, $newWave
setdimlabel 1,8, Frame9, $newWave
setdimlabel 1,9, Frame10, $newWave
variable i, nLayers
nLayers = dimsize($newWave, 2)
for (i=0; i < nLayers; i+=1)
string layerName = "episode"+ num2str(i)
setdimlabel 2, i, $layername, $newWave
endfor
print getdimlabel($newWave, 2, 1)
end
Function GetPeak(AvgWave)
//The function will take an averaged 2D wave of spectra and extract rows corresponding to CFP and YFP peaks.
wave avgWave
variable i, n=80, p=n+10
for (i = n; i<= p; i += 1)
Imagetransform/G=(i) getrow avgWave // getrow creates a 1D wave consisting of all points in a selected row,second command overwrites the existing 1D wave
string CFPName= "CFP_peak" + nameofwave(W_extractedRow) + num2str(i)
rename W_ExtractedRow $CFPName
endfor
variable j, m=240, q=m+10
for (j = m; j<= q; j += 1)
Imagetransform/G=(j) getrow avgWave // getrow creates a 1D wave consisting of all points in a selected row,second command overwrites the existing 1D wave
string YFPName= "YFP_peak" + nameofwave(W_extractedRow) + num2str(j)
rename W_ExtractedRow $YFPName
endfor
end
Function Plot(inwave)
//Takes the Row containing the fluorescence intensity and creates a 2D wave.
wave wName
variable colToTake
variable Nrows
variable Ncols
variable Nlayers
Nrows = DimSize(wName, 0)
Ncols = DimSize(wName, 1)
Nlayers = DimSize(wName, 2)
string NewWname = nameofwave(wName)+"_R"+num2str(coltotake)
make/O/N =(Nrows,Nlayers) wave2d
wave2d= wName[p][colToTake][q]
rename wave2d $NewWname
end
Function MakeMatrix(w2d, FramePerEpisode)
//The function will take a 2D wave containing all the frames as columns and redimension it into a 3d matrix
//with each frame of each episode aligned in the z-plane.
wave w2d
variable FramePerEpisode
variable Ncols
Variable NEps
variable Nrows
variable sizeOf1D
duplicate/O w2d, tempWave
Nrows = DimSize(w2d, 0)
Ncols = DimSize(w2d, 1)
NEps = Ncols / FramePerEpisode
sizeof1D = Nrows*Ncols
redimension/N=(sizeof1d) tempwave
redimension/N=(512, frameperepisode, Neps) tempwave
string newWave = nameofwave(w2d)+"Matrix"
rename tempwave $newWave
setdimlabel 1,0, Frame1, $newWave
setdimlabel 1,1, Frame2, $newWave
setdimlabel 1,2, Frame3, $newWave
setdimlabel 1,3, Frame4, $newWave
setdimlabel 1,4, Frame5, $newWave
setdimlabel 1,5, Frame6, $newWave
setdimlabel 1,6, Frame7, $newWave
setdimlabel 1,7, Frame8, $newWave
setdimlabel 1,8, Frame9, $newWave
setdimlabel 1,9, Frame10, $newWave
variable i, nLayers
nLayers = dimsize($newWave, 2)
for (i=0; i < nLayers; i+=1)
string layerName = "episode"+ num2str(i)
setdimlabel 2, i, $layername, $newWave
endfor
print getdimlabel($newWave, 2, 1)
end
Function GetPeak(AvgWave)
//The function will take an averaged 2D wave of spectra and extract rows corresponding to CFP and YFP peaks.
wave avgWave
variable i, n=80, p=n+10
for (i = n; i<= p; i += 1)
Imagetransform/G=(i) getrow avgWave // getrow creates a 1D wave consisting of all points in a selected row,second command overwrites the existing 1D wave
string CFPName= "CFP_peak" + nameofwave(W_extractedRow) + num2str(i)
rename W_ExtractedRow $CFPName
endfor
variable j, m=240, q=m+10
for (j = m; j<= q; j += 1)
Imagetransform/G=(j) getrow avgWave // getrow creates a 1D wave consisting of all points in a selected row,second command overwrites the existing 1D wave
string YFPName= "YFP_peak" + nameofwave(W_extractedRow) + num2str(j)
rename W_ExtractedRow $YFPName
endfor
end
Function Plot(inwave)
I'm stuck at the end of the MakeMatrix function.
Any suggestions will be greatly appreciated! Thanks in advance!
Start by figuring out how to rearrange your data in the best way. I can't say that I followed all the manipulations you made as they are unlikely to be necessary. I assume that you want to average either rows or columns of a layer. The easiest way to write the code is using MatrixOP. Specifically,
From here you could average rows or columns using MatrixOP with sumRows, sumCols, etc.
Another alternative is if your data are arranged as "beams" in which case you would use:
If you prefer you could also use either MatrixOP or ImageTransform with transposeVol to convert your 3D data packing to one that makes your subsequent data manipulation easier.
I hope this helps,
A.G.
WaveMetrics, Inc.
June 24, 2013 at 11:32 am - Permalink
wave Wave3D
variable NEp = dimsize(Wave3D, 2)
variable i,j, n=2, k=NEp/n
for (i=0; i<NEp; i+=1)
for (j=0; j<k; j+=1)
if(i==j*n)
print "even"
elseif(i==j*n-1)
print "odd"
endif
endfor
endfor
end
but it doesn't seem to work. Print cmd is just to check if the loops are working ok. For a matrix with 3 layers the function prints 3x"even".
Any suggestions would be greatly appreciated,
Chris
July 1, 2013 at 09:35 am - Permalink
Wave inWave
Variable layers=DimSize(inWave,2)
Variable nrows=DimSize(inWave,0)
Variable i
for(i=0;i<layers;i+=1)
MatrixOP/O/Free theLayer=inWave[][][i]
MatrixOP/O/Free avgCols=sumcols(theLayer)/nrows
if(i&1) // odd
// do something here with the avgCols
else // even
// do something here with the avgCols
endif
endfor
End
A.G.
WaveMetrics, Inc.
July 1, 2013 at 10:41 am - Permalink
Wave inWave
Variable layers=DimSize(inWave,2)
Variable i
for(i=0;i<layers;i+=1)
if(i&1) // odd
MatrixOP/O/Free OddLayers=inWave[][][i]
print "odd" // do something here with the avgCols
else // even
MatrixOP/O/Free EvenLayers=inWave[][][i]
print "even" // do something here with the avgCols
endif
endfor
MatrixOP/O/Free avgCols=sumbeams(oddLayers)/layers
MatrixOP/O/Free avgCols=sumbeams(evenLayers)/layers
End
<pre><code class="language-igor">
but when this runs the destwave is not created (the waves odd layers and evenLayers). Did I misunderstand MatrixOp?
July 3, 2013 at 04:29 am - Permalink
MatrixOP/O/Free avgCols=sumbeams(evenLayers)/layers
The flag "/Free" makes a free wave that disappears when, in this case, the function exits. It does not appear in the data browser as well. Try removing this flag.
July 3, 2013 at 05:02 am - Permalink
July 3, 2013 at 05:52 am - Permalink
No. The MatrixOP line will only extract one layer at a time and put it in a 2D wave. If you remove the "/O" flag, Igor will complain the second time through the loop that the waves already exist. You may need to do a running sum in the loop or create two 3D output waves each with half the layers of the original and store the odd and even layers there while the loop runs.
Alternatively, you might abandon MatrixOP and extract the odd and even layers using the following syntax
EvenLayers[][][] = InWave[p][q][2 * r]
Note that "p", "q" & "r" are the "special" functions that refer to row, column and layer indices for the destination wave.
July 3, 2013 at 07:11 am - Permalink
Well, if you are going to "abandon" MatrixOP why not use something more efficient than straight wave assignment. For example, ImageTransform with getPlane and setPlane.
AG
July 3, 2013 at 09:42 am - Permalink
If the OP needs a 3D data set is it more efficient to assemble it with a single wave assignment or a loop using MatrixOP or ImageTransform? My guess was that it would be more efficient without the loop. I really meant abandon the loop rather than abandon MatrixOP.
Looking back at my previous post, I see that the meaning was not clear, but, aparently, I knew what I meant.
July 3, 2013 at 10:59 am - Permalink
Best,
Chris
July 4, 2013 at 10:04 am - Permalink