Rescaling image contrast: excluding outlier pixel values

Hello,

I am working with some images and need to rescale the contrast on those image for visualisation. If I scale to the min-max for the image, it still appears quite dark because I have a handful of extremely bright pixels in the image. In ImageJ the way around this is a button called Auto which seems to trim the top (and bottom) 0.35% of pixels according to the shape of the image histogram. I have coded something similar in Igor* and it works, but it is slow. I'm wondering if there is a faster way to get these values. Is there a function to specify the 99.65th (or Xth) percentile of an image and use that for the scaling? Maybe a histogram function would work? Any help or other ideas welcome.

* I calculate the Min and Max like this and store them in a wave to then do the rescaling in the calling function. It's slow because I duplicate-redimension-sort for each channel and for each z-position.

STATIC Function MakeMinMaxWave(ImageMat,rCh,gCh,bCh)
    Wave ImageMat
    Variable rCh,gCh,bCh // index of chunk for each channel. -1 means blank
    String rgbList = num2str(rCh) + ";" + num2str(gCh) + ";" + num2str(bCh) + ";"
    Variable nZ = DimSize(ImageMat,2)
    Make/O/N=(nZ,2,3) minMaxW // rows for z, columns min and max, layers r g b
    Variable totalOrigPix = DimSize(ImageMat,0) * DimSize(ImageMat,1)
    Variable chSelect
   
    Variable i,j
   
    for(i = 0; i < nZ; i += 1)
        for(j = 0; j < 3; j += 1)
            chSelect = str2num(StringFromList(j,rgbList))
            if(chSelect == -1)
                continue
            endif
            Duplicate/O/FREE/RMD=[][][i][chSelect] ImageMat, tempW
            Redimension/N=(totalOrigPix) tempW
            Sort tempW,tempW
            minMaxW[i][0][j] = tempW[0]
            minMaxW[i][1][j] = tempW[floor(totalOrigPix-(totalOrigPix * 0.003))] // ImageJ does something like 0.35%
        endfor
    endfor
End

 

Are you mainly interested in thresholding for display within Igor? If so, it's not necessarily useful to change the values in the original image. As far as thresholding for display, Igor seems to have better transformations for display of single color images than it does three color images. For example, "ctab" and "cindex", which can respectively implement thresholding and non-linear transforms like gamma functions, don't appear to be available for RGB images.

A sort of work around that I have used is to split the RGB layers into 3 separate 2D matrices and then overlay them (with newImage then appendImage). This approach only became viable with the implementation transparency in Igor 7. One has to append them and then assign an appropriate 4-column custom color table to each, with all alpha values set to 1/3rd of maximum. 

This allows the following much faster thresholding for each image. For example, here is how you would threshold the red channel: 

ModifyImage originalImage_red ctab= {lowerThreshold,upperThreshold,myRedsWithTransparency,0}

With single color images, there are some helpful built in functions like the "Image Range Adjustment" dialog in the Image menu, which modifies the ctab thresholds.

 

If you really do want to modify the original data values, you can do the following (among other possibilities). Say you wanted to quickly apply thresholds set with the "Image Range Adjustment" Dialog. You can click "Save Range" in this dialog, save as "MIP_Range0" and then run:

MatrixOp thresholdedImage = clip(originalImage,MIP_Range0[0],MIP_range0[1])
//then transfer any dimension scaling from originalImage to thresholdedImage

 

In case you still want to determine the lower and upper thresholds as % pixel values: once you have each color component in separate 2D waves you could run the following (though there may be a less involved approach I am not aware of):

Variable lowerPercent=0.05 //say threshold pixels in bottom 5% of range
Variable upperPercent=0.9 //and above 9% of range

imagehistogram/i originalImageComponent
integrate w_imagehist/d=cumHist
Variable minLevel=cumHist[0]
Variable maxLevel=cumHist[dimsize(cumHist,0)-1]
Variable range=maxLevel-minLevel
FindLevel cumHist, minLevel + range*lowerPercent
Variable lowerThreshold=V_levelX
FindLevel cumHist,minLevel+range*upperPercent
Variable upperThreshold=V_levelX
//use lowerThreshold and upperThreshold in matrixop clip or modifyImage ctab, for example

You may also want to look at the "imageThreshold" function

Thanks aoa this is really helpful! The overlay approach is nice since I only want to change the contrast of the images as displayed. Possibly breaking each image into channels as 2D waves might not be needed. It looks like it is possible to display/append each layer using /G=1 and specify the trace instance to be a specific layer, but I need to investigate further. 

Hi sjr51,

 

Could you perhaps supply a sample image.

 

Best,

_sk

@_sk I am attaching a zipped tif file. This is one channel and one z-layer that I'm dealing with. This channel is the most difficult to scale. 

zipped tif file

sjr51,

I tried something like this on your image.

 

function fn_autoscale(s_graph)
    string s_graph // the name of the graph, i.e. graph0
   
    string s_list = imagenamelist(s_graph, ";") // list of images appended to graph
    variable i = itemsinlist(s_list)-1
   
    string s_img
    do
        s_img = stringfromlist(i, s_list)
        wave w_img = imagenametowaveref(s_graph, s_img)

        imagethreshold/q/m=2 w_img //produces v_threshold

        modifyimage/w=$s_graph $s_img ctab= {0,v_threshold,Grays,0}
       
        i -= 1
    while ( i != -1)
end

And idk if it works, but I think it does more or less what you wanted, unless I have misunderstood (which is quite possible).

 

Best,

_sk

 

edit: I am including some before and after images

before after

@_sk

Thank you for your help. I think this combined with the overlay approach suggested by aoa should work.

In reply to by sjr51

sjr51 wrote:

@_sk

Thank you for your help. I think this combined with the overlay approach suggested by aoa should work.

 

Yes, it seems like it will do what you wanted. Also note, that if you attach different images (i.e. different channels) to the same graph the above code will run through them, namely, the do.. while loop.

 

I thought I would leave another comment here to help anyone searching the forum for this topic.

aoa's code for thresholding worked great and was much faster than my code. I also looked at the way Igor does Image Range adjustment. I didn't know about this function when I wrote my question. It's similar to aoa's suggestion, the function is WMCont94ButtonProc in case anyone is interested.

While I liked the idea of using three planes and changing the display range, I couldn't get this to work. Three layers of one-third opacity did not sum correctly for me (see code below for an example). This left a greyish colour where there should have been black. Maybe this is because of the background of the graph window? In the end I was happy to make a copy of the images and adjust their values directly. So I didn't pursue further.

The second function below shows how to make the RGBA waves for display as mentioned above.

Function ThreeLayerDemo()
    Make/O/N=(100,100) aa=0
    MakeRGBAWavesForImageDisplay()
    // Make an overlay of three channels with on-third alpha
    NewImage/S=1 aa
    ModifyImage aa ctab= {0,65535,myRedW,0}
    AppendImage/T aa
    ModifyImage aa#1 ctab= {0,65535,myGreenW,0}
    AppendImage/T aa
    ModifyImage aa#2 ctab= {0,65535,myBlueW,0}
    // This is a pure black image
    NewImage/S=1 aa
    ModifyImage aa ctab= {0,65535,Grays,0}
End

Function MakeRGBAWavesForImageDisplay()
    Make/O/N=(256,4) myRedW=0,myGreenW=0,myBlueW=0
    myRedW[][3] = floor(65535 * 1/3)
    myGreenW[][3] = floor(65535 * 1/3)
    myBlueW[][3] = floor(65535 * 1/3)
    myRedW[][0] = p * 257
    myGreenW[][1] = p * 257
    myBlueW[][2] = p * 257
End

 

In reply to by sjr51

Looking at the three image overlay idea again, I agree that it's imperfect. I edited your code, sjr, to see what different color combinations look like. I think the main issue is that in this design is that each component is at most 1/3rd its maximum value.

I think it would be really nice if Igor would allow control of each color component with a dedicated wave. I don't think that's currently possible. 

 

 

Function ThreeLayerDemo()
    Make/O/N=(100,100) reds,greens,blues
    variable maxv=65535
    reds=0;greens=0;blues=0
   
    MakeRGBAWavesForImageDisplay()
    // Make an overlay of three channels with on-third alpha
    NewImage/S=1/k=1 reds
    ModifyImage reds ctab= {0,65535,myRedW,0}
    AppendImage/T greens
    ModifyImage greens ctab= {0,65535,myGreenW,0}
    AppendImage/T blues
    ModifyImage blues ctab= {0,65535,myBlueW,0}
    doupdate; sleep/s 1
   
    reds=maxv;greens=0;blues=0;doupdate;  sleep/s 1
    reds=maxv;greens=0;blues=0;doupdate;  sleep/s 1
    reds=0;greens=maxv;blues=0;doupdate;  sleep/s 1
    reds=0;greens=0;blues=maxv;doupdate;  sleep/s 1
    reds=maxv;greens=maxv;blues=0; doupdate;  sleep/s 1
    reds=0;greens=maxv;blues=maxv; doupdate; sleep/s 1
    reds=maxv;greens=0;blues=maxv; doupdate ;sleep/s 1
    reds=maxv;greens=maxv;blues=maxv ; doupdate    
End
Function MakeRGBAWavesForImageDisplay()
    Make/O/N=(256,4) myRedW=0,myGreenW=0,myBlueW=0
    myRedW[][3] = floor(65535 * 1/3)
    myGreenW[][3] = floor(65535 * 1/3)
    myBlueW[][3] = floor(65535 * 1/3)
    myRedW[][0] = p * 257
    myGreenW[][1] = p * 257
    myBlueW[][2] = p * 257
End