通过 ChooseMenuItem() 执行粒子分析 - 有什么问题

Performing Particle Analysis via ChooseMenuItem() - what is wrong

这是问题 here.
I am re-posting one of the answers to keep the topics cleaner.
The original question (and continued question) is from user6406828 的延续。


Trying to perform some (looped over) particle analysis, a couple of errors are thrown every now and then. What can be improved in this code?


这里有一些代码:

// $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()
    }
}

如果select "useBadImage"和"stopAtSrcImage",则显示第一个src图像。这对于粒子分析来说是非常糟糕的。 UI 操作将产生多个 "invalid index" 错误(尝试 select 直方图中前 15% 的强度)。如果 useBadImage==0,则图像表现更好。

我将代码扩展到 运行 通过图像堆栈的一些循环。对于具有 "good" 个图像的图像堆栈,手动粒子分析可以很好地完成所有层而没有问题,但循环几乎总是在某处生成 "Exceptions"(显示在结果 window 中)循环。似乎添加长延迟没有帮助。但毫不拖延肯定会使循环崩溃。 doThreshold(image histogram, number PctOffPeak, number AdditionalShift) 假设为那里的 ROI 范围的 "Find the max in histogram, start from there, go right to some percentage off the peak, add the AdditionalShift, and set the "lf" 值。但这并不总是按预期运行。

A few comments to the script as posted:


stack:=exprSize(x1,y1,z1)

不是创建大小为 x1/y1/z1 的 3D 图像,而是创建大小为 x1/y1 且每个像素的值为 z1 的 2D 图像。你可能想要

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
}

I actually did not run into an issue running the script after these changes using the "bad" image

唯一一次我做了运行下面抛出的错误

是当我有一个阈值时,它产生了一个二进制掩码,然后在移除边缘粒子步骤中,它被完全移除,因为没有粒子不接触边缘。 "find particle" 例程然后 正确地 抛出错误,因为找不到粒子(不再有二进制掩码)。

但是,您可以通过脚本访问阈值掩码,因此最好在调用 find-particle 之前检查一下。然后你就避免了那个错误。 (见其他答案)

While DM scirpting does not have access to the "particle analysis" itself - and thus the ChooseMenutItem() trick needs to be used - it does have script commands for threasholding and getting the threshold mask.

为清楚起见:"Thresholding" 在图像顶部创建绿色 "mask",这实际上只是每个像素的布尔数组:

粒子分析使用此(二进制)掩码来查找粒子。该例程从掩码生成多个粒子注释,以红色显示(带有黄色边框和白色编号):

第一步——创建和修改绿色遮罩——完全可以通过脚本来控制。

第二步 - 从绿色蒙版到 red/yellow 粒子 - 不能并且需要 ChooseMenuItem()

第三步 - 从 red/yellow 粒子蒙版计算值 - 也不能直接编写脚本,需要 ChooseMenuItem() 和一些创意编码。


DM 帮助文档(至少在最新版本中)在 RasterImageDisplay 对象

部分中有一个使用阈值和通过脚本获取掩码的示例

我只是简单地把示例脚本复制到这里:

// 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 )

上面的脚本适用于简单的阈值遮罩。

However, I've found that there is some strange behavior if
ChooseMenuItem( "Analysis", "Particles", "Remove Edge Particles" ) is called and removes parts of the mask. Then, one no longer can get the correct mask using the RasterImageDisplayAddThresholdToMask() command.

但是,可以通过使用 组件 的通用命令访问所需信息来解决该问题。以下脚本正确创建掩码图像:

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()
        }
    }
}

在尝试处理隐藏的错误时,我制定了自己的粒子分析版本。速度不是大问题。下面是一些核心功能,

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
}