Convert Triplet data from grid into Matrix

The snippet here pertains to a recent question about displaying Gizmo surface plots from triplet data (see Exchange node/3867). John Weeks suggested that often the triplet data is created from a matrix grid, and needs only rearrangement to be able to visualize it as a matrix surface. The function below, gridTripletToMatrix(triplet_wave), takes such triplet data and converts it to a matrix (without any interpolation). The function does not require any particular ordering of the triplet (xyz) sequences, although this might often occur in a matrix-like order. If the original data is in the form of three x, y, and z waves it should be converted to a triplet wave. An alternative is to modify the code below to accept three 1D waves directly (which it does internally from the triplet wave). The first function foo() tests the conversion function for a triplet wave which has a randomly altered sequence order. The conversion function should be treated more as a code example than a bullet-proof WaveMetrics-class utility.
#pragma rtGlobals=1     // Use modern global access method.
#include <MatrixToXYZ>

//  Generate simulated test data, and put it in a triplet wave with random row order.
//  Test the matrix rearrangement function, comparing it to the original matrix that
//  was used to create the simulated data. These can also be compared in Tables
//  or Gizmo plots.
function foo()    
    variable Nx0 = 14, Ny0 = 12
    make/O/N=(Nx0,Ny0) wave0    //    form a gridded original source
    setscale/P x, -0.1*Nx0/2, 0.1,"" wave0
    setscale/P y, -0.1*Ny0/2, 0.1,"" wave0
    wave0 = exp( - 4*(x^2+y^2)) + gnoise(0.025)    //    original 'data' quad surface wave
    Execute "MatrixToXYZTriplet(\"wave0\",\"triplet\",2)"
    wave triplet    //    in gridded sequence at this point

    variable Ntriplet = DimSize(triplet,0)
    Make/O/N=(Ntriplet) wx, wy, wz         //    convert to xyz 1D waves
    MatrixOP/O wx = col(triplet,0)
    MatrixOP/O wy = col(triplet,1)
    MatrixOP/O wz = col(triplet,2)

    shuffle(Ntriplet)                                 //    random permutation key
    wave wPerm
    Sort wPerm, wx, wy, wz                   //    randomize order of 1D components

    Make/O/N=((Ntriplet), 3 ) wRanTrip    //    re-assemble triplet wave, in random order
    wRanTrip[][0]   =   wx[p]
    wRanTrip[][1]   =   wy[p]
    wRanTrip[][2]   =   wz[p]                  //   simulated triplet data 

    gridTripletToMatrix(wRanTrip)    //    TEST THE SORTING AND GRIDDING FUNCTION
    wave wz2D                             //     the resulting 2D matrix wave
    MatrixOP/O wdiff = wave0 - wz2D
    MatrixOP/O var   = (1/(Nx0*Ny0)) * sum(wdiff*wdiff)
    printf "wdiff Variance = %g\r", var[0]  // error between original and rearranged test
end
//-----------------------------------------------------------------------------------------------------------------
//  The input must be an (N,3) triplet wave, whose x and y components are on a
//   rectangular grid, so each (x,y) pair is unique.    The dim-0 ordering of each triplet can
//   be arbitrary, since the function will sort the row order and find the grid dimensions.
//   The resulting 2D wave containing the 'z' values on a scaled x-y grid is 'wz2D'.   
//    Internal waves are created as /FREE; change this if you wish to inspect them later.
//    There is no error checking on the input wave to ensure the triplet criteria ! ! !
function gridTripletToMatrix(triplet_wave)    //    sort the triplet order, and grid the data
    wave triplet_wave                                 //    triplet input wave in arbitray order
    variable Npts = DimSize(triplet_wave,0)
    Make/O/N=(Npts)/FREE wxR, wyR, wzR   //    convert to x,y,z
    MatrixOP/O wxR = col(triplet_wave,0)
    MatrixOP/O wyR = col(triplet_wave,1)
    MatrixOP/O wzR = col(triplet_wave,2)
    sort wyR, wxR, wyR, wzR                   //    sort all waves by y values
    variable Nx, Ny, i
    variable tol = 1e-6
    for(i=0;i<Npts;i+=1)                            //    find the dimension boundary
        if ( wyR[i+1]> (wyR[i] + tol* abs(wyR[i])) )  //  allow for numeric error
            Nx = i+1                                   //    Nx is the x dimension size
            break
        endif
    endfor
    Ny = Npts / Nx                                       //    Ny is the y dimension size
    Make/O/N=0/FREE     wxB, wyB, wzB    //    initial 'big' target wave for concatenating chunks
    Make/O/N=(Nx)/FREE wxS, wyS, wzS   //    Nx chunks ('S' for 'small')
    for(i=0; i<Ny; i+=1)
        wxS = wxR[p+i*Nx];    wyS = wyR[p+i*Nx];    wzS = wzR[p+i*Nx]  
        Sort wxS, wxS, wzS                         //    sort x and z in each chunk wave by x values
        Concatenate/NP  {wxS}, wxB;  Concatenate/NP  {wyS}, wyB;  Concatenate/NP  {wzS}, wzB
    endfor
    Make/O/N=((Nx*Ny), 3 )/FREE wSortedTriplet    //    triplet wave, grid sorted
    wSortedTriplet[][0]  =  wxB[p]
    wSortedTriplet[][1]  =  wyB[p]
    wSortedTriplet[][2]  =  wzB[p]

    MatrixOP/O wz2D = col(wSortedTriplet,2)    //    sorted z-values in a 1D wave
    Redimension/N=(Nx,Ny) wz2D
//   Make/O/N=(Nx, Ny) wz2D   //  alternate method:set up the rectangular gridded wave
//   wz2D  =  wz1D[p+Nx*q]      //  for 'z' values, and"unfold" 'z' values from 1D wave into 2D wave
    variable dx = wSortedTriplet[1][0]   - wSortedTriplet[0][0]
    variable dy = wSortedTriplet[Nx][1] - wSortedTriplet[0][1]
    variable FirstX = wSortedTriplet[0][0] ,  LastX = wSortedTriplet[Nx-1][0]    + dx
    variable FirstY = wSortedTriplet[0][1] , LastY = wSortedTriplet[Npts-1][1] + dy
    SetScale x, FirstX, LastX, ""  wz2D  //  using SetScale/I or/P gives scale precision errors
    SetScale y, FirstY, LastY,""  wz2D 
end
//-----------------------------------------------------------------------------------------------------------------
function shuffle(Npts)    //    perform a random permutation of indices 0...Npts-1
    variable Npts
    make/O/N=(Npts) wPerm = p    //   will be used as sort key; initialize it
    variable i, j, temp
    for(i=Npts-1; i>0; i-=1)
        temp = wPerm[i]
        j = round(  i/2+enoise(i/2)  )   //   pick random index from remaining choices
        wPerm[i] = wPerm[j]             //  swap entries
        wPerm[j] = temp
    endfor                                     //   wPerm is the final random sort key wave
end
If the X or Y values are already in grid order, you can use the procedure that ships with Igor to convert separate x, y and z waves into a matrix:
#include <XYZtoMatrix>

//
// XYGridandZtoMatrix rearranges three waves containing X, Y, and Z values
// that comprise a regularly-spaced grid of X and Y values
// (already sorted in either column-major or row-major order)
// into a matrix of Z values that spans the min and max X and Y.
Macro XYGridandZtoMatrix(wx,mat,wy,mktbl,wz,mkimg)
...


--Jim Prouty
Software Engineer, WaveMetrics, Inc.

Forum

Support

Gallery

Igor Pro 8

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More