text:
The Characters in the String Go Round and Round
(Inspired by an early draft of PhiNotPi's fractal line challenge.)
You're given a width W > 1, a height H > 1 and string consisting of 2(W+H-2) printable ASCII characters. The task is to print this string wrapped around a rectangle of the given width and height, starting in the top left corner, in a clockwise sense. The rectangle's inside is filled with spaces. The test cases should hopefully make this very clear.
You may write a program or function, taking input via STDIN (or closest alternative), command-line argument or function argument, and either printing the result to STDOUT (or closest alternative) or returning it as a string.
There must be no leading or trailing spaces (apart from those that might be in the input string). You may optionally output a single trailing newline.
This is code golf, so the shortest submission (in bytes) wins.
Test Cases
Each test case is "String" W H followed by the expected output.
"Hello, World! "
5 4
Hello
    ,
!
dlroW

"+--+|||+--+|||"
4 5
+--+
|  |
|  |
|  |
+--+

">v<^"
2 2
>v
^<

"rock beats scissors beats paper beats "
11 10
rock beats
          s
s         c
t         i
a         s
e         s
b         o
          r
r         s
epap staeb

Note that the following string contains an escaped '"'.
"!\"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~"
46 3
!"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMN
~                                            O
}|{zyxwvutsrqponmlkjihgfedcba`_^]\[ZYXWVUTSRQP

Leaderboards
Here is a Stack Snippet to generate both a regular leaderboard and an overview of winners by language.
To make sure that your answer shows up, please start your answer with a headline, using the following Markdown template:
# Language Name, N bytes

where N is the size of your submission. If you improve your score, you can keep old scores in the headline, by striking them through. For instance:
# Ruby, <s>104</s> <s>101</s> 96 bytes

function answersUrl(e){return"http://api.stackexchange.com/2.2/questions/"+QUESTION_ID+"/answers?page="+e+"&pagesize=100&order=desc&sort=creation&site=codegolf&filter="+ANSWER_FILTER}function getAnswers(){$.ajax({url:answersUrl(page++),method:"get",dataType:"jsonp",crossDomain:true,success:function(e){answers.push.apply(answers,e.items);if(e.has_more)getAnswers();else process()}})}function shouldHaveHeading(e){var t=false;var n=e.body_markdown.split("\n");try{t|=/^#/.test(e.body_markdown);t|=["-","="].indexOf(n[1][0])>-1;t&=LANGUAGE_REG.test(e.body_markdown)}catch(r){}return t}function shouldHaveScore(e){var t=false;try{t|=SIZE_REG.test(e.body_markdown.split("\n")[0])}catch(n){}return t}function getAuthorName(e){return e.owner.display_name}function process(){answers=answers.filter(shouldHaveScore).filter(shouldHaveHeading);answers.sort(function(e,t){var n=+(e.body_markdown.split("\n")[0].match(SIZE_REG)||[Infinity])[0],r=+(t.body_markdown.split("\n")[0].match(SIZE_REG)||[Infinity])[0];return n-r});var e={};var t=0,c=0,p=-1;answers.forEach(function(n){var r=n.body_markdown.split("\n")[0];var i=$("#answer-template").html();var s=r.match(NUMBER_REG)[0];var o=(r.match(SIZE_REG)||[0])[0];var u=r.match(LANGUAGE_REG)[1];var a=getAuthorName(n);t++;c=p==o?c:t;i=i.replace("{{PLACE}}",c+".").replace("{{NAME}}",a).replace("{{LANGUAGE}}",u).replace("{{SIZE}}",o).replace("{{LINK}}",n.share_link);i=$(i);p=o;$("#answers").append(i);e[u]=e[u]||{lang:u,user:a,size:o,link:n.share_link}});var n=[];for(var r in e)if(e.hasOwnProperty(r))n.push(e[r]);n.sort(function(e,t){if(e.lang>t.lang)return 1;if(e.lang<t.lang)return-1;return 0});for(var i=0;i<n.length;++i){var s=$("#language-template").html();var r=n[i];s=s.replace("{{LANGUAGE}}",r.lang).replace("{{NAME}}",r.user).replace("{{SIZE}}",r.size).replace("{{LINK}}",r.link);s=$(s);$("#languages").append(s)}}var QUESTION_ID=47710;var ANSWER_FILTER="!t)IWYnsLAZle2tQ3KqrVveCRJfxcRLe";var answers=[],page=1;getAnswers();var SIZE_REG=/\d+(?=[^\d&]*(?:<(?:s>[^&]*<\/s>|[^&]+>)[^\d&]*)*$)/;var NUMBER_REG=/\d+/;var LANGUAGE_REG=/^#*\s*([^,]+)/
body{text-align:left!important}#answer-list,#language-list{padding:10px;width:290px;float:left}table thead{font-weight:700}table td{padding:5px}
<link rel=stylesheet type=text/css href="//cdn.sstatic.net/codegolf/all.css?v=83c949450c8b">

Leaderboard

	AuthorLanguageSize

Winners by Language

LanguageUserScore

{{PLACE}}	{{NAME}}{{LANGUAGE}}{{SIZE}}Link

{{LANGUAGE}}{{NAME}}{{SIZE}}Link

Your leaderboards are so cool.

Did you golf your leaderboard script?

@mbomb007 No, I ran the code through minifiers, so it wouldn't take up much space when expanded. (I guess removing line breaks would have been sufficient.) I've still got an unminified version on my hard drive.

If you renamed this to "The Chars in the String Go Round and Round", it fits the rhythm of the song better.

CJam, 27 bytes
Nl~:L/(os\2-{)L2-S*@(N@}*W%

I don't really CJam, but I think this beats Martin. The main difference is that we push a newline before reading input and print the first line immediately, negating the need to store the height.
Takes input in the order
H "String" W

Try it online.

Python 2, 95 bytes
s,m,n=input()
print s[:n]
for i in range(m-2):print s[~i]+' '*(n-2)+s[n+i]
print s[1-m::-1][:n]

Prints the first line, then the two vertical lines, then the last line.
There's got to be something shorter than writing print three times, but everything I've tried so far with saving to a variable and '\n'.join has been longer.

You could switch to Python 3 and store print in a variable...

@Omar That winds up being longer because you have to use eval on the input and bracket the print statements.

Oh, I hadn't taken eval into account! The parenthesis shouldn't be too much of a problem since print in Python 2 needs a space after it. Going from print blah to p(blah) still saves 3 characters.

Pyth, 47 46 45 40 37 36 bytes
This is the obvious approach implemented in Pyth. It prints the first line by indexing 0:width and then the middle, then the end.
Thanks to @Jakube for the tip with using z and Q for two inputs and using p.
AkYQ<zkV-Y2p*d-k2@zt_N@z+kN;<_<z-2Yk

Takes input from stdin as a string and as a tuple of dimensions, newline separated:
Hello, World!
5, 4

and writes to stdout.
Try it here.
A              Double assignment
 kY            The vars k and Y
 Q             The dimension tuple
<zk            Prints out first line by doing z[:width]
V-Y2           For N in height-2
 p             Print out everything
  *d           Repeat " "
   -k2         Width-2 times
  @z           Index z
   -_N1        At index -N-1
  @z           Index z
   +kN         At index k+N
;              Close out loop
<_<z-2Yk       Print last line

Using z for reading the string saves quite a lot of chars. Also t_N is the same thing as -_N1.

37 chars are possible with our approach.

@Jakube thanks for the tips!

One more char saving. Instead of ++ use p and switch the zt_N with *d-k2.

CJam, 31 30 bytes
At Optimizer's insistence, here is my own attempt. I'm not a fan of winning my own challenges, so I'm counting the APL family (or someone better at CJam) to beat this. ;)
l~:H;:V/(N@s{)V2-S*@(N@}H2-*W%

Takes input in the same order as given in the question:
"Hello, World! " 5 4

Test it here.
One byte saved thanks to Optimizer.
Explanation
Originally, I had a really nice idea for starting with the rectangle of spaces and then literally wrapping the string around it while rotating the entire grid four times. However, I couldn't seem to get that to work in case where width or height or both are 2. So I tried the naive approach (print top, loop over sides, print bottom), and surprisingly it turned out to be really short.
l~                             "Read and evaluate the input.";
  :H;                          "Store the height in H and discard it.";
     :V/                       "Store the width in V and split the input into chunks of size V.";
        (N                     "Slice off the first such chunk and push a newline.";
          @s                   "Pull up the other chunks and join them back together.";
            {          }H2-*   "Repeat this block H-2 times, printing the sides.";
             )                 "Slice off the last character of the string.";
              V2-S*            "Push V-2 spaces.";
                   @(          "Pull up the remaining string and slice off the first character.";
                     N@        "Push a newline and pull up the remaining string.";
                            W% "Reverse the remainder of the string, which is the bottom row.";

Because we can get the length of the string and we have V, there is no need to save H. Just repeat the block until only V characters remain. l~;:V/(N@s{)V2-S*@(N@_,V-}gW% saves 1 char.

@DocMax Unfortunately, that doesn't work for height 2. It's a good idea though, I'll look into whether I can make use of it in a different way somehow.

D'oh! You even mentioned to H=2 problem and I still forgot to guard against it.

J, 61 bytes
Method:
Starting from a (height-2)*(width-2) block of spaces we take the necessary amount of characters from the end of the string and add it to the current block. We repeat this 4 times. The total 5 states illustrated with the 'Hello, World! ' 5 4 example (spaces replaced with Xs for readability):
XXX   !_   orld   ,_W   Hello
XXX   XX   XXX!   XXo   _XXX,
      XX   XXX_   XXr   !XXX_
      XX          XXl   dlroW
                  _!d

The code:
4 :'|:>{:((}.~{:@$);({.~{:@$)|.@|:@,])&>/^:4(|.x);'' ''$~y-2'

Explicit function definition. The two-operand function takes a string as left argument and a list of two integers as right argument.
Example usage:
   wrap_on=.4 :'|:>{:((}.~{:@$);({.~{:@$)|.@|:@,])&>/^:4(|.x);'' ''$~y-2'

   'Hello, World! ' wrap_on 5 4
Hello
    ,
!
dlroW

   '>v<^' wrap_on 2 2
>v
^<

Try it online here.

Wow, I'm impressed that this works for width and height 2 in J.

Pyth, 38 37
AGHQ<zGFNC,_>z_ttH>zGj*dttGN)<>_zttHG

I originally had a different 38 solution, but it was basically golfed solution of Maltysen's answer. So I decided to go a little bit different.
Try it online.
              implicit: z=string from input, Q=pair of numbers from input
AGHQ          G=Q[0] (width), H=Q[1] (height)
<zG           print z[:G]
    _>z_ttH     last H-2 chars reversed
    >zG         all chars from the Gth position to end
  C,           zip these 2 strings to pairs
FN            for each pair N:
  j*dttGN       seperate the two chars by (G-2) spaces and print
)             end for
<>_zttHG     print last line z[::-1][H-2:][:G]

_>z_ttH is equivalent to <_zttH.

@isaacg Thanks, already seen something similar in Maltysen's answer.

JavaScript (ES6), 110 115
Function with 3 parameters, returning a string
F=(s,w,h,q=h+(w-=2),t=b='')=>
[for(c of s)q?t+=q<h?c+'\n'+s[w+h+w+q--]+' '.repeat(q&&w):(b+=s[w+q--],c):q]
&&t+b

Chrome version 119: no short format for functions, no default parameters. No reason to use for(of) even if it is supported
function F(s,w,h){
  for(q=h+(w-=2),t=b=i='';
      q;
      q<h?t+='\n'+s[w+h+w+q--]+' '.repeat(q&&w):b+=s[w+q--])
    t+=s[i++];
  return t+b
}

ES5 version 126: no for(of), no string.repeat
function F(s,w,h){
  for(q=h+(w-=2),t=b=i='';
      q;
      q<h?t+='\n'+s[w+h+w+q--]+Array(q&&-~w).join(' '):b+=s[w+q--])
    t+=s[i++];
  return t+b
}

Ungolfed
F=(s,w,h)=>
{
  var q = h+(w-=2), // middle length
      t = '', // top and body
      b = ''; // bottom row
  for(c of s)
    if (q > 0)
    {
      if (q < h)
      {
        t += c+'\n'; // right side, straight
        t += s[w+h+w+q]; // left side, backwards
        if (q > 1) // body fill, except for the last line
          t += ' '.repeat(w)
      }
      else
      {
        t+=c, // top, straight
        b+=s[w+q] // bottom, backwards
      }
      --q
    }
  return t+b

Test In Firefox/FireBug console
;[["Hello, World! ", 5, 4],["+--+|||+--+|||",4,5],[">v<^",2,2]
,["rock beats scissors beats paper beats ",11,10]
,["!\"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_`abcdefghijklmnopqrstuvwxyz{|}~",46,3]]
.forEach(test => console.log(F(...test)))

Output
Hello
    ,
!
dlroW

+--+
|  |
|  |
|  |
+--+

>v
^<

rock beats
          s
s         c
t         i
a         s
e         s
b         o
          r
r         s
epap staeb

!"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMN
~                                            O
}|{zyxwvutsrqponmlkjihgfedcba`_^]\[ZYXWVUTSRQP

Python 2, 97 bytes
def f(s,w,h):print s[:w];i=0;exec'print s[~i]+" "*(w-2)+s[w+i];i+=1;'*(h-2);print s[1-h:w+h-3:-1]

Taking the direct approach.

Haskell, 164 156 bytes
import Data.List
r=replicate
p w h(s:t)=unlines$fst$n$n$n$n(r h$r w ' ',(w,h,s:t++[s]))
n(_:b,(w,h,s))=(transpose$map reverse$(take w s):b,(h,w,drop(w-1)s))

The function p doesn't print the output, but returns it as a string, e.g p 4 5 "+--+|||+--+|||" -> "+--+\n|  |\n|  |\n|  |\n+--+\n". For better display use putStr:
putStr $ p 4 5 "+--+|||+--+|||"

+--+
|  |
|  |
|  |
+--+

How it works: I create a w x h block of spaces and replace the first line with the beginning of the input string. Then I rotate the block counterclockwise and repeat replacing the first line another three times.
To prevent cutting of the first character again after turn #4, I append it to the input string before starting.
"Hello World" example, 5 x 4

         |  Start               Turn #1          Turn #2     Turn #3   Turn #4
---------+--------------------------------------------------------------------
String   |  "Hello, World! H"   "o, World! H"    "World! H"  "d! H"    ""
left     |
         |
Block    |  <empty>             Hello            o, W        World     d! H
before   |                                       l                     l  e
rotating |                                       l           ,         r  l
         |                                       e           olleH     o  l
         |                                       H                     W ,o

Edit: found a better way to solve the cut-off-first-character-after-turn-#4 problem.

Ah nice... this is similar to what I tried in CJam, except it works. ;)

><>, 82 80 + 3 = 83 bytes
:2-&\
v!?:<oi-1
/?(0:i
\~ao{2-{~}
\{:?!v1-}o&:&
>:?v!~{{o}ao4.
^  >" "o1-
o;!?l<

Esolang page for ><> (Fish)
This turned out to be shorter than I expected. It uses the straightforward approach of printing the first line, then the columns padded with the central spaces, then the last line.
Input the string via STDIN and the height and width via command line with the -v flag, like so:
py -3 fish.py round.fish -v <height> <width>

Explanation
:2-&           Put W-2 in the register
:?!v1-io       Directly print the first W characters of the input
i:0(?/         Read the rest of the input
~ao{2-{~}      Pop a few leftovers 0s from above, decrement H by 2 and print a newline
               Stack now consists of H = H-2 at the bottom and the rest of the input reversed

[loop]

{:?!v          If H is 0...
  ~                Pop the 0
  l?!;o            Print the rest of the (reversed) input

               Otherwise...
  1-}              Decrement H
  o                Output the top of stack
  &:&              Copy I = W-2 from the register
  :?               If I is nonzero...
    " "o1-             Print a space and decrement I, then repeat from the previous line
  {{o}ao           Print the bottom input character and output a newline
  4.               Jump to the start of the loop (note that I = 0 is leftover from above)

Postscript, 62 bytes
This of course uses binary tokens, but it is equivalent to:
/Courier findfont setfont

0 h moveto

s [
    w {1 0} repeat pop pop
    h {0 -1} repeat pop pop
    w {-1 0} repeat pop pop
    h {0 1} repeat
] xyshow

Here is a hexdump of the file (xxd round.ps):
0000000: 91c7 9243 9295 3020 6892 6b73 5b77 7b31  ...C..0 h.ks[w{1
0000010: 2030 7d92 8392 7592 7568 7b30 202d 317d   0}...u.uh{0 -1}
0000020: 9283 9275 9275 777b 2d31 2030 7d92 8392  ...u.uw{-1 0}...
0000030: 7592 7568 7b30 2031 7d92 835d 92c3       u.uh{0 1}..]..

Run as:
gs -dw=11 -dh=10 -ss="rock beats scissors beats paper beats " round.ps

The output is really small (as a result of not scaling the font at all), so you need to zoom in a fair bit to see it.
This takes advantage of the xyshow operator to write out the string using custom character spacings. In this case, I use the negative vertical spacing to write down, then negative horizontal space to write backwards, then positive vertical space to write upward. Because of this, I don't need to use any sort of string manipulation.

Bash+coreutils, 124
A shell script to get you started:
echo "${3:0:$1}"
fold -1<<<"${3:$1*2+$2-2}"|tac|paste - <(fold -1<<<"${3:$1:$2-2}")|expand -t$[$1-1]
rev<<<"${3:$1+$2-2:$1}"

Pass input as command-line args:
$ ./roundnround.sh 5 4 "Hello, World! "
Hello
    ,
!
dlroW
$

JavaScript, 161 160 158 bytes
The method I've come up with turned out way too long, but oh well, it's been practice. (Also, I got it to spell out r+o[u]+'\n':d.)
function f(o,w,n){s=o.slice(0,w)+'\n';o=o.slice(w);n-=2;r='';for(u=w-2;u--;)r+=' ';for(u=d=0;d=o[2*n+w+~u],u<w+n;u++)s+=(u<n)?(d||' ')+r+o[u]+'\n':d;return s}

For input that doesn't make sense the output is undefined (literally, and a bunch of times), but it works for all the test cases.

slice is shorter than substr, it's not exactly the same but in this case you can use it

Groovy, 140
f={a,x,y->println a.substring(0,x);(1..y-2).each{println a[a.length()-it]+' '*(x-2)+a[it+x-1]}println a.substring(x+y-2,2*x+y-2).reverse()}

call:
f('rock beats scissors beats paper beats ',11,10)

output:
rock beats
          s
s         c
t         i
a         s
e         s
b         o
          r
r         s
epap staeb

K, 80 68 bytes
f:{[s;y;n]`0:(,n#s),({s[(#s)-x+2],((n-2)#" "),s@n+x}'!y-2),,n#|-4!s}

Shortened from 80 thanks to @JohnE.
Original:
f:{s:x;n:z;`0:(,s@!n),({s[(#s)+-2-x],({" "}'!n-2),s@n+x}'!y-2),,(|s@!-4+#s)@!n}

I barely even know how this thing works.
Example usage:
f["Hello, world! ";5;4]

There are some possible optimizations, but I keep making Kona segfault...

You can improve this a bit by using 'take' (dyadic #) and an explicit argument list:
f:{[s;y;n]\0:(,n#s),({s[(#s)-x+2],((n-2)#" "),s@n+x}'!y-2),,n#|-4!s}`. 68 characters by my count.

@JohnE Thanks! I knew about the explicit argument list, but it somehow slipped my mind; I had no clue about dyadic #, though.

K, 55 54 bytes
Using the same approach as randomra's J implementation; Start with a block of spaces and add from the tail of the string to the edge while rotating four times:
f:{`0:*4{((,r#*|x),|:'+*x;(r:-#*x)_*|x)}/((y-2)#" ";x)}

And some examples:
  f["Hello,_World!_";4 5]
Hello
_   ,
!   _
dlroW

  f[">v<^";2 2]
>v
^<

Breaking it down a little for readability,
Generate an NxM block:
  t:2 3#!6
(0 1 2
 3 4 5)

Rotate 90 degrees by using transpose(+) and reverse-each (|:'):
  |:'+t
(3 0
 4 1
 5 2)

So if we have a block of spaces t and a string s, we can prepend a slice of the tail of s to t:
  s: 12 18 17 8 9
12 18 17 8 9
  (,(-#t)#s),|:'+t
(8 9
 3 0
 4 1
 5 2)

We use the form 4 {[x] ... }/( ... ) to repeatedly apply a function to a tuple consisting of the string and the matrix we're building. Each time we do this rotate-and-concatenate step we also chop down the string.
edit:
Another idea is to try splitting the input string into the fragments we want at each rotation, which simplifies the main body of the program. Unfortunately this works out to be slightly longer at 56 bytes:
f:{`0:((y-2)#" "){|:'(,y),+x}/(+\(0,y[0 1 0]-2 1 1))_|x}

If there's a better way to calculate those split points I'm open to suggestions.
edit2:
Rearranging slightly lets me remove a pair of parentheses. 54 bytes!
f:{`0:((y-2)#" "){|:'(,y),+x}/(0,+\y[0 1 0]-2 1 1)_|x}

R, 178
This is an unnamed function taking s, w, h as parameters.  I wish there was a nicer way to split the string.
function(s,w,h){W=w+h-1;H=w+W-1;S=strsplit(s,'')[[1]];O=array(" ",c(h,w+1));O[,w+1]="\n";O[1:h,c(w,1)]=c(S[w:W],S[(W+W-1):H]);O=t(O);O[1:w,c(1,h)]=c(S[1:w],S[H:W]);cat(O,sep='')}

Ungolfed
W=w+h-1;                                 # additional index points
H=w+W-1;                                 # additional index points
S=strsplit(s,'')[[1]];                   # vectorize the string
O=array(" ",c(h,w+1));                   # create an array of spaces
O[,w+1]="\n";                            # set newlines
O[1:h,c(w,1)]=c(S[w:W],S[(W+W-1):H]);    # build middles lines
O=t(O);                                  # transpose array
O[1:w,c(1,h)]=c(S[1:w],S[H:W]);          # build top and bottom lines
cat(O,sep='')                            # cat out results

Test run
> (function(s,w,h){W=w+h-1;H=w+W-1;S=strsplit(s,'')[[1]];O=array(" ",c(h,w+1));O[,w+1]="\n";O[1:h,c(w,1)]=c(S[w:W],S[(W+W-1):H]);O=t(O);O[1:w,c(1,h)]=c(S[1:w],S[H:W]);cat(O,sep='')})("Hello, World! ",5,4)
Hello
    ,
!
dlroW
> (function(s,w,h){W=w+h-1;H=w+W-1;S=strsplit(s,'')[[1]];O=array(" ",c(h,w+1));O[,w+1]="\n";O[1:h,c(w,1)]=c(S[w:W],S[(W+W-1):H]);O=t(O);O[1:w,c(1,h)]=c(S[1:w],S[H:W]);cat(O,sep='')})("+--+|||+--+|||",4,5)
+--+
|  |
|  |
|  |
+--+
> (function(s,w,h){W=w+h-1;H=w+W-1;S=strsplit(s,'')[[1]];O=array(" ",c(h,w+1));O[,w+1]="\n";O[1:h,c(w,1)]=c(S[w:W],S[(W+W-1):H]);O=t(O);O[1:w,c(1,h)]=c(S[1:w],S[H:W]);cat(O,sep='')})(">v<^",2,2)
>v
^<
> (function(s,w,h){W=w+h-1;H=w+W-1;S=strsplit(s,'')[[1]];O=array(" ",c(h,w+1));O[,w+1]="\n";O[1:h,c(w,1)]=c(S[w:W],S[(W+W-1):H]);O=t(O);O[1:w,c(1,h)]=c(S[1:w],S[H:W]);cat(O,sep='')})("rock beats scissors beats paper beats ",11,10)
rock beats
          s
s         c
t         i
a         s
e         s
b         o
          r
r         s
epap staeb
> # Escaped the \ as well
> (function(s,w,h){W=w+h-1;H=w+W-1;S=strsplit(s,'')[[1]];O=array(" ",c(h,w+1));O[,w+1]="\n";O[1:h,c(w,1)]=c(S[w:W],S[(W+W-1):H]);O=t(O);O[1:w,c(1,h)]=c(S[1:w],S[H:W]);cat(O,sep='')})("!\"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_`abcdefghijklmnopqrstuvwxyz{|}~",46,3)
!"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMN
~                                            O
}|{zyxwvutsrqponmlkjihgfedcba`_^]\[ZYXWVUTSRQP
>

T-SQL, 307
While still horrendously long, this turned out to be quite a bit easier (and shorter) than I thought in a query.  Implemented as an inline table valued function for T-SQL.
CREATE FUNCTION f(@S VARCHAR(MAX),@ INT,@H INT)RETURNS TABLE RETURN WITH R AS(SELECT 2i,LEFT(@S,@)S,STUFF(@S,1,@,'')+'|'R UNION ALL SELECT i+1,CASE WHEN i<@H THEN LEFT(RIGHT(R,2),1)+REPLICATE(' ',@-2)+LEFT(R,1)ELSE REVERSE(LEFT(R,@))END,STUFF(STUFF(R,LEN(R)-1,1,''),1,1,'')FROM R WHERE i<=@H)SELECT S FROM R

This recurses through the string @h times.  First recursion clips @W characters from the string.  The middle recursions take the last and first from remaining string with string padding between.  The last recursion reverses what is left.  There is a few lost characters dealing with the way SQL Server treats trailing spaces on VARCHARS.
Test Run
WITH TestSet AS (
    SELECT *
    FROM (VALUES
        ('Hello, World! ',5,4),
        ('+--+|||+--+|||',4,5),
        ('>v<^',2,2),
        ('rock beats scissors beats paper beats ',11,10),
        ('!"#$%&''()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~andfoure',50,3)
        ) Test(S,W,H)
)
SELECT x.S
FROM TestSet
    CROSS APPLY (
        SELECT S FROM dbo.F(S,W,H)
        )x

S
----------------------------
Hello
    ,
!
dlroW
+--+
|  |
|  |
|  |
+--+
>v
^<
rock beats
          s
s         c
t         i
a         s
e         s
b         o
          r
r         s
epap staeb
!"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQR
e                                                S
ruofdna~}|{zyxwvutsrqponmlkjihgfedcba`_^]\[ZYXWVUT

(24 row(s) affected)

Pyth, 33 bytes
J-vw2<zQFkC,_z<>zQJj*-Q2dk;<>_zJQ

Demonstration.

MATLAB, 101
function f(H,W,S)
w=1:W;h=(0:H-3).';n=W+H-2;S(3*n)=' ';S([w;[2*n-h,3*n*ones(H-2,W-2),h+W+1];n-w+W+1])

C++, 398 bytes
Compiler used - GCC 4.9.2 with -std=c++14 flag
#include<bits/stdc++.h>
using namespace std;string s;vector<vector<char>> M;int w,h,p,i,j;void F(int x,int y){if(p<s.size()&&(((!y||y==h-1)&&x>=0&&x<w)||((!x||x==w-1)&&y>=0&&y<h))&&!M[y][x])M[y][x]=s[p++],F(x+1,y),F(x,y+1),F(x-1,y),F(x,y-1);}int main(){getline(cin,s);cin>>w>>h;M.resize(h,vector<char>(w,0));F(0,0);while(i<h){j=0;while(j<w){if(!M[i][j])M[i][j]=32;cout<<M[i][j++];}i++;cout<<endl;}}

Test it here.
Explanation
#include<bits/stdc++.h>
using namespace std;

string s; // input string
vector<vector<char>> M; // output matrix
int w, h, p, i, j;
// w = width
// h = height
// p = iterator over s
// i, j = iterators used later for printing answer

void F( int x, int y )
{
    // If the coordinates (x, y) are either on the first row/column or the last row/column and are not already populated with the input characters, populate them
    if ( p < s.size() && ( ( ( y == 0 || y == h - 1 ) && x >= 0 && x < w ) || ( ( x == 0 || x == w - 1 ) && y >= 0 && y < h ) ) && !M[y][x] )
    {
        M[y][x] = s[p++];
        F( x + 1, y );
        F( x, y + 1 );
        F( x - 1, y );
        F( x, y - 1 );
    }
}

int main()
{
    getline( cin, s );
    cin >> w >> h;
    // Input taken !!

    M.resize( h, vector<char>( w, 0 ) ); // Fill the matrix with null characters initially

    F( 0, 0 ); // This function does all the work

    // Now printing the matrix
    while ( i < h )
    {
        j = 0;
        while ( j < w )
        {
            if ( !M[i][j] )
            {
                M[i][j] = ' ';  // Replace '\0' with ' '
            }
            cout << M[i][j++];
        }
        i++;
        cout << endl;
    }

}

Couldn't you save characters by using char[][] instead?

No, vector<vector<char>> M;M.resize(h,vector<char>(w,0)); is slightly shorter than char** M;M=new char*[h];while(i<h)M[i++]=new char[w]();

Perl, 193 195 bytes
($s,$w,$h,$i,$y)=(@ARGV,0,2);
$o.=substr$s,$i,$w;
$i+=$w;
$o.=sprintf"\n%s%*s",substr($s,2*($w+$h)-$y++-3,1)||' ',$w-1,substr($s,$i++,1)while$y<$h;
print$o."\n".reverse(substr($s,$i,$w))."\n";

I'm sure this can be greatly improved. I'm a noob. >,<

Java 11, 180 bytes
(s,w,h)->{var r=s.substring(0,w)+"\n";int i=w;for(var S=s.split("");i<w+h-2;)r+=S[3*w+2*h-i-5]+" ".repeat(w-2)+S[i++]+"\n";return r+new StringBuffer(s.substring(i,i+w)).reverse();}

Try it online (NOTE: String.repeat(int) is emulated as repeat(String,int) for the same byte-count, because Java 11 isn't on TIO yet.)
Explanation:
(s,w,h)->{               // Method with String & 2 int parameters and String return-type
  var r=s.substring(0,w)+"\n";
                         //  Result-String, starting at the the first row of output,
                         //  which is a substring in the range [0, `w`)
  int i=w;               //  Index-integer, starting at `w`
  for(var S=s.split(""); //  Split the input-String into a String-array of characters
      i<w+h-2;)          //  Loop `i` in the range [`w`, `w+h-2`)
    r+=                  //   Append the result-String with:
       S[3*w+2*h-i-5]    //    The character at index `2*w+2*h-4 - i+w-1`
       +" ".repeat(w-2)  //    Then append `w-2` amount of spaces
       +S[i++]           //    Then append the character at index `i`
       +"\n";            //    And a trailing new-line
  return r               //  After the loop, return `r` as result
         +new StringBuffer(s.substring(i,i+w)).reverse();
                         //  Appended with the last row of output,
                         //  which is a substring in the range [`i`, `i+w`) reversed

Charcoal, 4 bytes
ＢＮＮＳ

Try it online (verbose) or try it online (pure).
Explanation:
Basically a builtin for this challenge. Take the first two integer-inputs as width and height, and prints a box with the third input-string as border:
Box(InputNumber(), InputNumber(), InputString())
ＢＮＮＳ


text:
Calculate the size of the moon
The size of the moon mystery
I am sure you have heard that the moon changes its size. When you're in love and you're lucky, the moon is almost twice in size compared to normal situations. Some people say the reason is the atmosphere which acts as a lens. Others think that it is only a matter of comparison to other objects such as trees nearby. Whatever explanation you read, it's quite subjective.
The size of the moon science
Ok, we're programmers, aren't we? We rely on facts, right? So here's the experiment:

Take a nice camera which supports setting time and aperture manually.
Set your camera to maximum zoom level.
Go out, take some photographs of the moon in order to detect the best settings so that the moon is sharp and lighting is just fine.
Remember the settings
Take photos of the moon with those settings every time you think the moon is large or small.
Calculate the size of the moon in pixels

The camera won't lie, would it? By counting the bright pixels we can effectively measure the size of the moon - at least in pixels.
If the size is the same across all photos, then it's a bug in our brain. If the size differs, then there's room for speculation

the moon really grows (but what does it eat?)
there's an atmospherical lens effect
the moon has an elliptical curve and is sometimes nearer, sometimes further away from earth
...

But I'll leave that open until your task is completed. Of course you want to know in advance if your software can calculate the moon size accurately.
The task
Given a few optimized pictures of the moon, please calculate the size of the moon. The optimization is: the pixels are either black or white. Nothing in between. No antialiasing. That makes it easy, doesn't it?
The caveat: the moon is not always full, you know ... it can be a sickle! But even in shape of a sickle, the size of the moon is larger. So you'll calculate the full size, please.

Your program takes a PNG as input, e.g. as file name command line argument, piped into stdin or as a Bitmap object (of a standard framework library) if you write a function instead of a program.
Your program works with any reasonable input bitmap size, not necessarily square. Minimum width and height of 150 pixels are guaranteed.
The full moon covers at least 25% of the picture.
Your program outputs the calculated size of the moon in pixels as if it were a full moon.
We assume that the moon is a perfect sphere.
The exact size is always an integer number, but you can output a decimal number if your calculation returns that.
The accuracy should be between 98% and 102%. (That's rather a guess than something I could guarantee to be achievable. If you think it's too hard to reach, please leave a comment.)

Update:

The center of the moon isn't necessarily in the middle of the picture.
The minimum visible area is 5% of the moon or 1.25% of the total number of pixels.
The picture is taken in a way that the whole moon would fit the image, i.e. the total number of pixels is an upper boundary for the moon size.
The moon will not be cropped / clipped.

The samples
You can generate your own samples using the blend file if you like. I have created the following pictures for you. You can count pixels in a PNG file using WhitePixelCounter.exe (needs .NET) to check whether the image contains black and white pixels only and how many of them.
The following 256x256 pixel images differ in the amount of white pixels, but should all result in a calculated moon size of 16416 pixels.

And these 177x177 pixel images should return 10241 pixels. The images are basically the same, but this time a camera with a different focal length was used.

Non-square and non-centered samples with a result of 9988:

Oh, I don't have a reference implementation for now and I even don't know whether I am able to implement something. But in my brain there is a strong belief that tells me it must be mathematically solvable.
The rules
This is Code Golf. The shortest code on 2015-03-30 gets accepted.

In all the examples, the center of the moon appears to be centered within the picture.  Can we assume the moon will always be centered?

your accuracy of +/-2% on area corresponds to +/-1% on diameter: example r=100 pixels, area = 10000pi; r=101 pixels, area = 10201pi. Your smaller image has r=72 therefore d=144 so it should just be possible. For images below d=100 I think the accuracy could not be met, however.

@DigitalTrauma: the center needn't be in the middle.

@MartinBüttner: the minimum visible percentage is 5% of the moon or 1,25% of the picture.

@MartinBüttner: ok, I've updated the question, updated the blend file to produce non-square, non-centered images by default. You can download all images here (*.png.zip). Updated pixel counter as well: outputs some more information and checks the 1.25% rule.

@steveverrill: I have written a kind of "reference implementation". It fails if the sickle is pointing straight upwards, but for all other pictures, the area is from -1.66% to +1.88%.

Mathematica 126 119 109 bytes
Mathematica can measure the elongation of a component in an image.  A full moon, being perfectly symmetric, has an elongation of 0, on a scale of 0 to 1.
A diminishing moon becomes increasingly elongated, to a maximum of roughly 0.8.
0.998 -0.788 x-0.578 x^2 was the empirically determined model (based on the large photos) for `predicting the fullness of the moon (by area), given its elongation.
I adjusted the model to 1- 0.788 x -0.578 x^2 so that with exactly zero elongation (full moon) the model will return 1 for the pixel scale factor.  It saves 4 bytes and still stays within the accuracy limits.
This model is used for any size images.  The moon image does not need to be centered. It also does not need to cover a fixed proportion of the photo.
Here are the data points (elongation, displayedMoonPixels/fullMoonPixels) for the large images and the parabolic model that was generated to fit the data. Linear models fit ok, but the quadratic model is dead on, within limits (see below).
Here the data are from the large pictures.  So is the model

Below, the data (the red points) are from the small pictures.  The model (the blue curve) is the one generated by the large pictures, the same one as displayed above.
The smallest crescent has 7.5% the area of a full moon.  (The smallest crescent among the large photos is 19% of a full moon.)  If the quadratic model had been based on the small photos the fit below would be better, only because it accommodated the small crescent.  A robust model, that would stand up under a wide range of conditions, including very small crescents, would be better made from a larger variety of pictures.
The closeness of fit shows that the model was not hard-coded for the given pictures. We can be fairly certain that the elongation of a moon is independent of the size of the photo, as one would expect.

f takes the image, i, as input and outputs the predicted size of the full moon, in pixels.  It works for off-center shots.
As the data below show, it all of the test cases except one.
The moons were arranged from full to most diminished.
i_~c~t_ := Max@ComponentMeasurements[i, t][[All, 2]];
f@i_ := i~c~"Count"/(1 - 0.788 x - 0.578 x^2 /. x -> i~c~"Elongation")

More than one image component may turn up in a photo.  Even a single pixel separated from the others will be considered a distinct component.  For this reason, it is necessary to search "all" components, to find the one that has the greater number of pixels.  (One of the small photos has more than one image component.)
Large pictures
Predictions of moon size made from the large photos were uniformly accurate.
{"predicted size of full moon", f[#] & /@ large}
{"accuracy", %[[2]]/16416}

{"predicted sizes of full moon", {16422., 16270.9, 16420.6, 16585.5, 16126.5, 16151.6}}
{"accuracy", {1.00037, 0.991161, 1.00028, 1.01033, 0.982367, 0.983891}}

Small pictures
Predictions of moon size made from the small photos were uniformly, with one great exception, the final picture.  I suspect the issue stems from the fact that the crescent is very narrow.
{"predicted sizes of full moon", f[#] & /@ small}
{"accuracy", %[[2]]/10241}

{"predicted sizes of full moon",{10247.3, 10161., 10265.6, 10391., 10058.9, 7045.91}}
{"accuracy", {1.00061, 0.992192, 1.0024, 1.01465, 0.982221, 0.68801}}

Seems like I should learn Mathematica one day. How long did it take you to solve it without golfing?

@Thomas W It took a 2-3 hours experimenting with various sorts of image processing features and other (linear) models until I obtained the graph you see posted. The coding was not very difficult.  And there is almost no golfing other than uniting separate functions into a single function.

104: i_~c~t_:=Max[#2&@@@i~ComponentMeasurements~t];f@i_:=i~c~"Count"/(1-0.788x-0.578x^2/.x->i~c~"Elongation")

For unknown reasons, the #2&@@@ suggestion does not work

Huh, I'll look into that later. Another way to shorten c is c=Max@ComponentMeasurements[##][[All,2]]&

I can't get c=Max@ComponentMeasurements[##][[All,2]]& to work either.

@ThomasW.  The code now calls ComponentMeasures twice, once for determining the size (Count) and again for obtaining the elongation of the displayed moon.  For that reason, what was a single function is now two functions, c (for ComponentMeasures) and f (which predicts the size of the full moon).

J, 227 207 bytes (maximal error 1.9%)
My main idea is that if we can find 3 points on the contour of the moon which are on the contour of the full moon too we can compute the circumcircle of these points. That circumcircle will be to full moon.
If we find two white points with maximal distance those will always be such points as they will be either a real diagonal in the full moon or the endpoints of the crescent. We can find the pair of points with the greatest distance in any graph by selecting the point furthest from any given starting point and then selecting the point furthest from the selected one.
We find a third point with a maximal value of the products of the distances from the previous points. This will always be on the contour and on the outer side of a crescent or the bigger side of a gibbous.
The diameter of the circumcircle is computed as the length of one side divided by the sinus of the opposite angle.
The time-complexity of this method is linear in the size of the input image.
Code
f=.3 :0
load'graphics/png'
i=.readpng y
p=.(,i=_1)#|:,"%.0 1|:,"0/&>/<@i."*$i
s=.%:+/|:*:(-1|.]) (([,],:m@(*&d))(m@d))(m=.p{~(i.>./)@])(d=.+/@:*:@((|:p)-])) 0{p
o.*:-:({.s)%0 o.((+/-2*{.)*:s)%2**/}.s
)

The function expects the input filename as a string.
(For a (little) more readable version check revision history.)
Code explanation

p is a list of white pixel coordinates (called points in the future)
function d computes distances between elements of p and a given point
the second part of the definition of s creates a 3-point list:

A is the furthest point from the first point in the list
B is the furthest point from A
C is a point with a maximal value of distance form A times distance from B

s is the side-lengths of triangle ABC
the last line computes the area of the circumcircle of ABC which is the full moon

Results
The largest error is 1.9%.
Images are in the same order as in the question.
Output  Accuracy
----------------
  16407 0.999453 NB. Large images
16375.3 0.997523
16223.9 0.988301
16241.5 0.989369
16262.6 0.990654
16322.1 0.994279
10235.3 0.999445 NB. Small images
10235.3 0.999444
10221.2 0.998067
10220.3 0.997978
  10212 0.997169
10229.6  0.99889
9960.42 0.997239 NB. Offset images
9872.22 0.988408
10161.8   1.0174
9874.95 0.988681
 9805.9 0.981768

+1 for participating and mentioning the approach. I'm sorry I didn't specify that the center needn't be in the middle. Accidentally the sample images are all centered. That's my fault.

@ThomasW. Temporally deleted my answer until I correct it.

Matlab 162 156 (not quite in the current error margin)
First of all: The accuracy is under 2% for all but one image in each of the two series, where it is greater (about 5% and 14%).
My approach was finding the two pixels of the moon that are the furthest away from each other, and then using this as an estimate for the diameter.
a=imread(input(''));                 %read input image
b=a(:,:,1)>0;                        %binarize red channel
s=size(b);                           %get size of the image
[x,y]=meshgrid(1:s(1),1:s(2));
z=(x+i*y).*b;z=z(z~=0);              %find the coordinates of all white pixels (as a list)
o=ones(size(z(:)))*z(:)';
disp(max(max(abs(o-o.').^2))*pi/4);  %calculate the maximum of the distance between each possible pair and evaluate area formula

These are the accuracy results (reltative deviation 1 - (predicted size / real size))
0.0006 0.0025 0.0169 0.0500 0.0521 0.0113 0.0006 0.0006 0.0026 0.0472 0.1383 0.0131

C# - 617
This solution does not work for all images, because on one of the images, the slope (m) becomes infinity.
The principle was mentioned before:

Find two points with maximum distance (red)
Imagine a line between them (red)
Imagine a line with rectangular angle in the middle (green)
Find white points on the green line
Use the one with maximum distance from other points (green)
Calculate the area of a circle from three points

The problematic case is this one, where the slope is infinity. It is possible to workaround by rotating the image 90° or in code, loop over the y axis instead of x.

double A(Bitmap b){var a=new List<P>();for(var y=0;y<b.Height;y++)for(var x=0;x<b.Width;x++)if(b.GetPixel(x,y).R>0)a.Add(new P{x=x,y=y});double c=0.0,d=0.0,e=0.0,f=0.0,g=0.0,n=double.MaxValue;foreach(var h in a)foreach(var i in a){var t=Math.Sqrt(Math.Pow(h.x-i.x,2)+Math.Pow(h.y-i.y,2));if(t>c){d=h.x;f=i.x;e=h.y;g=i.y;c=t;}}c=(f-d)/(e-g);for(int x=0;x<b.Width;x++){int y=(int)(c*x+((e+g)/2-c*(d+f)/2));if(y>=0&&y<b.Height&&b.GetPixel(x,y).R>0){var s=(g-e)/(f-d);var q=(y-g)/(x-f);var j=(s*q*(e-y)+q*(d+f)-s*(f+x))/(2*(q-s));var k=-(j-(d+f)/2)/s+(e+g)/2;var l=(j-d)*(j-d)+(k-e)*(k-e);if(l<n)n=l;}}return Math.PI*n;}

The minimum accuracy is

+1,89% for the 256 pixel images
-0.55% for the 177 pixel images
-1.66% for the non-square images

source: Stack Exchange

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

sort: votes


text:
Summing colored graphs
In some cases, often in physics, you have to sum graphs. Your challenge is to write, in a language of your choice, a program or function that takes multiple graphs as images, calculates all possible sums, and outputs the result.
Graphs
The graphs are images that contain a white (rgb(255, 255, 255)) background with a non-white pixel in each column. Examples:

The values of the script are represented as the Y positions of the colored pixels. The value at a certain X coordinate is equal to the Y position of the uppermost colored pixel in that column, with coordinates starting from 0 at the bottom left. There may or may not be additional colored pixels below those pixels for aesthetical reasons.
Task
Your task is to write, in a language of your choice, a program or function that takes multiple graphs as images, calculates all the possible 2^n - 1 sums, and outputs the result.
A sum of graphs is a graph where each column's value is equal to the sum of the values of the corresponding column in each of the input graphs.
The graphs will come in multiple colors. The result image must contain all possible sums of the graphs as other graphs, including the original graphs but excluding the zero sum.
The color of each sum is determined by the average of the colors of the graphs included, for example, graphs of colors rgb(255, 0, 255) and rgb(0, 255, 255) would produce a graph of rgb(128, 128, 255) (may also be rounded down).
The resulting image should be as high as required to fit all graphs. This means you may have to output an image larger than any of the inputs.
The order in which the resulting graphs are drawn to the resulting image does not matter, i.e. if result graphs overlap, you may choose which one is on top, but it has to be one of the graphs, not a combination of their colors.
You may assume that the input images are of equal width, that all columns of the images have at least one non-white pixel, and that the heights of the images (including output) are below 4096 pixels.
Example
Input A:

Input B:

Example output:

(In case someone is interested, I copy-pasted the data for these from random companies' stock charts. That was the first way I found to get realistic data as CSV.)
Rules

You may choose any bitmap image input file format.
You may choose any bitmap image output file format, which doesn't have to match the input.
You may use image processing libraries, however any functions to complete this task directly are banned.
Standard loopholes apply.
This is code-golf, so shortest code in bytes wins.

Graph generator script
Here is a Python 2 script that generates graphs. Input is given in lines, with the three first lines as the RGB color and the rest as data, terminated by EOF.
import PIL.Image as image
import sys

if len(sys.argv) < 2:
    sys.stderr.write("Usage: graphgen.py <outfile> [infile]")
    exit(1)
outfile = sys.argv[1]
if len(sys.argv) > 2:
    try:
        stream = open(sys.argv[2], "r")
        data = stream.read()
        stream.close()
    except IOError as err:
        if err.errno == 2:
            sys.stderr.write("File \"{0}\" not found".format(sys.argv[2]))
        else:
            sys.stderr.write("IO error {0}: {1}".format(err.errno, err.strerror))
        exit(1)
else:
    data = sys.stdin.read()

try:
    items = map(int, data.strip().split("\n"))
    red, green, blue = items[:3]
    items = items[3:]
    highest = max(items)
except (ValueError, TypeError, IndexError):
    sys.stderr.write("Invalid value(s) in input")

img = image.new("RGB", (len(items), highest + 1), (255, 255, 255))

prev = items[0]
img.putpixel((0, highest - items[0]), (red, green, blue))
for x, item in enumerate(items[1:]):
    img.putpixel((x + 1, highest - item), (red, green, blue))
    if item < prev:
        for i in range(item + 1, prev):
            img.putpixel((x, highest - i), (red, green, blue))
    else:
        for i in range(prev + 1, item):
            img.putpixel((x + 1, highest - i), (red, green, blue))
    prev = item

img.save(outfile, "png")

@MartinBüttner I'm currently making one for two graphs. I'm doing it by hand (no reference impl yet), so I don't know if I have the patience for 3. Also, the three I gave can't be summed since they are of different widths.

So if there are n input graphs, there will be 2^n - 1 lines in the output image?

@PeterTaylor Yes.

I take it the output doesn't actually need to contain any vertical lines? Only the top-most pixel in each column?

@MartinBüttner That is correct, since that data can still be parsed as a graph as defined in the first section.

If the graphs overlap, do we have to choose one of the colours, or could we, say, multiply them?

@MartinBüttner Let's say you have to choose one, as at least YouTube and Excel do.

MATLAB, 405
Call via: f('https://i.stack.imgur.com/ffCzR.png','https://i.stack.imgur.com/zHldg.png')
function f(varargin)
for k=1:nargin
i=im2double(imread(varargin{k}))
V(k,:)=size(i,1)-cellfun(@(V)find(any(V~=1,3),1),num2cell(i,[1,3]))
C(k,:)=i(find(any(i(:,1,:)~=1,3),1),1,:)
end
s=2^nargin-1
G=dec2bin(1:s)-'0'
C=bsxfun(@rdivide,G*C,sum(G,2))
V=G*V
m=max(V(:))
r=ones(m+1,size(V,2))
g=r
b=r
for i=1:s
M=bsxfun(@eq,(m:-1:0).',V(i,:))
r(M)=C(i,1)
g(M)=C(i,2)
b(M)=C(i,3)
end
imwrite(cat(3,r,g,b),'S.png')

Python, 422
Call from commandline python plotsum im1.png im2.png im3.png
import sys
from numpy import*
from scipy import misc as m
R=m.imread
r=range
a=array
N=sys.args[1:]
L=len(N)
P=[map(argmin,R(n,1).T)for n in N]               #converts image to list of heights, counting from the top
C=a([R(N[i])[P[i][0],0,:]for i in r(L)])         #finds and stores the colour
P=a([len(R(N[i]))-a(P[i])for i in r(L)])         #flips the numbers, measures actual heights from bottom
w=len(P[0])
h=max(sum(P,0))+1                                    #compute dimensions
G=ones((h,w,3))*255                                  #and make a white grid
for i in r(1,2**L):
 z=where(a(list(bin(i)[2:].zfill(L)))=='1');y=sum(P[z],0)    #sum the graphs
 for x in r(w):G[y[x],x,:]=average(C[z],0)                   #average the colours
m.imsave('S.png',G[::-1])                            #flip image vertically and save

Example output

Another example

That was a tricky one, high-level array operations and using arrays as indices helps a lot here. I dont expect to see solutions under 1000 bytes except in Mathematica and Matlab

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


text:
Golfing every ASCII character in 99
99 is a programming language I invented earlier this week for my challenge Write an interpreter for 99. (Invented but never needed to implement thanks to half a dozen of you. ;) ) The full language spec is in that challenge so I'm not bothering to repost it all here.
In 99 you can print individual ASCII characters to stdout, but due to the constraints of the language, it is not always clear how to print a particular character as concisely as possible.
For each of the 128 ASCII characters, write a 99 program that takes no input and outputs that sole character. You may code any or all of these answers by hand, or you can write another program (in any language you like) to generate them for you.
The sum of the characters in each of your 128 99 programs is your score. The lowest score wins. Newlines count as one character.
Remember, in 99, only even sized variables like 9999 output ASCII characters (odd sized variables output integers). Their value is divided by 9 and then taken mod 128, so values do not need to be in a certain range to map to ASCII characters. For example, the internal values 297, 1449, and -855 all correspond to the character ! because when they are divided by 9 and taken mod 128, they all become 33, which is the character code for !.
If you need an interpreter for 99, I'd suggest Mac's Python answer.
I know I said my next challenge would be more interactive but I'm still working on the stuff for that one.

One assignment, 2075 (optimal)
This should be the optimal value (unless I have a big error in reasoning or my testing sucks).
First of all. There are only 7 different numbers (mod 128) you can express in 99. All values of 7 or more 9 evaluate to the same number 71. (Because 10^8 mod 128 == 0, 10^9 mod 128 == 0, ...)
If one of the 4 values you can express with an even number of nines, than outputting this number is clearly the optimal solution.
Otherwise I try to reach the number with one assignment statement (assign to 99) and print 99. As it turns out the maximal program size with this approach is 22 char. Clearly using Goto takes definitely more than that. The only possibility that a one-assignment solution might be beaten is a solution with two assignments. I tested this (hopefully with no errors, the code for this is quite messy), and found no solution for any ASCII char.
Therefore only checking the 4 direct numbers and the one-assignment approach should be enough to find the optimal solution. The following Python (2 and 3 compatible) program generates all programs and sums up their lengths. It uses a simple IDA* approach.
from itertools import count

def nines_to_dec(nines):
    return ((10**nines - 1) // 9) % 128

def shortest_representation(ascii_value):
    # try simple output,
    # max code length is 8, (8 nines == 10 nines == 12 nines == ...)
    # if this works, than this is the shortest representation

    for nines in range(2, 9, 2):
        if nines_to_dec(nines) == ascii_value:
            return "9" * nines

    # otherwise try one assignment
    for length in count(1):
        result = assignment(ascii_value, length, [])
        if result:
            return "99 " + result + "\n99"

def assignment(value, left, nines_list):
    if left == 0:
        eval_numbers = [nines_to_dec(nines) for nines in nines_list]

        if (sum(eval_numbers[::2]) - sum(eval_numbers[1::2])) % 128 == value:
            return " ".join("9" * nines for nines in nines_list)
        else:
            return False

    for nines in range(1, 8):
        left2 = left - nines - 1 # -1 for space
        if left2 >= 0:
            result = assignment(value, left2, nines_list + [nines])
            if result:
                return result

    return False

lengths = []
for i in range(128):
    program =shortest_representation(i)
    lengths.append(len(program))
    print("ASCII-value: {}, ASCII-char: {}".format(i, chr(i)))
    print(program)

print(sorted(lengths))
print(sum(lengths))

The output is of the following form:
....
ASCII-value: 65, ASCII-char: A
99 9 999999 9999999
99
ASCII-value: 66, ASCII-char: B
99 9 99 9999 99
99
ASCII-value: 67, ASCII-char: C
99 9 99 9 99 9999
99
....

You can find the complete output at: http://pastebin.com/bKXLAArq
The char with the shortest program (2 char) is vertical tab - 11 with a program length of 2, the chars with the longest programs (22 char) are bell - 7 and A - 65.
The sum for all programs is 2075.
And by the way, I used the k/q interpreter from tmartin. I hat quite some troubles with the other ones (Python, Perl, CJam). Not sure if it was my fault.

It would help implementors of interpreters if you could describe what troubles you had. Great answer.

A Variety of Techniques, 42109
For the numbers, instead of calculating the large ASCII character, I just calculated the number's value. You only said to be able to output the character, so this should still work.
EDIT: Switched the numbers to use the ASCII characters, so disregard that. I left the original number code in the Java code but commented out in case anyone wanted to use it.
Some of these I did by hand, most I just wrote a program to type out.
These are made up of 1-4 lines each, so they're fairly friendly to just copy and paste into a program. It should be noted that they don't work in succession due to my generated code not preserving variable states.
The most common technique used here was the same as orlp's approach:

Keep subtracting 9 from 99, then output.

My version differs by using some custom cases and composing a lot of the math onto just one line.
Custom cases are just where the character can be represented with just a bunch of 9's and no math or my generation code could be shortened.
Programs
I've put the output on Pastebin for those of you who don't feel like running the program:
http://pastebin.com/Cs6WZUfb
Java Code I used:
public class AsciiGen99 {

  public static void main(String[] args) {
    long totalsize = 0;
    for (int i = 0; i < 128; i++) {
      System.out.println("\n The program for ASCII code " + i + " is as follows:\n");
      String yea = find(i);
      if (yea != null) {
        System.out.println(yea);
        totalsize += yea.length();
      } else {
        String v = "99 9 9\n9 99 9";
        if (i != 0) {
          v += "\n99";
          for (int j = 0; j < i; j++) {
            v += " 99 9";
          }
        }

        v += "\n99";

        System.out.println(v);
        totalsize += v.length();
      }
    }
    System.out.println(totalsize);
  }

  public static String find(int i) {
    switch (i) {
      case '\0':
        return "99 9 9\n99";
      case '\1':
        return "99 9\n99";
    }
//    if (48 <= i && i <= 57) {
//      switch (i) {
//        case '0':
//          return "9 9 9\n9";
//        case '1':
//          return "9";
//        case '2':
//          return "999 9 9\n9 999 9\n999 999 9 999 9\n999";
//        case '3':
//          return "999 9 9\n9 999 9\n999 999 9 999 9 999 9\n999";
//        case '4':
//          return "999 9 9\n9 999 9\n999 999 9 999 9 999 9 999 9\n999";
//        case '5':
//          return "999 9 9\n9 999 9\n999 999 9 999 9 999 9 999 9 999 9\n999";
//        case '6':
//          return "99 9 9\n9 99 9\n999 99 9 99 9 99 9 99 9 99 9 99 9\n999";
//        case '7':
//          return "99 9 9\n9 99 9\n999 99 9 99 9 99 9 99 9 99 9 99 9 99 9\n999";
//        case '8':
//          return "99 9 9\n9 99 9\n999 99 9 99 9 99 9 99 9 99 9 99 9 99 9 99 9\n999";
//        case '9'://ironic
//          return "99 9 9\n9 99 9\n999 99 9 99 9 99 9 99 9 99 9 99 9 99 9 99 9 99 9\n999";
//      }
//    }
    int x, a;
    for (x = 0; x < 100000; x++) {
      a = i + 128 * x;
      String s = "" + a*9;
      if (containsOnly9(s) && (s.length() & 1) == 0) {
        return ("" + (a * 9));
      }
    }

    return null;
  }
  public static boolean containsOnly9(String s) {
    for (char c : s.toCharArray()) {
      if (c != '9' && c != ' ' && c != '\n' && c != '\r' && c != '\t') {
        return false;
      }
    }
    return true;
  }
}

You do have to actually output the character, not just the number. So all the programs with 999 at the end need to be fixed.

Ah, alrighty then, I'll fix that in a moment.

Should be fixed now unless I missed something. I left the original code in but commented out in case someone wanted to use the numbers like that. Pastebin was also edited.

Great. Though for some I think you could have just added 99 999\n99 (to reassign 999 to 99 so it will print as a character).

Repeated subtraction, 65280
A trivial solution to compare against. Keep subtracting 9 from 99, then output. Example for ASCII character 10:
99 99 9
99

There are 128 programs. The first program is two characters long (99), each one after that is 8 characters (99 99 9\n) longer than the previous one.
Python program generating programs seperated by empty lines, and computing score:
score = 0
for n in range(128):
    program = "99 99 9\n" * n + "99"
    score += len(program)
    print(program + "\n")

print(score)


text:
Monte Carlo estimator of Pi
Happy Pi Day everyone! For no reason at all, I'm trying to construct a Monte Carlo estimator of Pi that is as short as possible. Can we construct one that can fit in a tweet?
To clarify, what I have in mind is the typical approach of drawing random points from the unit square and calculating the ratio that fall within the unit circle. The number of samples can be hard coded or not. If you hardcode them, you must use at least 1000 samples. The result may be returned or printed as a floating point, fixed point or rational number.
No trig functions or Pi constants, must be a Monte Carlo approach.
This is code golf, so the shortest submission (in bytes) wins.

are trig functions allowed? I suggest you explicitly ban them.

((0..4e9).map{rand**2+rand**2<1}.to_s.sub(/./,"$1.")

@JanDvorak How is that supposed to work? Doesn't the map give you an array of true and false?

@MartinBüttner Ah, oops, sorry. .filter{...}.size should work, though.

@JanDvorak Indeed. That's really neat :)

Good point, minimum of 1000, any return type.

What about types of programs? Full programs, functions. Do you require stdin->stdout?

@Maltysen we've got defaults for that (see the tag wiki).

I actually did post a π approximation  in a tweet: https://twitter.com/ERDekhayser/status/568170358664237057

80386 machine code, 40 38 bytes
Hexdump of the code:
60 33 db 51 0f c7 f0 f7 e0 52 0f c7 f0 f7 e0 58
03 d0 72 03 83 c3 04 e2 eb 53 db 04 24 58 db 04
24 58 de f9 61 c3

How to get this code (from assembly language):
    // ecx = n (number of iterations)
    pushad;
    xor ebx, ebx; // counter
    push ecx; // save n for later
myloop:
    rdrand eax; // make a random number x (range 0...2^32)
    mul eax; // calculate x^2 / 2^32
    push edx;
    rdrand eax; // make another random number y
    mul eax; // calculate y^2 / 2^32
    pop eax;
    add edx, eax; // calculate D = x^2+y^2 / 2^32 (range 0...2^33)
    jc skip; // skip the following if outside the circle
    add ebx, 4; // accumulate the result multiplied by 4
skip:
    loop myloop;
    push ebx; // convert the result
    fild dword ptr [esp]; // to floating-point
    pop eax;
    fild dword ptr [esp]; // convert n to floating-point
    pop eax;
    fdivp st(1), st; // divide

    popad;
    ret;

This is a function using MS fastcall calling convention (number of iterations is passed in register ecx). It returns the result in the st register.
Fun things about this code:

rdrand - just 3 bytes to generate a random number!
It uses (unsigned) integer arithmetic until the final division.
The comparison of squared distance (D) with squared radius (2^32) is performed automatically - the carry flag contains the result.
To multiply the count by 4, it counts the samples in steps of 4.

Comment should read "Calculate x^2 % 2^32"

@ColeJohnson No - the random number is in eax; the mul command multiples it by itself and puts the high part in edx; the low part in eax is discarded.

rdrand is not a 80386 instruction.

Matlab/Octave, 27 bytes
I know there already is a Matlab/Octave answer, but I tried my own approach. I used the fact that the integral of 4/(1+x^2) between 0 and 1 is pi.
mean(4./(1+rand(1,1e5).^2))

A different algorithm is always great! Also, more efficient!

CJam, 27 23 22 or 20 bytes
4rd__{{1dmr}2*mhi-}*//

2 bytes saved thanks to Runner112, 1 byte saved thanks to Sp3000
It takes the iteration count from STDIN as an input.
This is as straight forward as it gets. These are the major steps involved:

Read the input and run the Monte Carlo iterations that many times
In each iteration, get sum of square of two random floats from 0 to 1 and see if it is less than 1
Get the ratio of how many times we got less than 1 by total iterations and multiply it by 4 to get PI

Code expansion:
4rd                     "Put 4 on stack, read input and convert it to a double";
   __{            }*    "Take two copies, one of them determines the iteration"
                        "count for this code block";
      {1dmr}2*          "Generate 2 random doubles from 0 to 1 and put them on stack";
              mh        "Take hypot (sqrt(x^2 + y^2)) where x & y are the above two numbers";
                i       "Convert the hypot to 0 if its less than 1, 1 otherwise";
                 -      "Subtract it from the total sum of input (the first copy of input)";
                    //  "This is essentially taking the ratio of iterations where hypot";
                        "is less than 1 by total iterations and then multiplying by 4";

Try it online here

If average value of 1/(1+x^2) is also considered as Monte Carlo, then this can be done in 20 bytes:
Urd:K{4Xdmr_*)/+}*K/

Try it here

I tried a CJam answer as well, and managed to come in 2 bytes under your score. But my code came out so similarly to yours, I'd feel dirty posting it as a separate answer. Everything was the same except for the variable choice and these two optimizations: get a random number from 0 to 1 with 1dmr instead of KmrK/, and check if the sum of the squares is greater than 1 with i instead of 1> (I thought this one was particularly clever).

@Runer112 Thanks. the i trick is really neat! And damn the lack of documentation for 1dmr

R, 40 (or 28 or 24 using other methods)
mean(4*replicate(1e5,sum(runif(2)^2)<1))

mean(4*sqrt(1-runif(1e7)^2))

mean(4/(1+runif(1e7)^2))

Python 2, 56
Another Python one, if numpy is allowed, but pretty similar to Matlab/Octave:
import numpy;sum(sum(numpy.random.rand(2,8e5)**2)<1)/2e5

Mathematica, 42 40 39 bytes (or 31/29?)
I've got three solutions all at 42 bytes:
4Count[1~RandomReal~{#,2},p_/;Norm@p<1]/#&
4Tr@Ceiling[1-Norm/@1~RandomReal~{#,2}]/#&
4Tr@Round[1.5-Norm/@1~RandomReal~{#,2}]/#&

They are all unnamed functions that take the number of samples n andd return a rational approximating π. First they all generate n points in the unit square in the positive quadrant. Then they determine the number of those samples that lie within the unit circle, and then they divide by the number of samples and multiply by 4. The only difference is in how they determine the number of sampples inside the unit circle:

The first one uses Count with the condition that Norm[p] < 1.
The second one subtracts the norm of each point from 1 and then rounds up. This turns numbers inside the unit circle to 1 and those outside to 0. Afterwards I just sum them all up with Tr.
The third one does essentially the same, but subtracts the from 1.5, so I can use Round instead of Ceiling.

Aaaaaand while writing this up, it occurred to me that there is indeed a shorter solution, if I just subtract from 2 and then use Floor:
4Tr@Floor[2-Norm/@1~RandomReal~{#,2}]/#&

or saving another byte by using the Unicode flooring or ceiling operators:
4Tr@⌊2-Norm/@1~RandomReal~{#,2}⌋/#&
4Tr@⌈1-Norm/@1~RandomReal~{#,2}⌉/#&

Note that the three rounding-based solutions can also be written with Mean instead of Tr and without the /#, again for the same bytes.

If other Monte Carlo based approaches are fine (specifically, the one Peter has chosen), I can do 31 bytes by estimating the integral of √(1-x2) or 29 using the integral of 1/(1+x2), this time given as a floating point number:
4Mean@Sqrt[1-1~RandomReal~#^2]&
Mean[4/(1+1~RandomReal~#^2)]&

You've got three solutions to life, universe and everything and you decide to ruin it? Heresy.

@Sieg And I did it knowing full well what I was doing. ;)

Commodore 64 Basic, 45 bytes
1F┌I=1TO1E3:C=C-(R/(1)↑2+R/(1)↑2<1):N─:?C/250

PETSCII substitutions: ─ = SHIFT+E, / = SHIFT+N, ┌ = SHIFT+O
Generates 1000 points in the first quadrant; for each, adds the truthness of "x^2+y^2<1" to a running count, then divides the count by 250 to get pi.  (The presence of a minus sign is because on the C64, "true" = -1.)

What does (1) do?

@echristopherson, you're misreading it.  / isn't the division symbol, it's the character produced by typing SHIFT+N on a Commodore 64 keyboard.  R/(1) is the shortcut form for RND(1), ie. "produce a random number between 0 and 1 using the current RNG seed".

Oh, you're right! Good old PETSCII graphics characters.

Great answer. I found that if I seed the random numbers first like r=rnd(-1) then I consistently got the answer 3.14, otherwise I got an answer close to PI. My listing is optimised to be slightly more performant, but is (in lower-case to illustrate the keyword abbreviations) it is: 0dIi,c,x,y,r:x=1:r=rN(-x):y=2:fOi=xto1e3:c=c-(rN(x)^y+rN(x)^y<x):nE:?c/250

Correction, with r=rnd(-7) not -1 - so change r=rN(-x) with r=rN(-7)

Python 2, 77 75 bytes
from random import*;r=random;a=0;exec"a+=r()**2+r()**2<1;"*4000;print a/1e3

Uses 4000 samples to save bytes with 1e3.

You could get a little more accuracy at no cost with ...*8000;print a/2e3.

J, 17 bytes
Computes the mean value of 40000 sample values of the function 4*sqrt(1-sqr(x)) in the range [0,1].
Handily 0 o.x returns sqrt(1-sqr(x)).
   1e4%~+/0 o.?4e4$0
3.14915

Java, 108 bytes
double π(){double π=0,x,i=0;for(;i++<4e5;)π+=(x=Math.random())*x+(x=Math.random())*x<1?1e-5:0;return π;}

Four thousand iterations, adding 0.001 each time the point is inside the unit circle. Pretty basic stuff.
Note: Yes, I know I can shed four bytes by changing π to a single-byte character. I like it this way.

why not 9999 iterations ?

@Optimizer It makes the sum shorter. For 9999 iterations, I'd have to add a more precise number each time instead, which costs me digits.

You could save another byte and improve precision by using "4e5" and "1e-5" for the numbers.

@VilmantasBaranauskas Thanks! I always forget about that :) It's tempting to use 4e9 and 1e-9 instead, but that takes quite a while...

Protip: when golfing, you should actually reduce the bytes, not artificially increase them

Matlab or Octave 29 bytes (thanks to flawr!)
mean(sum(rand(2,4e6).^2)<1)*4

(I am no quite sure if <1 is OK. I read it should be <=1. But how big is the probability to draw exactly 1...)
Matlab or Octave 31 bytes
sum(sum(rand(2,4e3).^2)<=1)/1e3

Very nice idea! You could save two additional bytes with mean(sum(rand(2,4e6).^2)<1)*4.

><> (Fish), 114 bytes
:00[2>d1[   01v
1-:?!vr:@>x|  >r
c]~$~< |+!/$2*^.3
 .41~/?:-1r
|]:*!r$:*+! \
r+)*: *:*8 8/v?:-1
;n*4, $-{:~ /\r10.

Now, ><> doesn't have a built-in random number generator. It does however have a function that sends the pointer in a random direction. The random number generator in my code:
______d1[   01v
1-:?!vr:@>x|  >r
_]~$~< |+!/$2*^__
 __________
___________ _
_____ ____ _______
_____ ____~ ______

It basically generates random bits that make up a binary number and then converts that random binary number to decimal.
The rest is just the regular points in the square approach.
Usage: when you run the code you must make sure to prepopulate the stack (-v in python interpreter) with the number of samples, for example
pi.fish -v 1000

returns
3.164

Javascript: 62 Bytes
for(r=Math.random,t=c=8e4;t--;c-=r()**2+r()**2|0);alert(c/2e4)

I used the previous (now deleted) javascript answer and shaved 5 bytes.

You could link to cfern's answer to properly give credit.

Your answer appears to be a snippet which is disallowed I/O. Please either fix or delete your post.

Sorry, I'm new I didn't know how to put the link to the previous solution which now appears have been deleted. Regarding the snippet: I completely agree, but that was the code of the previous javascript solution which I also think it was invalid that reason.
I modified mine to be a program.

Yes; the previous answer was deleted as it was invalid -- I saw your answer before I recommended deletion, thus the comment. +1 for submitting a valid answer; welcome to PPCG!

Python 2, 90 85 81 bytes
from random import*;r=random;print sum(4.for i in[0]*9**7if r()**2+r()**2<1)/9**7

returns 3.14120037157 for example. The sample count is 4782969 (9^7). You can get a better pi with 9^9 but you will have to be patient.

You can save 3 by replacing range(9**7) with [0]*9**7 or something, since you don't use i. And the list isn't too long to run into memory issues.

Thanks. I wanted to get rid of range() but I had completely forgotten that trick.

I have a feeling [0]9**7 isn't valid syntax.

You are right. I have re-attached the lost asterix (it was under my desk).

GolfScript (34 chars)
0{^3?^rand.*^.*+/+}2000:^*`1/('.'@

Online demo
This uses fixed point because GS doesn't really have floating point. It slightly abuses the use of fixed point, so if you want to change the iteration count make sure that it's twice a power of ten.
Credit to xnor for the particular Monte Carlo method employed.

Ruby, 39 bytes
p (1..8e5).count{rand**2+rand**2<1}/2e5

One of the highlight is that this one is able to use 8e5 notation, which makes it extendable up to ~8e9 samples in the same program byte count.

Perl 6, 33 bytes
{4*sum((1>rand²+rand²)xx$_)/$_}

Try it online!
This is a function that takes the number of samples as an argument.

Scala, 87 77 66 bytes
def s=math.pow(math.random,2);Seq.fill(1000)(s+s).count(_<1)/250d

If you replace 1000 with 8000 and 250d with 2e4 you both save a byte and increase the number of samples by a factor of 8.

dc, 59 characters (whitespace is ignored)
[? 2^ ? 2^ + 1>i]su
[lx 1+ sx]si
[lu x lm 1+ d sm ln>z]sz

5k
?sn
lzx
lx ln / 4* p
q

I tested this on Plan 9 and OpenBSD, so I imagine it will work on Linux (GNU?) dc.
Explanation by line:

Stores code to [read and square two floats; execute register i if 1 is greater than the sum of their squares] in register u.
Stores code to [increment register x by 1] in register i.
Stores code to [execute register u, increment register m, and then execute register z if register m is greater than register n] in register z.

Set the scale to 5 decimal points.
Read the number of points to sample from the first line of input.
Execute register z.
Divide register x (the number of hits) by register n (the number of points), multiply the result by 4, and print.
Quit.

However, I cheated:
The program needs a supply of random floats between 0 and 1.
/* frand.c */
#include <u.h>
#include <libc.h>

void
main(void)
{
    srand(time(0));

    for(;;)
        print("%f\n", frand());
}

Usage:
#!/bin/rc
# runpi <number of samples>

{ echo $1; frand } | dc pi.dc

Test run:
% runpi 10000
3.14840

Now with less cheating (100 bytes)
Someone pointed out that I could include a simple prng.
http://en.wikipedia.org/wiki/RANDU
[lrx2^lrx2^+1>i]su[lx1+sx]si[luxlm1+dsmln>z]sz[0kls65539*2 31^%dsslkk2 31^/]sr?sn5dksk1sslzxlxlm/4*p

Ungolfed
[
Registers:
u - routine : execute i if sum of squares less than 1
i - routine : increment register x
z - routine : iterator - execute u while n > m++
r - routine : RANDU PRNG
m - variable: number of samples
x - variable: number of samples inside circle
s - variable: seed for r
k - variable: scale for division
n - variable: number of iterations (user input)
]c
[lrx 2^ lrx 2^ + 1>i]su
[lx 1+ sx]si
[lu x lm 1+ d sm ln>z]sz
[0k ls 65539 * 2 31^ % d ss lkk 2 31 ^ /]sr
? sn
5dksk
1 ss
lzx
lx lm / 4*
p

Test run:
$ echo 10000 | dc pigolf.dc
3.13640

Joe, 20 19 bytes
Note: this answer is non-competing, because version 0.1.2, which added randomness, was released after this challenge.
Named function F:
:%$,(4*/+1>/+*,?2~;

Unnamed function:
%$,(4*/+1>/+*,?2~;)

These both take the sample count as an argument and return the result. How do they work?
%$,(4*/+1>/+*,?2~;)
   (4*/+1>/+*,?2~;) defines a chain, where functions are called right-to-left
               2~;  appends 2 to the argument, giving [x, 2]
              ?     create a table of random values from 0 to 1 with that shape
            *,      take square of every value
          /+         sum rows, giving a list of (x**2+y**2) values
        1>           check if a value is less than 1, per atom
      /+             sum the results
    4*               multiply by four
%$,                  divide the result by the original parameter

Example runs:
   :%$,(4*/+1>/+*,?2~;
   F400000
3.14154
   F400000
3.14302

Pure Bash, 65 bytes
for((;i++<$1*4;a+=RANDOM**2+RANDOM**2<32767**2));{ :;}
echo $a/$1

Takes a single command-line parameter which is multiplied by 4 to give the number of samples.  Bash arithmetic is integer-only, so a rational is output.  This may be piped to bc -l for the final division:
$ ./montepi.sh 10000
31477/10000
$ ./montepi.sh 10000|bc -l
3.13410000000000000000
$

Pyth, 19
c*4sm<sm^OQ2 2*QQQQ

Give the desired number of iterations as input.
Demonstration
Since Pyth doesn't have a "Random floating number" function, I had to improvise. The program choose two random positive integers less than the input, squares, sums, and compared to the input squared. This performed a number of times equal to the input, then the result is multiplied by 4 and divided by the input.
In related news, I will be adding a random floating point number operation to Pyth shortly. This program does not use that feature, however.

If we interpret "The result may be returned or printed as a floating point, fixed point or rational number." liberally, then printing the numerator and denominator of the resulting fraction should be sufficient. In that case:
Pyth, 18
*4sm<sm^OQ2 2*QQQQ

This is an identical program, with the floating point division operation (c) removed.

Julia, 37 bytes
4mean(1-floor(sum(rand(4^8,2).^2,2)))

The number of samples is 65536 (=4^8).
A sligthly  longer variant: a function with number of samples s as the only argument:
s->4mean(1-floor(sum(rand(s,2).^2,2)))

C, 130 bytes
#include<stdlib.h>f(){double x,y,c=0;for(int i=0;i<8e6;++i)x=rand(),y=rand(),c+=x*x+y*y<1.0*RAND_MAX*RAND_MAX;printf("%f",c/2e6);}

Ungolfed:
#include <stdlib.h>
f(){
 double x,y,c=0;
 for(int i=0; i<8e6; ++i) x=rand(), y=rand(), c+=x*x+y*y<1.0*RAND_MAX*RAND_MAX;
 printf("%f",c/2e6);
}

of course, you should still probably post the version without the whitespace (keep the current version with the header "ungolfed/with whitespace" or something)

@DestructibleWatermelon done!

The solution doesn't work in GCC without a newline before f(). What compiler did you use?
See https://tio.run/##Pc49C4JAHIDx3U9xGMG9ZdYgwWkgtNbQ1BZ6L/UHO8M07hA/@2VDbb/p4ZGrm5QhLMDKZlA6f/WqgTq57w0mo2qHutHIcc9lkQrTdhhsj2A25DudCcaAuKKrrMKE@x8kKxx1zFOfb5KUnsvT4XosL3@IZzdnDI6XJuZyvdUZEVP4lh8VWPxuQRE0RgjNDyKawgc

101 bytes

Racket 63 bytes
Using method of R language answer by @Matt :
(/(for/sum((i n))(define a(/(random 11)10))(/ 4(+ 1(* a a))))n)

Ungolfed:
(define(f n)
   (/
    (for/sum ((i n))
      (define a (/(random 11)10))
      (/ 4(+ 1(* a a))))
    n))

Testing:
(f 10000)

Output (fraction):
3 31491308966059784/243801776017028125

As decimal:
(exact->inexact(f 10000))

3.13583200307849

Haskell, 116 114 110 96 bytes
d=8^9
g[a,b]=sum[4|a*a+b*b<d*d]
p n=(sum.take(floor n)$g<$>iterate((\x->mod(9*x+1)d)<$>)[0,6])/n

Because dealing with import System.Random; r=randoms(mkStdGen 2) would take too many precious bytes, I generate an infinite list of random numbers with the linear congruential generator that some say is almost cryptographically strong: x↦x*9+1 mod 8^9, which by the Hull-Dobell Theorem has the full period of 8^9.
g yields 4 if the random number point is inside the circle for pairs of random numbers in [0..8^9-1] because this eliminates a multiplication in the formula used.
Usage:
> p 100000
3.14208

Try it online!

Actually, 14 bytes
`G²G²+1>`nkæ4*

Try it online!
The number of samples is given as input (rather than hard-coded).
Explanation:
`G²G²+1>`nkæ4*
`G²G²+1>`n      do the following N times:
 G²G²+            rand()**2 + rand()**2
      1>          is 1 greater?
          kæ    mean of results
            4*  multiply by 4

Why the downvote?

Japt, 26 or 18 bytes
o r_+ÂMhMr p +Mr p <1Ã*4/U

Try it online!
Analogous to Optimizer's answer, mainly just trying to learn Japt.
Takes the number of iterations to run for as the implicit input U.
o                           Take the input and turn it into a range [0, U),
                            essentially a cheap way to get a large array.
  r_                        Reduce it with the default initial value of 0.
    +Â                      On each iteration, add one if
      MhMr p +Mr p          the hypotenuse of a random [0,1)x[0,1) right triangle
                   <1       is smaller than one.
                     Ã*4/U  Multiply the whole result by four and divide by input.

If 1/(1+x^2) is allowed (instead of two separate randoms), then we can achieve 18 bytes with the same logic.
o r_Ä/(1+Mr pÃ*4/U

You can save a few bytes by letting Mh calculate the hypotenuse rather than doing it yourself ;-) Also, you can use x to take the sum of an array, rather than reducing by addition: o x@MhMr Mr)<1Ã*4/U

@ETHproductions Neat, I didn't know you can use Mh like that, thanks! Your two-random answer is nearly as short as my answer with only one random, that's pretty cool. I'll keep x in mind, I tend to use reduction a lot when trying to golf stuff, so this will come in very handy.

F#, 149 bytes
open System;
let r=new Random()
let q()=
 let b=r.NextDouble()
 b*b
let m(s:float)=(s-Seq.sumBy(fun x->q()+q()|>Math.Sqrt|>Math.Floor)[1.0..s])*4.0/s

Try it online!
As far as I can make out, to do this kind of running total in F# it's shorter to create an array of numbers and use the Seq.sumBy method than use a for..to..do block.
What this code does it that it creates a collection of floating-point numbers from 1 to s, performs the function fun x->... for the number of elements in the collection, and sums the result. There are s elements in the collection, so the random test is done s times. The actual numbers in the collection are ignored (fun x->, but x is not used).
It also means that the application has to first create and fill the array, and then iterate over it. So it's probably twice as slow as a for..to..do loop. And with the array creation memory usage is in the region of O(f**k)!
For the actual test itself, instead of using an if then else statement what it does is it calculates the distance (q()+q()|>Math.Sqrt) and rounds it down with Math.Floor. If the distance is within the circle, it'll be rounded down to 0. If the distance is outside the circle it'll be rounded down to 1. The Seq.sumBy method will then total these results.
Note then that what Seq.sumBy has totalled is not the points inside the circle, but the points outside it. So for the result it takes s (our sample size) and subtracts the total from it.
It also appears that taking a sample size as a parameter is shorter than hard-coding the value. So I am cheating a little bit...

Perl 5, 34 bytes
$_=$a/map$a+=4/(1+(rand)**2),1..$_

The number of samples is taken from stdin. Requires -p.
Works because:

Try it online!

Fortran (GFortran), 84 83 bytes
CALL SRAND(0)
DO I=1,4E3
X=RAND()
Y=RAND()
IF(X*X+Y*Y<1)A=A+1E-3
ENDDO
PRINT*,A
END

Try it online!
This code is very bad wrote. It will fail if gfortran decide to initialize variable A with other value then 0 (which occurs 50% of the compilations, approximately) and, if A is initialized as 0, it will always generate the same random sequence for the given seed. Then, the same value for Pi is printed always.
This is a much better program:
Fortran (GFortran), 100 99 bytes
A=0
DO I=1,4E3
CALL RANDOM_NUMBER(X)
CALL RANDOM_NUMBER(Y)
IF(X*X+Y*Y<1)A=A+1E-3
ENDDO
PRINT*,A
END

Try it online!
(One byte saved in each version; thanks Penguino).

In each version you can save a byte by changing 'DO I=1,1E3' to  'DO I=1,4E3', changing 'A=A+1' to 'A=A+1E-3', and changing 'PRINT*,A/250' to 'PRINT*,A'

Yes, you are sure! Thank for the suggestion!

MathGolf, 11 bytes
0♫ô4ƒ²)/+♫/

Try it online!
Explanation
0            Push 0 (used for the sum)
 ♫           Push 10000
  ô          Start block of length 6
   4         Push 4
    ƒ²       Push random float in range [0,1) and square
      )      Increment by 1
       /     Divide
        +    Add result to the total sum (last instruction of block)
         ♫/  Divide by 10000

I'll add an explanation tomorrow. It performs 10000 iterations using the same integral as many other languages.
Assuming number of iterations as input:
MathGolf, 10 bytes
ô4ƒ²)/+\/(

Try it online!
Explanation
ô           Start block of length 6
 4          Push 4
  ƒ²        Push random float in range [0,1) and square
    )       Increment by 1
     /      Divide
      +     Add result to the total sum (last instruction of block)
            The first iteration, this adds the result to the implicit input, meaning that if the input is 10000 and the result is 3.1, 10003.1 is pushed on the stack. The end result is 10000 too large.
       \    Swap the top two elements on the stack (pushes the implicit input to the top)
        /   Divide (Since the sum is wonky, this is around 4.14159...)
         (  Decrement by one, resulting in 3.14159...)

APL(NARS), 79 chars,  158 bytes
r←q;f;i
   f←{1e8÷⍨?1e8}⋄i←r←0
A: i←i+1⋄r←r+1>(f*2)+f*2⋄→A×⍳i<1e5
   r←4×r÷1e5

%pi=3.1415926535897932385... test:
  q
3.14284
  q
3.14256
  q
3.13832

Go, 130 116 bytes
import."math/rand";func p(){r,c,i:=Float64,0,1e6;for;i>0;i--{a,b:=r(),r();if a*a+b*b<1{c+=4}};print(float64(c)/1e6)}

No real tricks here, other than aliasing rand.Float64 and float64 to save a few bytes.
Run it here: http://play.golang.org/p/F5YniIncor
Please note that random on play.golang.org cannot be seeded by time (since the date is always the same) so you need to change the seed manually.
1e6 is used since the run time is too long for play.golang.org if I set it higher.
Edit:
No longer takes an input, outputs the computed value itself.

Python 3, 69 bytes
Apparently shorter than the other python entries, save one that uses sympy numpy (misremembered)
from random import*
print(eval(('+(1>('+2*'+random()**2)')*2000)/500)

Try it online!
Wow, it's been a while since I posted a multiline answer on this site... I've been using Turtlèd too much for that
Explanation:
from random import*    We can use random.random now as random
print(                 Print the result of the following
      eval(            evaluate the value of the expression in the string
           ('+(1>('+2*'+random()**2)')*2000    the string construction
                                           )/500) divide by 500 and print

Explanation of the string expression:
                             *2000  repeat this 2000 times
    '+(1>('                         have each section joined by a +,
                                    and the first item have the unary +. (nop, essentially)
                                    check if the section is less than 1, if so,
                                    evaluate to 1, else 0
           +2*'+random()**2)'        sum two random numbers squared, close both parens.

                                   [Evaluate the sum]

This is about the max accuracy, because 3000 samples gives a recursion error.
(I'll improve this explanation real soon actually.)

https://docs.python.org/3/library/sys.html#sys.setrecursionlimit

many bytes for that!

Yeah, just saying it isn't the max accuracy

There seems to be a discrepancy between your TIO link's source and your post's source -- namely parentheses around the print.

Pyt, 22 bytes
Đ0⇹`⇹ṛ²ṛ²+1≤+⇹⁻łŕ⇹ᵮ/4*

Try it online!
Explanation:
          Implicit input [the number of random points to be sampled]
Đ         Duplicate input
0         Push zero
⇹         Swap the top two items on the stack
`...(ł)   Do ... (while the top of the stack is not 0)
⇹         Swap the top two items on the stack
ṛ²ṛ²+     Add the squares of two random floats (between 0 and 1)
1≤        Is the sum less than or equal to 1?
+         Add the top two items on the stack [auto-converts boolean to int]
⇹         Swap the top two items on the stack
⁻         Decrement the top of the stack by 1
(`)...ł   (Do ...) while the top of the stack is not 0
ŕ         Remove the top of the stack
⇹         Swap the top two items on the stack
ᵮ         Cast the top of the stack to a float
/         Divide
4*        Multiply by 4
          Implicit output

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


text:
Sort all the pixels of an image by the number of occurrences
Input
The name of a file in the raster graphics format of your choice. The chosen format must support at least 8 bits per channel and 3 channels.
Output
A file in the same format, with the same dimensions and pixels as the first, but whose pixels are grouped in descending order of the number of times they occur, sorted from left-to-right, top-to-bottom.

If certain colours of pixels appear the same number of times, their order is unspecified.
You must not overwrite the input file (use a different file name for the output).
Any and all third-party image processing libraries are permitted.

Example

Will give output similar to:

Especially in the lower parts of the image some variation may occur, due to different tie breaking between colours of equal frequency.

Can we use other file formats, like PNM, for the output?

@EMBLEM I'm not sure what you mean with “forego grayscale.”

@MartinBüttner Good point, I've changed it to input and output being the same formats.

@FUZxxl doesn't matter anymore, I've changed it to input and output being the same formats.

Can we provide input using the C preprocessor?

@MartinBüttner No, and by votes.

@EMBLEM I think he meant “ties” as in “two pixels occur with the same frequency.” Can I assume that their order is unspecified in this case?

@MartinBüttner Unspecified.

@FUZxxl Sure, C Preprocessor is fine.

To be clear on input/output... can we use a string (file name), a stream of bytes, or a complex File-type object...? That would make a large difference in some languages.

@Bigtoes limiting to filename.

@MartinBüttner Agreed. Sorry. :(

If using PPM, do we have to allo for comments or not?

@MartinBüttner Because it couldn't handle transparency.

@MartinBüttner Oops. I'd intended that to be a requirement, but it never got in there. Accepted.

J, 94 81 bytes
A function taking the name of a PNG file (with no transparency channel) and writing the result in the input filename prepended by "o".
f=.3 :0
load'graphics/png'
(($a)$#/|:\:~(#,{.)/.~,a=.readpng y)writepng'o',y
)

Method

Create a list of every color-number and the number of it's appearances
Sort list by number of it's appearances
Multiplicate every color-number equal it's appearances
Reshape the new list to the shape of the original image

you can save one character if you make the function dyadic and save the sorted PNG into x.

@FUZxxl It would actually save 3 chars but I don't think that's in line with the specification. I could also save the image to a 1-letter extensionless file but that's not too elegant.

@FUZxxl Well, turned out there was quite some extra chars in the body of the function...

I don't speak J, but in the specification there is a string as input. Here the input file name is hardcoded in the program.

@edc65 In J if you define a function it's argument is always accessed by the y character. (If you define a function with two arguments those are accessed by x and y and you cannot define functions with more arguments.)

Mathematica, 125 123 bytes
Export["a"<>#,i=ImageData@Import@#;Image@ArrayReshape[Join@@ConstantArray@@@SortBy[Tally[Join@@i],-Last@#&],Dimensions@i]]&

This defines an unnamed function which takes the file name in any common image format and writes the result to the a file with the same name but prepended with a. The result looks a bit different from the OP's, since Mathematica's SortBy breaks ties by default sort order, so the bottom bits where loads of ties occur look a bit neater:

The implementation itself is really straightforward:

ImageData to get a grid of colour values.
Join to flatten the array.
Tally to count the occurrence of each colour.
SortBy[...,-Last@#&] to sort by frequencies from highest to lowest.
ConstantArray and Join to expand the tallies again.
ArrayReshape to recover the original image shape (obtained with Dimensions).
Image to convert the data back to an image object.

FYI, 22 bytes are used on file I/O. An equivalent version that takes and returns an image object comes in at 103 bytes:
(i=ImageData@#;Image@ArrayReshape[Join@@ConstantArray@@@SortBy[Tally[Join@@i],-Last@#&],Dimensions@i])&

Python2/PIL, 244 226 225 223 222 202 186 182 170 159
from PIL.Image import*
s=raw_input();i=open(s);g=list(i.getdata());i.putdata(sum([[c[1]]*-c[0]for c in sorted((-g.count(r),r)for r in set(g))],[]));i.save(2*s)

Changelog

Edit 1: saved 18 bytes with new algorithm.

Edit 2: saved 1 byte with remembering that Python 2 has int division without from __future__ import division. Whoops.

Edit 3: saved 2 bytes by swapping cmp for key.

Edit 4: saved 1 byte by clever unpacking. w=i.size[0] -> w,h=i.size

Edit 5: moved list() to an earlier position, drastically improving runtime. No length change.

Edit 6: saved 20 bytes by remembering that there is also sequence sorting in Python. Probably slower, but also produces images more like the the Mathematica one.

Edit 7: saved 16 bytes by changing putpixel to a list and putdata. Under 200!

Edit 8: saved 4 bytes by changing range(...) to ...*[0] since we don't need the values.

Edit 9: saved 12 bytes by converting the for loop to a sum.

Edit 10: saved 11 bytes by noticing that the image size doesn't matter anymore (since edit 7). Whoops v2.0.

Shorter version by stokastic, 123
from PIL.Image import*
s=raw_input();i=open(s);d=list(i.getdata());i.putdata(sorted(d,key=lambda D:d.count(D)));i.save(2*s)

Well, let's at least try, even though it's already beaten.
It's extremely slow, Panda processed for several minutes on my laptop.
Saves with a filename with the original filename repeated twice.

you can save 32 bytes by using i=open(raw_input());d=list(i.getdata());i.putdata(sorted(d,key=lambda D:d.count(D)));i.save('o.png'), although it is very very slow for large images (calls list.count on each pixel).

@stokastic Well, I'll add that but credit you for it. I prefer reading the name to a variable (what if the user inputs an image called o.png?)

fair enough! :)

Python, 1197 bytes
# -*- coding: utf-8 -*-
from __future__ import division
from collections import Counter
import png
from sys import argv,exit
script,file_name,output_file_name=argv
def group(s,n):
    return zip(*[iter(s)]*n)
def flatten(list_of_lists):
    result=[]
    for lst in list_of_lists:
        result.extend(lst)
    return result
def group_list(old_list,tuples_per_list):
    new_list=[]
    i=1
    appended_row=[]
    for item in old_list:
        appended_row.extend(flatten([item]))
        if i==tuples_per_list:
            new_list.append(appended_row)
            i=1
            appended_row=[]
        else:
            i+=1
    return new_list
input_image=png.Reader(file_name)
image_data=input_image.read()
width=image_data[0]
height=image_data[1]
pixels=list(image_data[2])
if image_data[3]["alpha"]:
    ints_per_colour=4
elif image_data[3]["greyscale"]:
    ints_per_colour=2
else:
    ints_per_colour=3
colours=Counter(colour for row in pixels for colour in group(row,ints_per_colour)).most_common()
pixel_list=flatten([element]*count for element,count in colours)
ordered_rows=group_list(pixel_list,width)
f=open(output_file_name,"wb")
png_object=png.Writer(width,height,greyscale=image_data[3]["greyscale"],alpha=image_data[3]["alpha"])
png_object.write(f,ordered_rows)
f.close()

The png module I used.

Thank you for providing a reference solution! We usually put these into the challenge post itself so people can find them more easily, but posting them as an answer is okay, too.

C# 413
Complete program. Pass the file name on command line, the output is saved in the same format as file "o".
Not using some neat features of linq, like SelectMany and Enumerable.Range, as the program would be cleaner but longer.
using System.Collections.Generic;using System.Linq;using System.Drawing;
class P{
static void Main(string[]a){
var d=new Dictionary<Color,int>();var b=new Bitmap(a[0]);
int n,x,y,w=b.Width,h=b.Height;
for (x=w;x-->0;)for(y=h;y-->0;){var p=b.GetPixel(x,y);d[p]=d.ContainsKey(p)?d[p]+1:1;}
y=h;foreach(var q in d.OrderBy(v=>v.Value)){for(n=q.Value;n-->0;){
if(x<=0){x=w;--y;}b.SetPixel(--x, y, q.Key);}}b.Save("o");
}}

Readable formatting courtesy of VS2010
using System.Collections.Generic;
using System.Linq;
using System.Drawing;

class P
{
    static void Main(string[] a)
    {
        var d = new Dictionary<Color, int>();
        var b = new Bitmap(a[0]);
        int n,x,y,w = b.Width, h=b.Height;

        for (x = w; x-- > 0;)
            for (y = h; y-- > 0;)
            {
                var p = b.GetPixel(x, y);
                d[p] = d.ContainsKey(p) ? d[p]+1 : 1;
            }
        y = h;
        foreach (var q in d.OrderBy(v => v.Value))
        {
            for (n = q.Value; n-- > 0; )
            {
                if (x <= 0)
                {
                    x = w;
                    --y;
                 }
                 b.SetPixel(--x, y, q.Key);
            }
        }
        b.Save(a[0]+".");
    }
}

May be I am wrong but I think you can remove the newlines after { chars

And before the } chars

I quess something like https://tio.run/##VY/Ba8IwGMX/leApYU3R62IcOMHBNiYI20E8xDTaD2rSJrFtKP7tXaLbwVve9x7v9yIdlcaqcbw40Ce0Dc6rc/5qqkpJD0a7fK20siDZQ@ADdPN4WVnRRclkJZxDm8F54UGi1kCBPgVo7LyN/m4vyNAKiwquVYdWcKMIG@aRaWwG2i8wYSlxuCWW4M@ixmI33RMWXaSzPgtZxw/5DxS@zMr4elNwKj07Got73rGe0sWUkSQDL1m4yxu2jum18hvoVYVjEWHFrt7zIn5Z@zjTvauAa/KSrk@z5xm7pobYpIQscWpoEGhU5F@2UHYZcMsXbf4tqosiZEhEzZu7ZvqPC0fcz/mUDGkbpYFdD/n2fwOlfYZChpo8kgm7Jk@0Ck/MJKlx/AU would work!

Python 2 : 191 bytes
Here's my attempt. I figured I might be able to save some space by using Counter, but it didn't end up as small as Pietu1998's answer.
from collections import Counter
from PIL import Image
i=Image.open(raw_input())
s=reduce(lambda x,y:x+y,map(lambda(a,b):[a]*b,Counter(i.getdata()).most_common()))
i.putdata(s)
i.save("o.png")

Output from Panda