Square, Circle, Triangle, ...Gear? Using Algodoo and Paint I made these six 300×300 monochromatic images of four convenient shapes:
     
This class of images has the following properties:

They are always 300×300 pixels, monochromatic (black and white only), and have exactly four white regions that correspond to a square, a circle, a triangle, and a gear.
The shapes never overlap or touch each other, nor do they touch the image border or go out of bounds.
The shapes always have the same size, but they may be rotated and positioned in any way.

(The shapes also have equal areas, though when rastered like this their pixel counts are unlikely to be exactly equivalent.)
Challenge
Write the shortest program or function possible that takes in the filename of such an image and turns all the white pixels...

red (255, 0, 0) if they are in the square.
blue (0, 0, 255) if they are in the circle.
green (0, 255, 0) if they are in the triangle.
yellow (255, 255, 0) if they are in the gear.

e.g.

Details
Your program should work for effectively all possible input images. (Only valid 300×300 monochromatic images will be input.) The six images I have provided are merely examples, you may not hardcode their output into your program.
You may not use computer vision libraries or functions, built-in or external. The point is to do this using your own pixel-level operations. You may use image libraries that simply let you open and modify images (e.g. PIL for Python).
You may use any common lossless image file formats for input and output as long as you stick to the color scheme.
You can take in the image filename as a function argument, from stdin, or from the command line. The output image can be saved to a new file, the same file, or simply displayed.
Scoring
The submission with the fewest bytes wins. I may test submissions with additional images to determine their validity.
Mathematica, 459 392 bytes

f=(d=ImageData@Import@#/.{a_,_,_}:&gt;a;(For[a={};b={#&amp;@@d~Position~1},b!={},c=#&amp;@@b;b=Rest@b;d[[##&amp;@@c]]=0;a~AppendTo~c;If[Extract[d,c+#]==1,b=b⋃{c+#}]&amp;/@{e={1,0},-e,e={0,1},-e}];m=1.Mean@a;m=#-m&amp;/@a;n=Count[Partition[Norm/@SortBy[m,ArcTan@@#&amp;],300,1,1],l_/;l[[150]]==Max@l];(d[[##&amp;@@#]]=Round[n^.68])&amp;/@a)&amp;/@Range@4;Image[d/.n_Integer:&gt;{{0,0,0},,{0,1,0},{1,0,0},,,,{1,1,0},{0,0,1}}[[n+1]]])&amp;


Ungolfed:

f = (
 d = ImageData@Import@# /. {a_, _, _} :&gt; a;
 (
    For[a = {}; b = {# &amp; @@ d~Position~1},
     b != {},
     c = # &amp; @@ b;
     b = Rest@b;
     d[[## &amp; @@ c]] = 0;
     a~AppendTo~c;
     If[Extract[d, c + #] == 1, 
        b = b ⋃ {c + #}] &amp; /@ {e = {1, 0}, -e, e = {0, 1}, -e}
     ];
    m = 1. Mean@a; m = # - m &amp; /@ a;
    n = 
     Count[Partition[Norm /@ SortBy[m, ArcTan @@ # &amp;], 300, 1, 1], 
      l_ /; l[[150]] == Max@l];
    (d[[## &amp; @@ #]] = Round[n^.68]) &amp; /@ a
    ) &amp; /@ Range@4;
 Image[d /. 
   n_Integer :&gt; {{0, 0, 0}, , {0, 1, 0}, {1, 0, 0}, , , , {1, 1, 
       0}, {0, 0, 1}}[[n + 1]]]
) &amp;


I could save 6 more bytes by turning m=1.Mean@a;m=#-m&amp;/@a; into m=#-Mean@a&amp;/@a;, but that significantly blows up execution time, which is annoying for testing. (Note, that this is two optimisations: pulling out the computation of Mean@a out of loop and using exact symbolic types instead of floating point numbers. Interestingly, the use of exact types is a lot more significant than computing the mean in every iteration.)

So this is approach number three:


Detect areas by flood-fill.
Find approximate centre of each area by averaging all pixel coordinates.
Now for all pixels in the shape let's plot distance from vs angle to that centre:



The triangle has 3 clear maxima, the square 4, the gear 16, and the circle has tons, due to aliasing fluctuations about the constant radius.
We find the number of maxima by looking at slices of 300 pixels (ordered by angle), and count the slices where the pixel at position 150 is the maximum.
Then we just colour all pixels depending on the number of peaks (the circle is anything over 16, and usually yields around 20 peaks, due to the size of the slices).


Just for the record, if I use Ell's idea, and simply sort the regions by largest distance between any pixel and centre, I can do this in 342 bytes:

f=(d=ImageData@Import@#/.{a_,_,_}:&gt;a;MapIndexed[(d[[##&amp;@@#]]=#&amp;@@#2)&amp;,SortBy[(For[a={};b={#&amp;@@d~Position~1},b!={},c=#&amp;@@b;b=Rest@b;d[[##&amp;@@c]]=0;a~AppendTo~c;If[Extract[d,c+#]==1,b=b⋃{c+#}]&amp;/@{e={1,0},-e,e={0,1},-e}];a)&amp;/@Range@4,(m=Mean@#;Max[1.Norm[#-m]&amp;/@#])&amp;],{2}];Image[d/.n_Integer:&gt;{{0,0,0},{0,0,1},{1,1,0},{1,0,0},{0,1,0}}[[n+1]]])&amp;


But I do not intend to compete with that, as long as everyone else is using their own original algorithms, instead of golfing down those of others.