dm-script - 从线轮廓识别峰

Question

我想问是否可以在 DM 上定位强度分布的每个最大值和最小值的位置。

我如何想出一个简单的脚本来识别下面示例中峰的位置？

这是 STEM 图像沿 Y 方向的线强度分布的屏幕截图：

Answer 1

如果要过滤“严格”局部最大值，则可以使用图像表达式和条件“tert”运算符轻松完成此操作。以下示例说明了这一点：

image CreateTestSpec( number nChannels, number nPeaks, number amp, number back )
{
    image testImg := RealImage( "TestSpec", 4, nChannels )
    testImg = amp * cos( PI() +  2*PI() * nPeaks * icol/(iwidth-1) )
    testImg += back
    testImg = PoissonRandom( testImg )
    return testImg
}

image FilterLocalMaxima1D( image spectrumIn, number range )
{
    image spectrumOut := spectrumIn.ImageClone()
    for( number dx = -range; dx<=range; dx++  )
        spectrumout *= ( spectrumIn >= offset(spectrumIn,dx,0) ? 1 : 0 )

    spectrumout.SetName("Local maxima ("+range+") filtered")
    return spectrumOut
}

image test1 := CreateTestSpec(256,10,1000,5000)
image test2 := FilterLocalMaxima1D(test1,3)
test1.ShowImage()
test2.ShowImage()

但是，考虑到噪声（也在您的示例图像中），您可能希望围绕这些“局部最大值”执行拟合，以确保您真正得到您想要的。上面的数据只是这个的起点。

另外：有时您可以先对数据进行平均，然后找到局部最大值，而不是在每个峰值中进行真实数据拟合。如果您“非常清楚”实际峰的宽度，这尤其有效。

这将是一个例子：

image CreateTestSpec( number nChannels, number nPeaks, number amp, number back )
{
    image testImg := RealImage( "TestSpec", 4, nChannels )
    testImg = amp * cos( PI() +  2*PI() * nPeaks * icol/(iwidth-1) )
    testImg += back
    testImg = PoissonRandom( testImg )
    return testImg
}

image FilterLocalMaxima1D( image spectrumIn, number range )
{
    image spectrumOut := spectrumIn.ImageClone()
    for( number dx = -range; dx<=range; dx++  )
        spectrumout *= ( spectrumIn >= offset(spectrumIn,dx,0) ? 1 : 0 )

    spectrumout.SetName("Local maxima ("+range+") filtered")
    return spectrumOut
}

image FilterLocalMaxima1DAveraged( image spectrumIn, number range )
{
    image avSpectrum := spectrumIn.ImageClone()
    avSpectrum = 0
    for( number dx = -range; dx<=range; dx++  )
        avSpectrum += offset(spectrumIn,dx,0) 
    avSpectrum *= 1/(2*range+1)

    image spectrumOut := spectrumIn.ImageClone()
    for( number dx = -range; dx<=range; dx++  )
        spectrumout *= ( avSpectrum >= offset(avSpectrum,dx,0) ? 1 : 0 )

    spectrumout.SetName("Local maxima ("+range+") filtered, average")
    return spectrumOut
}

image test1 := CreateTestSpec(256,10,1000,5000)
image maxPeaks      := FilterLocalMaxima1D(test1,3)
image maxPeaksAv    := FilterLocalMaxima1DAveraged(test1,3)
test1.ShowImage()
test1.ImageGetImageDisplay(0).ImageDisplayAddImage( maxPeaks, "local max" )
test1.ImageGetImageDisplay(0).ImageDisplayAddImage( maxPeaksAv, "local max from Average" )

test1.ImageGetImageDisplay(0).LinePlotImageDisplaySetSliceComponentColor( 0, 1, 0.7, 0.7, 0.7 )

test1.ImageGetImageDisplay(0).LinePlotImageDisplaySetSliceDrawingStyle( 1, 2)
test1.ImageGetImageDisplay(0).LinePlotImageDisplaySetSliceComponentColor( 1, 1, 1, 0, 0 )
test1.ImageGetImageDisplay(0).LinePlotImageDisplaySetSliceTransparency( 1, 1, 0.7 )

test1.ImageGetImageDisplay(0).LinePlotImageDisplaySetSliceDrawingStyle( 2, 2)
test1.ImageGetImageDisplay(0).LinePlotImageDisplaySetSliceComponentColor( 2, 1, 0, 1, 0 )
test1.ImageGetImageDisplay(0).LinePlotImageDisplaySetSliceTransparency( 2, 1, 0.7 )

Answer 2

最简单的方法是使用 1 点（或 2 点）邻域来决定，中心是否最小（最大）。见下面的伪代码：

// assume 0 <= x <= maxX, y(x) is value at point x, radius is 1
x = 1;

while (x + 1 <= maxX)
{
    if (y(x) > y(x-1) and y(x) > y(x+1))
        // we found a maximum at x

    if (y(x) < y(x-1) and y(x) < y(x+1))
        // we found a minimum at x

    x = x + 1
}

对于 2 点邻域，最大条件可能是

if (y(x) > y(x-1) and y(x-1) >= y(x-2) and y(x) > y(x+1) and y(x+1) >= y(x+2))

注意 >= 这里。您可以使用 > 代替。

请注意，如果两个连续的 x 具有相同的值 y，则上述方法不会找到最大值，例如对于 y(0) = 0, y(1) = 1, y(2) = 1, y(3) = 0 它赢了'对于 x = 1 和 x = 2 都找不到最大值。