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

Question

_{这是问题here 的延续。
我重新发布one of the answers 以保持主题更清晰。
原始问题（和后续问题）来自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 图像。对于粒子分析来说，这将是一个非常糟糕的问题。 UI 操作会产生多个“invalid index”错误（尝试在直方图中选择前 15% 的强度）。如果 useBadImage==0，则图像表现更好。

我扩展代码以通过图像堆栈运行一些循环。对于具有“好”图像的图像堆栈，手动粒子分析可以顺利通过所有层而没有问题，但循环几乎总是在循环中的某处生成“异常”（显示在结果窗口中）。似乎增加长时间的延迟没有帮助。但是毫不拖延地肯定会使循环崩溃。 doThreshold(image histogram, number PctOffPeak, number AdditionalShift) 被假设为“在直方图中找到最大值，从那里开始，向右走峰值的某个百分比，添加 AdditionalShift，并在那里设置 ROI 范围的“lf”值。但这并不总是表现为预计。

Answer 1

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

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
}

Answer 2

虽然 DM 脚本无法访问“粒子分析”本身 - 因此需要使用 ChooseMenutItem() 技巧 - 它确实具有用于阈值化和获取阈值掩码的脚本命令。

为了清楚起见：“阈值”在图像顶部创建绿色“蒙版”，这实际上只是每个像素的布尔数组：

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

第一步 - 创建和修改绿色掩码 - 可以通过脚本完全控制。

第二步——从绿色蒙版到红色/黄色粒子——不能，需要ChooseMenuItem()。

第三步 - 从红色/黄色粒子蒙版计算值 - 也不能直接编写脚本，需要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 )

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

但是，我发现如果调用
ChooseMenuItem( "Analysis", "Particles", "Remove Edge Particles" ) 并删除部分掩码，会出现一些奇怪的行为。然后，使用RasterImageDisplayAddThresholdToMask() 命令将无法获得正确的掩码。

但是，可以通过使用 components 的通用命令访问所需信息来解决该问题。以下脚本正确创建蒙版图像：

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

Answer 3

发布的脚本中的一些 cmets：

---

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
}

---

在使用“坏”图像进行这些更改后，我实际上并没有遇到运行脚本的问题

我确实唯一一次遇到了以下抛出的错误

当我有一个阈值时，它产生了一个二进制蒙版，然后在 remove-edge-particle 步骤中，它被完全删除，因为没有粒子不接触边缘。 然后“查找粒子”例程正确抛出错误，因为找不到粒子（不再有二进制掩码）。

但是，您可以通过脚本访问阈值掩码，因此在调用 find-particle 之前 执行检查可能是个好主意。然后你就避免了这个错误。 （见其他答案）