Further example to translate to multi-thread?

Share icon

I appreciate the recent insight on using a thread safe function to avoid a loop. I'd like to ask whether the code below is amenable to being streamlined in the same way. I am using it to loop over all layers in a stack of grayscale images and remove a background.

Duplicate/FREE simg, fbimgSTACK
Redimension/D fbimgSTACK   
for (ic=0;ic<nimgs;ic+=1)
    MatrixOP/FREE fbimg = layer(fbimgSTACK,ic)
    ImageRemoveBackground/O/W/P=(po)/R=RoIMask fbimg
    WaveStats/Q/M=1 fbimg
        sf = 255/(V_max - V_min)
    MatrixOP/O fbimg = (sf*(fbimg - V_min))
    fbimgSTACK[][][ic] = fbimg[p][q]
endfor
MatrixOP/O wbimg = uint8(fbimgSTACK)

I have tried a few approaches. I cannot seem to wrap my head around how to pass the layer information into and back out of the ThreadSafe function.

 
threadsafe Function RemoveBackground(WAVE fbimgSTACK, variable ic)
    MatrixOP/FREE fbimg = layer(fbimgSTACK,ic)
    ImageRemoveBackground/O/W/P=(po)/R=RoIMask fbimg
    WaveStats/Q/M=1 fbimg
        sf = 255/(V_max - V_min)
    MatrixOP/O fbimg = (sf*(fbimg - V_min))
    fbimgSTACK[][][ic] = fbimg[p][q]
End

Duplicate/FREE simg, fbimgSTACK
Redimension/D fbimgSTACK

Make/N=(nimgs) junkWave
Multithread junkWave[] = RemoveBackground(fbimgSTACK, p)

MatrixOP/O wbimg = uint8(fbimgSTACK)

Something like the above could be done.But beware to only change the contents of fbimgSTACK in the threadsafe function and nothing else like its size.

 

April 16, 2020 at 07:49 am - Permalink

@Thomas -- Thanks. I have only started to learn about including the input designations in the function definition. I had not considered using a temporary wave.

Here is what works.

Make/N=(nimgs)/FREE tmpWave
Duplicate/FREE simg, fbimgSTACK
Redimension/D fbimgSTACK   
MultiThread tmpWave[] = TS_RemoveBackground(fbimgSTACK,ROIMask,po,p)
MatrixOP/O wbimg = uint8(fbimgSTACK)

ThreadSafe Function TS_RemoveBackground(WAVE fbimgSTACK, WAVE RoIMask, variable po, variable ic)

    variable sf
    MatrixOP/FREE fbimg = layer(fbimgSTACK,ic)
    ImageRemoveBackground/O/W/P=(po)/R=RoIMask fbimg
    WaveStats/Q/M=1 fbimg
    sf = 255/(V_max - V_min)
    MatrixOP/O fbimg = (sf*(fbimg - V_min))
    fbimgSTACK[][][ic] = fbimg[p][q]

    return 0

end

The above approach is faster by factor of 1.9x than the brute-force loop over an image stack with 15 images at 8-bit and 2448x2048.

I still getting a spinning wheel. I suspect, based on discussions with AG some time ago, that my RoIMask is over-indulgent (too many points for the plane fitting). That is yet another case to explore.

EDIT: Well, this worked once. But after that, it did not. I'll have to dig a bit further.

April 16, 2020 at 08:21 am - Permalink

Jeff,

You really want to minimize the number of pixels in the ROI mask.  Also, a minor improvement may come from using something along the lines of

Variable vmin=WaveMin(fbimg)
Variable vmax=WaveMax(fbimg)

which in IP9 you may be able to change to a single call.

Also, your layer assignment at the end is probably less efficient than it should be.  I'd try ImageTransform setPlane.

AG

April 16, 2020 at 09:23 am - Permalink

@jjweimer: Can you provide a MWE for playing around?

The last matrixop call in the function could also be ditched and integrated in the layer setting (the last line).

April 16, 2020 at 10:44 am - Permalink

@Thomas: I am incorporating changes in the Image Tools package (reducing the number of data points in the image background RoI mask for example). It will be easier for me to post the updated Image Tools package and then point to the location where the optimization would be useful. Even this effort may take a while between other obligations that I have.

@ AG ... I've not had luck with ImageTransform setPlane to replace a plane back into an image stack. I keep getting an error: Parameter is out of range.

This is the command that I am trying to replace.

wbimg[][][10]= fbimg[p][q]

If I understand ImageTransform setPlane, the /D flag is the DESTINATION stack and the /P flag sets the plane that will be REPLACED. So, this *should* be the equivalent:

ImageTransform/P=(10)/D=wbimg setPlane, fbimg

What is the correct approach here?

April 23, 2020 at 07:01 pm - Permalink