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

Question

Я хотел бы спросить, можно ли определить положение всех максимумов и минимумов профиля интенсивности на DM.

Как мне создать простой скрипт, который идентифицирует положения пиков в примере ниже?

Вот снимок экрана интенсивности линии профиля стволовой изображения вдоль Y-направлении:

Answer 1

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

Answer 2

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