Идентификация пиков из профиля линии

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
}

Для двухточечной окрестности максимальное условие может быть

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, которое он выиграл не найти максимум ни для х = 1, ни для х = 2.