The approach you outline is simple and useful, but suffers from terrible artifacts as shown. Avoid it. You need a parallel growth algorithm; for a single-threaded model, a round-robin approach follows:
- Randomly place various points in your map space. Normalise their distribution (avoids ugly clustering) using Gaussian distribution or by applying an iterative relaxation algorithm to move randomly placed points away from one another, as per Lloyd relaxation for Voronoi Diagrams. These points represent the centroids of your rooms-to-be.
- (Parallel / Repeat round-robin for all rooms) From each point, grow 4 vertices outward (a rectangular room), moving from the centre point per global iteration (you can grow different rooms at different rates in each axis rather than the same for all, and see how this turns out -- possibly for a more organic / varied result). At some point, some of your rectangles will begin to press against each other. At that point, limit growth in that axis, ensure the two rooms' adjacent edges are exactly touching, and move on.
- Repeat step 2, growing every room incrementally, till all growth is restricted by adjacent rooms or the bounds of the map.
- This will still leave some empty spaces. The problem now becomes one of locating and making rooms out of non-occupied spaces. In fact, if your underlying space is a (integer) grid (and every growth iteration snaps to that grid), then this is much easier to deal with, since you can maintain lists of occupied and unoccupied grid cells: Once you've placed and grown all your rooms, search through the unoccupied list for discrete spaces consisting of groups of adjacent cells. As many unused spaces will have non-rectangular shapes, you'll need to pick a cell at random from within that non-rectangular space, and grow it to its maximal size just as you did with rooms in step 2. This process may recurse within said non-rectangular space, till it is completely filled.
- Repeat step 4 till your map is 100% occupied.