Reduce the total internal border of a set of touching rectangles (using graphs)
Here's a naive attempt at a greedy solution.

I'm assuming the initial input is a polyomino (polygon made up of grid aligned squares) on a grid as in your "real world example". 

Let $b = (i, j, w, h)$ be the rectangle with lower left corner containing grid location $(i, j)$ with width $w$ and height $h$. A rectangle is valid if and only if all of its internal grid cells are part of the polyomino. 

For a valid rectangel $r$, the boundary edges of $r$ will either also be boundary edges of the polyomino itself, or internal edges. Denote by $B_r$ the number of boundary edges of the rectangle that are internal to the polyomino and $I_r$ the number of grid edges that are internal to both the boundary and the rectangle. In the teal rectangle in my example below, $I_r = 17$ and $B_r = 6$ (the red sketchy edges). Now define the quality of the rectangle $r$ to be $Q_r = I_r - B_r$. 

The thought here is that the quality of a rectangle measures how good it is at eliminating internal edges. 

Now, the greedy solution I have in mind is: 

Generate all of the $O(n^4)$ possible rectangles $r$. 
For each rectangle $r$, check if $r$ is valid. 
* If $r$ is valid, compute its quality score and maintain the rectangle with the maximum quality. 
End for loop. 
Add the rectangle that maximized the quality to the output set, and color the squares it covers white. Repeat until all squares are white. 

Unfortunately this is $O(n^5)$ where $n$ is the width and height of the grid. But if your problem size really is like the image, then that may not be an issue.

Also, I don't have a proof that it's optimal, but I'll think about it. 



Ok, I wrote some code for this. It's quick and dirty, but it returns exactly your optimal result on your original small example. It is quick and dirty, so don't judge ;). 

from functools import reduce

# Generates all valid rectangles that have a lower left hand corner in i and j
def validRectanglesFrom(grid, i, j):
    result = []
    w = 1
    while i + w &lt;= len(grid) and grid[i + w - 1][j] == 1:
        h = 1
        row_is_valid = True
        while j + h &lt;= len(grid[i + w - 1]) and row_is_valid:
            for k in range(i, i + w):
                if grid[k][j + h - 1] == 0: 
                    row_is_valid = False
                    break
            if row_is_valid: 
                result.append((i, j, w, h))
            h = h + 1
        w = w + 1
    return result

def left_edge_is_internal(grid, i, j):
    return i &gt; 0 and grid[i - 1][j] == 1 and grid[i][j] == 1

def right_edge_is_internal(grid, i, j):
    return i &lt; len(grid) - 1 and grid[i][j] == 1 and grid[i + 1][j] == 1

def bottom_edge_is_internal(grid, i, j):
    return j &gt; 0 and grid[i][j - 1] == 1 and grid[i][j] == 1

def top_edge_is_internal(grid, i, j):
    return j &lt; len(grid[i]) - 1 and grid[i][j] == 1 and grid[i][j + 1] == 1

def bottom_boundary_penalty(grid, rect): 
    return len(
        list(
            filter(
                lambda isInt: isInt, 
                map(
                    lambda i: bottom_edge_is_internal(grid, i, rect[1]), 
                    range(rect[0], rect[0] + rect[2])
                )
            )
        )
    )

def top_boundary_penalty(grid, rect): 
    return len(
        list(
            filter(
                lambda isInt: isInt, 
                map(
                    lambda i: top_edge_is_internal(grid, i, rect[1] + rect[3] - 1), 
                    range(rect[0], rect[0] + rect[2])
                )
            )
        )
    )

def left_boundary_penalty(grid, rect):
    return len(
        list(
            filter(
                lambda isInt: isInt, 
                map(
                    lambda j: left_edge_is_internal(grid, rect[0], j), 
                    range(rect[1], rect[1] + rect[3])
                )
            )
        )
    )

def right_boundary_penalty(grid, rect):
    return len(
        list(
            filter(
                lambda isInt: isInt, 
                map(
                    lambda j: right_edge_is_internal(grid, rect[0] + rect[2] - 1, j), 
                    range(rect[1], rect[1] + rect[3])
                )
            )
        )
    )

def boundary_penalty(grid, rect):
    return (
        bottom_boundary_penalty(grid, rect) + 
        top_boundary_penalty(grid, rect) +
        left_boundary_penalty(grid, rect) + 
        right_boundary_penalty(grid, rect)
    )

def internal_grid_edges(rect):
    return ((rect[2] - 1) * rect[3]) + (rect[2] * (rect[3] - 1))

def quality(grid, rect):
    return internal_grid_edges(rect) - boundary_penalty(grid, rect)


def get_rectangles(grid):
    result = []
    grid_locations = [(i, j) for i in range(0, len(grid)) for j in range(0, len(grid[i]))]
    unfilled = {cell for cell in grid_locations if grid[cell[0]][cell[1]] == 1}
    while len(unfilled) &gt; 0:
        rects = [rect for cell in unfilled for rect in validRectanglesFrom(grid, cell[0], cell[1])]
        qualityScores = list(map(lambda rect: quality(grid, rect), rects))
        bestRectIdx = reduce(lambda i, j: i if qualityScores[i] &gt; qualityScores[j] else j, range(len(qualityScores)))
        bestRect = rects[bestRectIdx]
        result.append(bestRect)
        for i in range(bestRect[0], bestRect[0] + bestRect[2]):
            for j in range(bestRect[1], bestRect[1] + bestRect[3]):
                unfilled.remove((i,j))
                grid[i][j] = 0
    return result

# Your initial example
grid = [[1, 1, 1], [0, 1, 0], [0, 1, 1], [0, 1, 0]] 
get_rectangles(grid)

# Your "real world" example
grid2 = [[0,0,0,0,0,0,0,0,1,1,1,0,0,0,0],
        [0,0,0,0,0,0,0,0,0,0,1,1,1,0,0],
        [0,0,0,0,0,0,0,0,0,1,1,1,0,0,0],
        [0,1,0,0,0,0,0,0,0,1,0,1,1,0,0],
        [1,1,0,0,0,0,0,0,0,0,0,1,1,1,0],
        [1,1,1,1,0,0,0,0,1,0,0,1,0,1,0],
        [0,0,0,1,1,0,0,1,1,1,1,1,0,1,0],
        [1,1,1,1,1,0,1,1,1,1,1,1,0,0,0],
        [1,1,1,0,1,0,1,1,1,1,1,0,0,0,0],
        [0,0,1,1,1,1,1,0,0,1,1,1,0,0,0],
        [0,0,0,1,0,0,0,0,0,0,1,0,0,0,0],
        [0,0,0,1,1,1,0,0,0,0,0,0,0,0,0],
        [0,1,1,1,0,1,1,0,0,0,0,0,0,0,0],
        [0,1,1,1,0,1,0,0,0,0,0,0,0,0,0],
        [0,1,0,1,1,1,1,1,0,0,0,0,0,0,0]]
get_rectangles(grid2)


Here's the result my code returns for your "real world" example.