title: How to find same-value rectangular areas of a given size in a matrix most efficiently?
body: <p>My problem is very simple but I haven't found an efficient implementation yet.</p>
  
  <p>Suppose there is a matrix A like this:</p>
  
  <pre><code>0 0 0 0 0 0 0
  4 4 2 2 2 0 0
  4 4 2 2 2 0 0
  0 0 2 2 2 1 1
  0 0 0 0 0 1 1
  </code></pre>
  
  <p>Now I want to find all starting positions of rectangular areas in this matrix which have a given size. An area is a subset of A where all numbers are the same.</p>
  
  <p>Let's say width=2 and height=3. There are 3 areas which have this size:</p>
  
  <pre><code>2 2   2 2   0 0
  2 2   2 2   0 0
  2 2   2 2   0 0
  </code></pre>
  
  <p>The result of the function call would be a list of starting positions (x,y starting with 0) of those areas. </p>
  
  <pre><code>List((2,1),(3,1),(5,0))
  </code></pre>
  
  <p>The following is my current implementation. "Areas" are called "surfaces" here.</p>
  
  <pre><code>case class Dimension2D(width: Int, height: Int)
  case class Position2D(x: Int, y: Int)
  
  def findFlatSurfaces(matrix: Array[Array[Int]], surfaceSize: Dimension2D): List[Position2D] = {
  
      val matrixWidth = matrix.length
      val matrixHeight = matrix(0).length
      var resultPositions: List[Position2D] = Nil
  
      for (y &lt;- 0 to matrixHeight - surfaceSize.height) {
          var x = 0
          while (x &lt;= matrixWidth - surfaceSize.width) {
              val topLeft = matrix(x)(y)
              val topRight = matrix(x + surfaceSize.width - 1)(y)
              val bottomLeft = matrix(x)(y + surfaceSize.height - 1)
              val bottomRight = matrix(x + surfaceSize.width - 1)(y + surfaceSize.height - 1)
              // investigate further if corners are equal
              if (topLeft == bottomLeft &amp;&amp; topLeft == topRight &amp;&amp; topLeft == bottomRight) {
                  breakable {
                      for (sx &lt;- x until x + surfaceSize.width;
                           sy &lt;- y until y + surfaceSize.height) {
                          if (matrix(sx)(sy) != topLeft) {
                              x = if (x == sx) sx + 1 else sx 
                              break
                          }
                      }
                      // found one!       
                      resultPositions ::= Position2D(x, y)
                      x += 1
                  }
              } else if (topRight != bottomRight) {
                  // can skip x a bit as there won't be a valid match in current row in this area
                  x += surfaceSize.width 
              } else {
                  x += 1
              }
          }   
      }
      return resultPositions
  }
  </code></pre>
  
  <p>I already tried to include some optimizations in it but I am sure that there are far better solutions. Is there a matlab function existing for it which I could port? I'm also wondering whether this problem has its own name as I didn't exactly know what to google for.</p>
  
  <p>Thanks for thinking about it! I'm excited to see your proposals or solutions :)</p>
  
  <p><strong>EDIT:</strong> Matrix dimensions in my application range from 300x300 to 3000x3000 approximately. Also, the algorithm will only be called <em>once</em> for the same matrix. The reason is that the matrix will always be changed afterwards (approx. 1-20% of it).</p>
  
  <h1>RESULTS</h1>
  
  <p>I implemented the algorithms of Kevin, Nikita and Daniel and benchmarked them in my application environment, i.e. no isolated synthetic benchmark here, but special care was taken to integrate all algorithms in their most performant way which was especially important for Kevin's approach as it uses generics (see below).</p>
  
  <p>First, the raw results, using Scala 2.8 and jdk 1.6.0_23. The algorithms were executed several hundred times as part of solving an application-specific problem. "Duration" denotes the total time needed until the application algorithm finished (of course without jvm startup etc.). My machine is a 2.8GHz Core 2 Duo with 2 cores and 2gig of memory, -Xmx800M were given to the JVM.</p>
  
  <p><strong>IMPORTANT NOTE:</strong> I think my benchmark setup is not really fair for parallelized algorithms like the one from Daniel. This is because the application is already calculating multi-threaded. So the results here probably only show an equivalent to single-threaded speed.</p>
  
  <p>Matrix size 233x587:</p>
  
  <pre><code>                  duration | JVM memory | avg CPU utilization
  original O(n^4) | 3000s      30M          100%  
  original/-server| 840s       270M         100%
  Nikita O(n^2)   | 5-6s       34M          70-80%
  Nikita/-server  | 1-2s       300M         100%
  Kevin/-server   | 7400s      800M         96-98%
  Kevin/-server** | 4900s      800M         96-99%
  Daniel/-server  | 240s       360M         96-99%
  </code></pre>
  
  <p>** with @specialized, to make <a href="http://lamp.epfl.ch/~dragos/files/scala-spec.pdf" rel="nofollow">generics faster</a> by avoiding type erasure</p>
  
  <p>Matrix size 2000x3000:</p>
  
  <pre><code>                  duration | JVM memory | avg CPU utilization
  original O(n^4) | too long   100M         100%  
  Nikita O(n^2)   | 150s       760M         70%
  Nikita/-server  | 295s (!)   780M         100%
  Kevin/-server   | too long, didn't try
  </code></pre>
  
  <p>First, a small note on memory. The -server JVM option uses considerably more memory at the advantage of more optimizations and in general faster execution. As you can see from the 2nd table Nikita's algorithm is slower with the -server option which is obviously due to hitting the memory limit. I assume that this also slows down Kevin's algorithm even for the small matrix as the functional approach is using much more memory anyways. To eliminate the memory factor I also tried it once with a 50x50 matrix and then Kevin's took 5secs and Nikita's 0secs (well, nearly 0). So in any case it's still slower and not just because of memory.</p>
  
  <p>As you can see from the numbers, I will obviously use Nikita's algorithm because it's damn fast and this is absolutely necessary in my case. It can also be parallelized easily as Daniel pointed out. The only downside is that it's not really the scala-way.</p>
  
  <p>At the moment Kevin's algorithm is probably in general a bit too complex and therefore slow, but I'm sure there are more optimizations possible (see last comments in his answer).</p>
  
  <p>With the goal of directly transforming Nikita's algorithm to functional style Daniel came up with a solution which is already quite fast and as he says would even be faster if he could use scanRight (see last comments in his answer).</p>
  
  <p><strong>What's next?</strong></p>
  
  <p>At the technological side: waiting for Scala 2.9, ScalaCL, and doing synthetic benchmarks to get raw speeds.</p>
  
  <p>My goal in all this is to have functional code, BUT only if it's not sacrificing too much speed.</p>
  
  <p><strong>Choice of answer:</strong></p>
  
  <p>As for choosing an answer, I would want to mark Nikita's and Daniel's algorithms as answers but I have to choose one. The title of my question included "most efficiently", and one is the fastest in imperative and the other in functional style. Although this question is tagged Scala I chose Nikita's imperative algorithm as 2s vs. 240s is still too much difference for me to accept. I'm sure the difference can still be pushed down a bit, any ideas?</p>
  
  <p>So, thank you all very very much! Although I won't use the functional algorithms <em>yet</em>, I got many new insights into Scala and I think I slowly get an understanding of all the functional crazyness and its potential. (of course, even without doing much functional programming, Scala is much more pleasing than Java... that's another reason to learn it)</p>