Message Boards Message Boards

Something About The Automated Stippling Drawing

1. Stippling drawing -- an introduction

results examples

Stippling is a kind of drawing style using only points to mimic lines, edges and grayscales. The entire drawing consists only of black points on a white background. The density of the points gives the impression of grayscale shading.

Back in 1510, stippling was first invented as an engraving technique, then became popular in many fields because of its requirement of just one color of ink. The style influence widely beyond engraving. We can see many technical illustrations done by it, especially in math and mechanics books. Even today many people are interested in it, and archaeologists, geologists and biologists are still using it in field works.

On the art side, stippling shows interesting relation with the divisionism and pointillism, which both belong to the neo-impressionism.

In the domain of so-called non-photorealistic computer graphics, stippling is one of the basic artistic painting style, and can be used as the foundation of many other styles such as mosaic, stained glass, hatching, etc. (Non-Photorealistic Computer Graphics, Thomas Strothotte, Stefan Schlechtweg; The Algorithms and Principles of Non-photorealistic Graphics, Weidong Geng)

2. The essential properties of a "good" stippling

2.1 Well-spaced points

When using points to approximate a grayscale, a random set of points with uniform distributuion is usually not good enough. To illustrate, both images below have 63024 points, but the one on the right is a good stippling (or usually called well-spaced). The one on the left has too many unwanted small artifacts -- clumps and voids do not exist in the original image. (And please be patient, we'll talk about how to generate the "good" one in the next section:)

comparison between random and stippling - non-uniform case

It can be seen that there are too many unwanted small artifacts in the former one -- clumps and voids that do not exist in the original image.

The phenomenon can be seen even more clearly in the uniform grayscale case:

comparison between random and stippling - uniform case

I hope it's plain to see, that in the stippling graphics, the distance between any point and the one nearest to it is nearly a constant. So it looks like we found the essential property of a "good" stippling.

But wait... What about the Lena example? Clearly it's "good" ( if you allow me to say so) but the distances between closest points are various.

It turns out we just need to take one more step to notice that the constant distance implies equilateral triangles:

comparison between random and stippling - uniform DelaunayMesh

So the so-called good or well-spaced property can be described more precisely as:

The cells of the Delaunay mesh induced by the points are as equilateral triangles as possible.

The statement can be confirmed by the statistics on the set of interior angles of all the cells

interior angles statistics

2.2 Another (equivalent) definition

There is still another view point to judge if the positions of the points are good or not.

Consider the same blur taken on all the three Lena images: the original one, the random points one and the well-spaced one:

comparison among blurred images

The so-called projection algorithm looks for a point distribution, when convolved with a certain kernel h, which is as the same as the one obtained by convolving the original image with the same h. We are not going to talk about the details here. Readers who are interested in it can refer to the mentioned paper.

3. Stippling of a uniform area -- centroidal Voronoi diagram approach

Now let's do the real fun :)

For an easy start we consider the stippling technique for a uniform area (i.e. an area with constant grayscale). There is a de facto standard method for it called centroidal Voronoi diagram (CVD) which is usually generated by the Lloyd's algorithm.

Basically, the algorithm acts as the following relaxation steps:

  1. Generate $n$ random points inside the interested region

  2. Generate the Voronoi diagram of the $n$ points

  3. Find the centroid (i.e. center of mass) of each Voronoi cell

  4. Use the $n$ centroids as the resulting points

  5. If satisfied, algorithm stop, else goto step 1

Here the key techniques are the Voronoi diagram generation and the centroid finding. The former one is a perfect match for the bounded version of the built-in function VoronoiMesh. The latter one, as the Voronoi cells for a closed region are always closed convex polygons, has a simple formula.

Suppose the cell is defined by $n$ vertices ordered (counter-)clockwise $\{\boldsymbol{P}_1=(x_1,y_1),\boldsymbol{P}_2=(x_2,y_2),\dots,\boldsymbol{P}_n=(x_n,y_n),\boldsymbol{P}_{n+1}=(x_1,y_1)\}$, then its centroid $\boldsymbol{C}$ can be determined as follwing Eq. 1:

$$\left\{\begin{eqnarray} \boldsymbol{C}&=&\frac{1}{6A}\sum_{k=1}^n (\boldsymbol{P}_i+\boldsymbol{P}_{i+1})\,\left|\begin{matrix}x_i&y_i\\x_{i+1}&y_{i+1}\end{matrix}\right|\\ A&=&\frac{1}{2}\sum_{k=1}^n \left|\begin{matrix}x_i&y_i\\x_{i+1}&y_{i+1}\end{matrix}\right| \end{eqnarray}\right.$$

As a test we generate 2500 uniformly distributed random points in a square region $[-1,1]\times[-1,1]$:

initPts = RandomPoint[Rectangle[{-1, -1}, {1, 1}], 2500];

whose Voronoi diagram is:

VoronoiMesh[initPts, {{-1, 1}, {-1, 1}}] //
   Graphics[{
               GrayLevel[.6],
               Line[Identity @@@ MeshPrimitives[#, 1]],
               Red, AbsolutePointSize[2],
               Point[initPts]
            }] &

Voronoi of random points

The Lloyd's algorithm can be expressed as the following findCentroid function:

Clear[findCentroid]
findCentroid[p : Polygon[{{__Real} ..}]] :=
        Module[{pts = p[[1]], pairs, dets, area},
               pairs = Partition[pts, 2, 1, {1, 1}];
               dets = Det /@ pairs;
               area = 1/2 Plus @@ dets;
               dets.(Plus @@@ pairs)/(6 area)
              ]

The cell Polygon-s of the Voronoi mesh can be extracted with MeshPrimitives[...,2], which then should be piped to the findCentroid to complete one iteration:

Module[{vm},
        Function[pts,
                 vm = VoronoiMesh[pts, {{-1, 1}, {-1, 1}}];
                 findCentroid /@ MeshPrimitives[vm, 2]
                ]@initPts
      ];

Now we are ready to animate the first 50 iteration results to give a intuitive feeling about the CVD:

intermRes =
            Module[{vm},
                NestList[
                    Function[pts,
                        vm = VoronoiMesh[pts, {{-1, 1}, {-1, 1}}];
                        findCentroid /@ MeshPrimitives[vm, 2]
                        ],
                    initPts,
                    500
                    ]
                ]; // AbsoluteTiming
refinedPts = intermRes[[-1]];
allframes =
        Function[pts,
                VoronoiMesh[pts, {{-1, 1}, {-1, 1}}] //
                    Graphics[{
                                GrayLevel[.6],
                                Line[Identity @@@ MeshPrimitives[#, 1]],
                                Red, AbsolutePointSize[2],
                                Point[pts]
                                }, ImageSize -> 600] &
                ] /@ intermRes[[;; 50]];
ListAnimate[allframes, 24, DisplayAllSteps -> True, AnimationRepetitions -> 1]

{197.984, Null}

uniform CVD

There are various ways to show the difference between the points distributions before and after the process.

For example we can use NearestNeighborGraph to illustrate the connectivities of them, which will highlight the unwanted voids in the former case:

MapThread[
        Function[{pts, label},
            Graphics[
                    Line[EdgeList[NearestNeighborGraph[pts, 3]] /. 
                            UndirectedEdge -> List], ImageSize -> 400] //
                {Style[label, 20, FontFamily -> "Constantia"], #} &
            ],
        {{initPts, refinedPts}, {"initial points", "refined points"}}
        ] // Grid[# // Transpose, Alignment -> Left, Frame -> All] &

Comparison of connectivity between random and stippling

Or as having been shown in the previous section, we can compare the statistics of the interior angles:

Plot[
    Evaluate[MapThread[
            Legended[PDF[KernelMixtureDistribution[#1], ?], 
                    Style[#2, 15]] &,
            {
                (
                            # // DelaunayMesh // MeshPrimitives[#, 2] & //
                                    Function[pts, 
                                                VectorAngle[#2 - #1, #3 - #1]/Degree & @@@ 
                                                    Partition[pts, 3, 1, {1, 1}]] @@@ # & // Flatten
                            ) & /@ {initPts, refinedPts},
                {"initial points", "refined points"}
                }]],
    {?, 0, 180},
    PlotRange -> All, GridLines -> {{60}, None}, Frame -> True, 
    FrameLabel -> (Style[#, Bold, 15] & /@ {"angle in °", "PDF"})
    ]

interior angles statistics

To give another intuitive impression on the "well-spaced" property from a different point of view, we compare the discrete Fourier transformation of the initial points with the one of the refined points:

MapThread[
        Function[{pts, label},
            Graphics[{AbsolutePointSize[0], Point@pts}] //
                                                            Rasterize[#, ImageSize -> {701, 701}] & // Binarize // 
                                                    ColorNegate //

                                                Function[img, 
                                                    ImageMultiply[img, 
                                                        DiskMatrix[.4 ImageDimensions[img], 
                                                                ImageDimensions[img]] // Image]] //
                                            ImageData // # - Mean[Flatten@#] & //
                                    Fourier // Abs // Rescale // 
                        RotateRight[#, Floor[Dimensions[#]/2]] & //
                    Image[1.2 - (1 - #)^5, "Real", ImageSize -> 400] & //
                {Style[label, 20, FontFamily -> "Constantia"], #} &
            ],
        {{initPts, refinedPts}, {"initial points", "refined points"}}
        ] // Grid[#//Transpose, Alignment -> Left, Frame -> All] &

comparison of FFT between random and stippling

It can be clearly seen that in the refined points' case, we effectively achieved both isotropism and low-stop filter, which is another facet of view to understand the term well-spaced. Note that the property is also related to the concept of "colors" of noise.

4. Stippling of a non-uniform area

So far we only talked about the method for stippling a uniform area. But how about a non-uniform area, like the Lena example?

The Lloyd's algorithm described above can not be directly adapted here without modification, as it always smoothes out the distribution, results a uniform CVD. So how can we generalize it to non-uniform case?

4.1 Approach using the conformal mapping

Recall that we were looking for a Delaunay mesh with cells that resembles equilateral triangles as much as possible, one thing that immediately came to our mind is the conformal map which transforms between two spaces while preserving angles locally. Thus a good stippling on a uniform area is guaranteed to be mapped to a good one on a specified non-uniform area.

A simple example of conformal map is a complex holomorphic function $f$, say $f(z)=(z-z^2/2)\,\exp\left[-(z-1)^2/5\right]$:

cmfunc = Compile[{{z, _Complex}},
                 Module[{w = 1. I},
                        w = (z - 1/2 z^2) Exp[-((z - 1)^2/5)];
                        {Re@w, Im@w}
                       ],
                 RuntimeOptions -> "Speed",
                 RuntimeAttributes -> {Listable},
                 Parallelization -> True
                ];

which transforms the stippling points refinedPts in the uniform square $[-1,1]\times[-1,1]$ to a new points distribution transPts:

transPts = refinedPts.{1, I} // cmfunc;
transPts // Graphics[{AbsolutePointSize[2], Point@#}] &

conformal transformed points

The result looks very nice in the sense of "well-spaced", which can also be confirmed with the Delaunay mesh, Voronoi mesh and its connectivity:

transPts // DelaunayMesh

DelaunayMesh of the transPts

Module[{vm = VoronoiMesh[refinedPts, {{-1, 1}, {-1, 1}}]},
    Graphics[GraphicsComplex[
            (Identity @@@ MeshPrimitives[vm, 0]).{1, I} // cmfunc,
            vm // MeshCells[#, 1] & // Line[Identity @@@ #] &
            ]]
    ]

VoronoiMesh of the transPts

transPts // NearestNeighborGraph[#, 3] & // 
                EdgeList // # /. UndirectedEdge -> List & // Line // Graphics

connectivity of the transPts

So far so good. However, this theoretically elegant approach is not easily generalizable. To find the right $f$ for an arbitrary target distribution will need sophisticated mathematical skill. Again we are not going to talk about the details here. For readers who are interested in it, there is a whole dedicated research field called computing conformal mapping.

4.2 Approach using the weighted Voronoi diagram

Despite the elegance of the conformal mapping, the popular method for stippling non-uniform area comes from a modification of the CVD, which is called the weighted centroidal Voronoi diagram (A. Secord. Weighted Voronoi stippling. In Proceedings of the second international symposium on Non-photorealistic animation and rendering, pages 37-43. ACM Press, 2002.) (A. Secord. Random Marks on Paper, Non-Photorealistic Rendering with Small Primitives, Thesis for the Degree of Master of Science, 2002).

The algorithm is similar to the Lloyd's algorithm, only in step 3, when looking for the centroid of the Voronoi cell, a variable areal density $\rho(\boldsymbol{P})$ is considered, which is proportional to the grayscale at the location $\boldsymbol{P}$. Thus instead of using Eq. 1, we shall calculate the centroid according to the following definition (Eq. 2):

$$\boldsymbol{C}_{\text{Cell}(\boldsymbol{P})} = \frac{\int_{\boldsymbol{x}\in\text{Cell}(\boldsymbol{P})}\boldsymbol{x}\rho(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}}{\int_{\boldsymbol{x}\in\text{Cell}(\boldsymbol{P})}\boldsymbol{x}\,\mathrm{d}\boldsymbol{x}}$$

Clearly the integrations are much more time-consuming than Eq. 1. In his paper and master thesis Secord presented an efficient way to compute them, which involves a precomputation of certain integrations. However, I noticed (without theoretical proof) that the weighted CVD can be sufficiently approximated in a much cheaper way if we accept a compromise not to stick to the exact formula (i.e. Eq. 2) of centroid but only to emphasis the core idea of choosing $\boldsymbol{C}$ closer to points with larger weights.

The new idea is simple. For a cell of $n$ vertices $\{\boldsymbol{P}_1,\boldsymbol{P}_2,\dots,\boldsymbol{P}_n\}$, the algorithm acts as follows:

approx algorithm sketch

  1. Compute the geometric centroid $\boldsymbol{C}$

  2. Compute the weights of the vertices as $\{w_1,\dots,w_n\}$

  3. Compute normalized weights as $W_k = \frac{w_k}{\max(w_1,\dots,w_n)}$

  4. For every vertex $\boldsymbol{P}_k$, move it along $\overrightarrow{\boldsymbol{C}\boldsymbol{P}_k}$ with factor of $W_k$ to new position $\boldsymbol{P}_k'$ (So vertex with largest weight does not move, vertex with smallest weight moves most)

  5. Compute the geometric centroid $\boldsymbol{C}'$ of the new cell defined by $\{\boldsymbol{P}_1',\dots,\boldsymbol{P}_n'\}$ as the final result

Note that the convergency of our simplified algorithm is not obvious, so it might be wise to do an "early stop" during iteration.

Written in the Wolfram Language, it's this findCentroidW function:

Clear[findCentroidW]
findCentroidW[p : Polygon[{{__Real} ..}], densityFunc_] :=
    Module[{cent, pts = p[[1]], wlst},
           wlst = #/Max[#] &[densityFunc @@@ pts];
           cent = findCentroid[p];
           cent + findCentroid@Polygon@MapThread[#2 (#1 - cent) &, {pts, wlst}]
          ]

Numerical experiments show that this approximation gives fairly good results. In the next section we'll demonstrate its application on the Lena image.

4.3 Numerical experiment on Lena

In this section we test our non-uniform stippling method on the famous Lena's photo.

First let's import the original image. We'll keep both the grayscale and color versions for later use in the artistic rendering section:

imgOrigColored = ExampleData[{"TestImage", "Lena"}];
imgOrig = imgOrigColored // ColorConvert[#, "Grayscale"] & // ImageAdjust;

For Interpolation's convenience, we'll use ColorNegate. For a better visual feeling, we enhance the edges with large gradient:

img = Function[img,
                 ImageAdd[ImageMultiply[img, # // ColorNegate], #] &[
                             GradientFilter[img, 1] // ImageAdjust[#, {2, 2, 1}] &
                            ]
              ][imgOrig // ColorNegate]

Lena for density image

The image coordinate is rescaled to the rectangle region $[-1,1]\times[-1,1]$ for personal convenience and Interpolation-ed to get a smooth description of the grayscale field:

densityFunc = 
        img // Reverse[ImageData[#, "Real"]] & // Transpose // Function[array,
                        Module[{dim = Dimensions@array}, 
                            MapIndexed[{(1 + dim - 2 #2)/(1 - dim), #1} &, array, {-1}]]
                        ] // Flatten[#, 1] & // Interpolation;

To have a good initial points distribution, we would like to sample points so the local density is proportional to the local grayscale (though we don't need to have this precisely, as the weighted Voronoi process will anyway smooth the distribution.) So taking advantage of ContourPlot, we generate a few regions according to the level set of the densityFunc:

levelRegions =
            Module[{ctplot, pts, polys},
                ctplot = 
                    ContourPlot[densityFunc[x, y], {x, -1, 1}, {y, -1, 1}, 
                        Contours -> 10];
                {pts, polys} = 
                    Cases[ctplot, GraphicsComplex[pts_, e_, ___] :> {pts,
                                    Cases[e, GraphicsGroup[r : {__Polygon}] :> r, ?]
                                    }, ?][[1]];
                BoundaryDiscretizeGraphics[GraphicsComplex[pts, #]] & /@ polys
                ]; // AbsoluteTiming

{111.862, Null}

level sets

For each region in levelRegions, we sample points inside it on regular grid, with area of the grid cell inversely proportional to the level counter of the region. Notice the regular grid can be a steady state of our algorithm, to ensure a isotropic result, initial randomness is needed. For that purpose we add dithering effect on the points, with its strength specified by a parameter $0\leq\kappa\leq 1$ where $0$ gives no dithering while $1$ gives a globally random distribution:

Clear[ditherFunc]
ditherFunc[pts_, width_, ?_: .5] := 
    pts + ? width/2 RandomReal[{-1, 1}, {Length@pts, 2}]

levelPts = 
            Module[{baseGridWidth = 1/50., width, num, min = -1, 
                    max = 1, ? = 1},
                MapIndexed[
                    Function[{region, idx},
                        width = baseGridWidth/Sqrt[idx[[1]]];
                        num = Ceiling[(max - min)/width];
                        Tuples[Range@num, 2] // 
                                    Rescale[#, {1, num}, {min + width/2, max - width/2}] & // N //

                            ditherFunc[Select[#, RegionMember[region]], width, ?] &
                        ],
                    levelRegions,
                    {1}]
                ]; // AbsoluteTiming
initPts = levelPts // Join @@ # &;

{15.3046, Null}

The total number of points we sampled is huge:

initPts // Length

63024

But their quality is poor:

Graphics[{AbsolutePointSize[0], Point@initPts}]

Lena of random sample

Now we perform the weighted Voronoi relaxation process, which is similar to the one in the uniform case, though the computation is a bit slow due to the amount of points:

refinedPts = Module[{vm},
                Nest[
                    Function[pts,
                        vm = VoronoiMesh[pts, {{-1, 1}, {-1, 1}}];
                        findCentroidW[#, densityFunc] & /@ MeshPrimitives[vm, 2]
                        ],
                    initPts,
                    30
                    ]
                ]; // AbsoluteTiming

{1106.06, Null}

In spite of only 30 iterations, the result is in my opinion fairly good (Note that some visual artifacts in the following image is due to the rasterizing process during display, try right-click and "open image in new tab"):

Graphics[{AbsolutePointSize[0], Point@refinedPts}]

Lena of stippling

The pattern of the connectivity rather interestingly forms some kind of self-similar multi-resolution tiles:

refinedPts // NearestNeighborGraph[#, 2] & // 
                EdgeList // # /. UndirectedEdge -> List & // Line // Graphics

Lena of connectivity

Statistics on the interior angles of the corresponding DelaunayMesh indicates we have indeed achieved the well-spaced distribution:

Plot[
    Evaluate[MapThread[
            Legended[PDF[KernelMixtureDistribution[#1], ?], 
                    Style[#2, 15]] &,
            {
                (
                            # // DelaunayMesh // MeshPrimitives[#, 2] & // 
                                    Function[pts, 
                                                VectorAngle[#2 - #1, #3 - #1]/Degree & @@@ 
                                                    Partition[pts, 3, 1, {1, 1}]] @@@ # & // Flatten
                            ) & /@ {initPts, refinedPts},
                {"initial points", "refined points"}
                }]],
    {?, 0, 180},
    PlotRange -> All, GridLines -> {{60}, None}, Frame -> True, 
    FrameLabel -> (Style[#, Bold, 15] & /@ {"angle in °", "PDF"})
    ]

interior angles statistics of Lena

5. Artistic styles based on stippling

5.1 Hatching

As we mentioned in the introduction, some artistic rendering effects can be simulated based on the stippling result. One of them is hatching, for our example which will act like a pencil sketch. The lengths and orientations of the strokes are controled by densityFunc i.e. local grayscale, while the positions of the them are controled by the stippling points. Note there are artifacts in the dark regions due to our naïve approach:

Module[{strokeLength, strokeOrient, dens, stroke},
    Function[{x, y},
                    dens = densityFunc[x, y];
                    strokeLength = .03 dens^.3 (1 + .1 RandomReal[{-1, 1}]);
                    strokeOrient = ?/
                            3 + ?/6 (1 - (2 dens - 1)^2)^20 RandomReal[{-1, 1}];
                    stroke = strokeLength Through[{Cos, Sin}[strokeOrient]];
                    {{x, y} - stroke, {x, y} + stroke}
                    ] @@@ refinedPts // Line // Graphics[{AbsoluteThickness[0], #}] &
    ]

Lena sketch

5.2 Pointillism

Another styling is the simulation of the pointillism.

For this simulation, we need to take care of the colors:

imgColorChannels = imgOrigColored // ColorSeparate[#, "RGB"] &;
densityFuncColorChannels = 
        Function[img, 
                img // Reverse[ImageData[#, "Real"]] & // Transpose // 
                            Function[array,
                                Module[{dim = Dimensions@array}, 
                                    MapIndexed[{(1 + dim - 2 #2)/(1 - dim), #1} &, array, {-1}]]
                                ] // Flatten[#, 1] & // Interpolation] /@ imgColorChannels;

We add a random deviation of the color on each points:

randomColoredImg =
    {
                                    RGBColor @@ (

                                            Through[
                                                    densityFuncColorChannels[##]] (1 + .3 RandomReal[{-1, 1},
                                                                    3])
                                            ),
                                    Point[{##}]
                                    } & @@@ refinedPts //
                        Graphics[{AbsolutePointSize[0], #} // Flatten, 
                                ImageSize -> 1000] & //
                    Rasterize[#, ImageSize -> 1000] & //
                ImageConvolve[#, GaussianMatrix[{10, 1.3}]] & // ImageAdjust // 
        ImageCrop

Lena color stippling

We finalize the process with the morphological Opening operation and match the histogram with the original color image using HistogramTransform:

Opening[randomColoredImg, DiskMatrix[6]] // HistogramTransform[#, imgOrigColored] &

Lena pointillism

5.3 Pointillism -- another attempt

The rendering in the last section uses strokes with identical size, but in practice it can often vary with some local properties of the targets being drawn. So let's try to include that characteristic here.

For the special property to be reflected, we choose the ImageSaliencyFilter, but it is of course totally fine to choose other ones:

ImageAdjust[
                ImageSaliencyFilter[imgOrigColored, Method -> #]] & /@ {"Itti", 
            "IttiColor", "IttiIntensity", "IttiOrientation"};
coarseMask = ImageAdd @@ (Image[#, "Byte"] & /@ %) // Blur[#, 20] &

coarseMask

Like the densityFunc, we interpolation the coarseMask and generate some level regions:

coarseFunc = 
        coarseMask // Reverse[ImageData[#, "Real"]] & // Transpose // 
                    Function[array,
                        Module[{dim = Dimensions@array}, 
                            MapIndexed[{(1 + dim - 2 #2)/(1 - dim), #1} &, array, {-1}]]
                        ] // Flatten[#, 1] & // Interpolation;
coarseLevelRegions =
            Module[{ctplot, pts, polys},
                ctplot = 
                    ContourPlot[coarseFunc[x, y], {x, -1, 1}, {y, -1, 1}, 
                        Contours -> 3];
                {pts, polys} = 
                    Cases[ctplot, GraphicsComplex[pts_, e_, ___] :> {pts,
                                    Cases[e, GraphicsGroup[r : {__Polygon}] :> r, ?]
                                    }, ?][[1]];
                BoundaryDiscretizeGraphics[GraphicsComplex[pts, #]] & /@ polys
                ]; // AbsoluteTiming

{11.1717, Null}

So now the stippling points can be grouped according to the coarseLevelRegions:

coarseLevelPts = Select[refinedPts, RegionMember[#]] & /@ coarseLevelRegions;

The stroke template is composed with squares centered at the stippling points, with their color determined and randomized in a way similar to that in the last section:

strokeFunc = Function[{pt, size},
            GeometricTransformation[
                    Rectangle[{-1, -1}, {1, 1}],
                    RotationMatrix[
                            RandomReal[?/
                                2]].ScalingMatrix[(1 + .2 RandomReal[{-1, 1}]) size {1, 1}]
                    ] // GeometricTransformation[#, TranslationTransform[pt]] &
            ];

Now we are ready to paint the picture layer by layer. NOTE THAT the following code can be very time-consuming.

coarseLayers =
    Function[{pts, size},
                    {
                                        FaceForm[RGBColor @@ (
                                                    Append[
                                                        (1 + .5 RandomReal[{-1, 1}, 3]) Through[
                                                                densityFuncColorChannels @@ #],
                                                        RandomReal[{.2, 1}]
                                                        ]
                                                    )],
                                        strokeFunc[#, pts]} & /@ pts //
                            Graphics[Flatten@{EdgeForm[], #}, 
                                    PlotRange -> 1.5 {{-1, 1}, {-1, 1}}] & //
                        Rasterize[#, ImageSize -> 1000, Background -> None] &
                    ] @@ # & /@ ({coarseLevelPts, {.2, .1, .05, .02}}//Transpose)

paint layers

Composing the layers, we get our final "neo-impressionism" result:

composedImg = FoldList[ImageCompose, coarseLayers[[1]], coarseLayers[[2 ;;]]]

paint stages

composedImg[[-1]] // ImageCrop

Lena pointillism 2

Note the local vivid color blocks especially those at the feather and hat regions, and how the rendering still preserved an accurate scene at the macro-scale. :)

6. Epilogue

In this small writings, we discussed a morden automated method for simulating an ancient art skill called stippling drawing. We also demonstrated how its power can be extended to more advanced artistic renderings. But readers should be reminded that there are other methods for generating the well-spaced point distribution, and there are still many fields related to non-photorealistic computer graphics which are not necessarily related to it. Finally, although not mentioned in this article, the applications of the so-called well-spaced point distribution goes way beyond the examples shown here. The idea can be generalized to high dimensions, and its nice properties makes it an import tool in many scientific feilds, like re-meshing in the finite element method, sparse representation and compressed sensing.

POSTED BY: Silvia Hao
9 Replies
Posted 3 years ago

Excellent post! (And I would imagine that what you've shown would have implications for statistical sampling: completely random vs something a bit more structured in estimating quantities in landscapes from GIS images.)

I do have one issue however: Lena needs to be retired. I don't believe you (and zillions of others) think that there is an issue and certainly intend no offense. But I hope you can understand that it is an issue for many (and, of course, I think it should be an issue for most). Certainly there are alternative images to use.

For whatever it's worth I'm very grateful that you participate so much with this and other Mathematica forums. There does need to be a greater diversity in Mathematica users.

Jim

POSTED BY: Jim Baldwin

Nice work Sylvia. I am exploring your code for some of my work. I have most of it running on Mathematica 12 on the Mac. I am having trouble with the Hatching block of code here.

Module[{strokeLength, strokeOrient, dens, stroke},
Function[{x, y},
dens = densityFunc[x, y];
strokeLength = .03 dens^.3 (1 + .1 RandomReal[{-1, 1}]);
strokeOrient = ?/
3 + ?/6 (1 - (2 dens - 1)^2)^20 RandomReal[{-1, 1}];
stroke = strokeLength Through[{Cos, Sin}[strokeOrient]];
{{x, y} - stroke, {x, y} + stroke}
] @@@ refinedPts // Line // Graphics[{AbsoluteThickness[0], #}] &
]

For some reason there are two ? (Question marks) showing up on my browser at Line 5 and Line 6!! Could you help identify what variables these question marks correspond to? I suspect they are the source of my issue. Thanks Murali

Does anyone know what is the fix for this line:

Line[Identity @@@ MeshPrimitives[#, 1]]

My version (Windows V10.1) says:

Identity::argx: Identity called with 2 arguments; 1 argument is expected. >>

Is this related to my version of Mathematica?

Thanks in advance for any assistance.

POSTED BY: Aeyoss Antelope

enter image description here - another post of yours has been selected for the Staff Picks group, congratulations !

We are happy to see you at the top of the "Featured Contributor" board. Thank you for your wonderful contributions, and please keep them coming!

POSTED BY: EDITORIAL BOARD
Posted 9 years ago

Robot can make a color point on paper, usually a round disk. The disk diameter can be changed in a range, and the color can be changed also.

It looks like ImageEffect, "Oil Painting", but ImageEffect is not cool enough, and it has no processing output. enter image description here

I am looking for a program has a painting processing (like your last example) to guide a robot. Meanwhile, the painting points can overlap, but less overlapping area to save oil or color material. The whole painting can be made just two or three loops, then finish all the works. What's you opinion?

POSTED BY: Frederick Wu
Posted 9 years ago

Wow ! This post is almost same quality level as The Mathematica Journal.

Hi Silvia, Your last image style looks like many different squares overlapping each other with big area. I am working on a small project, we try to use a robot to make or fake a painting. If the processing like your last image, which have many small graphics overlapping each other, it's take very long time and waste too many painting material. Will you suggest some simple algorithms or methods for this kind of robot application? Maybe, the small key graphics is a disk (like a painting point). Those layers are not too many (maybe 3-4 layers), but the color points are more average distributed on a whole paper?

POSTED BY: Frederick Wu

Hi @Frederick Wu . Thank you very much for your replying!

I myself was very much interested in making a painting machine, too!

The last example in my post is not optimized for real painting purpose. A region overlapping check can be done to avoid shadowed strokes. Also, instead of square shape strokes, other shapes, such as disk, ellipse, strip, etc., can be considered.

I'm not sure how is your machine painting. Is it using real oil brushes? How will you ensure it replicate the desired stroke shape? How will you precisely mix the colors? Will you use computer vision to guide the process or is it a blind process?

I would really love to have more detailed conversation with you on the topic!

POSTED BY: Silvia Hao
Posted 9 years ago

Hi @Silvia Hao,

I am very glad to know you on Wolfram community. It's very proud, we have so good female Mathematica player in Beijing. ^_^

I find a few videos below for you, so you can just watch and get the idea of robot painting. From my view point, the robot is nothing too new, but I believe, there are still many space for AI and algorithm.

e-David Robot ( use camera with robot, possible camera feedback for color control )

http://v.youku.com/v_show/id_XNTgzNTIyNjQ4.html?from=s1.8-1-1.2

e-David Robot ( points and lines make a drawing )

http://v.youku.com/v_show/id_XNTgzNTIzNjQw.html?from=s1.8-1-1.2

7-Bot (small robot but low cost, for personal desktop )

http://v.youku.com/v_show/id_XMTMzNzIyMzcwMA==.html?from=s1.8-1-1.2

POSTED BY: Frederick Wu

Hi Fred, thank you very much! (Though I'd prefer to hide the gender behind my code :)

And thank you for your videos! Very interesting! Are you making / planning to make a painting machine yourself? It would be great if we can have some words on your thought and design shared here!

POSTED BY: Silvia Hao
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract