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

Если вы хотите отфильтровать «строгие» локальные максимумы, вы можете легко сделать это с помощью выражений изображения и условного оператора «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 )

Самый простой способ - использовать 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.