[Question]
  [
Imagine a bunch of rectangles drawn in the plane, each rectangle with its vertices at integer coordinates and its sides parallel to the axes:
![enter image description here](https://i.stack.imgur.com/aSLrum.jpg)
The rectangles partition the plane into a number of disjoint regions, coloured red and blue below:
![enter image description here](https://i.stack.imgur.com/5mlACm.jpg)
Your goal is to find the number of such regions which are perfect squares. In the example above, there are three:
![enter image description here](https://i.stack.imgur.com/6RpmKm.jpg)
Note that the big square in the middle is not counted as it is not a single region, but is instead made up of several smaller disjoint regions.
## Input
You may write a function or a full program for this challenge.
Input will be `4n` nonnegative integers defining `n` rectangles in the plane. Each rectangle is represented by two opposing vertices, e.g. `4 9 7 8` represents the rectangle with opposing vertices `(4, 9)` and `(7, 8)`. Note that this rectangle could also be represented as `7 8 4 9` or `4 8 7 9`.
The exact input format is flexible (e.g. space-separated string, comma-separated string, single array of integers, list of coordinate tuples, and so on) but please [be reasonable](http://meta.codegolf.stackexchange.com/questions/1061/standard-loopholes-which-are-no-longer-funny) and give an example of how to run your code in your post. You may not reorder the input.
For simplicity, you can assume that no two edges will be overlapping — this includes overlapping at a vertex. In particular, this implies that no two rectangles will be touching edge-to-edge or corner-to-corner, and the rectangles will have nonzero area.
## Output
Your program should print or return a single integer, which is the number of square regions.
## Scoring
This is code golf, so the code in the fewest bytes wins.
---
## Test cases
**Input:**
```
0 0 5 5
6 8 10 4
14 16 11 13
19 1 18 2
```
**Output:**
```
4
```
This is simply four disjoint squares:
![enter image description here](https://i.stack.imgur.com/geLR3m.jpg)
---
**Input:**
```
2 1 3 11
1 10 5 19
6 10 11 3
8 8 15 15
13 13 9 5
15 1 19 7
17 19 19 17
```
**Output:**
```
3
```
This is the example test case at the start of the post.
---
**Input:**
```
0 9 15 12
6 3 18 15
9 6 12 20
13 4 17 8
```
**Output:**
```
7
```
![enter image description here](https://i.stack.imgur.com/PT1eXm.jpg)
---
**Input:**
```
5 9 11 10
5 12 11 13
6 8 7 14
9 8 10 14
13 8 14 9
13 10 14 14
```
**Output:**
```
14
```
![enter image description here](https://i.stack.imgur.com/qGYAIm.jpg)
---
**Input:**
```
0 99999 100000 0
```
**Output:**
```
0
```
This is just one big rectangle.
---
**Input:**
```
0 99999 100000 0
2 1 142857 285714
```
**Output:**
```
1
```
Two big rectangles which overlap.
  
[Answer]
# Python 2, ~~480 436 386~~ 352 bytes
```
exec u"""s=sorted;H=[];V=[]
FRIinput():
 S=2*map(s,zip(*R))
 FiI0,1,2,3:
    c=S[i][i/2];a,b=S[~i]
    FeIs(H):
     C,(A,B)=e
     if a<C<b&A<c<B:e[:]=C,(A,c);H+=[C,(c,B)],;V+=[c,(a,C)],;a=C
    V+=[c,(a,b)],;H,V=V,H
print sum(a==A==(d,D)&c==C==(b,B)&B-b==D-d&1-any(d<X[0]<D&b<y<B Fy,XIH)Fb,aIH FB,AIH Fd,cIV FD,CIV)""".translate({70:u"for ",73:u" in ",38:u" and "})
```
Takes a list of coordinate pairs through STDIN in the format:
```
[  [(x, y), (x, y)],  [(x, y), (x, y)],  ...  ]
```
and prints the result to STDOUT.
  
The actual program, after string replacement, is:
```
s=sorted;H=[];V=[]
for R in input():
 S=2*map(s,zip(*R))
 for i in 0,1,2,3:
    c=S[i][i/2];a,b=S[~i]
    for e in s(H):
     C,(A,B)=e
     if a<C<b and A<c<B:e[:]=C,(A,c);H+=[C,(c,B)],;V+=[c,(a,C)],;a=C
    V+=[c,(a,b)],;H,V=V,H
print sum(a==A==(d,D) and c==C==(b,B) and B-b==D-d and 1-any(d<X[0]<D and b<y<B for y,X in H)for b,a in H for B,A in H for d,c in V for D,C in V)
```
## Explanation
Instead of fiddling with complex polygons, this program deals with simple line segments.
For each input rectangle, we add each of its four edges to a collective segment list, individually.
Adding a segment to the list goes as follows:
we test each of the existing segments for intersection with the new segment;
if we find an intersection, we divide both segments at the point of intersection and continue.
To make things easier, we actually keep two separate segment lists: a horizonal one and a vertical one.
Since segments don't overlap, horizontal segments can only intersect vertical segments and vice versa.
Better yet, it means that all intersections (not considering the edges of the same rectangle) are "proper," i.e., we don't have T-shaped intersections, so "both sides" of each segment are truly divided.
Once we've constructed the segment list(s), we start counting squares.
For each combination of four segments (particularly, two horizontal segments and two vertical ones,) we test if they form a square.
Furthermore, we verify that no vertex lies within this square (which can happen if, for example, we have a small square inside a bigger one.)
This gives us the desired quantity. Note that even though the program tests each combination four times in different orders, the particular ordering of the segment coordinates guarantees that we count each square only once.
[Answer]
## SQL (POSTGIS), ~~286~~ ~~269~~ ~~261~~ ~~240~~ ~~226~~ ~~218~~ 216
This is a query for the PostGIS extension to PostgreSQL.
I haven't counted the input values in the total.
```
SELECT SUM(1)FROM(SELECT(ST_Dump(ST_Polygonize(g))).geom d FROM(SELECT ST_Union(ST_Boundary(ST_MakeEnvelope(a,b,c,d)))g FROM(VALUES
-- Coordinate input
(2, 1, 3, 11)
,(1, 10, 5, 19)
,(6, 10, 11, 3)
,(8, 8, 15, 15)
,(13, 13, 9, 5)
,(15, 1, 19, 7)
,(17, 19, 19, 17)
)i(a,b,c,d))i)a WHERE(ST_XMax(d)-ST_XMin(d))^2+(ST_YMax(d)-ST_YMin(d))^2=ST_Area(d)*2
```
**Explanation**
The query builds geometries for each coordinate pair. Unions the exterior rings to properly node the lines. Turns the results into polygons and tests ~~the width against height and~~ the area doubled against the sum of each side squared.
It will run as a standalone query on any PostgreSQL database with the PostGIS Extension.
**Edit** Found a couple more.
[Answer]
# Haskell, ~~276 266 250 237 225 222~~ 217 bytes
It keeps getting shorter... and more obfuscated.
```
(x#i)l=mapM id[[min x i..max x i-1],l]
(y!j)l f=and[p l==p(f[y,j])|p<-map elem$f[y..j]]
s[h,v]=sum[1|[x,j]<-h,[y,i]<-v,x<i,i-x==j-y,(y!j)h$x#i,(x!i)v$y#j]
n=s.foldr(\(x,y,i,j)->zipWith(++)[x#i$[y,j],y#j$[x,i]])[[],[]]
```
Evaluate `n [(0,0,5,5),(6,8,10,4),(14,16,11,13),(19,1,18,2)]` for the first test case. I think I'm getting close to the limits of golfing this algorithm on Haskell.
This function is so slow (at least *O(n3)* where *n* is the total perimeter of all rectangles in the input) that I cannot evaluate it on the last two test cases. When I compiled it with optimizations turned on and run it on the 400-times shrunk version `[(0,249,250,0),(2,1,357,714)]` of the last test, it finished in a little over 12 seconds. Based on this, the actual test case would finish in about 25 years.
### Explanation (partial, I will expand this when I have time)
We first build two lists `h` and `v` as follows. For each rectangle in the input, we split its border into segments of length 1. The west endpoints of horizontal segments are stored in `h`, and the south endpoints of vertical segments in `v`, as lists `[x,y]` of length 2. The coordinates in `v` are stored in a reversed form as `[y,x]` for golfing reasons. Then we just loop over both lists and search for horizontal edge `[x,j]` and vertical edge `[i,y]` such that `x < i` and `i-x == j-y` (so they are the northwest and southeast corners of a square), and check that the borders of the square are in the correct lists `h` and `v`, while the interior coordinates are not. The number of the positive instances of the search is the output.
]