Describe a 2D Bitmap Using Filled Rectangle Mapping
Rust

Try it online

Average case num rectangles: 57.9054 (10000 iterations of 20x20 matrices, with 1/7 probability of "0")

Solution to example input:
{10,0,13,3,8,1,14,3,4,2,15,2,5,3,15,3,0,5,6,7,17,5,21,5,0,6,10,7,20,6,21,6,0,7,15,7}

The algorithm is pretty simple: for each horizontal line of 1s, use that as the top side of a rectangle, extend this rectangle down while possible, and if it covers some previously uncovered area, add it to the list. This assumes rectangles can overlap. I have not verified whether the output is always correct, but the example input looks ok. The output was verified against the input for 10000 iterations, and it looks correct.

The code includes a benchmark function to calculate the score, and a function to draw the result. The input can be changed in the main function.

extern crate rand;
use rand::{Rng, thread_rng};

#[derive(Copy, Clone, Debug)]
struct Rect {
    x0: usize,
    y0: usize,
    x1: usize,
    y1: usize
}

impl Rect {
    fn contains_xslice(&amp;self, xa: usize, xb: usize) -&gt; bool {
        self.x0 &lt;= xa &amp;&amp; xb &lt;= self.x1
    }
    fn contains(&amp;self, r: &amp;Rect) -&gt; bool {
        // r is completely inside of self
        self.x0 &lt;= r.x0 &amp;&amp; r.x1 &lt;= self.x1 &amp;&amp; self.y0 &lt;= r.y0 &amp;&amp; r.y1 &lt;= self.y1
    }
}

fn display_rects(r: &amp;[Rect]) -&gt; String {
    let s = r.iter()
        .map(|z| format!("{},{},{},{}", z.x0, z.y0, z.x1, z.y1))
        .collect::&lt;Vec&lt;_&gt;&gt;()
        .join(",");

    // transform 1,1 into {1,1}
    format!("{{{}}}", s)
}

fn draw_rects(r: &amp;[Rect], w: usize, h: usize) -&gt; String {
    let mut s = format!("");
    for y in 0..h {
        for x in 0..w {
            let p = Rect { x0: x, x1: x, y0: y, y1: y };
            if r.iter().any(|a| a.contains(&amp;p)) {
                s.push_str("1");
            } else {
                s.push_str("0");
            }
        }
        s.push_str("\n");
    }
    s
}

#[derive(Clone, Debug, PartialEq, Eq)]
struct Matrix {
    v: Vec&lt;u8&gt;,
    size: (usize, usize),
}

impl Matrix {
    fn new(w: usize, h: usize) -&gt; Self {
        Self {
            v: vec![0; w*h],
            size: (w, h),
        }
    }
    fn from_str(s: &amp;str) -&gt; Self {
        let mut v = vec![];
        let mut w = None;
        let mut h = 0;
        for l in s.lines() {
            h += 1;
            let mut tw = 0;
            for c in l.chars() {
                tw += 1;
                let x = match c {
                    '0' =&gt; 0,
                    '1' =&gt; 1,
                    a =&gt; panic!("Invalid input! {:?} at {},{}", a, tw, h),
                };
                v.push(x);
            }
            if let Some(ww) = w {
                if tw != ww {
                    panic!("All lines must have the same width! {} != {}", tw, ww);
                }
            } else {
                w = Some(tw);
            }
        }

        let w = w.unwrap_or(0);

        Matrix {
            v,
            size: (w, h),
        }
    }
    fn random(w: usize, h: usize) -&gt; Self {
        let mut rng = thread_rng();
        let mut v = vec![];

        for _ in 0..w*h {
            let n = rng.gen_range(0, 6+1);
            let x = if n == 0 { 0 } else { 1 };
            v.push(x);
        }

        Matrix {
            v,
            size: (w, h),
        }
    }
    fn get(&amp;self, x: usize, y: usize) -&gt; u8 {
        self.v[y * self.size.0 + x]
    }
    fn hline(&amp;self, y: usize) -&gt; &amp;[u8] {
        &amp;self.v[y * self.size.0 .. (y + 1) * self.size.0]
    }
    fn hline_1_pieces(&amp;self, y: usize) -&gt; Vec&lt;Rect&gt; {
        let mut r = vec![];
        let line = self.hline(y);
        let xmax = self.size.0 - 1;
        let y0 = y;
        let y1 = y;
        let mut lineit = line.iter();
        let mut xoff = 0;
        while lineit.len() &gt; 0 {
            // Find first non-zero pixel
            if let Some(x0) = lineit.position(|&amp;x| x != 0) {
                let x0 = xoff + x0;
                // x0..x1 is a line of 11111
                let x1 = x0 + lineit.position(|&amp;x| x == 0).unwrap_or(xmax - x0);
                r.push(Rect { x0, x1, y0, y1 });
                // Remember offset because the iterator is being consumed
                xoff = x1 + 2;
            }
        }

        r
    }
    fn rectangles(&amp;self) -&gt; Vec&lt;Rect&gt; {
        let mut r: Vec&lt;Rect&gt; = vec![];
        let ymax = self.size.1 - 1;
        // For each hline, add rectangles covering the maximum width possible
        let hpieces: Vec&lt;Vec&lt;Rect&gt;&gt; = (0..self.size.1).map(|i| self.hline_1_pieces(i)).collect();
        // Extend each rectangle down while possible
        for (y0, layer) in hpieces.iter().enumerate() {
            for rect in layer {
                let mut y = y0 + 1;
                while y &lt;= ymax {
                    if hpieces[y].iter().any(|&amp;r| r.contains_xslice(rect.x0, rect.x1)) {
                        y += 1;
                    } else {
                        break;
                    }
                }
                let y1 = y - 1;
                let rn = Rect { y1, ..*rect };
                if r.iter().any(|&amp;r1| r1.contains(&amp;rn)) {
                    // An equivalent rectangle is already in the vec
                } else {
                    // This rectangle covers some new area, push it
                    r.push(rn);
                }
            }
        }

        r
    }
}

fn benchmark(iterations: usize) -&gt; f64 {
    let mut sum = 0;
    for _ in 0..iterations {
        let m = Matrix::random(20, 20);
        let r = m.rectangles();
        let output = draw_rects(&amp;m.rectangles(), 20, 20);
        assert_eq!(m, Matrix::from_str(&amp;output));
        sum += r.len();
    }
    sum as f64 / iterations as f64
}

fn main() {
    let input = "\
0000000000111100000000
0000000011111110000000
0000111111111111000000
0000011111111111000000
0000000000000000000000
1111111000000000011111
1111111111100000000011
1111111111111111000000
";
    let m = Matrix::from_str(input);
    assert_eq!(input, draw_rects(&amp;m.rectangles(), 22, 8));
    println!("{}", display_rects(&amp;m.rectangles()));
    println!("{}", benchmark(10000));
}