Find peaks in 2D-Spectra and calculate peak-volume?

Share icon

Fri, 09/04/2015 - 02:23 am

Dear all,

I have some three-dimensional peaks on a plane (two-dimensional spectrum) and wondered whether there is an existing macro to automatically find these peaks and calculate their peak volumes?

Maybe someone stumbled over this in the past or is actively working with something like this?

If not, what would be the best general way to approach this problem?

Especially I wondered what would be the best (mathematically sound) way to calculate peak volumes?

Thanks for any help!

Best regards,
Peter
Well, here is a function I wrote for another customer that uses ImageAnalyzeParticles to find spots in a 2D wave. It then fits a two-dimensional Gaussian to each spot it finds. The coefficients for the 2D Gaussian are stored in a multi-column wave, one row per spot. You would then need to figure out how to compute the volume for each spot.

Note that this won't do the right thing if the spots overlap. If there are overlapping spots, the problem gets much more complicated.

Function DoParticleAnalysis(imageW)
    Wave imageW
   
    // First, find the spots
    ImageThreshold/I/M=0/T=800/Q imageW
    Wave M_ImageThresh
    ImageAnalyzeParticles /E/W/Q/M=3/A=3 stats, M_ImageThresh
   
    Wave W_ImageObjArea
    Variable nSpots = DimSize(W_ImageObjArea, 0)
   
    // Now we have an estimate of the spot locations- we can now loop through the spots and do fits
    Wave W_xmin, W_xmax, W_ymin, W_ymax
    WaveStats/Q W_ImageObjArea
    Variable maxArea = V_avg+3*V_sdev   // try to eliminate obviously bad "spots"
    Variable i
    Variable V_fitOptions=4
    Variable V_fitError
    Make/N=(nSpots, 7)/O $(nameOfWave(imageW)+"fits")
    Wave results = $(nameOfWave(imageW)+"fits")
    results = nan
    for (i = 0; i < nSpots; i += 1)
        if (W_ImageObjArea[i] > maxArea)
            continue        // this spot looks like a bad one- go around for another
        endif
       
        V_fitError = 0
        // Fit a 2D gaussian to just one spot. The -3, +3 bit is to include sufficient background to get a reliable fit
        CurveFit/N/Q Gauss2d imageW[W_xmin[i]-3, W_xmax[i]+3][W_ymin[i]-3, W_ymax[i]+3]
        Wave W_coef
        if (V_fitError == 0)
            results[i][] = W_coef[q]
        endif
    endfor
end


John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

September 8, 2015 at 11:46 am - Permalink

The image analysis expert has pointed out to me that the function has a hard-coded threshold value of 800. You will need to adjust the code for your data.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

September 9, 2015 at 10:21 am - Permalink