Search algorithm for specific values in a 3D array/matrix being super slow

Question

I'm trying to create code that reads a text file (containing (double) numbers between 0 and 1) and filling up a 3D array/matrix (calling it matrix from now on) with those values. After that is done I need to figure out the current max-value inside this matrix and then look at the values around it (in every possible direction) to check if those values are> than a set threshold. Since you can't get the index of the maxvalue-element in a 3D matrix I'm checking this by using 3 for-loops to find the coordinates of the maxvalue-element.

After having found the current maxvalue I set the element containing it to 0 and get a new maxvalue after everything else is done. The whole process is being repeated until the BubbleFrame Matrix only consists of values that are < Threshold.

Now this code is working, but it's super slow. I used it in Matlab and it only took about 1:30h to find everything (27370 matches aka bubbles) compared to C# where I stopped the code at 6137 bubbles after about 2 hours.

How do I improve my code to improve the operation time?

class Blasensuche
{
 public static Int32 BlasenSuche(ref int Sensor, ref int n_fr, ref string BubbleFile, ref double Threshold)
 {
 bool SearchDone = false;
 //double Threshold = 0.12;
 double[,,] BubbleFrame = new double[Sensor,Sensor,n_fr];
 int[,,] BubbleCollection = new int[Sensor, Sensor, n_fr];
 int BubbleCounter = 0;
 //int BubbleNumber = 0;
 double Maxvalue = 0;
 //int Index = 0;
 int CPlus1 = 0;
 int CMinus1 = 0;
 using (StreamReader BubbleReader = new StreamReader(BubbleFile))
 {
 string calibline = "";
 string[] numbers;
 do
 {
 for (int c = 0; c < n_fr; c++)
 {
 for (int b = 0; b < Sensor; b++)
 {
 calibline= BubbleReader.ReadLine();
 numbers = calibline.Split(' ');
 for (int a = 0; a < Sensor; a++)
 {
 BubbleFrame[b, a, c] = Convert.ToDouble(numbers[a]);
 }
 }
 }
 } while (!BubbleReader.EndOfStream);
 }
 while (SearchDone == false)
 {
 Maxvalue = BubbleFrame.Cast<double>().Max();
 if (Maxvalue < Threshold)
 {
 SearchDone = true;
 break;
 }
 BubbleCounter++;
 //Blasenanzahl.SetLabel(BubbleCounter);
 for (int c = 0; c < n_fr; c++)
 {
 for (int b = 0; b < Sensor; b++)
 {
 for (int a = 0; a < Sensor; a++)
 {
 if (BubbleFrame[b, a, c] == Maxvalue)
 {
 BubbleFrame[b, a, c] = 0D;
 BubbleCollection[b, a, c] = BubbleCounter;
 if (c > 0)
 {
 CMinus1 = c - 1;
 BubbleSearchCMinus1(ref BubbleFrame, ref BubbleCollection, ref CMinus1, ref b, ref a, ref Threshold, ref BubbleCounter, ref Sensor);
 }
 BubbleSearchC(ref BubbleFrame, ref BubbleCollection, ref c, ref b, ref a, ref Threshold, ref BubbleCounter, ref Sensor);
 CPlus1 = c + 1;
 if (CPlus1 < n_fr)
 {
 BubbleSearchCPlus1(ref BubbleFrame, ref BubbleCollection, ref CPlus1, ref b, ref a, ref Threshold, ref BubbleCounter, ref Sensor, ref n_fr);
 }
 }
 }
 }
 }
 //Index = Array.IndexOf(BubbleFrame, Maxvalue);
 }
 return BubbleCounter;
 }
 public static void BubbleSearchCMinus1(ref double[, ,] BubbleFrame,ref int[, ,] BubbleCollection, ref int CMinus1, ref int b, ref int a, ref double Threshold, ref int BubbleCounter, ref int Sensor)
 {
 //C-1 Ebene
 int BPlus1 = b + 1; int BMinus1 = b - 1;
 int APlus1 = a + 1; int AMinus1 = a - 1;
 if (b > 0 && a > 0 && BubbleFrame[BMinus1, AMinus1, CMinus1] > Threshold) //1
 {
 BubbleFrame[BMinus1, AMinus1, CMinus1] = 0D;
 BubbleCollection[BMinus1, AMinus1, CMinus1] = BubbleCounter;
 }
 if (b > 0 && BubbleFrame[BMinus1, a, CMinus1] > Threshold) //2
 {
 BubbleFrame[BMinus1, a, CMinus1] = 0D;
 BubbleCollection[BMinus1, a, CMinus1] = BubbleCounter;
 }
 if (b > 0 && (APlus1) < Sensor && BubbleFrame[BMinus1, APlus1, CMinus1] > Threshold)//3
 {
 BubbleFrame[BMinus1, APlus1, CMinus1] = 0D;
 BubbleCollection[BMinus1, APlus1, CMinus1] = BubbleCounter;
 }
 if (a > 0 && BubbleFrame[b, AMinus1, CMinus1] > Threshold)//4
 {
 BubbleFrame[b, AMinus1, CMinus1] = 0D;
 BubbleCollection[b, AMinus1, CMinus1] = BubbleCounter;
 }
 if (BubbleFrame[b, a, CMinus1] > Threshold)//5
 {
 BubbleFrame[b, a, CMinus1] = 0D;
 BubbleCollection[b, a, CMinus1] = BubbleCounter;
 }
 if ((APlus1) < Sensor && BubbleFrame[b, APlus1, CMinus1] > Threshold)//6
 {
 BubbleFrame[b, APlus1, CMinus1] = 0D;
 BubbleCollection[b, APlus1, CMinus1] = BubbleCounter;
 }
 if ((BPlus1) < Sensor && a > 0 && BubbleFrame[BPlus1, AMinus1, CMinus1] > Threshold)//7
 {
 BubbleFrame[BPlus1, AMinus1, CMinus1] = 0D;
 BubbleCollection[BPlus1, AMinus1, CMinus1] = BubbleCounter;
 }
 if ((BPlus1) < Sensor && BubbleFrame[BPlus1, a, CMinus1] > Threshold)//8
 {
 BubbleFrame[BPlus1, a, CMinus1] = 0D;
 BubbleCollection[BPlus1, a, CMinus1] = BubbleCounter;
 }
 if ((BPlus1) < Sensor && (APlus1) < Sensor && BubbleFrame[BPlus1, APlus1, CMinus1] > Threshold)//9
 {
 BubbleFrame[BPlus1, APlus1, CMinus1] = 0D;
 BubbleCollection[BPlus1, APlus1, CMinus1] = BubbleCounter;
 }
 }
 public static void BubbleSearchC(ref double[, ,] BubbleFrame, ref int[, ,] BubbleCollection, ref int c, ref int b, ref int a, ref double Threshold, ref int BubbleCounter, ref int Sensor)
 {
 int BPlus1 = b + 1; int BMinus1 = b - 1;
 int APlus1 = a + 1; int AMinus1 = a - 1;
 //C Ebene
 if (b > 0 && a > 0 && BubbleFrame[BMinus1, AMinus1, c] > Threshold) //1
 {
 BubbleFrame[BMinus1, AMinus1, c] = 0D;
 BubbleCollection[BMinus1, AMinus1, c] = BubbleCounter;
 }
 if (b > 0 && BubbleFrame[BMinus1, a, c] > Threshold) //2
 {
 BubbleFrame[BMinus1, a, c] = 0D;
 BubbleCollection[BMinus1, a, c] = BubbleCounter;
 }
 if (b > 0 && (APlus1) < Sensor && BubbleFrame[BMinus1, APlus1, c] > Threshold)//3
 {
 BubbleFrame[BMinus1, APlus1, c] = 0D;
 BubbleCollection[BMinus1, APlus1, c] = BubbleCounter;
 }
 if (a > 0 && BubbleFrame[b, AMinus1, c] > Threshold)//4
 {
 BubbleFrame[b, AMinus1, c] = 0D;
 BubbleCollection[b, AMinus1, c] = BubbleCounter;
 }
 /*if (BubbleFrame[b, a, c] > Threshold)//5 Entfällt!
 {
 BubbleFrame[b, a, c] = 0D;
 BubbleCollection[b, a, c] = BubbleCounter;
 }*/
 if ((APlus1) < Sensor && BubbleFrame[b, APlus1, c] > Threshold)//6
 {
 BubbleFrame[b, APlus1, c] = 0D;
 BubbleCollection[b, APlus1, c] = BubbleCounter;
 }
 if ((BPlus1) < Sensor && a > 0 && BubbleFrame[BPlus1, AMinus1, c] > Threshold)//7
 {
 BubbleFrame[BPlus1, AMinus1, c] = 0D;
 BubbleCollection[BPlus1, AMinus1, c] = BubbleCounter;
 }
 if ((BPlus1) < Sensor && BubbleFrame[BPlus1, a, c] > Threshold)//8
 {
 BubbleFrame[BPlus1, a, c] = 0D;
 BubbleCollection[BPlus1, a, c] = BubbleCounter;
 }
 if ((BPlus1) < Sensor && (APlus1) < Sensor && BubbleFrame[BPlus1, APlus1, c] > Threshold)//9
 {
 BubbleFrame[BPlus1, APlus1, c] = 0D;
 BubbleCollection[BPlus1, APlus1, c] = BubbleCounter;
 }
 }
 public static void BubbleSearchCPlus1(ref double[, ,] BubbleFrame, ref int[, ,] BubbleCollection, ref int CPlus1, ref int b, ref int a, ref double Threshold, ref int BubbleCounter, ref int Sensor, ref int n_fr)
 {
 //C+1 Ebene
 int BPlus1 = b + 1; int BMinus1 = b - 1;
 int APlus1 = a + 1; int AMinus1 = a - 1;
 if (b > 0 && a > 0 && BubbleFrame[BMinus1, AMinus1, CPlus1] > Threshold) //1
 {
 BubbleFrame[BMinus1, AMinus1, CPlus1] = 0D;
 BubbleCollection[BMinus1, AMinus1, CPlus1] = BubbleCounter;
 }
 if (b > 0 && BubbleFrame[BMinus1, a, CPlus1] > Threshold) //2
 {
 BubbleFrame[BMinus1, a, CPlus1] = 0D;
 BubbleCollection[BMinus1, a, CPlus1] = BubbleCounter;
 }
 if (b > 0 && (APlus1) < Sensor && BubbleFrame[BMinus1, APlus1, CPlus1] > Threshold)//3
 {
 BubbleFrame[BMinus1, APlus1, CPlus1] = 0D;
 BubbleCollection[BMinus1, APlus1, CPlus1] = BubbleCounter;
 }
 if (a > 0 && BubbleFrame[b, AMinus1, CPlus1] > Threshold)//4
 {
 BubbleFrame[b, AMinus1, CPlus1] = 0D;
 BubbleCollection[b, AMinus1, CPlus1] = BubbleCounter;
 }
 if (BubbleFrame[b, a, CPlus1] > Threshold)//5
 {
 BubbleFrame[b, a, CPlus1] = 0D;
 BubbleCollection[b, a, CPlus1] = BubbleCounter;
 }
 if ((APlus1) < Sensor && BubbleFrame[b, APlus1, CPlus1] > Threshold)//6
 {
 BubbleFrame[b, APlus1, CPlus1] = 0D;
 BubbleCollection[b, APlus1, CPlus1] = BubbleCounter;
 }
 if ((BPlus1) < Sensor && a > 0 && BubbleFrame[BPlus1, AMinus1, CPlus1] > Threshold)//7
 {
 BubbleFrame[BPlus1, AMinus1, CPlus1] = 0D;
 BubbleCollection[BPlus1, AMinus1, CPlus1] = BubbleCounter;
 }
 if ((BPlus1) < Sensor && BubbleFrame[BPlus1, a, CPlus1] > Threshold)//8
 {
 BubbleFrame[BPlus1, a, CPlus1] = 0D;
 BubbleCollection[BPlus1, a, CPlus1] = BubbleCounter;
 }
 if ((BPlus1) < Sensor && (APlus1) < Sensor && BubbleFrame[BPlus1, APlus1, CPlus1] > Threshold)//9
 {
 BubbleFrame[BPlus1, APlus1, CPlus1] = 0D;
 BubbleCollection[BPlus1, APlus1, CPlus1] = BubbleCounter;
 }
 }
}

• Simon Forsberg: Welcome to Code Review! It's an interesting question you have here, I feel that it would be helpful if you could provide some sample input/output data that can be used to test your code. If you also could share what use-case this have it would be very interesting to know! (Although that's not completely necessary if your project is top-secret).
• Stef: Here is a text file you can use as an input (be aware that its very large (120MB)!): dropbox.com/s/aoykvwuxms9hl6m/AllFrames.txt . Set Sensor = 32, n_fr=2000*5 and Threshold = 0.12 . I hope it works and helps! I'm afraid I'm not allowed to tell you to much about what this code is for but : basically it's to evaluate data from a sensor with 32x32 (or 64x64) measurement-points collected over n_fr Frames (time*frequency). This code is supposed to find bubbles (values between 0.12 (default) and 1) and count those bubbles. What's missing is the part that unites adjacent bubbles.
• mjolka: None of the BubbleSearch* methods make use of n_fr, is that a bug?
• Falco: Ok - I wont give you a proper code Review - just the way I would solve this problem for good performance: Use the right data-Structure for the right job! When reading the values and saving them into your 3D-Array fill a second datastructure: An ordered Stack with all the values and their coordinates. Then just always pop the first item of the stack and you got the right coordinates! You don't have to search at all! This will likely triple the memory footprint, but should be blazing fast!
• Ben Aaronson: It looks like you've edited some of mjolka's suggestions into your question. As explained here, it's preferred that you don't do this as it goes against the question and answer nature of the site.

Answer 1

The ref keyword

You don't need to use the ref keyword anywhere. Please see Method Parameters and Types on MSDN.

Naming

Please try to follow C# naming and capitalization conventions.

Globalization

Use the version of Convert.ToDouble that takes an IFormatProvider. The sample file you gave has numbers in the form 0,0421216848673947, which will be converted differently depending on the current thread's culture. On my machine, that string is converted to 421216848673947.

You will want something like this:

var provider = CultureInfo.GetCultureInfo("de-DE").NumberFormat;
...
bubbleFrame[b, a, c] = Convert.ToDouble(numbers[a], provider);

Repetition

There are 26 if statements, one for each surrounding element. That's a lot of repetition, and a lot of chances to make mistakes.

You can write this more succinctly by storing the offsets of the surrounding elements in an array.

private static readonly int[][] Offsets =
 (from i in Enumerable.Range(-1, 3)
 from j in Enumerable.Range(-1, 3)
 from k in Enumerable.Range(-1, 3)
 where !(i == 0 && j == 0 && k == 0)
 select new[] { i, j, k }).ToArray();
...
foreach (var offset in Offsets)
{
 var offsetA = a + offset[0];
 var offsetB = b + offset[1];
 var offsetC = c + offset[2];
 if (IsInRange(offsetA, sensor)
 && IsInRange(offsetB, sensor)
 && IsInRange(offsetC, frames)
 && bubbleFrame[offsetB, offsetA, offsetC] > threshold)
 {
 bubbleFrame[offsetB, offsetA, offsetC] = 0;
 bubbleCollection[offsetB, offsetA, offsetC] = bubbleCounter;
 }
}
...
private static bool IsInRange(int i, int max)
{
 return i >= 0 && i < max;
}

Algorithm

You can get a list of candidate bubbles in one pass through the matrix.

private static IEnumerable<Point> GetCandidateBubbles(int sensor, int frames, double threshold, double[,,] bubbleFrame)
{
 var candidates = new List<Point>();
 for (var c = 0; c < frames; c++)
 {
 for (var b = 0; b < sensor; b++)
 {
 for (var a = 0; a < sensor; a++)
 {
 var value = bubbleFrame[b, a, c];
 if (value >= threshold)
 {
 candidates.Add(new Point(a, b, c, value));
 }
 }
 }
 }
 return candidates;
}

Then process each candidate in order:

foreach (var candidate in GetCandidateBubbles(sensor, frames, threshold, bubbleFrame).OrderByDescending(point => point.Value))
{
 var a = candidate.A;
 var b = candidate.B;
 var c = candidate.C;
 var value = bubbleFrame[b, a, c];
 // We might have already popped the candidate.
 if (value < threshold)
 {
 continue;
 }
 if (!previousBubble.HasValue || previousBubble.Value != value)
 {
 bubbleCounter++;
 }
 ...

Counting bubbles in your 120MB sample file using this approach takes ~15s on my computer, ~13s of which is just reading the input.

(You might also get better performance with jagged arrays.)

• radarbob: Re: ref. Yeah, not needed; but as it does not cause boxing of value-types, I don't see any performance issue. Is there?
• mjolka: @radarbob Using ref is not a performance issue, it's a style issue. Performance is only addressed in my last point.
• Falco: Great Answer! Clear and precise style. But I think there may be a little problem. The original algorithm processed the maximum value in the whole array first, then the next maximum which wasn't popped. In your case a smaller value may be processed first and popp a bigger value? I think you need to sort your Candidates List once biggest-first before processing it!
• Falco: And createion of the Candidates List would be even faster if done when reading in all the values, so just add all read in values > threshold on creation into the Candidates List, sort it and you will probably get another performance gain of a ew seconds :-)
• mjolka: @Stef the second parameter is the number of elements, so Enumerable.Range(-1, 3) will give { -1, 0, 1 }.

Answer 2

Cast-Aways

Make your constants double at compile time.

BubbleFrame[b, a, c] = 0D;, 0.12D, BubbleFrame[b + 1, a - 1, c - 1] = 0D etc. Well, 0.12 is already a double, but what the heck.

BubbleCounter is used as a value in the matrix. Declare it as double and pass it as a double.

Guessing At Performance

Not this:

 calibline= BubbleReader.ReadLine();
 numbers = calibline.Split(' ');
 for (int a = 0; a < Sensor; a++)
 {
 BubbleFrame[b, a, c] = Convert.ToDouble(numbers[a]);
 }

but this:

BubbleFrame[b,a] = Array.ConvertAll(BubbleReader.ReadLine().Split(' '), double.Parse);

I honestly don't know that this is more performant. And I hope the new array does not have > Sensor values!

---

Why this?

Maxvalue = BubbleFrame.Cast<double>().Max();

Aren't BubbleFrame values already double?

---

Calculate once; and avoid unnecessary call:

Cminus1 = c-1;
if (Cminus1 > 0)
 BubbleSearchCMinus1(ref BubbleFrame, ref BubbleCollection, ref Cminus1, ref b, ref a, ref Threshold, ref BubbleCounter, ref Sensor, ref n_fr);

then inside:

 BubbleSearchCMinus1(ref double[, ,] BubbleFrame,ref double[, ,] BubbleCollection, ref int c, ref int b, ref int a, ref double Threshold, ref int BubbleCounter, ref int Sensor, ref int n_fr) {
 Bminus1 = b-1; Bplus1 = b+1;
}

• Snowbody: The Cast<double> is a trick to flatten out a multi-dimensional array.
• Stef: Thanks for all the hints and suggestions! I made use of almost all of it (only thing I didnt use was your version of reading the data, simply because that this operation is only done once and the real search is slowing the whole operation down, I will try it in the near future when I got more time to focus on this part of the code) and updated my Code in the Startpost. Like Snowbody said, the cast is necessary to use .Max() on a multi-dimensional array. God knows why they didn't make .Max() available for that kind of arrays.
• radarbob: @Stef, didn't know that. I'd have thought it would have been BubbleFrame.Max().Cast<double>(); if anything. As I searched (briefly, admittedly) I did not run across this trick.
• mjolka: @radarbob I believe the implicit conversion from int to double in the statement BubbleFrame[b, a, c] = 0 will be performed at compile-time.

Answer 3

You may be using the wrong data structure.

While it is easy to understand a 3-dimensional array, that data structure doesn't make it easy or efficient to solve the problem.

You're basically repeatedly finding the maximum, and then setting it to 0, otherwise known as "extracting" it.

This is a job for a priority queue, with the dominating function being "max". Also known as a MaxHeap. See https://stackoverflow.com/a/13776636/1108056 for one set of code that does it.

Also, since you know the range of the values, why not use a fixed-point type like Decimal or scale them up to longs?

Answer 4

Your current approach is \$O(k*N)\$ where k is the number of max values found/zeroed and \$N\$ is the total number of values. If you add all the values to a max-heap or priority queue, then the algorithm would be \$O(N + k log N)\$. Though, it might take a bit of work to ensure that the value remains linked to it's position in the matrix when it's added to the data structure.

Answer 5

You should keep the max (with its position) in a variable while you fill the matrix. Then once it is filled, you already have the max and don't have to go through the whole matrix.

• Stef: My bad, I should have mentioned that once I found the maxvalue I set the element containing it to 0 and repeat everything (searching for maxvalue and checking for values > threshold next to it) until every element in my matrix is < Threshold.