Algorithm to slice into square matrix a matrix

i'm searching for a algorithm that take a matrix (in fact, a double entry array) and return an array of matrix that:
is square (WIDTH = HEIGHT)
all of the element in the matrix has the same value.
I don't know if that is clear, so imagine that you have a image made of pixels that is red, blue or green and i want to get an array that contained the least possible squares. Like the pictures shows

EDIT:

Ok, maybe it's not clear: I've a grid of element that can have some values like that:

0011121

0111122

2211122

0010221

0012221

That was my input, and i want in output somethings like that:

|0|0|111|2|1|

|0|1|111|22|

|2|2|111|22|

|00|1|0|22|1|

|00|1|2|22|1|

When each |X| is an array that is a piece of the input array.
My goal is to minimize the number of output array

This problem does not seem to have an efficient solution.

Consider a subset of instances of your problem defined as follows:

• There are only 2 values of matrix elements, say 0 and 1.
• Consider only matrix elements with value 0.
• Identify each matrix element m_ij with a unit square in a rectangular 2D grid whose lower left corner has the coordinates (i, n-j).
• The set of unit squares SU chosen this way must be 'connected' and must not have 'holes'; formally, for each pair of units squares (m_ij, m_kl) \in SU^2: (i, j) != (k, l) there is a sequence <m_ij = m_i(0)j(0), m_i(1)j(1), ..., m_i(q)j(q) = m_kl> of q+1 unit squares such that (|i(r)-i(r+1)| = 1 _and_ j(r)=j(r+1)) _or_ (i(r)=i(r+1) _and_ |j(r)-j(r+1)| = 1 ); r=0...q (unit squares adjacent in the sequence share one side), and the set SUALL of all unit squares with lower left corner coordinates from the integers minus SU is also 'connected'.

Slicing matrices that admit for this construction into a minimal number of square submatrices is equivalent to tiling the smallest orthogonal polygon enclosing SU ( which is the union of all elements of SU ) into the minimum number of squares.

This SE.CS post gives the references (and one proof) that show that this problem is NP-complete for integer side lengths of the squares of the tiling set.

Note that according to the same post, a tiling into rectangles runs in polynomial time.

Some hints may be useful.<br>
<br> For representation of reduced matrix, maybe a vector is better because it's needed to be stored (start_x,start_y,value ... not sure if another matrix very useful).<br>
<br>Step 1: loop on x for n occurrences (start with y=0)
<br>Step 2: loop on y for/untill n occurrences. Most of cases here will be m lees then n.
(case m greater then n excluded since cannot do a square) Fine, just keep the min value[m]
<br>Step 3: mark on vector (start_x,start_y, value)
<br> Repeat Step 1-3 from x=m until end x
<br>Step 4: End x, adjust y starting from most left_x found(m-in vector, reiterate vector).
...
<br>keep going till end matrix.<br>
<br>Need to be very careful of how boundary are made(squares) in order to include in result full cover of initial matrix.
<br> Reformulate full-initial matrix can be recomposed exactly from result vector.
<br>(need to find gaps and place it on vector derived from step_4)
<br><br> Note ! This is not a full solution, maybe it's how to start and figure out on each steps what is to be adjusted.