text:
Determine the Dimensions of a Rotated Rectangle
This Stack Snippet draws an aliased white rectangle on a black background given parameters for its dimensions, position, angle, and the grid dimensions:

grid w:<input id='gw' type='text' value='60'> grid h:<input id='gh' type='text' value='34'> w:<input id='w' type='text' value='40'> h:<input id='h' type='text' value='24'> x:<input id='x' type='text' value='0'> y:<input id='y' type='text' value='0'> θ:<input id='t' type='text' value='12'>° <button type='button' onclick='go()'>Go</button>
Image
<canvas id='c'>Canvas not supported</canvas>
Text
<textarea id='o' rows='36' cols='128'></textarea>

(JSFiddle Version)
The text representation has XX wherever there is a black pixel in the image and .. wherever there is a white pixel. (It looks squished if they are X and ..)
Write a program that takes the text representation of a rectangle produced by the Snippet and outputs the approximate width and height of the rectangle, both to within ± 7% of the actual width and height.
Your program should work for effectively all possible rectangles that can be drawn by the snippet, with the constraints that:

The rectangle width and height are 24 at minimum.
The grid width and height are 26 at minimum.
The rectangle never touches nor goes out of the grid bounds.

So the input rectangle may have any rotation, position, and dimensions, and the grid may have any dimensions, as long as the three constraints above are met. Note that except for the grid dimensions, the Snippet parameters can be floats.
Details

Take the raw text rectangle as input or take the filename of a file that contains the raw text rectangle (via stdin or command line). You may assume the text rectangle has a trailing newline.
You may assume the text rectangle is made from any two distinct printable ASCII characters other than X and . if desired. (Newlines must stay newlines.)
Output the measured width and height as integers or floats to stdout in any order (since there's no way to determine which one actually went with which parameter). Any format that clearly shows the two dimensions is fine, e.g. D1 D2, D1,D2, D1\nD2, (D1, D2), etc.
Instead of a program, you may write a function that takes the text rectangle as a string or the filename to it and prints the result normally or returns it as a string or list/tuple with two elements.
Remember that XX or .. is one 'pixel' of the rectangle, not two.

Examples
Ex. 1
Parameters: grid w:60 grid h:34 w:40 h:24 x:0 y:0 θ:12 (Snippet defaults)
Input
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX....XXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX............XXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX........................XXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX..................................XXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX............................................XXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX....................................................XXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX..............................................................XXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXX..........................................................................XXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXX................................................................................XXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXX................................................................................XXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXX..................................................................................XXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXX..........................................................................XXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXX..............................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXX....................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXX............................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXX..................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXX........................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXX............XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXX....XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

Example Outputs

40 24
24 40
[40.0, 24.0]
42.8, 25.68 (+7%)
37.2, 22.32 (-7%)

Ex. 2
Parameters: grid w:55 grid h:40 w:24.5 h:24 x:-10.1 y:2 θ:38.5
Input
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX..XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX......XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX............XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX..............XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXX..................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXX......................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXX............................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXX..............................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXX..................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXX......................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXX............................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXX................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXX..................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXX......................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XX............................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XX..............................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XX................................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXX..............................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXX..............................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXX............................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXX......................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXX....................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXX................................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXX............................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXX......................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXX..................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXX................................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXX..........................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXX......................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXX..................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXX................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXX..........XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXX......XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX..XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

Example Outputs

24.0 24.5
25.68 26.215 (+7%)
22.32 22.785 (-7%)

Scoring
The shortest code in bytes wins. Tiebreaker is highest voted post.

Shouldn't a solution meet the precision requirements to be accepted? The one you accepted is far off for certain input values.

Matlab, 226 bytes
The idea is simple: First I try to find out how much the rectangle was turned, then turn the image accordingly such that the rectangle is upright. Then I just 'sum' up all the pixels in the row columns sperately and try count how many of the sums are above than average (simple thresholding) for determining the width and height. This simple method works surprisingly reliably.
How can I detect the angle?
I just try each step (one degree each) and sum along the columns and get a vector of sums. When the rectangle is upright, I should ideally get only two sudden changes in this vector of sums. If the square is on the tip, the changes are very gradual. So I just use the first derivative, and try to minimize the number of 'jumps'. Here you can see a plot of the criterion we are trying to minimize. Notice that you can see the four minima which corresponds to the fours possible upright orientations.

Further thoughts: I am not sure how much it could be golfed as the exhausive angle search takes a lot of characters and I doubt you could achieve that so well with built in optimization methods, because as you can see there are a lot of local minima that we are not looking for. You can easily improve the accuracy (for big pictures) by choosing a smaller step size for the angle and only search 90° instead of 360° so you could replace 0:360 with 0:.1:90 or somehting like that. But anyways, for me the challenge was more finding a robust algorithm rather than golfing and I am sure the entries of the golfing languages will leave my submission far behind=)
PS: Someone should really derive a golfing language from Matlab/Octave.
Outputs
Example 1:
 25    39

Example 2:
 25    24

Code
Golfed:
s=input('');r=sum(s=='n');S=reshape(s',nnz(s)/r,r)';S=S(:,1:2:end-2)=='.';m=Inf;a=0;for d=0:360;v=sum(1-~diff(sum(imrotate(S,d))));if v<m;m=v;a=d;end;end;S=imrotate(S,a);x=sum(S);y=sum(S');disp([sum(x>mean(x)),sum(y>mean(y))])

Ungolfed:
s=input('');
r=sum(s=='n');
S=reshape(s',nnz(s)/r,r)';
S=S(:,1:2:end-2)=='.';
m=Inf;a=0;
for d=0:360;
    v=sum(1-~diff(sum(imrotate(S,d))));
    if v<m;
        m=v;a=d;
    end;
end;
S=imrotate(S,a);
x=sum(S);y=sum(S');
disp([sum(x>mean(x)),sum(y>mean(y))])

CJam, 68 65 64 bytes
This can be golfed a little more..
qN/2f%{{:QN*'.#Qz,)mdQ2$>2<".X"f#_~>@@0=?Qz}2*;@@-@@-mhSQWf%}2*;

How it works
The logic is pretty simple, if you think about it.
All we need from the input X. combinations is 3 coordinates of two adjacent sides. Here is how we get them:
First
In any orientation of the rectangle, the first . in the whole of the input is going to be one of the corners. For example..
XXXXXXXXXXXXXX
XXXXXXX...XXXX
XXXX.......XXX
X............X
XX.........XXX
XXXX...XXXXXXX
XXXXXXXXXXXXXX

Here, the first . is in the 2nd line, 8th  column.
But that's not it, we have to do some adjustment and add the width of the . run on that line to the coordinates to get the coordinate of the right end.
Second
If we transpose the above rectangle (pivoted on newlines), then the bottom left corner takes the above step's place. But here, we do not compensate for the . run length as we would have wished to get the bottom left coordinate of the edge anyways (which in transposed form will still be the first encountered .)
Rest two
For rest two coordinates, we simply flip horizontally, the rectangle and perform the above two steps. One of the corners here will be common from the first two.
After getting all 4, we simply do some simple math to get the distances.
Now this is not the most accurate method, but it works well within the error margin and well for all possible orientations of the rectangle.
Code expansion (bit outdated)
qN/2f%{{:QN*'.#Q0=,)md}:A~1$Q='.e=+QzA@@-@@-mhSQWf%}2*;
qN/2f%                               e# Read the input, split on newlines and squish it
      {   ...   }2*                  e# Run the code block two times, one for each side
{:QN*'.#Q0=,)md}:A~                  e# Store the code block in variable A and execute it
 :QN*                                e# Store the rows in Q variable and join by newlines
     '.#                             e# Get the location of the first '.'
        Q0=,)                        e# Get length + 1 of the first row
             md                      e# Take in X and Y and leave out X/Y and X%Y on stack
1$Q=                                 e# Get the row in which the first '.' appeared
    '.e=+                            e# Get number of '.' in that row and add it to X%Y
         QzA                         e# Transpose the rows and apply function A to get
                                     e# the second coordinate
            @@-@@-                   e# Subtract resp. x and y coordinates of the two corners
                  mh                 e# Calculate (diff_x**2 + diff_y**2)**0.5 to get 1 side
                    SQWF%            e# Put a space on stack and put the horizontally flipped
                                     e# version of the rows/rectangle all ready for next two
                                     e# coordinates and thus, the second side

Try it online here

Try a grid size of 50x50, rectangle size of 45x45, and angle -2. The error is about 28%. I tried a similar approach (it was my initial idea, before seeing yours), and getting it precise enough turns out to be trickier than expected, particularly if the sides are close to horizontal/vertical. Works great if they are closer to diagonal. I think this requires either more logic (e.g. also searching for extremes in the diagonal direction), or an entirely different approach.

@RetoKoradi Oh. That's just because all negative angles need the . width adjustment on the second coordinate, instead of first. Will fix. Should be short fix.

@RetoKoradi should be fixed now.

Try the 40x24 rectangle with angle 0.

@RetoKoradi Good points. Unaccepted for now.

It's been a while and still no fix. I'm accepting the next shortest  answer.

Python 3, 347 337 bytes
This turned out harder than I expected. Work in progress...
def f(s):
 l=s.split('\n');r=range;v=sorted;w=len(l[0]);h=len(l);p=[[x,y]for x in r(w)for y in r(h)if'X'>l[y][x]];x,y=[sum(k)/w/h for k in zip(*p)];g=[[x/2,y]];d=lambda a:((a[0]/2-a[2]/2)**2+(a[1]-a[3])**2)**.5
 for i in r(3):g+=v(p,key=lambda c:~-(c in g)*sum(d(j+c)for j in g))[:1]
 print(v(map(d,[g[1]+g[2],g[2]+g[3],g[1]+g[3]]))[:2])

Try it online! Comes with a test harness that tries 100 random valid rectangles. Note: The test harness seems to produce some errors and prints the corresponding grids. It seems like these cases all have rather round corners - an issue which seems to be caused by the grid generation code, which I translated from the original snippet.
Defines a function f taking the string as an argument and printing the result to STDOUT.
Pyth, 87 84 82 81 75 72 71 bytes (broken)
Km%2d.zJf<@@KeThTG*UhKUKPSm.adfqlT2ytu+G]ho*t}NGsm.a,kNGJ3]mccsklhKlKCJ

This seems to have been broken in 2015. Since then, something in Pyth has changed so that the old versions now crash and the new ones return [0, 0]. Right now I have neither the slightest idea what this monstrosity does nor the time to figure it out. I'm leaving it here as crossed-out for historical reasons.
Way Still too long. Basically a port of the previous. Loving Pyth's .a Euclidean distance. Takes input via STDIN and gives output via STDOUT. Expects the non-rectangle character to be lower-case x (well, anything with ASCII value 98 or more).
Algorithm
Both of these use the same algorithm. I basically start with an array containing the center of mass of the rectangle area. I then add three points to the array of all points in the rectangle, always choosing the one with maximum sum of distances to the points already in the array. The result is always three points in different corners of the rectangle. I then just calculate all the three distances between those three points and take the two shortest ones.

The Pyth solution doesn't work at all. The two examples from the OP give the results [33.0, 59.0] instead of  [40, 24]and [39.0, 54.0] instead of [24.0, 24.5].

@Jakube Weird. I'll investigate once I get home. Unfortunately I'm on a class trip to Lapland until June 9th.

I wouldn't call a trip to Lapland unfortunately ;-)

Have you come home yet?

@pppery Yes, lol. Edited answer.

Python 2, 342 bytes
import sys
r=[]
h=.0
for l in sys.stdin:w=len(l);r+=[[x*.5,h]for x in range(0,w,2)if l[x:x+2]=='..'];h+=1
x,y=.0,.0
for p in r:x+=p[0];y+=p[1]
n=len(r)
x/=n
y/=n
m=.0
for p in r:
 p[0]-=x;p[1]-=y;d=p[0]**2+p[1]**2
 if d>m:m=d;u,v=p
m=.0
for p in r:
 d=p[0]*v-p[1]*u
 if d>m:m=d;s,t=p
print ((u-s)**2+(v-t)**2)**.5+1,((u+s)**2+(v+t)**2)**.5+1

This took inspiration from @Pietu1998's algorithm. It takes the idea of determining one corner as the point farthest from the center, but differs from there:

I determine the second corner as the point with the largest cross product with the vector from center to first corner. This gives the point with the largest distance from the line from center to first corner.
There's no need to search for a third corner, since it's just the mirror image of the second corner relative to the center.

So the code follows this sequence:

First loop is over the lines in the input, and builds a list r of rectangle points.
Second loop calculates the average of all rectangle points, giving the center of the rectangle.
Third loop finds the point farthest from the center. This is the first corner. At the same time, it subtracts the center from the points in the list, so that the point coordinates are relative to the center for the remaining calculation.
Fourth loop finds the point with the largest cross product with the vector to the first corner. This is the second corner.
Prints out the distance between first corner and second corner, and distance between first corner and mirror image of second corner.
1.0 is added to the distances because the original distance calculations use pixel indices. For example, if you have 5 pixels, the difference between the index of the last and first pixel was only 4, which needs compensation in the final result.

Precision is quite good. For the two examples:
$ cat rect1.txt | python Golf.py
24.5372045919 39.8329756779
$ cat rect2.txt | python Golf.py
23.803508502 24.5095563412

Python 2, 272 bytes
Posting this as a separate answer since it's an entirely different algorithm than my previous one:
import sys,math
y,a,r=0,0,0
l,t=[1<<99]*2
for s in sys.stdin:
 c=s.count('..')
 if c:a+=c;x=s.find('.')/2;l=min(l,x);r=max(r,x+c);t=min(t,y);b=y+1
 y+=1
r-=l
b-=t
p=.0
w,h=r,b
while w*h>a:c=math.cos(p);s=math.sin(p);d=c*c-s*s;w=(r*c-b*s)/d;h=(b*c-r*s)/d;p+=.001
print w,h

This approach does not identify corners at all. It is based on the observation that the size (width and height) of the bounding box and the area of the rotated rectangle are sufficient to determine the width and height of the rectangle.
If you look at a sketch, it is fairly easy to calculate the width (wb) and height (hb) of the bounding box with w/h the size of the rectangle, and p the rotation angle:
wb = w * cos(p) + h * sin(p)
hb = w * sin(p) + h * cos(p)

wb and hb can be extracted directly from the image. We can also quickly extract the total area a of the rectangle by counting the number of .. pixels. Since we're dealing with a rectangle, this gives us the additional equation:
a = w * h

So we have 3 equations with 3 unknowns (w, h and p), which is enough to determine the unknowns. The only bummer is that the equations contain trigonometric functions, and at least with my patience and math skills the system cannot easily be solved analytically.
What I implemented is a brute force search for the angle p. Once p is given, the first two equations above become a system of two linear equations, which can be resolved for w and h:
w = (wb * cos(p) - hb * sin(p)) / (cos(p) * cos(p) - sin(p) * sin(p))
h = (hb * cos(p) - wb * sin(p)) / (cos(p) * cos(p) - sin(p) * sin(p))

With these values, we can then compare w * h with the measured area of the rectangle. The two values would ideally be equal at some point. This is of course not going to happen in floating point math.
The value of w * h decreases as the angle increases. So we start at angle 0.0, and then increment the angle by small steps until the first time w * h is less than the measured area.
The code only has two main steps:

Extract size of bounding box and rectangle area from input.
Loop over candidate angles until termination criterion is reached.

The precision of the output is good for rectangles where the width and height are significantly different. It gets somewhat iffy with rectangles that are almost square and rotated close to 45 degrees, just barely clearing the 7% error hurdle for test example 2.
The bitmap for example 2 actually looks slightly odd. The left corner looks suspiciously dull. If I add one more pixel at the left corner, it both looks better (to me), and gives a lot better precision for this algorithm.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Character counts in source code
Write a program that outputs a list of the number of occurrences of each unique character in its source code.
For example, this hypothetical program {Source_Print_1}; should produce this output:
; 1
P 1
S 1
_ 2
c 1
e 1
i 1
n 1
o 1
p 1
r 2
t 1
u 1
{ 1
} 1

The formatting should match this example. No extraneous whitespace is allowed, except an optional final newline.
Your program may not read its own source code from the source file.
The characters listed must be in one of two orders. Either the order of character values in the character encoding used by your language (probably ASCII), or the order the characters appear in your source.
This question inspired by this comment by Jan Dvorak.

Similar: http://codegolf.stackexchange.com/questions/34431/introspective-programming-code-that-analyzes-its-source-and-its-output

A zero-length program would work in quite a few languages.  Does this count as a standard loophole?

Let's go with... yes.

@DigitalTrauma: Added to the list.

Can the code contain newlines?

@user23013 good question. I did not consider newlines. I guess if you include them, I'd accept an answer that prints them out literally, so there would be one double-newline in the file somewhere.

CJam, 14 bytes
{S2N`/}`{S2N}/

Try it here.
Output is in the order they firstly appears:
{ 2
S 2
2 2
N 2
` 2
/ 2
} 2

It simply appends <SP>2<NL> to each character in {S2N`/}.

///, 12 bytes
  4
4 4

 4

A big thank you goes to @user23013, who suggested this improvement over my CJam code, outgolfing his own, highest-scoring answer in the process.
The characters are sorted by appearance. This code works in any language that just prints its own source code under the given circumstances (PHP, ASP, etc.).

CJam, 20 bytes
''S5N'5S5N'NS5N'SS5N

This approach doesn't use any built-in character counting.
Try it online in the CJam interpreter.
How it works
''S5N e# Push a single quote, a space, the integer 5 and a linefeed.
'5S5N e# Push the character 5, a space, the integer 5 and a linefeed.
'NS5N e# Push the character N, a space, the integer 5 and a linefeed.
'SS5N e# Push the character S, a space, the integer 5 and a linefeed.

+1 for not using standard quine techniques.

I really really hope this one remains tied for the lead. I'll happily give it the checkmark over its quine counterpart.

Now newlines are allowed. I think this answer should be better merged into yours.

@user23013: That's even shorter than your CJam answer. Thanks!

CJam, 20 bytes
{`"_~"+$e`{)S@N}%}_~

How it works
We first start off with one of the standard quine in CJam
{`"_~"}_~

which pushes the first block on stack, copies it, and runs the copy, which makes it print the source code itself finally.
Then we add the logic to compute the character count from the source code:
{`"_~"+                         e# At this point, we have the full source code with us
       $e`                      e# Sort to get similar characters together and run RLE to
                                e# get count of each character as [count char] array
          {    }%               e# Run each array element through this loop
           )S@N                 e# Pop the character, put a space, rotate the count after
                                e# space and then finally put a newline after the trio
                 }_~            e# Second half of the standard quine explained above

Try it online here

Python 3.5.0b1, 107 73 bytes
s="t='s=%r;exec(s)'%s;[print(c,t.count(c))for c in sorted({*t})]";exec(s)

Rather than the usual string replacement quine, which requires writing everything twice, here's an exec quine.

Mathematica, 101 bytes
Apply[Print[#1, " ", #2] &, Tally[Characters[StringJoin[ToString[#0, InputForm], "[];"]]], {1}] & [];

Unfortunately, I can't use any of the normal golfing tricks like removing whitespace, <> for StringJoin, # instead of #1, @ for prefix function calls or @@@ instead of Apply[...,{1}], because ToString[...,InputForm] thinks it has to pretty print everything...
This prints the characters in the order they first appear in the code. If I can assume that this isn't run in a REPL environment (which is rather unusual for Mathematica) I can save two bytes by omitting the two ;.

InputForm is annoying... OutputForm is better but it doesn't quote strings.

Haskell, 178 bytes
main=putStr(unlines[s:' ':show t|(s,t)<-zip" \"'(),-0123456789:<=S[\\]aehilmnoprstuwz|"[3,3,3,3,3,41,4,1,6,19,12,5,5,2,2,2,2,3,2,2,2,3,3,3,2,2,2,4,2,2,4,2,3,2,5,5,3,2,2,2]])--178

Nothing fancy. All characters of the program are in a literal list (String). So are the frequencies. Zip both lists and print. Output:
  3
" 3
' 3
( 3
) 3
, 41
- 4
0 1
1 6
2 19
3 12
4 5
5 5
6 2
7 2
8 2
9 2
: 3
< 2
= 2
S 2
[ 3
\ 3
] 3
a 2
e 2
h 2
i 4
l 2
m 2
n 4
o 2
p 3
r 2
s 5
t 5
u 3
w 2
z 2
| 2

Dart - 214 127
A direct version:
main(){print("  22\n\" 3\n( 3\n) 3\n1 3\n2 15\n3 8\n4 1\n5 2\n8 2\n; 2\n\\ 23\na 2\ni 3\nm 2\nn 23\np 2\nr 2\nt 2\n{ 2\n} 2");}

The "4" is just a fiddling factor to make the numbers add up.
See/run on DartPad.
Original: Standard quine tactic, and Dart's function names are a little too long for good golfing.
main({m,v,q:r'''main({m,v,q:r''}'')''{m={};for(v in q.split(''))m[v]=m[v]==null?2:m[v]+2;m.forEach((k,v)=>print("$k $v"));}'''}){m={};for(v in q.split(''))m[v]=m[v]==null?2:m[v]+2;m.forEach((k,v)=>print("$k $v"));}

See/run it on DartPad.

Vyxal D, 12 bytes
`:qpĊ⁋`:qpĊ⁋

Try it Online!

Haskell, 146 bytes
main=mapM putStrLn[a:" "++show s|a<-[' '..],s<-[sum[2|b<-show"main=mapM putStrLn[a: ++show s|a<-[' '..],s<-[sum[2|b<-show,a==b]],s>0]",a==b]],s>0]

Try it online!
Output:
  8
" 4
' 4
+ 4
, 6
- 6
. 4
0 2
2 2
: 2
< 6
= 6
> 2
L 2
M 2
S 2
[ 8
] 8
a 10
b 4
h 4
i 2
m 6
n 4
o 4
p 4
r 2
s 12
t 4
u 4
w 4
| 4

(Plus an additional newline)
Explanation:
The code is
main=mapM putStrLn[a:" "++show s|a<-[' '..],s<-[sum[2|b<-show"<code>",a==b]],s>0]

where "<code>" is a string of the program code without the ".
a goes through the ascii characters starting with a space. sum[2|b<-show"<code>",a==b] counts how often the character appears in the string, with each occurrence counted twice. a:" "++show s builds a string of the current character, a space and the character count. Finally mapM putStrLn prints each string in the list with a trailing newline.
The hardest part was getting the count of " right. Using just b<-"<code>" would count zero quotation marks because there are none in the string. Using show"<code>" adds a " to the front and end to the string, resulting in a count of four. So I had to put two additional quotation marks in the code, so instead of the (shorter) a:' ':show s I used a:" "++show s.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
Double-slit Quine
Objective:

Code two programs where each of the programs outputs both source codes
  interlaced per character like a zipper, a Double-slit Quine. The output from a Double-slit Quine
  is starting with the first character from the source code of the first
  program. If the source code of one of the programs is shorter in length
  than the other, then the rest of the output must be filled with the
  rest of the longer source code.

Rules:

You can use any programming language for both programs, not necessary the same programming language for both.
Your programs should not take any input from a file, file name, network, the other program or anything else.

Mandatory criteria:

There is a catch, somewhere in the output the word QUINE in capital letters must exist, uninterrupted.
You need to state what two programming languages you are using. If you are using the same programming language for both, then you only need to state one programming language.
Both programs should be able to execute, or be interpreted, respectively independent of the other program.

Example:
Having this example source code of program one:
"QIE"

Having this example source code of program two:
"UN"

Then the valid Double-slit Quine output from both programs must be:
""QUINE""

This is code-golf, the least bytes when summing the length of both source codes, and obviously also the length of each programs output, wins!

Each individual program doesn't have to be a quine also, correct? And must each program contain at least one character?

@mbomb007 That is correct. Each individual program should output both source code characters like a zipper from start to end. Well, I guess you'll need at least 1 byte in a programing language to output QUINE? The output from both programs must be identical.

I'm not sure I understand the rest of the output must be filled with the rest of the longer source code correctly. How would A and XYZ be interleaved? AXYZ?

@Dennis Just like .+ or z works

@Dennis Correct. What is left from the longer source code when the shorter runs out of bytes must be appended to the output from both programs.

Must the two programs have a Levenshtein distance of at least 1?

@quintopia I don't think this task is possible with two identical programs, since the output must contain the word QUINE.

@Dennis I'm thinking you are right.

CJam, 49 47 bytes
{`:_"__~~e#QUINE"}_~

and
{`:_"__~~e#QUINE"}_~e#QUINE

both print
{{``""__~~""++::__""ee##QQUUIINNEE""}}__~~e#QUINE

Try it online: program 1, program 2, proof of validity
How they work
{                }   e# Push a code block.
                  _~ e# Execute a copy.
 `                   e# Convert the code block into a string.
  :_                 e# Duplicate each character of the resulting string.
    "__~~e#QUINE"    e# Push this string.

Program 1 finishes here and has a string representation of its entire source code (which each character repeated twice) as well as the string e#QUINE on the stack.
Program 2 additionally parses e#QUINE, which is a comment.
In both cases, CJam prints the two strings automatically, resulting in the aforementioned output.

Seriously, 46 56 42 bytes
QUNX;RZtdXεj.ƒ£REIQ

Hex Dump:
51554e583b525a746458ee6a2e7f9f9c524549510a

The second program is this exact program reversed. It does contain an invisible character, but the byte count is correct. Both programs output this palindromic string:

QQUINEXR;£RƒZ⌂t.djXεεXjd.t⌂ZƒR£;RXENIUQQ

(For some reason it displays the invisible character when it outputs it in my terminal. I don't really understand that 7F byte.)
How it works:
Q                     Push the source to the stack
 U                    Nop
  NX                  Push and pop lyrics to 99 bottles of beer
    ;R                Make a copy of the source reversed.
      Z               Zip them together.
       t              Completely flatten the list.
        dX            Delete the trailing newline
          εj          Join it into a string.
            .         Print it with a trailing newline.
             ⌂        (This is that 7F byte that won't show up. It halts.)

The rest of the program is not executed.
Then the other direction:
\n                   Nop
  Q                  Push the source code.
   IE                Nops
     R               Reverse source code to get source of first program
      £              Eval it into a function
       ƒ             Call it.

At this point, we are now running the first program above, so execution continues as described there.
I am heavily exploiting the fact that several commands do nothing to strings or empty stacks here. I'm also exploiting the unexpected behavior of t when there is only one item on the stack. Do not expect this program to work in future versions of Seriously.
As the online version of the Seriously interpreter continues to be broken, you'll have to download the interpreter and run it on your own machine to test it.

@Dennis did you test it on your computer? That may be an artifact of my own build environment.

@Dennis nvm you're right. All of Seriously's output methods add a linefeed. I can fix it at the cost of 8 bytes. Looks like your method will be shorter (which is sad because this method is cooler IMO).

GolfScript, 46 bytes
{`{.}%"..~~QUINE"}.~

and
{`{.}%"..~~QUINE"}.~QUINE

both print
{{``{{..}}%%""....~~~~QQUUIINNEE""}}..~~QUINE

Try it online: program 1, program 2, proof of validity
How it works
{`{.}%"..~~QUINE"}.~

{                }    # Push a code block.
                  .~  # Execute a copy.
 `                    # Convert the code block into a string.
  {.}%                # Duplicate each character of the resulting string.
      "..~~QUINE"     # Push that string.

The source code of program 2 additionally parses QUINE\n, which are two undefined tokens.

Perl, 61+60=121 bytes
Program 1:
$x=q<print"\$x=q<$x>;eval\$x#"=~s/./$&$&/gr,QUINE>;eval$x#QIE

Program 2:
$x=q<print"\$x=q<$x>;eval\$x#"=~s/./$&$&/gr,QUINE>;eval$x#UN

I thought I'd give this a go in a non-golfing language. This is basically just a universal quine constructor in Perl modified to double every character before printing it, and append QUINE to the end. Then we just have to stick a comment at the end of the code to compensate for the added text.
(I wrote this without really looking at the other answers. Turns out it'd be possible to save a byte by putting the entire comment in one program, but I'm not sure if I should just be blatantly copying algorithms like that.)

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Improving a Rubik's Cube scrambler's byte count
Okay, I recently wrote a Javascript Rubik's Cube scrambler that outputs a 25-move scramble. It's 135 bytes. The Gist is here, but for reference, I have the formatted code:
function(){
  for(a=s=y=r=[],m='RLUDFB',x=Math.random;++a<26;y=r,s.push(r+["'",2,''][0|x()*2]))
  for(;r==y;r=m[0|x()*5]);
  return s.join(' ')
}

Any tips?
EDIT: Code is updated, with @Optimizer's edits.

I am assuming that you need to run this in browser, so ES5 and below only ?

Some of the answers are likely to be complete rewrites rather than tips for lowering the byte count of your algorithm/implementation.

@Optimizer Yeah, you have to run in a browser console.

140 135 127 bytes
Still a big one, but here's for starters:
function(){for(a=s=y=r=[],x=Math.random;a<25;s[a++]=r+["'",2,''][0|x(y=r)*3])for(;r==y;r='RLUDFB'[0|x()*6]);return s.join(' ')}

Some of the few tips used:

Loose the var. There is no point in using var when your aim is to golf!
When a = [], a++ gives 0. We use this fact to remove a 0 and ,
Combine the a<25 and a++
You don't need to split a string. You can access its characters anyways via [i]
Get rid of {} in for loops

More to come.

Just wondering, why are the parentheses after the function not there?

@molarmanful Whoops


text:
Chinese checkerboard
The game of Chinese checkers is played on a board with spaces in the shape of a six-pointed star:

Image from Wikipedia
We can create an ASCII-art representation of this board, using . for empty spots and the letters GYORPB for the six colored starting locations:
            G
           G G
          G G G
         G G G G
B B B B . . . . . Y Y Y Y
 B B B . . . . . . Y Y Y
  B B . . . . . . . Y Y
   B . . . . . . . . Y
    . . . . . . . . .
   P . . . . . . . . O
  P P . . . . . . . O O
 P P P . . . . . . O O O
P P P P . . . . . O O O O
         R R R R
          R R R
           R R
            R

To make it more interesting, we can also change the size. We'll measure the size of a board by the side length of its triangular starting locations: the board above is size 4.
Since it's really a pain to type all that by hand, let's write a program (or function) to do it!
Details
Your code should take a positive integer representing the size of the board, via STDIN, ARGV, or function argument. Output the checkerboard pattern to STDOUT (you may alternately return it as a string if your submission is a function).
Output must either

have no trailing spaces at all, or
have exactly enough trailing spaces to fill out the pattern to a perfect rectangle of width 6 * N + 1.

Output may optionally have a trailing newline. No other extra (leading, trailing) whitespace is permitted.
Examples
Size 1:
   G
B . . Y
 . . .
P . . O
   R

Size 2:
      G
     G G
B B . . . Y Y
 B . . . . Y
  . . . . .
 P . . . . O
P P . . . O O
     R R
      R

Size 4:
            G
           G G
          G G G
         G G G G
B B B B . . . . . Y Y Y Y
 B B B . . . . . . Y Y Y
  B B . . . . . . . Y Y
   B . . . . . . . . Y
    . . . . . . . . .
   P . . . . . . . . O
  P P . . . . . . . O O
 P P P . . . . . . O O O
P P P P . . . . . O O O O
         R R R R
          R R R
           R R
            R

Scoring
This is code-golf: the shortest code in bytes wins.

Can the output have empty lines of spaces before and after?

I'm going to say no.

You mentioned trailing spaces, but what about leading spaces? Does the image need to be flush left, or can it have an equal amount of leading spaces on each line?

Flush left, as shown in the sample output.

Can there be spaces beyond the right edge but still forming a rectangle?

Where were you guys when this was in the sandbox?  ;)  No.  Question edited to clarify whitespace policy.

Ruby, 141 127
Returns a rectangular string
->n{(-2*n..2*n).map{|i|j=i.abs
k=j>n ?0:j
(([i>0??P:?B]*k+[j>n ?i>0??R:?G:?.]*(2*n+1-j)+[i>0??O:?Y]*k)*" ").center(6*n+1)}*$/}

Ungolfed in test program
f=->n{
  (-2*n..2*n).map{|i|                    #Iterate rows from -2*n to 2*n
    j=i.abs                              #Absolute value of i
    k=j>n ?0:j                           #Value of j up to n: for PBYO
    (                                    #An array of characters forming one line
      ([i>0??P:?B]*k+                    #B or P * (k=j or 0 as appropriate)
       [j>n ?i>0??R:?G:?.]*(2*n+1-j)+    #R,G or . * (2*n+1-j) to form centre diamond
       [i>0??O:?Y]*k                     #O or Y * (k=j or 0 as appropriate)
      )*" "                              #Concatenate the array of characters into a string separated by spaces.
    ).center(6*n+1)                      #pad the string to the full width of the image, adding spaces as necessary.
  }*$/                                   #Concatenate the array of lines into a string separated by newlines.
}

puts f[gets.to_i]

Python 2, 152
n=input();x=N=2*n
while~N<x:s='';y=n*3;exec"a=x+y;q=[0,a>N,x-y>N,-x>n,-a>N,y-x>N,x>n,1];s+=' BYROPG.'[q.index(sum(q)<~a%2*3)];y-=1;"*(y-~y);print s;x-=1

This is, in retrospect, the wrong approach for Python, but I'm posting it here in case someone can make use of it. Rather than explaining this mess of code, I'll try to say the idea behind it.
The idea is to use triangular coordinates, in which the triangular lattice corresponds to integer triples (a,b,c) with a+b+c=0.

(Here, the lattice points are drawn as hexagons.)
We can convert Cartesian coordinates to triangular ones as
a = (x+y)/2
b = (x-y)/2
c = -x

noting that x and y must have the same parity, or otherwise it's off-checkerboard and we should print a space.
In triangular coordinates, the bounding lines of the six-sided star have equations: a==n, b==n, c==n, a==-n, b==-n, c==-n.
So, we can determine what region we're in by which of [a,b,c,-a,-b,-c] are greater than n.

If none are, we're in the center and print a dot.
If exactly one is, we're in one of the six outer triangles, and print the letter corresponding to the index.
If two or more are, we're outside the board, and print a space.

The bounding rectangle requires that we do this for x in the closed interval [-2*n,2*n] and y in the closed interval [-3*n,3*n].

The code does not work for me.

@BadAtGeometry It works for me.

What version are you using?

@BadAtGeometry TIO is using 2.7.15. What happens when you run it?

Python 2, 140 bytes
n=input()
for k in range(4*n+1):x=abs(k-2*n);y=2*n-x;p,q,r=" BP G..R YO "[(k-~k)/(n-~n)::4];print(" "*y+" ".join(p*x+q*-~y+r*x)+" "*y)[n:-n]

Not great, but here's my initial bid.
The whitespace rules added a lot of bytes. For comparison, here's a 120 byte Python 3 program which is only correct visually, and doesn't follow the whitespace rules:
def f(n):
 for k in range(4*n+1):x=abs(k-2*n);y=2*n-x;p,q,r=" BP G..R YO "[(k-~k)//(n-~n)::4];print(" "*y,*p*x+q*-~y+r*x)

And here's my slightly longer recursive 149 byte Python 3 attempt:
def f(n,k=0):x=2*n-k;s=" ".join(["B"*x+"."*-~k+"Y"*x,"G"*-~k][k<n]).center(6*n+1);print(s);k<n*2and[f(n,k+1),print(s.translate({71:82,66:80,89:79}))]

Retina, 234 bytes
.
P
.+
iP$0$0x$0j$0x$0Px$0kqw
P(?=P*xP*j)
s
P(?=P*j)
R
P(?=P*xP*k)
c
P(?=P*k)
O
x

+`i(s+R+)R
is$1#$1R
+`(s*)P(P*c*)(O*)O(?=k)
$0#s$1$2c$3
j|k
#
s

+`([^#]+#)q(.*)
q$1$2$1
R(?=.*w)
G
P(?=.*w)
B
O(?=.*w)
Y
w[^#]*#|q|i

\w
$0
c
.
 #
#

Takes input in unary.
Each line should go to its own file and # should be changed to newline in the file. This is impractical but you can run the code as is as one file with the -s flag, keeping the # markers and maybe changing them to newlines in the output for readability if you wish.
The code has minimal regex-complexity. The main steps in the generation are the followings:

Create the last G line and the first B.Y line (delimited by markers ijk and actual used letetrs are RPO).
Duplicate the topmost G line with a plus space, minus a G until there is only one G.
Duplicate the bottom B.Y line with a plus space and dot, minus a B and Y until there are no B and Y left.
Copy all the lines in reverse order after the current string (with the help of the marker q). We keep a marker (w) in the middle.
We change the letters RPO to GBY if they are before the marker.
Add the missing in-between spaces.

The results after each of the above points (delimited by ='s) for the input 1111 (unary 4):
1111
==============================
isssssssssRRRRjPPPPcccccOOOOkqw
==============================
issssssssssssR
sssssssssssRR
ssssssssssRRR
sssssssssRRRRjPPPPcccccOOOOkqw
==============================
issssssssssssR
sssssssssssRR
ssssssssssRRR
sssssssssRRRRjPPPPcccccOOOO
sPPPccccccOOO
ssPPcccccccOO
sssPccccccccO
ssssccccccccckqw
==============================
qi            R
           RR
          RRR
         RRRR
PPPPcccccOOOO
 PPPccccccOOO
  PPcccccccOO
   PccccccccO
    ccccccccc
w    ccccccccc
   PccccccccO
  PPcccccccOO
 PPPccccccOOO
PPPPcccccOOOO
         RRRR
          RRR
           RR
i            R
==============================
qi            G
           GG
          GGG
         GGGG
BBBBcccccYYYY
 BBBccccccYYY
  BBcccccccYY
   BccccccccY
    ccccccccc
w    ccccccccc
   PccccccccO
  PPcccccccOO
 PPPccccccOOO
PPPPcccccOOOO
         RRRR
          RRR
           RR
i            R
==============================
            G
           G G
          G G G
         G G G G
B B B B . . . . . Y Y Y Y
 B B B . . . . . . Y Y Y
  B B . . . . . . . Y Y
   B . . . . . . . . Y
    . . . . . . . . .
   P . . . . . . . . O
  P P . . . . . . . O O
 P P P . . . . . . O O O
P P P P . . . . . O O O O
         R R R R
          R R R
           R R
            R

JavaScript (ES6) 228
Construction line by line. Incredibly long compared to @Sp3000 that does the same.
Using template string to save 3 more bytes for newlines. All newlines are significant and counted.

f=w=>(i=>{r=(n,s=b=' ')=>s.repeat(n),l=c=>(c='GBYPOR'[c])+r(i,b+c),t=n=>r(w*3-i)+l(n)+`
`,s=n=>r(w-1-i)+l(n)+b+r(w+w-i,'. ')+l(n+1)+`
`;for(o='',q=r(w)+r(w+w,'. ')+`.
`;++i<w;o+=t(0))q+=s(3);for(;i--;o+=s(1))q+=t(5)})(-1)||o+q

// LESS GOLFED

u=w=>{
  r =(n,s=b=' ') => s.repeat(n),
  l = c => (c='GBYPOR'[c])+r(i, b+c),
  t = n => r(w*3-i) + l(n) + '\n',
  s = n => r(w-1-i) + l(n) + b + r(w+w-i,'. ') + l(n+1) + '\n',
  o = '',
  q = r(w) + r(w+w,'. ') + '.\n';
  for(i=0; i<w; i++)
    o += t(0), q += s(3);
  for(;i--;)
    o += s(1), q += t(5);
  return o+q
}

go=()=> O.innerHTML=f(I.value|0)

go()
<input id=I value=5><button onclick='go()'>-></button>


text:
Map of Islands (and a river)
Introduction
For many centuries, there has been a certain river that has never been mapped. The Guild of Cartographers want to produce a map of the river, however, they have never managed to succeed -- for some reason, all the cartographers they have sent to map the river have been eaten by wild animals in the area. A different approach is required.
Input Description
The area is a rectangular grid of cells of length m and width n. The cell in the bottom left would be 0,0, and the cell in the top right would be m-1,n-1. m and n are provided in the input as a tuple m,n.
By using long distance geographical sounding techniques the location of islands around the river have been identified. The size of the islands (i.e. how many cells the island occupies) have also been identified but the shape has not. We supply this information in a tuple s,x,y where s is the size of the island, and x and y are the x and y positions of one particular cell of that island. Each tuple in the input is space separated, so here is an example input:
7,7 2,0,0 2,3,1 2,6,1 2,4,3 2,2,4 8,0,6 1,2,6 3,4,6

To illustrate more clearly, here are is the input on a graph:
 y 6|8 1 3
   5|
   4|  2
   3|    2
   2|
   1|   2  2
   0|2
     =======
     0123456
     x

Output Description
Output a map using ASCII characters to represent parts of the area. Each cell will either be # (land) or . (water). The map should follow these rules:

Definition. An island is a orthogonally contiguous group of land cells that is bounded entirely by river cells and/or the border of the area.
Definition. A river is an orthogonally contiguous group of water cells that is bounded entirely by land cells and/or the border of the area, and does not contain "lakes" (2x2 areas of water cells).
Rule. The map shall contain exactly one river.
Rule. Each numbered cell in the input shall be part of an island containing exactly s cells.
Rule. Every island in the map shall contain exactly one of the numbered cells in the input.
Rule. There exists a single unique map for every input.

Here is the output of the example input:
#.#.##.
#....#.
#.##...
##..##.
###....
...##.#
##....#

Here is another input and output.
Input:
5,5 3,0,1 1,4,1 2,0,4 2,2,4 2,4,4

Output:
#.#.#
#.#.#
.....
###.#
.....

Note: this is the same question as a Nurikabe solver.

Can we take input in any convenient format, or should we stick to the one in the question?

this is also problem 4 from the 2012 Dyalog competition

@ngn Since when is "post a cryptographic hash"... usual? (but I suppose it's allowed when a challenge explicitly allow it)

@user202729 "usual" in the sense that I've allowed that in previous bounties I've posted. Nobody ever took advantage of it, but I still think it's a good idea and should be allowed by default.

A river is defined as not having a "lakes", 2x2 regions, but are "lakes" allowed to exist on there own? I am presuming not.

@ngn will a straightforward backtracking search out from each number be too slow to solve 15x15 puzzles?

@TheMatt no, there can't be isolated 2x2 regions - "Rule. 3 The map shall contain exactly one river."

@Jonah I don't know but it should be easy to test

here's a bookmarklet for https://puzzle-nurikabe.com/ - it converts the current puzzle to a valid input for this challenge and shows it in red just below the board: javascript:(_=>{var t=Game.nurikabe().task,m=t.length,n=t[0].length,s=[m,n];for(var i=0;i<m;i++)for(var j=0;j<n;j++)if(t[i][j]>=0)s+=' '+[t[i][j],i,j];puzzleContainerDiv.insertAdjacentHTML('beforeend','
<tt style=color:red>'+s+'</tt>
')})();void(0)

Well, I dont think any golfing language would be able to run that under a min

@Chromium Having solved some puzzles by hand, I think it should be possible. Most of the time you can find a square that's obviously water (where "obviously" is efficiently computable) or obviously land, so there would be very little branching in a backtracking algorithim with the right optimisations.

C + PicoSAT, 2345 995 952 bytes
#include<picosat.h>
#define f(i,a)for(i=a;i;i--)
#define g(a)picosat_add(E,a)
#define b calloc(z+1,sizeof z)
#define e(a,q)if(a)A[q]^A[p]?l[q]++||(j[++k]=q):s[q]||(i[q]=p,u(q));
z,F,v,k,n,h,p,q,r,C,*x,*A,*i,*l,*s,*j,*m;u(p){s[m[++n]=p]=1;e(p%F-1,p-1)e(p%F,p+1)e(p>F,p-F)e(p<=F*v-F,p+F)}t(){f(q,k)l[j[q]]=0;f(q,n)s[m[q]]=0;k=n=0;i[p]=-1;u(p);}main(){void*E=picosat_init();if(scanf("%d,%d",&F,&v)-2)abort();z=F*v;for(x=b;scanf("%d,%d,%d",&r,&p,&q)==3;g(p),g(0))x[p=F-p+q*F]=r;f(p,F*v-F)if(p%F)g(p),g(p+1),g(p+F),g(p+F+1),g(0);for(A=b,i=b,l=b,s=b,j=b,m=b;!C;){picosat_sat(E,C=h=-1);f(p,F*v)A[p]=picosat_deref(E,p)>0,i[p]=0;f(p,F*v)if(x[p])if(i[q=p]){for(g(-q);i[q]+1;)q=i[q],g(-q);g(C=0);}else if(t(),r=n-x[p]){f(q,r<0?k:n)g(r<0?j[q]:-m[q]);g(C=0);}f(p,F*v)if(!i[p])if(t(),A[p]){g(-++z);f(q,k)g(j[q]);g(C=0);f(q,n)g(-m[q]),g(z),g(0);}else{C&=h++;f(q,k)g(-j[q]);g(++z);g(++z);g(0);f(q,F*v)g(s[q]-z),g(q),g(0);}}f(p,F*v)putchar(A[p]?35:46),p%F-1||puts("");}

Try it online!
(Warning: this TIO link is a 30 kilobyte URL that contains a minified copy of PicoSAT 965, so you might not be able to load it in some browsers, but it loads in at least Firefox and Chrome.)
How it works
We initialize the SAT solver with a variable for each cell (land or water), and only the following constraints:

Every numbered cell is land.
Every 2×2 rectangle has at least one land.

The rest of the constraints are difficult to encode directly into SAT, so instead we run the solver to get a model, run a sequence of depth-first searches to find the connected regions of this model, and add additional constraints as follows:

For every numbered cell in a land region that’s too big, add a constraint that there should be at least one water cell among the current land cells in that region.
For every numbered cell in a land region that’s too small, add a constraint that there should be at least one land cell among the current water cells bordering that region.
For every numbered cell in the same land region as another numbered cell, add a constraint that there should be at least one water cell along the path of current land cells between them (found by walking the parent pointers left over from the depth-first search).
For every land region including no numbered cells, add constraints that either

all of those current land cells should be water, or
at least one of the current water cells bordering that region should be land.

For every water region, add constraints that either

all of those current water cells should be land, or
every cell other than those current water cells should be land, or
at least one of the current land cells bordering that region should be water.

Taking advantage of the incremental interface to the PicoSAT library, we can immediately rerun the solver including the added constraints, preserving all the previous inferences made by the solver.  PicoSAT gives us a new model, and we continue iterating the above steps until the solution is valid.
This is remarkably effective; it solves 15×15 and 20×20 instances in a tiny fraction of a second.
(Thanks to Lopsy for suggesting this idea of interactively constraining connected regions in an incremental SAT solver, a while back.)
Some larger test cases from puzzle-nurikabe.com
A random page of 15×15 hard puzzles (5057541, 5122197, 5383030, 6275294, 6646970, 6944232):
15,15 1,5,1 3,9,1 5,4,2 1,6,2 2,11,2 2,2,3 3,9,3 2,4,4 1,10,4 5,12,4 3,1,5 1,3,5 3,8,5 1,13,5 5,5,6 1,12,6 1,2,8 2,9,8 1,1,9 2,6,9 6,11,9 3,13,9 5,2,10 2,4,10 4,10,10 1,5,11 2,12,11 2,3,12 2,8,12 5,10,12 1,5,13 1,9,13 1,6,14 1,8,14
15,15 4,2,0 2,5,0 1,3,1 2,14,2 1,3,3 2,11,3 1,13,3 1,5,4 11,7,4 1,9,4 1,4,5 1,8,5 2,10,5 12,14,5 3,5,6 1,4,7 2,10,7 3,9,8 4,0,9 1,4,9 1,6,9 3,10,9 1,5,10 1,7,10 8,9,10 1,1,11 10,3,11 2,11,11 6,0,12 1,11,13 2,9,14 1,12,14
15,15 2,2,0 8,10,0 2,3,1 2,14,2 2,3,3 3,5,3 3,9,3 2,11,3 5,13,3 6,0,4 3,7,4 3,3,5 2,11,5 2,6,6 1,8,6 1,4,7 2,10,7 1,6,8 2,8,8 5,3,9 2,11,9 2,7,10 7,14,10 2,1,11 4,3,11 2,5,11 1,9,11 2,11,11 2,0,12 4,6,13 1,11,13 3,4,14 1,12,14
15,15 2,0,0 2,4,0 3,6,1 2,10,1 1,13,1 2,5,2 2,12,2 3,0,3 2,2,3 4,7,3 2,9,3 1,14,3 1,4,4 1,8,4 2,12,5 4,2,6 3,4,6 1,14,6 7,7,7 1,10,8 2,12,8 3,2,9 2,14,9 2,0,10 2,6,10 1,10,10 2,5,11 4,7,11 2,12,11 1,14,11 3,2,12 3,9,12 1,1,13 2,4,13 3,8,13 2,10,14 5,14,14
15,15 1,3,0 1,14,0 3,7,1 3,10,1 2,13,1 3,1,2 4,5,2 2,12,3 3,3,4 1,8,4 1,1,5 3,5,5 1,9,5 5,13,5 3,3,6 1,8,6 2,2,7 2,12,7 1,6,8 1,8,8 2,11,8 2,1,9 4,5,9 2,9,9 2,13,9 2,6,10 4,11,10 1,2,11 3,9,12 2,13,12 3,1,13 2,4,13 3,7,13 1,0,14
15,15 2,8,0 2,4,1 2,7,1 1,10,1 6,4,3 1,1,4 12,5,4 3,11,4 5,13,4 3,10,5 3,0,6 1,6,6 2,8,6 4,13,7 2,3,8 1,6,8 3,8,8 2,14,8 2,4,9 5,1,10 4,3,10 1,9,10 6,13,10 3,8,11 1,10,11 3,4,13 2,7,13 3,10,13 1,6,14 1,14,14

A random page of 20×20 normal puzzles (536628, 3757659):
20,20 1,0,0 3,2,0 2,6,0 1,13,0 3,9,1 3,15,1 2,7,2 3,13,2 3,0,3 2,3,3 3,18,3 3,5,4 2,9,4 2,11,4 2,16,4 1,0,5 2,7,5 1,10,5 1,19,5 3,2,6 1,11,6 2,17,6 2,0,7 3,4,7 5,6,7 2,9,7 4,13,7 3,15,7 1,3,8 1,10,8 1,14,9 2,18,9 3,1,10 2,4,10 1,8,10 1,10,10 3,12,10 3,16,10 1,9,11 1,17,11 2,19,11 2,0,12 2,2,12 1,4,12 4,6,12 2,13,12 2,15,12 1,14,13 2,17,13 1,3,14 2,5,14 4,7,14 2,15,14 3,0,15 1,2,15 2,13,15 3,18,15 3,7,16 7,10,16 1,17,16 2,0,17 2,3,17 2,5,17 3,11,17 3,15,17 1,0,19 1,2,19 1,4,19 2,6,19 5,8,19 1,11,19 1,13,19 3,15,19 2,18,19
20,20 1,0,0 1,4,0 5,8,0 1,17,0 1,19,0 2,17,2 3,6,3 2,10,3 2,12,3 4,14,3 6,0,4 3,4,4 4,7,4 1,11,4 1,18,4 1,6,5 3,12,5 4,15,5 4,4,6 2,16,6 2,19,6 6,0,7 3,10,7 2,12,8 2,17,8 3,3,9 2,5,9 4,8,9 2,10,9 3,0,10 1,2,10 5,14,10 2,16,10 2,19,10 7,7,11 3,12,12 2,17,12 2,2,13 4,4,13 3,6,13 4,14,13 3,0,14 1,3,14 1,5,14 3,16,14 1,2,15 1,9,15 2,11,15 5,13,15 3,19,15 1,4,16 3,6,16 1,3,17 1,12,17 1,14,17 1,16,17 6,0,19 2,2,19 3,5,19 2,7,19 5,9,19 1,11,19 2,13,19 1,15,19 4,17,19

GLPK/GNU MathProg, 1286 1227 1086 998 947 bytes
The main idea is to encode the "contiguous islands" constraint by introducing a preserved flow variable that has a pre-specified source at the hint location. The decision variable is_island then plays the role of placeable sinks. By minimizing the total sum of this flow, the islands are forced to remain connected in the optimum. The other constraints rather directly implement the various rules. What puzzles me that I still seem to need island_fields_have_at_least_one_neighbor.
Performance of this formulation is not great, though. By directly eating all the complexity with a large amount of constraints, the solver takes close to 15 seconds for the 15 x 15 example.
Code (golfed)
param M;param N;param P;set X:=0..N-1;set Y:=0..M-1;set Z:=0..P-1;set V:=X cross Y;set E:=setof{(p,q)in V,(r,s)in V:((abs(p-r)=1)&&(q=s))||((p=r)&&(abs(q-s)=1))}(p,q,r,s);param I{X,Y,Z},default 0;param t{z in Z}:=sum{x in X,y in Y}I[x,y,z];var u{X,Y},binary;var v{X,Y,Z},integer;var w{E,Z}>=0;minimize o:sum{(p,q,r,s)in E,z in Z}w[p,q,r,s,z];a{(x,y)in V,z in Z}:I[x,y,z]+sum{(p,q,r,s)in E:x=r&&y=s}w[p,q,x,y,z]-sum{(p,q,r,s)in E:x=p&&y=q}w[x,y,r,s,z]=v[x,y,z];b{(x,y)in V,z in Z}:v[x,y,z]*99>=I[x,y,z];c{(x,y)in V}:u[x,y]+sum{z in Z}v[x,y,z]=1;d{z in Z}:sum{(x,y)in V}v[x,y,z]=t[z];e{(x,y)in V:x>0&&y>0}:3>=u[x,y]+u[x-1,y-1]+u[x-1,y]+u[x,y-1];f{(x,y)in V}:sum{(p,q,r,s)in E:x=r&&y=s}u[p,q]>=u[x,y];g{(x,y)in V,z in Z:t[z]>1}:sum{(p,q,r,s)in E:x=r&&y=s}v[p,q,z]>=v[x,y,z];h{(x,y)in V,z in Z}:sum{(p,q,r,s)in E,oz in Z:x=p&&y=q&&z!=oz}v[r,s,oz]<=4*P*(1-v[x,y,z]);solve;for{y in M-1..0 by -1}{for{x in X}printf if u[x,y]=1 then"."else"#";printf"\n";}

Code (ungolfed)
# SETS
param M > 0, integer; # height
param N > 0, integer; # width
param P > 0, integer; # no. of islands

set X := 0..N-1;  # set of x coords
set Y := 0..M-1;  # set of y coords
set Z := 0..P-1;  # set of I

set V := X cross Y;
set E within V cross V := setof{
  (sx, sy) in V, (tx, ty) in V :

  ((abs(sx - tx) = 1) and (sy = ty)) or
  ((sx = tx) and (abs(sy - ty) = 1))
}
  (sx, sy, tx, ty);

# PARAMETERS
param I{x in X, y in Y, z in Z}, integer, default 0;
param island_area{z in Z} := sum{x in X, y in Y} I[x, y, z];

# VARIABLES
var is_river{x in X, y in Y}, binary;
var is_island{x in X, y in Y, z in Z}, binary;
var flow{(sx, sy, tx, ty) in E, z in Z} >= 0;

# OBJECTIVE
minimize obj: sum{(sx, sy, tx, ty) in E, z in Z} flow[sx, sy, tx, ty, z];

# CONSTRAINTS
s.t. I_are_connected{(x, y) in V, z in Z}:

    I[x, y, z]
  - is_island[x, y, z]
  + sum{(sx, sy, tx, ty) in E: x = tx and y = ty} flow[sx, sy, x, y, z]
  - sum{(sx, sy, tx, ty) in E: x = sx and y = sy} flow[x, y, tx, ty, z]
  = 0;

s.t. island_contains_hint_location{(x, y) in V, z in Z}:

    is_island[x, y, z] >= if I[x, y, z] > 0 then 1 else 0;

s.t. each_square_is_river_or_island{(x, y) in V}:

    is_river[x, y] + sum{z in Z} is_island[x, y, z] = 1;

s.t. I_match_hint_area{z in Z}:

    sum{(x, y) in V} is_island[x, y, z] = island_area[z];

s.t. river_has_no_lakes{(x,y) in V: x > 0 and y > 0}:

  3 >= is_river[x, y] + is_river[x - 1, y - 1]
     + is_river[x - 1, y] + is_river[x, y - 1];

s.t. river_squares_have_at_least_one_neighbor{(x, y) in V}:

    sum{(sx, sy, tx, ty) in E: x = tx and y = ty} is_river[sx, sy]
 >= is_river[x, y];

s.t. island_fields_have_at_least_one_neighbor{(x, y) in V, z in Z: island_area[z] > 1}:

    sum{(sx, sy, tx, ty) in E: x = tx and y = ty} is_island[sx, sy, z]
 >= is_island[x, y, z];

s.t. I_are_separated_by_water{(x, y) in V, z in Z}:

    sum{(sx, sy, tx, ty) in E, oz in Z: x = sx and y = sy and z != oz}
        is_island[tx, ty, oz]
 <= 4 * P * (1 - is_island[x, y, z]);

solve;

for{y in M-1..0 by -1}
{
    for {x in X}
    {
        printf "%s", if is_river[x, y] = 1 then "." else "#";
    }
    printf "\n";
}

Problem data
5 x 5
data;
param M := 5;
param N := 5;
param P := 5;
param I:=
# x,y,z,area
  0,1,0,3
  4,1,1,1
  0,4,2,2
  2,4,3,2
  4,4,4,2;
end;

7 x 7
data;
param M := 7;
param N := 7;
param P := 8;
param I:=
# x,y,z,area
  0,0,0,2
  3,1,1,2
  6,1,2,2
  4,3,3,2
  2,4,4,2
  0,6,5,8
  2,6,6,1
  4,6,7,3;
end;

15 x 15
data;
param M := 15;
param N := 15;
param P := 34;
param I:=
# x,y,   z,area
5,  1,   0, 1
9,  1,   1, 3
4,  2,   2, 5
6,  2,   3, 1
11, 2,   4, 2
2,  3,   5, 2
9,  3,   6, 3
4,  4,   7, 2
10, 4,   8, 1
12, 4,   9, 5
1,  5,  10, 3
3,  5,  11, 1
8,  5,  12, 3
13, 5,  13, 1
5,  6,  14, 5
12, 6,  15, 1
2,  8,  16, 1
9,  8,  17, 2
1,  9,  18, 1
6,  9,  19, 2
11, 9,  20, 6
13, 9,  21, 3
2,  10, 22, 5
4,  10, 23, 2
10, 10, 24, 4
5,  11, 25, 1
12, 11, 26, 2
3,  12, 27, 2
8,  12, 28, 2
10, 12, 29, 5
5,  13, 30, 1
9,  13, 31, 1
6,  14, 32, 1
8,  14  33, 1;
end;

Usage
Have glpsol installed, model as one file (e.g. river.mod), input data in separate file(s) (e.g. 7x7.mod). Then:
glpsol -m river.mod -d 7x7.mod

This plus some patience will print the solution in the specified output format (together with "some" diagnostic output).

I think this answer should be considered to be competing as it is possible for other people to verify that it works. The differences in IO format should be covered by the assumption that any reasonable IO format should be accepted.

@fəˈnɛtɪk Agreed. It wasn’t eligible for ngn’s bounty that just ended, which specifically required an answer runnable on TIO, but that isn’t a requirement of the question itself.

Given that I haven't started golfing it, I will not consider it competing (yet). Once I'm sure that I have pruned all redundant constraints, I'll one-char all declarations.

Skipped the pruning, went straight to golfing.... after >1 year delay.

And now that I'm below 1000, I want to get below 952... that will need changes to the actual formulation to make it less redundant.

Ok, that was easier than expected.

Python 3, 1295 bytes
This is a python only solution. It uses no libraries and loads the board through standard input. Further explanation to come. This is my first attempt at such a large golf. There is a link to the commented and none-golfed code at the bottom.
L,S,O,R,F=len,set,None,range,frozenset
U,N,J,D,I=S.update,F.union,F.isdisjoint,F.difference,F.intersection
def r(n,a,c):
 U(c,P)
 if L(I(N(Q[n],C[n]),a))<2:return 1
 w=D(P,N(a,[n]));e=S();u=S([next(iter(w))])
 while u:n=I(Q[u.pop()],w);U(u,D(n,e));U(e,n)
 return L(e)==L(w)
def T(a,o,i,c,k):
 s,p,m=a
 for _ in o:
  t=s,p,N(m,[_]);e=D(o,[_])
  if t[2] in c:o=e;continue
  c.add(t[2]);n=D(Q[_],m);U(k,n)
  if not J(i,n)or not r(_,N(m,i),k):o=e
  elif s==L(t[2]):yield t
  else:yield from T(t,N(e,n),i,c,k)
s,*p=input().split()
X,Y=eval(s)
A=[]
l=1,-1,0,0
P=F((x,y)for y in R(Y)for x in R(X))
exec("Q%sl,l[::-1]%s;C%s(1,1,-1,-1),l[:2]*2%s"%(('={(x,y):F((x+i,y+j)for i,j in zip(',')if X>x+i>-1<y+j<Y)for x,y in P}')*2))
for a in p:a,x,y=eval(a);k=x,y;A+=[(a,k,F([k]))]
A.sort(reverse=1)
k=F(a[1]for a in A)
p=[O]*L([a for a in A if a[0]!=1])
g,h=p[:],p[:]
i=0
while 1:
 if g[i]is O:h[i]=S();f=O;g[i]=T(A[i],Q[A[i][1]],D(N(k,*p[:i]),[A[i][1]]),S(),h[i])
 try:p[i]=g[i].send(f)[2]
 except:
  f=I(N(k,*p[:i]),h[i]);g[i]=p[i]=O;i-=1
  while J(p[i],f):g[i]=p[i]=O;i-=1
 else:
  i+=1
  if i==L(p):
   z=N(k,*p)
   if not any(J(z,F(zip([x,x+1]*2,[y,y,y+1,y+1])))for x in R(X-1)for y in R(Y-1)):break
   for c in h:U(c,z)
b=[X*['.']for i in R(Y)]
for x,y in z:b[y][x]='#'
for l in b[::-1]:print(''.join(l))

Try it online!
Un-golfed code here.