Erosion Script

davek604

Hello

I'm using the code below on a distance map image of a cell. The idea is that the script makes 5 ~equal areas from the distance map. I do this by calculating the max pixel value which should be the middle of the distance map and then gradually decrease the threshold to select more and more of the image until a defined area is measured. I then save the selection to the features manager for use on other images. The script works but it is painfully slow (around 10mins per run) Are there any other commands to achieve this more quickly, I'm guessing the IQL commands would do it. All help/suggestions welcome.

For z = 1 To 5
            Stopper = Stopper + Incrementer

            'Delete any measurements already present
            With Measure.MeasurementsCommands.DeleteAll(Nothing)
                .Run(ThisApplication.ActiveImage)
            End With

            Do
                With Application.DocumentCommands.ActiveImage(Nothing)
                    .Run(ThisApplication.ActiveImage)
                End With

                With Measure.MeasurementsCommands.DeleteAll(Nothing)
                    .Run(ThisApplication.ActiveImage)
                End With

                With Measure.MeasurementsCommands.Options(Nothing)
                    .Segmentation.ApplyFilters = False
                    .Run(ThisApplication.ActiveImage)
                End With

                '** Sets Count Size to Manual   **
                With Measure.MeasurementsCommands.Options(Nothing)
                    .Segmentation.AutoFindPhase = MediaCy.IQL.Features.mcFindPhase.mcfpManual
                    .Segmentation.AutoThresholding = eAutoThresholdMethod.MinimumVariance
                    .Segmentation.SegmentationType = McMMOptions.mcmmSegmentationType.mcmmstThresholdSegmentation
                    .Run(ThisApplication.ActiveImage)
                End With

                With Measure.ThresholdToolCommands.Thresholds(Nothing)
                    .AllowOverlap = False
                    .Interpretation = eInterpretation.Mono
                    .Classes = New System.Collections.Generic.List(Of SegmentationClass)
                    .Classes.Add(New SegmentationClass("Class 1",System.Drawing.Color.Yellow,New Double(){LowVal,HighVal}))
                    .Run(ThisApplication.ActiveImage)
                End With

                 '** Section counts with whatever count parameters are set **
                With Measure.MeasurementsCommands.ExecuteCount(Nothing)
                    .Run(ThisApplication.ActiveImage)
                End With

                 '** Section to get Max value from Distance map **
                With measure.MeasurementsCommands.GetData(Nothing)
                    .Run(ThisApplication.ActiveImage, data)
                End With

                With ThisApplication.Output
                    For Each sf As McMMSubFeature In data.SubFeatures
                        MaxArea = sf.Value(eMeasures.RgnArea)
                    Next
                End With

                LowVal = LowVal -2

            Loop Until MaxArea >= Stopper

YuriG

Hello David,

You can  use just one Count operation to speed it up. Add multiple ranges to thresholds and use it with Count (note that the ranges are overlapping), it will generate objects with increasing areas. Here is a sample macro:

  Public Sub MeasureDistMap
        Dim HighVal As Double=500'some high value
        Dim LowVal As Double=200'first threshold
        Dim ringStep As Double=5
        Dim MaxArea As Double

        With Measure.ThresholdToolCommands.Thresholds(Nothing)
            .AllowOverlap = True'specifically ALLOW overlap!!!
            .Interpretation = eInterpretation.Mono
            .Classes = New System.Collections.Generic.List(Of SegmentationClass)
            'create multiple segmentation ranges
            For z As Integer = 1 To 5
                .Classes.Add(New SegmentationClass("Class " & z,System.Drawing.Color.Yellow,New Double(){LowVal-ringStep*z,HighVal}))
            Next z
            .Run(ThisApplication.ActiveImage)
        End With
        '** Sets Count Size to Manual   **
        With Measure.MeasurementsCommands.Options(Nothing)
            .Segmentation.AutoFindPhase = MediaCy.IQL.Features.mcFindPhase.mcfpManual
            .Segmentation.SegmentationType = McMMOptions.mcmmSegmentationType.mcmmstThresholdSegmentation
            .Run(ThisApplication.ActiveImage)
        End With
        'execute count
        With Measure.MeasurementsCommands.ExecuteCount(Nothing)
            .Run(ThisApplication.ActiveImage)
        End With
        'print areas
        With ThisApplication.Output
            .Show
            Dim data As McMMData=ThisApplication.ActiveImage.MeasurementsData
            For Each sf As McMMSubFeature In data.SubFeatures
                MaxArea = sf.Value(eMeasures.RgnArea)
                .PrintMessage (MaxArea.ToString)
            Next
        End With
End Sub

Regards,

Yuri  

davek604

Hello Yuri

Thanks for this, it works but not in the way the user wants. I did it a similar way to this but it gives areas which start high and decrease as you get closer to the centre as this method gives rings which are equal diameter but the area is very different. The user wants the areas the same so you end up with a very thin outer ring and a quite thick innermost ring. The only way I could think of doing it was dividing the area of the whole cell by 5 and then gradually increasing the threshold (measuring area with each iteration) until the area is correct. Its this gradual iteration that takes an age.

Regards

David

YuriG

Hello David,

I missed that you want to have equal ring areas. 
If you want to have rings of equal area the best way would be to set threshold ranges after analyzing Histogram of the distance map, turning it to accumulated (to report area) and then calculate ranges based on that. It can still be done in one Count operation.

Here is the complete macro that will do that (note that you have to add reference to MediaCy.IQL.Operations.dll!):

  Public Sub SplitDistanceMapToEqualAreaRings
        Dim NumberOfRings As Integer=5'define number of rings
        Dim Ranges() As Double
        ReDim Ranges(NumberOfRings)
        Dim MaxArea As Double
        Dim i As Integer

        ThisApplication.Output.Show
        ThisApplication.Output.Clear

        Dim hist As MediaCy.IQL.Operations.McHistogram=ThisApplication.ActiveImage.Histogram
        'calculate accumulated area starting from 1
        Dim TotalArea As Double=0
        For i=1 To hist.BinCount-1
            TotalArea+=hist.Values(i)
        Next
        Dim targetRingArea As Double=TotalArea/NumberOfRings
        ThisApplication.Output.PrintMessage String.Format("Total distance map area = {0} Ring area = {1}",TotalArea,targetRingArea)
        'calculate threshold ranges for rings of (approximatelly) the same area
        Ranges(0)=1'start from 1
        Dim sum As Double=0
        Dim bin As Integer=1
        For i=1 To hist.BinCount-1
            sum+=hist.Values(i)
            If sum>=targetRingArea*bin Then
                Ranges(bin)=hist.BinXValues(i)'set threshold
                ThisApplication.Output.PrintMessage String.Format("Threshold {0} = {1}",bin,Ranges(bin))
                bin+=1'increment bin
            End If
        Next

        With Measure.ThresholdToolCommands.Thresholds(Nothing)
            .AllowOverlap = False
            .Interpretation = eInterpretation.Mono
            .Classes = New System.Collections.Generic.List(Of SegmentationClass)
            'create multiple segmentation ranges, RINGS based on the threshold ranges calculated from histogram
            For z As Integer = 1 To NumberOfRings
                .Classes.Add(New SegmentationClass("Class " & z,System.Drawing.Color.Yellow,New Double(){Ranges(z-1),Ranges(z)}))
            Next z
            .Run(ThisApplication.ActiveImage)
        End With
        '** Sets Count Size to Manual   **
        With Measure.MeasurementsCommands.Options(Nothing)
            .Segmentation.AutoFindPhase = MediaCy.IQL.Features.mcFindPhase.mcfpManual
            .Segmentation.SegmentationType = McMMOptions.mcmmSegmentationType.mcmmstThresholdSegmentation
            .Run(ThisApplication.ActiveImage)
        End With
        'execute count
        With Measure.MeasurementsCommands.ExecuteCount(Nothing)
            .Run(ThisApplication.ActiveImage)
        End With
        'print areas
        With ThisApplication.Output
            Dim data As McMMData=ThisApplication.ActiveImage.MeasurementsData
            For Each sf As McMMSubFeature In data.SubFeatures
                .PrintMessage String.Format("Ring {0} area = {1}",sf.Name, sf.Value(eMeasures.RgnArea))
            Next
        End With
End Sub

I've also attached the screenshot that shows rings with different thickness but the same area.

Regards,

Yuri

davek604

Hello Yuri

Still does the same thing on my computer. The areas are very different. I also need access to the regions to use them on other images. I put them in the features manager. I think I'll send you the whole thing off list so you can see exactly what its doing.

Regards

David

YuriG

David,

You can just try my project (attached) on the demo images I used and then check what is the difference between your image and mine causing different areas. (areas in my test are the same as you can see from previous screenshot).

You can also send me your macro and image. Maybe the problem is in other part of your macro.

Regards,

Yuri

YuriG

One note: the macro posted previously works properly for 3D distance map, which is generated as floating point image. 2D distance map is generated as 16-bit image and the default number of bins in histogram (256) is not enough to accurately split the distance map to equal area rings. The number of bins of the histogram must be set to maximum for 16-bit image (65536), so the underlined code has to be inserted to make the macro work properly with 2D distance maps:
...

        Dim hist As MediaCy.IQL.Operations.McHistogram=ThisApplication.ActiveImage.Histogram
        'calculate accumulated area starting from 1
        If ThisApplication.ActiveImage.Type.BitsPerChannel=16 Then
            'for 2D distance map
            'set maximal number of bins
            'to avoid merging bins (default is 256)
            hist.BinCount=65536
        End If
        Dim TotalArea As Double=0

...

I've also attached the modified project.

Yuri 

davek604

Thanks Yuri, I'll put all this into my macro.

Regards

David