Reproduce an image using lines Write a program that takes in a true-color RGB image I, the maximum number of lines to draw L, and the minimum m and maximum M length of each line. Output an image O that looks as much as possible like I and is drawn using L or fewer straight lines, all of which have Euclidean length between m and M.

Each line must be one solid color, have both endpoints in the bounds of O, and be drawn using Bresenham's line algorithm (which most graphics libraries will already do for you). Individual lines can only be 1-pixel thick.

All lines, even those of length 0, should take up at least one pixel. Lines may be drawn on top of one another.

Before drawing any lines you may initialize the background of O to any solid color (that may depend on I).

Details


O should have the same dimensions as I.
L will always be a nonnegative integer. It may be greater than the area of I.
m and M are nonnegative floating point numbers with M >= m. The distance between two pixels is the Euclidean distance between their centers. If this distance is less than m or greater than M, then a line between those pixels is not allowed.
The lines should not be antialiased.
Opacity and alpha should not be used.
Your program should not take more than an hour to run on a decent modern computer on images with less than a million pixels and L less than 10,000.


Test Images

You should of course show us your most accurate or interesting output images (which I expect will occur when L is between 5% and 25% of the number of pixels in I, and m and M are around one tenth of the diagonal size).

Here are some test images (click for originals). You may also post your own.

   


Simpler images:

 

This is a popularity contest. The highest voted submission wins.

Notes


It may be helpful to let L be derived from a percentage of total pixels in I as well as an absolute value. e.g. &gt;&gt;&gt; imageliner I=img.png L=50% m=10 M=20 would be the same thing as &gt;&gt;&gt; imageliner I=img.png L=32 m=10 M=20 if img.png were an 8 by 8 pixel image. Something similar could be done for m and M. This is not required.
Since lines cannot go out of bounds, the longest lines possible will be the diagonal length of I. Having M higher than this should not break anything though.
Naturally, if m is 0 and L is greater than or equal to the number of pixels in I, O could be identical to I by having length 0 "lines" at each pixel location. This behavior is not required.
Arguably, reproducing the shape of I is more important than reproducing the color. You may want to look into edge detection.
C++ - somewhat random lines and then some

First some random lines

The first step of the algorithm randomly generates lines, takes for the target image an average of the pixels along this, and then calculates if the summed up square of rgb space distances of all pixels would be lower if we would paint the new line (and only paint it, if it is). The new lines color for this is chosen as the channel wise average of the rgb values, with a -15/+15 random addition. 

Things that I noticed and influenced the implementation:


The initial color is the average of the complete image. This is to counter funny effects like when making it white, and the area is black, then already something as a bright green line is seen better, as it is nearer to black than the already white one.
Taking the pure average color for the line is not so good as it turns out to be unable to generate highlights by being overwritten by later lines. Doing a little random deviation helps a bit, but if you look at starry night, it fails if the local contrast is high at many places.


I was experimenting with some numbers, and chose L=0.3*pixel_count(I) and left m=10 and M=50. It will produce nice results starting at around 0.25 to 0.26 for the number of lines, but I chose 0.3 to have more room for accurate details.

For the full sized golden gate image, this resulted in 235929 lines to paint (for which it took a whopping 13 seconds here). Note that all images here are displayed in reduced size and you need to open them in a new tab/download them to view the full resolution.



Erase the unworthy

The next step is rather expensive (for the 235k lines it took about an hour, but that should be well within the "an hour for 10k lines on 1 megapixel" time requirement), but it is also a bit surprising. I go through all the previously painted lines, and remove those that do not make the image better. This leaves me in this run with only 97347 lines that produce the following image:



You probably need to download and compare them in an appropriate image viewer to spot most differences.

and start over again

Now I have a lot of lines that I can paint again to have a total of 235929 again. Not much to say, so here is the image:



short analysis

The whole procedure seems to work like a blurring filter that is sensitive to local contrast and object sizes. But it is also interesting to see where the lines are painted, so the program records these too (For each line, the pixel color will be made one step whiter, at the end the contrast is maximized). Here are the corresponding ones to the three colored above.





animations

And since we all love animations, here are some animated gifs of the whole process for the smaller golden gate image. Note that there is significant dithering due to gif format (and since creators of true color animation file formats and browser manufacturers are in war over their egos, there is no standard format for true color animations, otherwise I could have added an .mng or similar).




Some more

As requested, here are some results of the other images (again you might need to open them in a new tab to not have them downscaled)






Future thoughts

Playing around with the code can give some intresting variations.


Chose the color of the lines by random instead of being based on the average. You might need more than two cycles.
The code in the pastebin also contains some idea of a genetic algorithm, but the image is probably already so good that it would take too many generations, and this code is also too slow to fit into the "one hour" rule.
Do another round of erase/repaint, or even two...
Change the limit of where lines can be erased (e.g. "must make the image at lest N better")


The code

These are just the two main useful functions, the whole code doesn't fit in here and can be found at http://ideone.com/Z2P6Ls

The bmp classes raw and raw_line function do access pixels and lines respectively in an object that can be written to bmp format (It was just some hack lying around and I thought that makes this somewhat independent from any library).

The input file format is PPM

std::pair&lt;bmp,std::vector&lt;line&gt;&gt;  paint_useful( const bmp&amp; orig, bmp&amp; clone, std::vector&lt;line&gt;&amp; retlines, bmp&amp; layer, const std::string&amp; outprefix, size_t x, size_t y )
{
        const size_t pixels = (x*y);
        const size_t lines = 0.3*pixels;
//      const size_t lines = 10000;

//      const size_t start_accurate_color = lines/4;

        std::random_device rnd;

        std::uniform_int_distribution&lt;size_t&gt; distx(0,x-1);
        std::uniform_int_distribution&lt;size_t&gt; disty(0,y-1);
        std::uniform_int_distribution&lt;size_t&gt; col(-15,15);
        std::uniform_int_distribution&lt;size_t&gt; acol(0,255);

        const ssize_t m = 1*1;
        const ssize_t M = 50*50;

        retlines.reserve( lines );

        for (size_t i = retlines.size(); i &lt; lines; ++i)
        {
                size_t x0;
                size_t x1;

                size_t y0;
                size_t y1;

                size_t dist = 0;
                do
                {
                        x0 = distx(rnd);
                        x1 = distx(rnd);

                        y0 = disty(rnd);
                        y1 = disty(rnd);

                        dist = distance(x0,x1,y0,y1);
                }
                while( dist &gt; M || dist &lt; m );

                std::vector&lt;std::pair&lt;int32_t,int32_t&gt;&gt; points = clone.raw_line_pixels(x0,y0,x1,y1);

                ssize_t r = 0;
                ssize_t g = 0;
                ssize_t b = 0;

                for (size_t i = 0; i &lt; points.size(); ++i)
                {
                        r += orig.raw(points[i].first,points[i].second).r;
                        g += orig.raw(points[i].first,points[i].second).g;
                        b += orig.raw(points[i].first,points[i].second).b;
                }

                r += col(rnd);
                g += col(rnd);
                b += col(rnd);

                r /= points.size();
                g /= points.size();
                b /= points.size();

                r %= 255;
                g %= 255;
                b %= 255;

                r = std::max(ssize_t(0),r);
                g = std::max(ssize_t(0),g);
                b = std::max(ssize_t(0),b);

//              r = acol(rnd);
//              g = acol(rnd);
//              b = acol(rnd);

//              if( i &gt; start_accurate_color )
                {
                        ssize_t dp = 0; // accumulated distance of new color to original
                        ssize_t dn = 0; // accumulated distance of current reproduced to original
                        for (size_t i = 0; i &lt; points.size(); ++i)
                        {
                                dp += rgb_distance(
                                                                                orig.raw(points[i].first,points[i].second).r,r,
                                                                                orig.raw(points[i].first,points[i].second).g,g,
                                                                                orig.raw(points[i].first,points[i].second).b,b
                                                                        );

                                dn += rgb_distance(
                                                                                clone.raw(points[i].first,points[i].second).r,orig.raw(points[i].first,points[i].second).r,
                                                                                clone.raw(points[i].first,points[i].second).g,orig.raw(points[i].first,points[i].second).g,
                                                                                clone.raw(points[i].first,points[i].second).b,orig.raw(points[i].first,points[i].second).b
                                                                        );

                        }

                        if( dp &gt; dn ) // the distance to original is bigger, use the new one
                        {
                                --i;
                                continue;
                        }
                        // also abandon if already too bad
//                      if( dp &gt; 100000 )
//                      {
//                              --i;
//                              continue;
//                      }
                }

                layer.raw_line_add(x0,y0,x1,y1,{1u,1u,1u});
                clone.raw_line(x0,y0,x1,y1,{(uint32_t)r,(uint32_t)g,(uint32_t)b});
                retlines.push_back({ (int)x0,(int)y0,(int)x1,(int)y1,(int)r,(int)g,(int)b});

                static time_t last = 0;
                time_t now = time(0);
                if( i % (lines/100) == 0 )
                {
                        std::ostringstream fn;
                        fn &lt;&lt; outprefix + "perc_" &lt;&lt; std::setw(3) &lt;&lt; std::setfill('0') &lt;&lt; (i/(lines/100)) &lt;&lt; ".bmp"; 
                        clone.write(fn.str());
                        bmp lc(layer);
                        lc.max_contrast_all();
                        lc.write(outprefix + "layer_" + fn.str());
                }

                if( (now-last) &gt; 10 )
                {
                        last = now;
                        static int st = 0;
                        std::ostringstream fn;
                        fn &lt;&lt; outprefix + "inter_" &lt;&lt; std::setw(8) &lt;&lt; std::setfill('0') &lt;&lt; i &lt;&lt; ".bmp";
                        clone.write(fn.str());

                        ++st;
                }
        }
        clone.write(outprefix + "clone.bmp");
        return { clone, retlines };
}


void erase_bad( std::vector&lt;line&gt;&amp; lines, const bmp&amp; orig )
{
        ssize_t current_score = evaluate(lines,orig);

        std::vector&lt;line&gt; newlines(lines);

        uint32_t deactivated = 0;
        std::cout &lt;&lt; "current_score = " &lt;&lt; current_score &lt;&lt; "\n";
        for (size_t i = 0; i &lt; newlines.size(); ++i)
        {
                newlines[i].active = false;
                ssize_t score = evaluate(newlines,orig);
                if( score &gt; current_score )
                {
                        newlines[i].active = true;
                }
                else
                {
                        current_score = score;
                        ++deactivated;
                }
                if( i % 1000 == 0 )
                {
                        std::ostringstream fn;
                        fn &lt;&lt; "erase_" &lt;&lt; std::setw(6) &lt;&lt; std::setfill('0') &lt;&lt; i &lt;&lt; ".bmp";
                        bmp tmp(orig);
                        paint(newlines,tmp);
                        tmp.write(fn.str());
                        paint_layers(newlines,tmp);
                        tmp.max_contrast_all();
                        tmp.write("layers_" + fn.str());
                        std::cout &lt;&lt; "\r i = " &lt;&lt; i &lt;&lt; std::flush;
                }
        }
        std::cout &lt;&lt; "\n";
        std::cout &lt;&lt; "current_score = " &lt;&lt; current_score &lt;&lt; "\n";
        std::cout &lt;&lt; "deactivated = " &lt;&lt; deactivated &lt;&lt; "\n";


        bmp tmp(orig);

        paint(newlines,tmp);
        tmp.write("newlines.bmp");
        lines.clear();
        for (size_t i = 0; i &lt; newlines.size(); ++i)
        {
                if( newlines[i].is_active() )
                {
                        lines.push_back(newlines[i]);
                }
        }
}