[Question]
  [
Using [Algodoo](http://www.algodoo.com/) and Paint I made these six 300×300 monochromatic images of four convenient shapes:
![Image 1](https://i.stack.imgur.com/Oa2O1.png) ![Image 2](https://i.stack.imgur.com/C6IKR.png) ![Image 3](https://i.stack.imgur.com/YZfc6.png) ![Image 4](https://i.stack.imgur.com/F1ZDM.png) ![Image 5](https://i.stack.imgur.com/chJFi.png) ![Image 6](https://i.stack.imgur.com/DLM3y.png)
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.
![Image 1 colored in](https://i.stack.imgur.com/EIRiN.png)
# 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](https://mothereff.in/byte-counter) wins. I may test submissions with additional images to determine their validity.
  
[Answer]
## Mathematica, ~~459~~ 392 bytes
```
f=(d=ImageData@Import@#/.{a_,_,_}:>a;(For[a={};b={#&@@d~Position~1},b!={},c=#&@@b;b=Rest@b;d[[##&@@c]]=0;a~AppendTo~c;If[Extract[d,c+#]==1,b=b⋃{c+#}]&/@{e={1,0},-e,e={0,1},-e}];m=1.Mean@a;m=#-m&/@a;n=Count[Partition[Norm/@SortBy[m,ArcTan@@#&],300,1,1],l_/;l[[150]]==Max@l];(d[[##&@@#]]=Round[n^.68])&/@a)&/@Range@4;Image[d/.n_Integer:>{{0,0,0},,{0,1,0},{1,0,0},,,,{1,1,0},{0,0,1}}[[n+1]]])&
```
Ungolfed:
```
f = (
 d = ImageData@Import@# /. {a_, _, _} :> a;
 (
    For[a = {}; b = {# & @@ d~Position~1},
     b != {},
     c = # & @@ b;
     b = Rest@b;
     d[[## & @@ c]] = 0;
     a~AppendTo~c;
     If[Extract[d, c + #] == 1, 
        b = b ⋃ {c + #}] & /@ {e = {1, 0}, -e, e = {0, 1}, -e}
     ];
    m = 1. Mean@a; m = # - m & /@ a;
    n = 
     Count[Partition[Norm /@ SortBy[m, ArcTan @@ # &], 300, 1, 1], 
      l_ /; l[[150]] == Max@l];
    (d[[## & @@ #]] = Round[n^.68]) & /@ a
    ) & /@ Range@4;
 Image[d /. 
   n_Integer :> {{0, 0, 0}, , {0, 1, 0}, {1, 0, 0}, , , , {1, 1, 
       0}, {0, 0, 1}}[[n + 1]]]
) &
```
I could save 6 more bytes by turning `m=1.Mean@a;m=#-m&/@a;` into `m=#-Mean@a&/@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:
![enter image description here](https://i.stack.imgur.com/DHpbr.png)
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_,_,_}:>a;MapIndexed[(d[[##&@@#]]=#&@@#2)&,SortBy[(For[a={};b={#&@@d~Position~1},b!={},c=#&@@b;b=Rest@b;d[[##&@@c]]=0;a~AppendTo~c;If[Extract[d,c+#]==1,b=b⋃{c+#}]&/@{e={1,0},-e,e={0,1},-e}];a)&/@Range@4,(m=Mean@#;Max[1.Norm[#-m]&/@#])&],{2}];Image[d/.n_Integer:>{{0,0,0},{0,0,1},{1,1,0},{1,0,0},{0,1,0}}[[n+1]]])&
```
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.
[Answer]
# Java, ~~1204~~ ~~1132~~ ~~1087~~ 1076
Just to prove to myself that ***I*** can do this.
I included imports right next to the function declarations; these would have to be outside the class for this to work:
```
import java.awt.*;import java.awt.image.*;import java.io.*;import java.util.*;import javax.imageio.*;
BufferedImage i;Set<Point>Q;void p(String a)throws Exception{i=new BufferedImage(302,302,1);i.getGraphics().drawImage(ImageIO.read(new File(a)),1,1,null);Set<Set<Point>>S=new HashSet<>();for(int y=0;y<300;y++){for(int x=0;x<300;x++){if(!G(x,y)){Point p=new Point(x,y);Q=new HashSet<>();if(!S.stream().anyMatch(s->s.contains(p)))S.add(f(x,y));}}}Object[]o=S.stream().sorted((p,P)->c(p)-c(P)).toArray();s(o[0],255);s(o[1],255<<16);s(o[2],0xFF00);s(o[3],0xFFFF00);ImageIO.write(i.getSubimage(1,1,300,300),"png",new File(a));}boolean G(int x,int y){return i.getRGB(x,y)!=-1;}Set<Point>f(int x,int y){Point p=new Point(x,y);if(!Q.contains(p)&&!G(x,y)){Q.add(p);f(x-1,y);f(x+1,y);f(x,y-1);f(x,y+1);}return Q;}int c(Set<Point>s){return(int)s.stream().filter(p->G(p.x-2,p.y-1)||G(p.x-2,p.y+1)||G(p.x+1,p.y-2)||G(p.x-1,p.y-2)||G(p.x+2,p.y-1)||G(p.x+2,p.y+1)||G(p.x+1,p.y+2)||G(p.x-1,p.y+2)).count();}void s(Object o,int c){((Set<Point>)o).stream().forEach(p->{i.setRGB(p.x,p.y,c);});}
```
Ungolfed (and runnable; ie boilerplate added):
```
import java.awt.Point;
import java.awt.image.BufferedImage;
import java.io.File;
import java.util.HashSet;
import java.util.Set;
import javax.imageio.ImageIO;
public class SquareCircleTriangleGear {
    public static void main(String[]args){
        try {
            new SquareCircleTriangleGear().p("filepath");
        } catch (Exception ex) {
        }
    }
    BufferedImage i;
    Set<Point>Q;
    void p(String a)throws Exception{
        i = new BufferedImage(302,302,BufferedImage.TYPE_INT_RGB);
        i.getGraphics().drawImage(ImageIO.read(new File(a)),1,1,null);
        Set<Set<Point>>set=new HashSet<>();
        for(int y=0;y<300;y++){
            for(int x = 0;x<300;x++){
                if(i.getRGB(x,y)==-1){
                    Point p = new Point(x,y);
                    Q=new HashSet<>();
                    if(!set.stream().anyMatch((s)->s.contains(p))){
                        set.add(fill(x,y));
                    }
                }
            }
        }
        Object[]o=set.stream().sorted((p,P)->c(p)-c(P)).toArray();
        s(o[0],0x0000FF);
        s(o[1],0xFF0000);
        s(o[2],0x00FF00);
        s(o[3],0xFFFF00);
        ImageIO.write(i.getSubImage(1,1,300,300), "png", new File(a));
    }
    Set<Point>fill(int x, int y){
        Point p=new Point(x,y);
        if(!Q.contains(p)&&!i.getRGB(x,y)!=-1) {
        Q.add(p);
            fill(x-1,y);
            fill(x+1,y);
            fill(x,y-1);
            fill(x,y+1);
        }
        return Q;
    }
    int c(Set<Point>s){return (int)s.stream().filter(p->isBoundary(p.x,p.y)).count();}
    boolean isBoundary(int x, int y){
        return i.getRGB(x-2,y-1)!=-1||i.getRGB(x-2,y+1)!=-1||i.getRGB(x+1,y-2)!=-1||
               i.getRGB(x-1,y-2)!=-1||i.getRGB(x+2,y-1)!=-1||i.getRGB(x+2,y+1)!=-1||
               i.getRGB(x+1,y+2)!=-1||i.getRGB(x-1,y+2)!=-1;
    }
    void s(Object o,int c){
        ((Set<Point>)o).stream().forEach(p->{i.setRGB(p.x,p.y,c);});
    }
}
```
This works by iterating over every pixel of the image and flood-filling each time we reach a "hole". We add each flood-fill result as a `Set<Point>` to a `Set`. Then we determine which shape is which. This is done by looking at the number of boundary pixels of the shape. I defined boundary as a knight's move away from a black tile, since that would stay more constant between rotations and such. When we do this, it becomes clear that the shapes can be sorted by that value: Circle, Square, Triangle, Gear. So I sort and set all pixels of that shape to the correct color.
Note that the image I'm writing to is not directly taken from the file, because if I were to do that, Java would treat the image as black and white and filling with colors wouldn't work. So I have to create my own image with `TYPE_INT_RGB` (which is `1`). Also note that the image I'm doing work on is `302` by `302`; this is so that the Knight's distance algorithm doesn't need to worry about attempting to read out-of-bounds on the image. I fix this discrepancy in size by calling `i.getSubImage(1,1,300,300)`. Note: I may have forgotten to fix this when I uploaded the images, in which cases the images are 2 pixels too wide, but except for this fact, they should be correct
The function will overwrite the file whose path is passed in. Outputs:
![enter image description here](https://i.stack.imgur.com/zFLKA.png) ![enter image description here](https://i.stack.imgur.com/DVY87.png) ![enter image description here](https://i.stack.imgur.com/sJKPM.png) ![enter image description here](https://i.stack.imgur.com/4Uvqv.png)
[Answer]
# Python, ~~571 567~~ 528 bytes
Similarly to Quincunx's solution, it starts by flood-filling each shape with an index from 1 to 4. It then determines the identity of the shapes by the radius of their bounding circle. A color-palette is constructed accordingly and the image is saved as an indexed-color image.
**EDIT:** Missed the fact the shapes are guaranteed not to touch the image border. Shorter it is, then!
```
from PIL.Image import*;from numpy import*
I=open(sys.argv[1]).convert("P")
D=list(I.getdata())
W=300;R=range(W*W);N=range(5)
O=[[0,i,array([0,0])]for i in N];n=0
for i in R:
 if D[i]>4:
    n+=1;S=[i]
    while S:
     j=S.pop()
     if D[j]>4:D[j]=n;O[n][0]+=1;O[n][2]+=j%W,j/W;S+=[j+1,j-1,j+W,j-W]
for o in O[1:]:o[2]/=o[0];o[0]=0
for i in R:
 if D[i]:o=O[D[i]];v=(i%W,i/W)-o[2];o[0]=max(o[0],dot(v,v))
O.sort()
C=[0]*5+[255]*3+[0,255,0,0]*2;P=C[:]
for i in N:j=3*O[i][1];P[j:j+3]=C[3*i:3*i+3]
I.putdata(D);I.putpalette(P);I.save("o.png")
```
Takes an input filename on the command line and writes the output to `o.png`.
[Answer]
## J - 246,224 185 bytes
```
load'viewmat'
(viewmat~0,(255*4 3$_2|.#:3720)/:/:@(<([:}.(>./%+/%#)@:(+/&:*:@(-"1)+/%#)@(4$.$.@:=)&>)<"0@~.@,))@(~.@,i.])@(>./**@{.)@((0,(,-)#:>:i.3)&|.)^:_@(*i.@:$)@(0<readimg_jqtide_)
```
This was a fun one!
I reused the connected-components part I used for the ["Am I in the biggest room" challenge](https://codegolf.stackexchange.com/questions/32015/are-you-in-the-biggest-room/40205#40205), and used the ratio between average and maximum distance of all points to each component's center. I settled for this, as it's both scale and rotation invariant, and apparently good enough to distinguish between the shapes as given. Ranking this value from low to high gives me the order circle, gear, square and triangle, used to permute the colormap.
Displays the result using the viewmap addon. No toolboxes used except for file reading and output.
Robustness does not seem to be a requirement, this takes off 18 bytes. 2 more unnecessary spaces, replaced `&.>` by `&>` in `ratio` and `&.:` by `&:` in dcent for another 2 bytes.
Huge gain in both shortness and performance of `comp` using shifting instead of `cut` (`;.`). This way the image is replicated and shifted in all 8 directions instead of scanning it with a 3x3 window.
The `id` function was ridiculously complex for what it needed to do. Now it assigns id's to pixels in objects by multiplying the image with an array of unique numbers, hence setting the BG to zero.
Code a bit more explained:
```
load'viewmat'                                 NB. display only
imnames =: < ;. _2 (0 : 0)
C6IKR.png
DLM3y.png
F1ZDM.png
Oa2O1.png
YZfc6.png
chJFi.png
)
images =: (0<readimg_jqtide_) each imnames    NB. read all images in boxed array
id =: *i.@:$                                  NB. NB. assign one number to each non-background (non-zero) pixel
comp =: (>./ * *@{.)@shift^:_@id              NB. 8 connected neighbor using shift
  shift =: (>,{,~<0 _1 1)&|.                  NB. generate the original, and 8 shifted versions (automatically padding and cropping).
result =: comp each images                    NB. Execute comp verb for each image
col =: (~.@, i. ])                            NB. Color: give each component and BG a separate color.
NB. BG in 0, 0 Get all max distance to center % mean distance to center ratios
ratio  =: (< ([:}.rat@:dcent@getInd &>)  <"0@~.@,)
  getInd =: 4 $. $.@:=                        NB. get indices for component y in array x
  dcent  =: +/&.:*:@(-"1) +/%#                NB. distence from center each point
  rat    =: >./ % +/%#                        NB. ratio from distances
cm=: (255*4 3$_2|.#:3720)                     NB. colormap (except black).
(viewmat~ 0,cm /: /:@ratio )@col each result  NB. for each image, show the result, permuting the colormap according to ratio's
NB. almostgolf this
P1 =: (>./**@{.)@((0,(,-)#:>:i.3)&|.)^:_@(*i.@:$)@(0<readimg_jqtide_) NB. reading till components
P2 =: (<([:}.(>./%+/%#)@:(+/&:*:@(-"1)+/%#)@(4$.$.@:=)&>)<"0@~.@,) NB. recognition: get fraction mean vs max distance to center per component, toss BG.     
P3 =: (viewmat~0,(255*4 3$_2|.#:3720)/:/:@P2)@(~.@,i.])@P1    NB. piece together : permute colormap, display components
NB. seriousgolf
load'viewmat'
f =:(viewmat~0,(255*4 3$_2|.#:3720)/:/:@(<([:}.(>./%+/%#)@:(+/&:*:@(-"1)+/%#)@(4$.$.@:=)&>)<"0@~.@,))@(~.@,i.])@((>./**@{.)@shift^:_)@(*i.@:$)@(0<readimg_jqtide_)
NB. example usage:
f&> imnames NB. do for all images
```
This one is a bit long to explain in detail, but will do if there is interest.
[Answer]
# Mathematica 225
---
**Update**:
The OP decided that this approach uses computer vision functions, so it's no longer in the running. I'll leave it posted however. Perhaps someone may find it of interest.
---
```
f@i_ := (m = MorphologicalComponents[ImageData@i];
Image@Partition[Flatten[(m)] /. 
   Append[ ReplacePart[SortBy[ComponentMeasurements[m, "Circularity"], Last], 
   {{1, 2} -> Yellow, {2, 2} -> Green, {3, 2} -> Red, {4, 2} -> Blue}], 0 -> Black], 
Dimensions[m][[2]]])
```
`ImageData` returns the image as a matrix of 0's and 1's.
`Flatten` converts that matrix to a list.
`Morphological Components` finds the 4 clusters of pixels and assigns a distinct integer,
1, 2, 3, 4 to each pixel according to the cluster. 0 is reserved for the (black) background.
`ComponentMeasurements` tests the circularity of the clusters.
From most to least circular will always be: the circle, the square, the triangle, and the gear.
`ReplacePart` replaces each component integer with the respective RGB color, using the circularity sort.
`Partition...Dimensions[m][[2]]` takes the list of pixel colors, and returns a matrix
the same dimensions as the input image.
`Image` converts the matrix of pixel colors to a colored image.
![inputs](https://i.stack.imgur.com/si2y1.png)
```
{f[img1],f[img2],f[img3],f[img4]}
```
![outputs](https://i.stack.imgur.com/lr0HQ.png)
[Answer]
# Mathematica, ~~354~~ ~~345~~ ~~314~~ ~~291~~ 288
Still golfing, could be shortened by a few more chars, but performance becomes unbearable. Uses the Variance to identify shapes:
```
f=(w=Position[z=ImageData@Import@#,1];r=Nearest;v@x_:=Variance@N[Norm[Mean@x-#]&/@x];Image[Plus@@(ReplacePart[0z/. 0->{0,0,0},#->r[{108,124,196,115}->List@@@{Blue,Red,Green,Yellow},v@#][[1]]]&/@Rest@NestList[(m=r[w=w~Complement~#];FixedPoint[Union@@(m[#,{8,2}]&/@#)&,{#&@@w}])&,{},4])])&
```
With spacing:
```
f = (w = Position[z = ImageData@Import@#, 1];
     r = Nearest; 
     v@x_ := Variance@N[Norm[Mean@x - #] & /@ x];
     Image[Plus @@ (ReplacePart[ 0 z /. 0 -> {0, 0, 0}, # -> r[{108, 124, 196, 115} -> 
                                              List @@@ {Blue, Red, Green, Yellow}, v@#][[1]]] & /@
     Rest@NestList[(m = r[w = w~ Complement~#];
                   FixedPoint[Union @@ (m[#, {8, 2}] & /@ #) &, {# & @@ w}]) &
                   , {}, 4])]) &
```
Testing:
```
s = {"https://i.stack.imgur.com/Oa2O1.png", "https://i.stack.imgur.com/C6IKR.png", 
     "https://i.stack.imgur.com/YZfc6.png", "https://i.stack.imgur.com/F1ZDM.png", 
     "https://i.stack.imgur.com/chJFi.png", "https://i.stack.imgur.com/DLM3y.png"};
Partition[f /@ s, 3] // Grid
```
![Mathematica graphics](https://i.stack.imgur.com/8kjeZ.png)
Here it is completely ungolfed. Will add explanations later:
```
findOneZone[{universe_List, lastZone_List}] :=
 Module[{newUniverse, proximityFindFunc, seedElement},
  newUniverse = Complement[universe, lastZone];
  proximityFindFunc = Nearest@newUniverse;
  seedElement = {First@newUniverse};
  {newUniverse, FixedPoint[Union @@ (proximityFindFunc[#, {8, 2}] & /@ #) &, seedElement]}]
colorAssign[zone_List] :=
 Module[{
   vlist = {108, 124, 196, 115},
   cols = List @@@ {Blue, Red, Green, Yellow},
   centerVariance},
  centerVariance[x_List] := Variance@N[Norm[Mean@x - #] & /@ x];
  First@Nearest[vlist -> cols, centerVariance@zone]]
colorRules[zones_List] := (# -> colorAssign[#] & /@ zones)
main[urlName_String] := 
 Module[{pixels, FgPixelPositions, rawZones, zones},
  pixels = ImageData@Import@urlName;
  FgPixelPositions = Position[pixels, 1];
  (*fill and separate the regions*)
  rawZones = NestList[findOneZone[#] &, {FgPixelPositions, {}}, 4];
  zones = Rest[rawZones][[All, 2]];
  (*Identify,colorize and render*)
  Image@ReplacePart[ConstantArray[{0, 0, 0}, Dimensions@pixels], 
    colorRules[zones]]]
s = {"https://i.stack.imgur.com/Oa2O1.png"};
main /@ s
```
[Answer]
# Python, ~~579~~ ~~577~~ ~~554~~ ~~514~~ ~~502~~ 501 bytes
For each shape, flood-fills it, then compute the distance between the centroid and the farthest point.
then the real surface of the shape is compared to the surface of a triangle, square, disc or wheel that would have the same size.
```
import math;from PIL.Image import*;A,R,_,I=abs,range(300),255,open(sys.argv[1]).convert('P');Q=I.load()
for j in R:
 for i in R:
  if Q[i,j]==_:
   X,Y,s,z,p=0,0,0,[],[(i,j)]
   while p:
    a,b=n=p.pop()
    if not(Q[n]!=_ or n in z):
     X+=a;Y+=b;z+=[n];p+=[(a,b-1),(a+1,b),(a,b+1),(a-1,b)];s+=1
   r=max([math.hypot(X/s-x,Y/s-y) for x,y in z]);C={1:A(s-(1.4*r)**2),2:A(s-r*r/3),3:A(s-math.pi*r*r),4:A(s-2.5*r*r)}
   for p in z:
    Q[p]=min(C,key=C.get)
I.putpalette([0,0,0,_]*3+[_,_,0])
I.show()
```
[Answer]
**C# 1086 bytes**
Yet another floodfill-solution, just for the record since there is no C# version here.
Like [Quincunx](https://codegolf.stackexchange.com/a/40806/22846) i wanted to prove myself that I can do it and there is not much differnce to his approach in Java.
* This solution does not use any recursion (but a stack) because i kept
running into StackOverflows.
* The detection of borderpixels is
simplified by looking at the 4 next pixels, if any of those is black,
the current is a border pixel.
It accepts every imageformat.
* Parameter 1 = InputPath
* Parameter 2 = OutputPath
It probably can be stripped down a few charachters by removing all the static stuff and creating an instance of Program.
Readable Version:
```
using System;
using System.Collections.Generic;
using System.Drawing;
using System.Drawing.Imaging;
using System.Linq;
class Program
{
    static Bitmap img;
    static int w, h;
    static ISet<Point> pointsDone = new HashSet<Point>();
    static void Main(string[] a)
    {
        img = new Bitmap(a[0]);
        w = img.Width;
        h = img.Height;
        Bitmap clone = new Bitmap(w,h, PixelFormat.Format32bppArgb);
        Graphics.FromImage(clone).DrawImage(img, 0, 0, w, h);
        img = clone;
        Color[] colors = new[] { Color.Blue, Color.Red, Color.Green, Color.Yellow };
        var shapes = new List<ISet<Tuple<bool, Point>>>();
        for(int x=0;x<w;x++)
            for (int y = 0; y < h; y++)
            {
                Point p = new Point(x, y);
                if (pointsDone.Add(p) && _isWhitePixel(p))
                    shapes.Add(_detectShape(p));
            }
        int index = 0;
        foreach (var shp in shapes.OrderBy(shp => shp.Count(item => item.Item1)))
        {
            foreach (var pixel in shp)
                img.SetPixel(pixel.Item2.X, pixel.Item2.Y, colors[index]);
            index++;
        }
        img.Save(a[1]);
    }
    private static ISet<Tuple<bool, Point>> _detectShape(Point p)
    {
        var todo = new Stack<Point>(new[] { p });
        var shape = new HashSet<Tuple<bool, Point>>();
        do
        {
            p = todo.Pop();
            var isBorderPixel = false;
            foreach (var n in new[] { new Point(p.X + 1, p.Y), new Point(p.X - 1, p.Y), new Point(p.X, p.Y + 1), new Point(p.X, p.Y - 1) })
                if (_isWhitePixel(n))
                {
                    if (pointsDone.Add(n))
                        todo.Push(n);
                }
                else isBorderPixel = true; // We know we are at the border of the shape
            shape.Add(Tuple.Create(isBorderPixel, p));
        } while (todo.Count > 0);
        return shape;
    }
    static bool _isWhitePixel(Point p)
    {
        return img.GetPixel(p.X, p.Y).ToArgb() == Color.White.ToArgb();
    }
}
```
Golfed:
```
using System;using System.Collections.Generic;using System.Drawing;using System.Drawing.Imaging;using System.Linq;class P{static Bitmap a;static int w,h;static ISet<Point> d=new HashSet<Point>();static void Main(string[] q){a=new Bitmap(q[0]);w=a.Width;h=a.Height;var c=new Bitmap(w,h,PixelFormat.Format32bppArgb);Graphics.FromImage(c).DrawImage(a,0,0,w,h);a=c;var e=new[]{Color.Blue,Color.Red,Color.Green,Color.Yellow};var f=new List<ISet<dynamic>>();for(int x=0;x<w;x++)for(int y=0;y<h;y++){Point p=new Point(x,y);if (d.Add(p)&&v(p))f.Add(u(p));}int i=0;foreach(var s in f.OrderBy(s=>s.Count(item=>item.b))){foreach(var x in s)a.SetPixel(x.p.X,x.p.Y,e[i]);i++;}a.Save(q[1]);}private static ISet<dynamic> u(Point p){var t=new Stack<Point>(new[]{p});var s=new HashSet<dynamic>();do{p=t.Pop();var b=false;foreach(var n in new[]{new Point(p.X+1,p.Y),new Point(p.X-1,p.Y),new Point(p.X,p.Y+1),new Point(p.X,p.Y-1)})if(v(n)){if (d.Add(n))t.Push(n);}else b=true;s.Add(new{b,p});}while (t.Count>0);return s;}static bool v(Point p){return a.GetPixel(p.X,p.Y).ToArgb()==Color.White.ToArgb();}}
```
]