Clustering with shape prior EDIT STARTS
After seeing the comments and answers, I believe I started in a wrong direction. I have a set of rectangles, which I want to cluster as shown below. . The approach I took was to consider central points of each rectangle as a data point in $R^2$ and cluster them using Euclidean distance (K-means, K-mediods approach or any other method). I am not trying to discover "shapes" in the data, as I know the best shape would be rectangles. This is the reason I was trying to know if there is any way to include "shape priors" in the clustering methods. There might be algorithmic formulations of this problem, I would be glad if someone can point me to. 
I hope it is clear now what I want to do, I will be glad if someone can point me to right directions.
EDIT ENDS
Suppose, we have a set of points in $R^d$ which we want to cluster. We do have some idea about the shape of the clusters, i.e., we know the clusters should be spherical, ellipsoidal or rectangular. Is there a way to include this "shape prior" in the clustering process? (it would be more useful for me if this has been done in some well known algorithms like k-means)
I did find some papers like this one and this one which talks about imposing some constraint into clustering, but the constraints are primarily in the form of "pairwise constraints" i.e. specifying whether two instances should be in same or different clusters. 
 A: There is the frequent claim that k-means "prefers" spherical clusters. Mathematically, it produces Voronoi cells, but there exists a close relationship between Voronoi cells, nearest neighbors and euclidean spheres.
In R, when you look at the Mclust function, you can specify which model to use. If you allow variances, it will model the data using spheres of different size, if you add in covariance they can also be rotated.
I'm not aware of an approach preferring rectangular clusters in any meaningful way.
You could however use an R*-tree index, and use the index pages as clusters.
The question you didn't answer is: what are you trying to do?
If you have prerequirements such as rectangles, you are most likely not trying to discover structure, but instead you want to squeeze your data into a predefined model. At which point you are in the domain of constraint optimization and data modelling, not structure discovery (clustering).
When specifying such constraints, you should also talk about how to handle model overlaps and these things.
A: Okay, things are a bit clearer with your updated question. Thoughts that occur:


*

*How dense is the data within each cluster compared to the rest of the space?


In your diagrams the square regions all look separated by fairly wide gaps. If you have enough samples within each cluster, and you have few/zero outliers, then you might be able to just ignore the shape prior and cluster based on pairwise sample distances. A graph cut type approach, or spectral clustering, for example.


*

*Have you tried simple clustering algorithms but with different distance metrics?


You rectangular regions all look axis aligned, yet you talk about euclidean distance clustering. Have you tried tweaking an algorithm to use $L_1$ distances instead? It won't be ideal (some of the above shapes are cuboids, not cubes), but if your regions are suitable spread out, it might be good enough.


*

*Have you tried parametric modelling?


It's going to be fairly easy to define a parametric model for a cuboid, and an associated likelihood for samples based on this (e.g. 1/(vol. of cuboid) for samples within it, zero otherwise). You could do this for a mixture model of cuboids. From here you could take an optimisation route (throw the log likelihood, plus any log prior or regularisation terms you desire, into your favourite optimisation algorithm) or you could do the whole bayesian model selection thing, marginalising over the parameters. My biggest worry here is that if your dimensionality is high, you may end up with a lot of parameters to deal with.
Any of that sound doable, or are there issues with it I've missed?

EDIT: As requested, more details on parametric modelling.
Wikipedia isn't a bad place to start (http://en.wikipedia.org/wiki/Parametric_model). 
Let's talk about this particular problem though. It's easiest to think about it in 1D at first. Let's say we have a uniform distribution in a 1d space, specified by a start point $a$ and an end point $b$. These are our parameters. If we observe a sample $x$ in our 1D space, the probability of that sample given $a$ and $b$ is
$$
p(x|a,b) = \begin{cases} \frac{1}{|b-a|} & \text{if  } a \leq x < b \\ 0 & \text{ otherwise}\end{cases}
$$
which is also known as the likelihood of the parameters.
We can extend this beyond one dimension very easily. In that case, $a$ and $b$ become opposite corner vertices in our space, the test $a \leq x < b$ is performed in each dimension separately,  and $|b-a|$ become the volume of the cuboid they describe.
E.G. for a 2D space maybe $a = (2,2)$ and $b = (4,5)$. Any $x$ which lies within the cuboid they specify would have a $p(x|a,b)$ equal to $1/((4-2)*(5-2)) = 1/6$. A  diagram is probably makes this clearer...
      _____________.(4,5)
      |            |
      | p(x|a,b) = |
      | 1/6        |   p(x|a,b) = 0
(2,2).|____________|

You, however, don't have one cuboid cluster. You have a number of them. So it makes sense to use a mixture model (https://en.wikipedia.org/wiki/Mixture_model). A mixture model is where we take a group of simple distirbutions and add them up to make a more complicated one - just like the way the distribution you're interested in is equal to the sum of lots of cuboid distributions. Let's start with a two component model. In this situation you have two cuboid clusters in your space, so we need to specify four corner parameters $a_1, b_1, a_2, b_2$ (two for each cuboid) and two weighting parameter $w_1$ and $w_2$ which are always positive and sum to 1 (to ensure our resulting distribution will integrate to 1). Then we have
$$
p(x|a_1, b_1, a_2, b_2, w) = w_1*p(x|a_1,b_1) + w_2*p(x|a_2,b_2)
$$
where $p(x|a_1,b_1)$ is just the likelihood for a single cuboid we defined earlier. The $w$ parameters control how common it is for a sample to be part of a particular cluster, and the $a$ and $b$ parameters control where that particular cluster is. If we add another cluster to our earlier example and set $w_1 = w_2 = 0.5$ we'd have something which looks like this:
      ______________.(4,5)
      |             |
      |w1*p(x|a1,b1)|   p(x|a1,b1, w1, a2, b2, w2) = 0
      | =1/12       |  
(2,2).|_____________|    ________________________________. (2,10)
                        |   w2*p(x|a2,b2) = 1/8          |
                  (1,6).|________________________________|

which is already starting to look a bit like your examples.


*

*It's easy to extend this from two components to $K$ components:


$$
p(x|a_1, b_1, w_1 \ldots a_K, b_K, w_K) = \sum_{k=1}^K w_k*p(x|a_k, b_k)
$$
where $w_k\geq 0$ and $\sum_{k=1}^K w_k = 1$. 


*

*It's also easy to extend this to many observed sample instead of just 1. If you have $N$ samples, then assuming they are iid (https://en.wikipedia.org/wiki/Iid) given the parameters, we have


$$
p(x_1 \ldots x_N|a_1,b_1, w_1, \ldots a_K, b_K, w_K) = \prod_{i=1}^N p(x_i|a_1,b_1, w_1, \ldots a_K, b_K, w_K).
$$
That, then, gives you the likelihood of a $K$ component cuboid mixture model, evaluated on a particular dataset of $N$ samples. You could try optimising that directly to get the maximum likelihood result (protip - optimise the log likelihood, and include some regularisation to stop the cuboids getting arbitrarily small). You could also include a prior over the parameters to get the MAP result.
You may run into optimisation problems due to the discontinuities at the edge of each cuboid. You're going to need to model outliers somehow too,if they occur. You're also going to have to find a way to select the number of clusters, $K$. And actual performance is going to be heavily dependant on how you initialise your parameters. But these are problems to think about once you've decided if this is actually a route you want to go down (as opposed to something more out of the box), and have checked how much of a problem they cause in practise.
The only reference I can find right now to prior work (though I'm sure people have done this many times before): 
Parameter estimation for finite mixtures of uniform distributions
Craigmile, PF (Craigmile, PF); Titterington, DM (Titterington, DM)
Source: COMMUNICATIONS IN STATISTICS-THEORY AND METHODS  Volume: 26   Issue: 8   Pages: 1981-1995   DOI: 10.1080/03610929708832026   Published: 1997
Also that took way more time to explain than I anticipated, which makes me wonder if writing the code might be too big an investment of your time too if it's all new to you :P. Try R-trees first, as other people have suggested.
