3D matrix rearrangement

Share icon

Fri, 06/21/2013 - 02:03 am

Hi all,

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)


I'm stuck at the end of the MakeMatrix function.

Any suggestions will be greatly appreciated! Thanks in advance!
Hello Chris,

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,
MatrixOP/O sampleLayer=my3DWave[][][layerNumber]

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:
MatrixOP/O myBeam=beam(my3DWave,row,col)

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

Actually I want to average the columns of all the even-numbered layers and all the odd-numbered layers. So I'm trying to write a function that separates my initial 3D matrix into 2 3D matrices: one with odd numbered layers and one with the even-numbered layers. I was thinking of making a loop like this:

function Matrix(wave3D)

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

I suggest a slightly different approach:

Function average3d(inWave)
    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

Thanks a lot for the suggestion! I changed the code a bit to this:
Function average3d(inWave)
    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

Thanks! Another question: shouldn't MatrixOP/O OddLayers=inWave[][][i] create a 3D output wave, combining all the odd-numbered layers? At this point I get a 1D wave, I'm guessing with date from the last layer it cycled through. Is there a way to overcome this? Removing the /O flag?

July 3, 2013 at 05:52 am - Permalink

L.K.Chris wrote:
Thanks! Another question: shouldn't MatrixOP/O OddLayers=inWave[][][i] create a 3D output wave, combining all the odd-numbered layers? At this point I get a 1D wave, I'm guessing with date from the last layer it cycled through. Is there a way to overcome this? Removing the /O flag?


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

OddLayers[][][] = InWave[p][q][2 * r + 1]
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

jtigor wrote:

Alternatively, you might abandon MatrixOP and extract the odd and even layers using the following syntax


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

Igor wrote:
jtigor wrote:

Alternatively, you might abandon MatrixOP and extract the odd and even layers using the following syntax


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


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