Performing Particle Analysis via ChooseMenuItem()

Question

_{This is a continuation of the question 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?

---

Here are some 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()
    }
}

If select "useBadImage" and "stopAtSrcImage", the first src image is shown. That would be a very bad one for particle analysis. The UI operation will generate multiple "invalid index" error (try select top 15% intensity in histogram). If useBadImage==0, then the image behaves better.

I extend the code to run some loop through an image stack. For an image stack with "good" images, the manual particle analysis can go through well with all the layers without problem, but the loop almost always generate "Exceptions" (shown in Result window) somewhere in the loop. Seems adding long delay does not help. But without delay definitely crashes the loop. The doThreshold(image histogram, number PctOffPeak, number AdditionalShift) was suppose to "Find the max in histogram, start from there, go right to some percentage off the peak, add the AdditionalShift, and set the "lf" value for ROI range there. But this does not always behave as expected.

Answer 1

While trying to deal with the hidden error, I worked out my own version of particle analysis. The speed is not a big issue. Here are some core functions,

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

A few comments to the script as posted:

---

stack:=exprSize(x1,y1,z1)

does not create a 3D image of size x1/y1/z1 but a 2D image of size x1/y1 with value z1 in each pixel. You possibly want

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

---

_FloatingModelessDialog( message, buttontext )

is a custom-script command not available for others. However, I've guess it should be a modeless dialog with two buttons, so I've used the following script code:

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

The only times I did run into the following thrown error

was when I had a treshold value which produced a binary mask which then, in the remove-edge-particle step, was completely removed because there was no particle not touching an edge. The "find particle" routine then rightfully throws the error, because no particle could be found (there was no longer a binary mask).

However, you can access the threashold-mask by scripting, so it might be a good idea to perform a check on that before calling find-particle. Then you avoid that error. (See other answer)

Answer 3

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.

For clarity: "Thresholding" creates the green "mask" on top of an image, which really is just a Boolean array for each pixel:

This (binary) mask is used by the particle analysis to find the particles. The routine generates multiple particle-annotations from the mask, which are shown in red (with yellow border and white numbering):

The first step - creating and modifying the green mask - can be completely crontrolled by scripting.

The second step - getting from the green mask to the red/yellow particles - can not and requires ChooseMenuItem().

The third step - computing values from the red/yellow particles mask - also can not be scripted directly and requires ChooseMenuItem() and some creative coding.

---

The DM help documentation (in recent versions at least) has an example of using thresholding and getting the mask by script in the section on the RasterImageDisplay Object

I simply copy the example script here:

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

The script above works for simple tresholding 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.

However, one can work around that issue by using the general commands for components to access the needed information. The following script creates the mask image correctly:

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

With the ability to get the thresholdmask, one can check if that mask is actually containing any pixels before calling the Particle Analysis command. This avoids the errors to be thrown.

---

The original script modified to avoid the problem:

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