Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am Mathematica novice and am trying to find a fast way to replace labels in some rather large matrices (512x512 or 1024x1024) that have resulted from segmenting images with WatershedComponents or MorphologicalComponents. At present I have written the function below that takes a "label matrix" (matrix that results from segmentation) and a list of "keeper" labels (labels that I want to remain in place in the returned matrix) and returns a new matrix with only those labels remaining in place that are in the listOfKeepersIn. All other labels are replaced in the returned matrix with a single "background" label, here designated as 100.

I suspect there are many more efficient ways to achieve this than what I have written here, so I welcome your suggestions. I have tried to convert this to a Compiled Function using Compile but haven't had any success (get an error message about "tensors of different ranks"). I also tried to use ParallelTable instead of Table but didn't notice much improvement. Thanks in advance for any help you can provide. -GR

ComponentKeeperFunction[matrixIn_, listOfKeepersIn_] :=

Block[{matrixToReturn},

matrixToReturn = 

Table[

 Which[

      MemberQ[listOfKeepersIn, matrixIn[[r, c]]] == True,

    (* If label value is on the keeper list then just keep that
    label in its place, (i.e. make the value in the table the same as    
    matrixIn[[r,c]])*)

    matrixIn[[r, c]], 


    MemberQ[listOfKeepersIn, matrixIn[[r, c]]] == False,

    (* If label value is not on the keeper list then replace that label with a          
     designated background label  *)        

    100  (*here I designate the background label as having a value of 100 *)

 ] (* close Which *),


{r, 1, Dimensions[matrixIn][[1]]}, {c, 1, Dimensions[matrixIn][[2]]}

 ] (*close Table *);

  Return[matrixToReturn]

 ](*close Block*)
share|improve this question
    
You can remove matrixToReturn = and Return[matrixToReturn. If you do not put a semicolon at the end of Table it will return its result automatically. Also Module is recommended rather than Block, which you can read about here. Overall, using Table like this is typical for someone with a background in procedural languages, but Mathematica is faster with functional programming so you should try to use corresponding functional alternatives, such as Map. –  Pickett Jul 19 at 22:52
    
Thanks for the tips. Good to know about not needing to explicitly Return since then I don't need to use Block or Module as much. -GR –  user13999 Jul 20 at 3:22

3 Answers 3

In order to apply a function to every element we can use Map with the level specification:

Map[If[MemberQ[keep, #], #, 100] &, matrix, {2}]

Another option using the Listable attribute:

Function[{x}, If[MemberQ[keep, x], x, 100], Listable]@matrix;

This turned out to be a lot slower though. I ran these on Michael E2's test case and the Listable approach took 13.87 seconds to complete, while the Map version took only 0.34 seconds. But Michael E2's compiled Listable approach is faster than both.

share|improve this answer

SelectComponents is pretty fast but it labels the background with 0, not 100. You might be able to work with that.

SelectComponents[mat, "Label", MemberQ[keep, #] &]

but this is a bit faster:

sel = Compile[{{label, _Integer}, {keep, _Integer, 1}},
  If[MemberQ[keep, label], label, 0],  (* or 100 if necessary *)
  RuntimeAttributes -> {Listable}, Parallelization -> True]

Example

mat = MorphologicalComponents @ Binarize[ExampleData[{"TestImage", "U2"}], 0.2];
Dimensions@mat
Max[mat]
(*
  {1024, 1024}  -- dimensions
  1923          -- number of components
*)

SeedRandom[1];
keep = RandomSample[Range@Max[mat], Round[Max[mat]/10]]
(*
  {1675, 29, 17, 1552, 1697, 1826, 1022, 1568, 138, 1203, 1085, 1539, 
   ... 713, 114, 567, 656, 1747, 134, 94}
*)
SelectComponents[mat, "Label", MemberQ[keep, #] &]; // AbsoluteTiming
(*
  {0.739382, Null}
*)

sel[mat, keep]; // AbsoluteTiming
(*
  {0.282210, Null}
*)

If I compile to C (CompilationTarget -> "C"), I save another 20% in timing.

If keep is about half of the components, then the timings rise to 2.5 sec. and 0.9 sec. respectively. If keep is larger, it would be more efficient to take its complement and delete.

share|improve this answer
    
Thanks very much! SelectComponents is quite a bit faster and seeing a Compile example I can make work is helpful. Can you suggest any similar methods to efficiently join one or more components to another component? -GR –  user13999 Jul 20 at 3:26
    
@user13999 I think you just have to give them the same label (integer), somewhat like the reverse of the compiled function sel. In a sense, sel joins all the non-keep components to the background. Just switch the logic, If[MemberQ[join, label], Min[join], label], where join is a list of components to join. –  Michael E2 Jul 20 at 3:56
    
Thanks again for the suggestion, I have a joinComponents function working along these lines now and it is much faster than what I was doing before. –  user13999 Jul 20 at 19:43
    
Btw, one thing I am still wondering is if there is some way to use SparseArray or something similar to store the final segmentation mask (i.e. store the resulting label matrix after finishing using joinComponents and SelectComponents or your "sel" function) so it does not take up so much memory (I am doing this for several hundred timelapse images). The actual components of interest typically take up between 5% to 25% of the whole mask so there's a lot of background in each mask. -GR –  user13999 Jul 20 at 19:53

I shall not attempt to replicate the exact function of your code but rather to address the problem posed in text of your Question.

As a starting point I suggest you build a Dispatch table of the replacements you wish to make and then apply it with Replace.

First some sample data:

SeedRandom[0]
m = RandomInteger[66, {1024, 1024}];

keep = Array[Prime, 18];

keep is my label list; I chose prime numbers merely as an example.

Now build the Dispatch table:

rules = Append[Thread[keep -> keep], _ -> 100] // Dispatch;

And apply it:

(result = Replace[m, rules, {2}];) // Timing // First
0.187201

Check the result with a Tally:

result // Flatten // Tally
{{100, 766948}, {43, 15736}, {17, 15664}, {53, 15823}, {19, 15503}, {31, 15726}, {5, 
  15609}, {37, 15492}, {11, 15713}, {13, 15687}, {59, 15828}, {41, 15837}, {29, 
  15722}, {47, 15552}, {2, 15519}, {23, 15648}, {61, 15590}, {7, 15435}, {3, 15544}}

If all of your keep labels are in a sequential band there are much faster numeric methods easily applicable. See: Replace values which obey certain criteria

share|improve this answer
    
Wow, yet another way to do this! I'll give it a try, never came across the Dispatch function before. Also, thanks for the link. -GR –  user13999 Jul 20 at 19:55

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.