Выполнение анализа частиц с помощью ChooseMenuItem() - что не так

Это продолжение вопроса здесь.
Я публикую один из ответов, чтобы темы были чище.
Оригинальный вопрос (и продолжение вопроса) от пользователя user6406828.

Пытаясь выполнить некоторый (зацикленный) анализ частиц, время от времени возникает пара ошибок. Что можно улучшить в этом коде?

Вот некоторый код:

// $BACKGROUND$
number useBadImage, stopAtSrcImage // flags for test
image histogram, src, stack
taggroup imgLst
number i, imgCt,imgID,x0,y0,x1,y1, z1
if (OkCancelDialog("Generate a bad image?")) useBadImage=1
else useBadImage=0  
if (OkCancelDialog("Stop at testing image generation?")) stopAtSrcImage=1
else stopAtSrcImage=0  

x0=128;y0=128
x1=512;y1=512
image img:=exprSize(x0,y0,0)
z1=10
stack:=exprSize(x1,y1,z1)

img=random()*100
src:=exprSize(x1,y1,0)
if (useBadImage) {
    src=img.warp(icol*x0/x1,irow*y0/y1)
    for (i=0;i<z1;i++) {
        stack[0,0,i,x1,y1,i+1]=src+10
    }
}
else {
    src=sin(icol/pi())+cos(irow/pi())
    for (i=0;i<z1;i++) {
        stack[0,0,i,x1,y1,i+1]=src+0.1
    }
}
if (stopAtSrcImage) {
    src.showImage()
    exit(0)
}
void doThreshold(image histogram, number PctOffPeak, number AdditionalShift) {
    number max, lf,rt,maxX,maxY, i, val, threshold,d0,d1
    ROI r = NewROI(); // foreground ROI
    histogram.getSize(d0,d1)
    //histogram[0,0,1,5]=0 //remove zero peak
    max=histogram.max(maxX,maxY)
    threshold=max*PctOffPeak/100
    //okdialog("After trim zero peak, maxX=" +maxX)
    lf=0;rt=d0
    for (i=maxX;i<d0;i++) {
        val=histogram.getPixel(i,0)
        if (val<=threshold) {
            lf=i
            i=d0
        }
    }
    lf=lf+AdditionalShift
    r.ROISetRange(lf,rt);
    histogram.ImageGetImageDisplay(0).ImageDisplayAddRoi(r);
}
//main loop
for (i=0;i<z1;i++) {
    src:=stack[0,0,i,x1,y1,i+1].imageclone()
    src.showImage()
    ChooseMenuItem( "Analysis", "Particles", "Start Threshold" )
    histogram := GetFrontImage();
    doThreshold(histogram, 20,2)
    if( !_FloatingModelessDialog( "Adjust threshold level","Proceed!" ) ) {; histogram.deleteimage(); exit(0); };
    else {
        src.showImage()
        delay(60)
        ChooseMenuItem( "Analysis", "Particles", "Remove Edge Particles" )
        delay(60)
        ChooseMenuItem( "Analysis", "Particles", "Find Particles" )
        delay(60)
        ChooseMenuItem( "Analysis", "Particles", "Analyze Particles" )
        delay(60)
        histogram.deleteImage()
    }
}

Если выбрать "useBadImage" и "stopAtSrcImage", будет показано первое изображение src. Это было бы очень плохо для анализа частиц. Операция пользовательского интерфейса вызовет многократную ошибку "неверный индекс" (попробуйте выбрать верхние 15% интенсивности в гистограмме). Если useBadImage==0, то изображение ведет себя лучше.

Я расширяю код для запуска некоторого цикла через стек изображений. Для стека изображений с "хорошими" изображениями ручной анализ частиц может пройти без проблем все слои, но цикл почти всегда генерирует "Исключения" (показанные в окне результатов) где-то в цикле. Кажется, добавление большой задержки не помогает. Но без промедления однозначно вылетает петля. doThreshold(image histogram, number PctOffPeak, number AdditionalShift) Предполагалось, что "Найти максимум в гистограмме, начать оттуда, перейти вправо к некоторому проценту от пика, добавить AdditionalShift и установить там значение" lf "для диапазона ROI. Но это не всегда ведет себя так, как ожидалось.

Источник

3 ответа

Решение

Пытаясь справиться со скрытой ошибкой, я разработал собственную версию анализа частиц. Скорость не является большой проблемой. Вот некоторые основные функции,

image dxImg,dyImg
void getSearchImg() {
    dxImg:=[8,1]:
    {
        {-1, 0, 1, 1, 1, 0,-1,-1}
    }
    dyImg:=[8,1]:
    {
        {-1,-1,-1, 0, 1, 1, 1, 0}
    }
}
number pop(taggroup &tg, number &x, number &y) {
    number n, mode
    n=tg.TagGroupCountTags()
    if (n==0) return -1
    else { //stack, last in, first out
        tg.TagGroupGetIndexedTagAsLongPoint(n-1,x,y)
        tg.TagGroupDeleteTagWithIndex(n-1)
    }
    return 1
}
void push(taggroup &tg, number x, number y) {
    tg.TagGroupInsertTagAsLongPoint(infinity(),x,y)  
}
taggroup dfs(image img) {//Depth-first search, Using stack
    number isOn, i,xInc, yInc,d0,d1
    number status, x0, y0,x1,y1, x, y, val
    taggroup tgStack=newtaglist()
    taggroup tgClusters=newtaglist()
    taggroup tgSingleCluster=newTaglist() //this is global
    roi r
    image imgTemp:=img.ImageClone()
    imgTemp.getSize(d0,d1)
    getSearchImg()
    imgTemp=tert(icol==0||irow==0||icol==d0-1||irow==d1-1, 0, imgTemp)
    while (1) {
        status=1
        val = imgTemp.max(x, y)
        if (val<=0) break
        else {
            tgSingleCluster.TagGroupInsertTagAsLongPoint(infinity(),x,y)
            tgStack.push(x,y)
            imgTemp[x,y]=0
            x0=x;y0=y
            while (status==1) {
                for (i=0;i<8;i++) {
                    xInc=dxImg.getPixel(i,0)
                    yInc=dyImg.getPixel(i,0)
                    x=x0+xInc;y=y0+yInc
                    isOn=imgTemp.getPixel(x,y)
                    if (isOn) {
                        imgTemp[x,y]=0
                        tgSingleCluster.TagGroupInsertTagAsLongPoint(infinity(),x,y)
                        tgStack.push(x,y)
                    }
                }
                status=tgStack.pop(x0,y0)
            }
            tgClusters.TagGroupInsertTagAsTagGroup(infinity(), tgSingleCluster) 
            tgSingleCluster=newTaglist()
            tgStack=newtaglist()
        }
    }
    return tgClusters
}
Источник

В то время как DM scirpting не имеет доступа к самому "анализу частиц" - и, следовательно, ChooseMenutItem() Нужно использовать трюк - в нем есть команды сценариев для установки и получения пороговой маски.

Для ясности: "Thresholding" создает зеленую "маску" поверх изображения, которая на самом деле является просто логическим массивом для каждого пикселя:

Изображение с пороговой маской

Эта (двоичная) маска используется анализом частиц, чтобы найти частицы. Процедура генерирует несколько аннотаций частиц из маски, которые показаны красным (с желтой рамкой и белой нумерацией):

Изображение с маской частиц

Первый шаг - создание и изменение зеленой маски - может полностью контролироваться сценариями.

Второй шаг - добраться от зеленой маски до красных / желтых частиц - не может и требует ChooseMenuItem(),

Третий шаг - вычисление значений по красной / желтой маске частиц - также не может быть записан непосредственно в сценарии и требует ChooseMenuItem() и немного творческого кодирования.

Справочная документация DM (по крайней мере, в последних версиях) содержит пример использования порогового значения и получения маски по сценарию в разделе объекта RasterImageDisplay.

Помощь DM

Я просто копирую пример сценария здесь:

// open image
image myImage := GetFrontImage()
ImageDisplay imageDisp = myImage.ImageGetImageDisplay( 0 )

number low, high
imageDisp.ImageDisplayGetContrastLimits( low, high )

number width = myImage.ImageGetDimensionSize( 0 )
number height = myImage.ImageGetDimensionSize( 1 )

// trun thresholding on
imageDisp.RasterImageDisplaySetThresholdOn( 1 ) 

// set limits
imageDisp.RasterImageDisplaySetThresholdLimits( low, (low + high)/2 ) 

// create mask image, should be binary or 8 bit signed or unsigned
image mask := IntegerImage("Mask", 1, 0, width, height )
imageDisp.RasterImageDisplayAddThresholdToMask( mask, 0, 0, height, width ) 

// turn thresholding off
imageDisp.RasterImageDisplaySetThresholdOn( 0 ) 

// display the mask
ShowImage( mask )

Сценарий выше работает для простой маски трешолдинга.

Тем не менее, я обнаружил, что есть странное поведение, если
ChooseMenuItem( "Analysis", "Particles", "Remove Edge Particles" ) вызывается и удаляет части маски. Тогда человек больше не сможет получить правильную маску, используя RasterImageDisplayAddThresholdToMask() команда.

Однако можно обойти эту проблему, используя общие команды для компонентов для доступа к необходимой информации. Следующий скрипт правильно создает изображение маски:

image GetThresholdMaskFromDisplay( imageDisplay disp )
{
    number kRasterMaskType = 3633   // DM defined constant
    if ( !disp.ImageDisplayIsValid() ) Throw( "Invalid display")
    number sx = disp.ImageDisplayGetImage().ImageGetDimensionSize(0)
    number sy = disp.ImageDisplayGetImage().ImageGetDimensionSize(1)
    image mask := IntegerImage( "Mask", 1, 0, sx, sy )

    component rasterMask = disp.ComponentGetNthChildOfType( kRasterMaskType, 0 )
    if ( rasterMask.ComponentIsValid() )
    {
        TagGroup compTgs = NewTagGroup()
        rasterMask.ComponentExternalizeProperties( compTgs )
        compTgs.TagGroupGetTagAsArray( "MaskData", mask )
    }
    return mask 
}

GetFrontImage().ImageGetImageDisplay(0).GetThresholdMaskFromDisplay().ShowImage()

Имея возможность получить пороговую маску, можно проверить, содержит ли эта маска какие-либо пиксели, прежде чем вызывать команду анализа частиц. Это позволяет избежать ошибок, которые будут выброшены.

Исходный скрипт изменен, чтобы избежать проблемы:

// $BACKGROUND$
number _fModelessDialog(string prompt, string buttonName){

    number sem = NewSemaphore()
    ModelessDialog(prompt, buttonName, sem)
    try GrabSemaphore(sem)
    catch return 0
    return 1
}

image GetThresholdMaskFromDisplay( imageDisplay disp )
{
    number kRasterMaskType = 3633   // DM defined constant
    if ( !disp.ImageDisplayIsValid() ) Throw( "Invalid display")
    number sx = disp.ImageDisplayGetImage().ImageGetDimensionSize(0)
    number sy = disp.ImageDisplayGetImage().ImageGetDimensionSize(1)
    image mask := IntegerImage( "Mask", 1, 0, sx, sy )

    component rasterMask = disp.ComponentGetNthChildOfType( kRasterMaskType, 0 )
    if ( rasterMask.ComponentIsValid() )
    {
        TagGroup compTgs = NewTagGroup()
        rasterMask.ComponentExternalizeProperties( compTgs )
        compTgs.TagGroupGetTagAsArray( "MaskData", mask )
    }
    return mask 
}

number useBadImage, stopAtSrcImage // flags for test
image histogram, src, stack
taggroup imgLst
number i, imgCt,imgID,x0,y0,x1,y1, z1
if (OkCancelDialog("Generate a bad image?")) useBadImage=1
else useBadImage=0  
if (OkCancelDialog("Stop at testing image generation?")) stopAtSrcImage=1
else stopAtSrcImage=0  

x0=128;y0=128
x1=512;y1=512
image img:=exprSize(x0,y0,0)
z1=10
stack:=exprSize(x1,y1,z1,0)

img=random()*100
src:=exprSize(x1,y1,0)
if (useBadImage) {
    src=img.warp(icol*x0/x1,irow*y0/y1)
    for (i=0;i<z1;i++) {
        stack[0,0,i,x1,y1,i+1]=src+10
    }
}
else {
    src=sin(icol/pi())+cos(irow/pi())
    for (i=0;i<z1;i++) {
        stack[0,0,i,x1,y1,i+1]=src+0.1
    }
}
if (stopAtSrcImage) {
    src.showImage()
    exit(0)
}
void doThreshold(image histogram, number PctOffPeak, number AdditionalShift) {
    number max, lf,rt,maxX,maxY, i, val, threshold,d0,d1
    ROI r = NewROI(); // foreground ROI
    histogram.getSize(d0,d1)
    //histogram[0,0,1,5]=0 //remove zero peak
    max=histogram.max(maxX,maxY)
    threshold=max*PctOffPeak/100
    //okdialog("After trim zero peak, maxX=" +maxX)
    lf=0;rt=d0
    for (i=maxX;i<d0;i++) {
        val=histogram.getPixel(i,0)
        if (val<=threshold) {
            lf=i
            i=d0
        }
    }
    lf=lf+AdditionalShift
    r.ROISetRange(lf,rt);
    histogram.ImageGetImageDisplay(0).ImageDisplayAddRoi(r);
}
//main loop
for (i=0;i<z1;i++) {
    src:=stack[0,0,i,x1,y1,i+1].imageclone()
    src.showImage()
    ChooseMenuItem( "Analysis", "Particles", "Start Threshold" )
    histogram := GetFrontImage();
    doThreshold(histogram, 20,2)
    if( !_fModelessDialog( "Adjust threshold level", "Proceed" ) ) {; histogram.deleteimage(); exit(0); };
    else {
        src.showImage()
        delay(60)
        ChooseMenuItem( "Analysis", "Particles", "Remove Edge Particles" )
        delay(60)
        // Check if there even is a particle mask now!
        number validPixelsInMask = sum( src.ImageGetImageDisplay(0).GetThresholdMaskFromDisplay() )

        if ( !validPixelsInMask )
            Result("\n Skipping plane, no particles found." )
        else
        {
            ChooseMenuItem( "Analysis", "Particles", "Find Particles" )
            delay(60)
            ChooseMenuItem( "Analysis", "Particles", "Analyze Particles" )
            delay(60)
            histogram.deleteImage()
        }
    }
}
Источник

Несколько комментариев к сценарию как опубликовано:

stack:=exprSize(x1,y1,z1)

создает не трехмерное изображение размером x1/y1/z1, а двумерное изображение размером x1 / y1 со значением z1 в каждом пикселе. Вы возможно хотите

stack:=exprSize(x1,y1,z1,0)

_FloatingModelessDialog( message, buttontext )

команда нестандартного сценария, недоступная для других Тем не менее, я предполагаю, что это должен быть немодальный диалог с двумя кнопками, поэтому я использовал следующий код скрипта:

number _fModelessDialog(string prompt, string buttonName){

    number sem = NewSemaphore()
    ModelessDialog(prompt, buttonName, sem)
    try GrabSemaphore(sem)
    catch return 0
    return 1
}

Я на самом деле не столкнулся с проблемой запуска скрипта после этих изменений, используя "плохой" образ

Единственный раз, когда я столкнулся со следующей ошибкой

Это было, когда у меня было пороговое значение, которое создавало двоичную маску, которая затем на этапе удаления частицы была полностью удалена, потому что не было частицы, не касающейся края. Затем процедура "найти частицу" по праву выдает ошибку, потому что никакая частица не может быть найдена (больше не было двоичной маски).

Тем не менее, вы можете получить доступ к маске treashold при помощи сценариев, поэтому было бы неплохо проверить это перед вызовом find-частиц. Тогда вы избежите этой ошибки. (См. Другой ответ)

Источник
Другие вопросы по тегам