Resampling a matrix with irregularly-spaced X and Y coords
Nasser
Wed, 10/26/2011 - 10:28 am
I'd like to resample my data, consisting of a 2D matrix. the stamps (X waves) are known for each of the dimensions.
ImageInterpolate does not help me because I can't specify the X waves in each dimension.Resample allows me to specify /W={X_wave1, X_wave2}, but it seems to works only for bilinear interpolation. Thus, I'm trying to resample using in both dimension, succesively. Here is for the 1st dimension (it fails)
Function Resample_Spectrum(spectre_source, str_Resampled_spectrum, str_Xwave)//ex:Resample_spectrum(SigPos0, "Resampled_spectrum", "RT")
WAVE spectre_source
string str_Resampled_spectrum
string str_Xwave
WAVE Xwave = $str_Xwave
variable i
Duplicate/O spectre_source, $str_Resampled_spectrum
WAVE Resampled_spectrum = :$str_Resampled_spectrum
for(i=0; i<DimSize(spectre_source,1); i+=1)
//We extract the profile that is going to be re-sampled
MatrixOP /O /S temp_profile = col(spectre_source, i)
Interpolate2 /T=2 /N=(474) /E=2 /Y=interp_temp_profile $str_Xwave, temp_profile
Resampled_spectrum[][i] = interp_temp_profile
endfor
End
WAVE spectre_source
string str_Resampled_spectrum
string str_Xwave
WAVE Xwave = $str_Xwave
variable i
Duplicate/O spectre_source, $str_Resampled_spectrum
WAVE Resampled_spectrum = :$str_Resampled_spectrum
for(i=0; i<DimSize(spectre_source,1); i+=1)
//We extract the profile that is going to be re-sampled
MatrixOP /O /S temp_profile = col(spectre_source, i)
Interpolate2 /T=2 /N=(474) /E=2 /Y=interp_temp_profile $str_Xwave, temp_profile
Resampled_spectrum[][i] = interp_temp_profile
endfor
End
Only the first column seems to be correctly interpolated. Then, all rows of the other columns are filled with the value of the last row of the column.
Is there a smarter way to resample a matrix according to X_waves ?
Thank you in advance for your support.
It seems like ImageInterpolate with the XYWaves keyword and /W={xWave, yWave} flag should do the trick.
Note that an Igor 2D wave has an X dimension and a Y dimension, not two X dimensions.
October 26, 2011 at 11:26 am - Permalink
Thank you for your support.
ImageInterpolatedoes not enable to specify the number of samples in the interpolated matrix, like withInterpolate2and the /N flag.So I tried with /F={1,1} and there was no difference between the raw and the interpolated matrices.
Then, with /F={2,2}, the data seems to be interpolated, but it is impossible to display the matric as an image (Table is OK). Even if I redefine th scaling, the problm is still the same. It seems like the index are not integers (!) Would you have an idea please ?
I've tried this :
//ex: Resample_Spectrum(SigPos0, "Interpolated_spectrum", "RT", "CV")
WAVE spectre_source
string str_Interpolated_spectrum
string str_Xwave, str_Ywave
WAVE Interpolated_spectrum = $str_Interpolated_spectrum
WAVE Xwave = $str_Xwave
WAVE Ywave = $str_Ywave
Duplicate/O spectre_source, Interpolated_spectrum
ImageInterpolate /DEST=Interpolated_spectrum /W={$str_XWave, $str_YWave} /f={2,2}/D=3 /FDS spline spectre_source
SetScale /I y, 0, 500, Interpolated_spectrum
SetScale /I x, -30, 10, Interpolated_spectrum
End
October 26, 2011 at 02:26 pm - Permalink
If you want me to investigate further I will need some sample data to play around with including your matrix and X and Y waves. Post an experiment containing your waves and function and the command you tried or send it to Wavemetrics support.
On another topic, your function takes three input wave names and one output wave name. Two of the input waves are specified as string parameters. The three input waves should all be passed as wave references, not as strings. Since the waves must exist, you can define the function to take wave references. The output wave name parameter must be passed as a string because the output wave may not exist when the function is called.
Also, this statement is problematic:
The problem is that the output wave, whose name is in str_Interpolated_spectrum, may not exist at this point. In that case, the Wave statement will fail and you will have a null wave reference. If you have WAVE,NVAR,SVAR Checking on (Procedure menu), which is a good idea, Igor will drop into the debugger at that point.
I would rewrite your function like this:
Function Resample_Spectrum(spectre_source, xWave, yWave, str_Interpolated_spectrum)
WAVE spectre_source // Input 2D wave
Wave xWave, yWave // Input XY 1D waves
String str_Interpolated_spectrum // Contains name of output wave
ImageInterpolate /DEST=$str_Interpolated_spectrum /W={xWave,yWave} /F={2,2} /FDS XYWaves, spectre_source
Wave output = $str_Interpolated_spectrum // Create wave reference after the wave exists
SetScale /I y, 0, 500, output
SetScale /I x, -30, 10, output
End
The SetScale calls at the end looks suspicious. Is that really the range of the output matrix? I would need to play with the data to be sure.
October 26, 2011 at 08:27 pm - Permalink
But the spline method is more accurate than the linear interpolation. I am going to test anyway. Hopefully the difference will not be noticeable.
Thank you. If I'm still stuck, I will attach an experiment with those waves.
Prior to calling the Resampling function, the waves are created. So is it really prohibited to use strings instead of waves ? Is there unexpected and unnoticeable results if I do so ?
Actually I wrote this in many different ways, beacause I suspected that the direct using of strings in
ImageInterpolateorInterpolate2led to syntax problems...Anyway, I was wrong.That's not really the range, I took round values. Actually, the range is defined by the min and max values of the X and Y 1D waves (containing irregularly spaced scalars) associated to the matrix that we want to resample. This range may also be different from a dataset to the other (that's another good reason to resample : I think it will facilitate image registration (alignment)).
Thanks a lot. I will try again, and keep you informed.
Best Regards.
October 27, 2011 at 05:51 am - Permalink
In order to resample the matrix
SigPos0, I've interpolated it using (inResample_Spectrum2(...)):ImageInterpolate /DEST=$str_Interpolated_spectrum /S={0,1,474,0,1,100} /W={xWave,yWave} XYWaves spectre_source
as you suggested. To do so, I also had to add a NaN point in the 1D X and Y waves (since they must have one more point than the data).
As you can see, the result in a truncated matrix...
Ideally, I would like to keep the same number of samples (with
ImageInterpolate, it may be done with/RESL={DimSize(spectre_source,0),DimSize(spectre_source,1)}), since the next step will be to interpolate the resampled matrix using splines.This future interpolation (
Interpolate_Spectrum(...)) will correctly reconstruct the signal, since the samples will be evenly spaced.I have also tried this :
ImageInterpolate /DEST=$str_Interpolated_spectrum /U=1 /W={xWave,yWave} /FDS /RESL={DimSize(spectre_source,0),DimSize(spectre_source,1)} splineBut the
/W={xWave,yWave}flag seems to be ignored, which is actually the most important parameter.Your help will be really appreciated !
Thanks in advance.
October 27, 2011 at 06:38 am - Permalink
Do you mean that the output wave must be created before calling the function? If so then you can pass the output wave as a wave reference. However, it is generally best to allow the function to create the output wave since the calling routine, in general, may not know how to create it (how many points, what data type, etc.).
It's not prohibited and it will not create a problem but it is also unnecessary and inelegant. You could pass numeric arguments as strings also and convert them to numbers using str2num. It would work but you would not be using the most appropriate parameter type.
October 27, 2011 at 08:15 am - Permalink
Don't add a NaN. Add a reasonable value. Interpolation needs to know the XY position of the start and end of each interval. If you have a 1x1 matrix, you have two X positions, the start of row 0 and the end of row 0. So you need two X positions for a 1 row matrix. In general you need n+1 X positions for an n-row matrix.
An alternative is to not add any points to your XY waves but instead delete one row and one column from your matrix.
As I said in an earlier post, the /W flag is supported for the XYWaves method only.
I have written a procedure that does what I think you want to do. It is in the attached experiment.
The procedure, ImageInterpXY, takes wave parameters specifying your image wave and your original X and Y waves.
The output image wave name is determined algorithmically by the function. For input SigPos0, the output image wave will be named SigPos0Interp.
ImageInterpXY creates tempXWave and tempYWave with one extra point. The extra point is a linear extrapolation of the last two points of the original waves. It kills the temporary waves when it is finished with them.
ImageInterpXY calculate the appropriate values for the /S flag.
It returns a wave reference which you can use in the calling function, if any.
October 27, 2011 at 10:55 am - Permalink
First of all, thank you for the time you spend on this topic.
Since the aim is to convert a matrix (with unevenly spaced samples on the X and Y waves), into a scaled matrix (with evenly spaced samples),
I would be tempted to scale the matrix like that :
SetScale/P y (y0),(dy),"", wOut
However, it is written in the manual that :
"The data domain is defined between the centers of the first and last pixels"
xmin=(xWave[0]+xWave[1])/2
xmax=(xWave[last]+xWave[last-1])/2
Is x0 still xWave[0] ?
Is endX still xWave[numXPoints-1] ?
Is dx stil (endX - x0) / numXPoints ?
I don't think so. I think it would be correct to write :
variable xlast_scale = (xWave[*]+xWave[*-1])/2
variable y0_scale = (yWave[0]+yWave[1])/2
variable ylast_scale = (yWave[*]+yWave[*-1])/2
SetScale/I x (x0_scale),(xlast_scale),"", wOut
SetScale/I y (y0_scale),(ylast_scale),"", wOut
To simplify, let's consider profiles, ignoring the other dimension.
To validate this, it would be nice to superimpose :
- a profile extracted from the original image vs. X and Y
- and a profile extracted at the same column, without any X and Y, but scaled
Unfortunately, this validation seems to be tricky, because the 1st index of the interpolated profile corresponds to the "1st 1/2" index on the original profile
-->No way to superimpose them.
Moreover, (still for validation purpose) one may not perform any interpolation one the other dimension, in order to ensure that the extracted columns a representing the same slice.
It is funny how it can be possible to go crazy with very simple mathematical operations !
When I tried
Interpolate2(that I actually discovered from the "XY Pair to Waveform" menu), I was really happy because it did exactly what I needed, without any problems, but only on 1D. Unfortunatelly, the "XYZ to Matrix" doesn't fit to my issue. That's why my first attempt was to perform Interpolate2 in both directions.Would you have some ideas about that please ?
Thank you
Your code was :
Function/WAVE CreatePlusOneWave(wIn, outputWaveName)
Wave wIn
String outputWaveName
Duplicate/O wIn, $outputWaveName
Wave wOut = $outputWaveName
Variable numPointsIn = numpnts(wOut)
InsertPoints numPointsIn, 1, wOut
// Set value of extra point assuming based on preceeding points
Variable lastPoint = numPointsIn
wOut[lastPoint] = wOut[lastPoint-1] + (wOut[lastPoint-1] - wOut[lastPoint-2])
return wOut
End
//-------------------------------------
// Example:
// ImageInterpXY(SigPos0, RT, CV) // Creates SigPos0Interp
Function/WAVE ImageInterpXY(imageWave, xWave, yWave)
Wave imageWave
Wave xWave // Assumed to have the same number of points as imageWave has rows
Wave yWave // Assumed to have the same number of points as imageWave has columns
Wave tempXWave = CreatePlusOneWave(xWave, "ImageInterpTempXWave")
Wave tempYWave = CreatePlusOneWave(yWave, "ImageInterpTempYWave")
String outputWaveName = NameOfWave(imageWave) + "Interp"
Variable numXPoints = DimSize(imageWave, 0)
Variable x0 = xWave[0]
Variable endX = xWave[numXPoints-1]
Variable dx = (endX - x0) / numXPoints
Variable numYPoints = DimSize(imageWave, 1)
Variable y0 = yWave[0]
Variable endY = yWave[numYPoints-1]
Variable dy = (endY - y0) / numYPoints
ImageInterpolate /DEST=$outputWaveName /S={x0,dx,endX,y0,dy,endY} /W={tempXWave,tempYWave} XYWaves imageWave
KillWaves/Z tempXWave, tempYWave
Wave wOut = $outputWaveName
return wOut
End
October 28, 2011 at 07:34 am - Permalink
It makes sense that it would be the same as what you specified using the /S flag as I believe /S tells ImageInterpolate at what XY locations to do each interpolation.
Here is my procedure with the SetScales added:
Wave wIn
String outputWaveName
Duplicate/O wIn, $outputWaveName
Wave wOut = $outputWaveName
Variable numPointsIn = numpnts(wOut)
InsertPoints numPointsIn, 1, wOut
// Set value of extra point assuming based on preceeding points
Variable lastPoint = numPointsIn
wOut[lastPoint] = wOut[lastPoint-1] + (wOut[lastPoint-1] - wOut[lastPoint-2])
return wOut
End
// Example:
// ImageInterpXY(SigPos0, RT, CV) // Creates SigPos0Interp
Function/WAVE ImageInterpXY(imageWave, xWave, yWave)
Wave imageWave
Wave xWave // Assumed to have the same number of points as imageWave has rows
Wave yWave // Assumed to have the same number of points as imageWave has columns
Wave tempXWave = CreatePlusOneWave(xWave, "ImageInterpTempXWave")
Wave tempYWave = CreatePlusOneWave(yWave, "ImageInterpTempYWave")
String outputWaveName = NameOfWave(imageWave) + "Interp"
Variable numXPoints = DimSize(imageWave, 0)
Variable x0 = xWave[0]
Variable endX = xWave[numXPoints-1]
Variable dx = (endX - x0) / numXPoints
Variable numYPoints = DimSize(imageWave, 1)
Variable y0 = yWave[0]
Variable endY = yWave[numYPoints-1]
Variable dy = (endY - y0) / numYPoints
ImageInterpolate /DEST=$outputWaveName /S={x0,dx,endX,y0,dy,endY} /W={tempXWave,tempYWave} XYWaves imageWave
KillWaves/Z tempXWave, tempYWave
Wave wOut = $outputWaveName
SetScale/P x, x0, dx, "", wOut
SetScale/P y, y0, dy, "", wOut
return wOut
End
October 28, 2011 at 07:56 am - Permalink
I agree, it makes sense to do the scaling using the same values thant those used for the /S flag. But in this case, don't you think that x0, dx, yo and dy should be changed ? Because the coordinate of the first sample in the interpolated 1D waveform is not the same as the coordinate of the first sample in the raw 1D wave. The first new sample is located between the two first samples of the raw 1D wave.
If so, everything will be shifted, and the result would be worse than before interpolation (provided that I'm right).
If interpolate2 could be used, it would simply be great ! It works perfectly on 1D waves (I've checked by superimposing profiles), and the scaling is automatically assigned thanks to the /N flag.
October 28, 2011 at 08:26 am - Permalink
It depends on how your XY data is to be interpreted. If your X wave represents the position of the left edge of each pixel then the way I did it is right. If your X wave represents the position of the the center of each pixel then your way is correct.
But it also depends on how ImageInterpolate interprets the /S flag - whether it represents the edges or center of the pixels. I'm don't know the answer to that.
I do know that the x0 parameter that you supply to SetScale/P represents the left edge of the first pixel.
October 28, 2011 at 08:59 am - Permalink
In this case, the interpolated samples are located between the original ones. Therefore, the linear interpolation is less accurate than if the pixels themselves, and not their edges, were considerated. (BTW These data represent chemical compounds at low concentrations)
Apparently, Interpolate2 allows to perform the same thing more easily, since the scaling step is included in this function.
Thank you for the info regarding Scale /P, I didn't know.
Could you please have a look to my attempt to implement interpolate2 ? Then it is the week-end, I promise I'll stop boring you ;)
Thank you
Best Regards
October 28, 2011 at 01:48 pm - Permalink
ImageInterpolate, like the display of images in Igor, needs to know the edges of each pixel, not the center. To see why, consider a 1-pixel image with the pixel center at x=0. How wide is the pixel? You can not say - it could be any width. The center does not determine the width. The same is true for a 2-pixel image or for an n-pixel image - specifying the centers does not determine where the edges are.
I thought that perhaps knowing the center of each pixel AND the width of the first pixel would permit determining the edges of each pixel. I think this might be true. However, I tried it with your X data and assuming that the width of the first pixel RT[1] - RT[0] where RT is your X wave. This led to a non-monotonic edge wave, meaning that the right edge of at least one pixel was to the left of that pixel.
My conclusion is that I don't have enough information to create an edge wave from your center wave and so I can't solve the problem.
For the reasoning explained above, to get an accurate output you need to know the edges, assuming that the width of the pixels has physical significance.
What do the X and Y waves represent?
The main problem with your function was this:
It needed to be this:
Without the added [p], meaning "current row", the assignment statement was reading from the current row and column. Since the right side has only one column, Igor clipped the index and returned the last point in interp_temp_profile.
As part of understanding your function to debug it, I changed some of your names, changed the order of the parameters so that the inputs come before the output, and changed your string parameter for the X wave to a Wave parameter. The result is this:
WAVE inputMatrix
Wave xCenterWave
String outputMatrixName
Duplicate/O inputMatrix, $outputMatrixName
WAVE outputMatrix = :$outputMatrixName
Variable numRows = DimSize(inputMatrix,0)
Variable numColumns = DimSize(inputMatrix,1)
Variable column
for(column=0; column<numColumns; column+=1)
//We extract the profile that is going to be re-sampled
MatrixOP /O /S /FREE temp_profile = col(inputMatrix, column)
Interpolate2 /T=2 /N=(numRows) /E=2 /Y=interp_temp_profile xCenterWave, temp_profile
outputMatrix[][column] = interp_temp_profile[p]
endfor
KillWaves/Z interp_temp_profile
End
Note that this solution takes the X positions of the pixels into account but ignores their widths and their Y positions and heights. I suspect that the X widths and Y positions and heights are irrelevant for your purposes. In this case, you don't really have an image, you have a set of 1D XY pairs. If that is correct then using Interpolate2 is the right solution.
October 30, 2011 at 10:20 am - Permalink
I agree. The problem cannot be properly resolved in this way. But fortunatelly, as you guessed, this is not really an image...
Y is a voltage, X is the elapsed retention time of a GC column, and the matrix is a voltage readout (corresponding to a concentration) at the coords (X,Y). Actually, while the column elutes, a DC voltage sweeps from, let's say -40 to +10 V, between each time stamp.
Thank you very much for your support. I think that this topic can now be marked as "Solved" !
Best Regards
October 31, 2011 at 02:46 am - Permalink