Newton-Raphson Weeds: This is why your trial root sometimes goes off into the rough.

-fill

You have an image and wish to fill it with something. If you are new to G'MIC and coming from a paint program like GIMP, then your first thought after encountering this command might possibly be: “Ah HA! the -fill command! At last I can fill circles and squares and things with red!”

Well yes. And no. If paint buckets come to mind, then go over to -fill_color, next aisle down, the one labeled “Color Manipulation.” That’s OK. Lots of people get lost in here first time through. G’MIC is a big place.

That isn’t to say you can’t use -fill to make an image red. It’s just that you will have do it with a little finesse:

gmic -input 2,2,1,3 -fill. 220,220,220,220,20,20,20,20,60,60,60,60

The command’s numeric parameters inject a carefully arranged data stream into a freshly minted 2 X 2 pixel, single slice, three channel image. The data are so arranged as to fill this tiny image with the X-11 color Crimson: 0xdc143c, or decimal: 220, 20, 60. The individual numeric elements take their seats in the target image by first filling channel zero, row zero, starting at column zero, then going to column one, then two and so on. You have to arrange your data stream in a like manner; -fill doesn’t know about any other ordering schemes.

G'MIC geeks may find this excessively verbose, preferring:

gmic '(220,220;220,220^20,20;20,20^60,60;60,60)'

which is terser, but dispenses with -fill entirely. Instead, the second example implicitly invokes the -input command. The input specification to which that command adheres uses commas, semicolons, carets, slashes and parentheses to describe image topology. This makes the terser example somewhat more informative at a glance.

However, since -fill is the topic of this article and we maintain a sentimental attachment to our topics, we’ll just send those interested in terse efficiency over to the aisle marked “Input-Output”. They’re having a blue-plate special on -input streams today. -fill itself does not use -input's file specification. Its world-view is limited to ordered, one dimensional streams.

You may gather from this that -fill is a fairly low-level command. Indeed. It fills images with pels. That’s it.

‘Pel’, you ask? It is shorthand for Picture Element, an indivisible quantum of an image. Could be a byte long. Could be longer. Could even be shorter. Might be a float, or an unsigned integer – or even a signed integer. The particulars of a pel’s organization stems from the image format. But a pel is not a pixel. Pixels are vector-like; they have components. Pels are scalars, Pixels are made of pels. Pels are components of pixels. Careless scriveners may put ‘pixel’ when they mean ‘pel’ but we make the distinction here because the essential item that -fill manipulates are pels. An RGBA pixel has four pels, and, we should note, -fill could touch upon that particular kind of pixel as many as four distinct times while filling an image. At its most basic, then, -fill pours pels into images from some notional stream.

The notional stream can be specified in three ways:

A literal stream of numbers – pels – that is the form we have discussed so far.
Another image already on the list – same approach, really, but the pels have been packaged as an image instead of a literal sequence.
A formula. Specifically, a processor written as a G'MIC math expression – this is really quite an interesting way. Instead of draining pels from a notional input stream, a G'MIC math expression executes, producing the requisite pel. This processor has access to the pel’s former value, or the pixel it inhabits, or the image containing the pixel, or the G'MIC image list holding the image, and the processor can access any of these data sources in any order – or none of it at all. With a math expression, you can do pretty much do what you want.

Here are these choices expressed syntactically:

-fill

pel, pel, pel,… |
[image selection] |
‘math expression’

Pel Streams

The easiest streaming method to grasp is probably not the most practical – though there are exceptions to this. You can just stream comma-separated pels into an image, an open-ended argument list.

Notice that, in contrast to the image specification streams used by -input, there are no structure markers – No opening or closing parentheses, commas do not denote columns – they just separate pels, the semicolons that separate successive columns into rows, the carets which separate a plane of rows into channels, the slashes which separate a series of channels into slices – none of those appear in -fill arguments. Nada. Ziltch. They aren’t needed. Images already have a built-in structure of slices, channels, rows and columns. -fill just – well – fills them, following the topology of the target image.

In light of this, order your data stream according to the dimensions of the target image. The largest (implicit) divisions in the pel stream are image slices. The -fill command fills the first slice of a selected image entirely before starting on the second and subsequent slices, should they exist. Within a slice, -fill populates channel zero completely before beginning on channel one and subsequent channels, if they exist. Within a channel, -fill populates the top image row completely before the second and subsequent rows. Finally, -fill populates rows from column zero through to column w – 1, with w, the width of the image. So you block out your data first by slices, within slices by channels, within channels, by rows. As a consequence of this, the pels that make up a multi-channel pixel are separated by a stride. For an image in the RGB color space, the red component is in channel zero. w×h pels downstream comes the green component. w×h pels further downstream comes the blue component. Revisit our initial example, above, for a tiny stride example.

Data streams need not be properly sized, which gives rise to a class of pattern generation. If -fill runs out of a source stream while there is still space left in the target image, it starts at the beginning of the source stream again. -fill loops over streams as many times as is necessary to fill a particular image. Once the target image is filled, the command discards the remaining excess.

So to make any one of a bazillion or so herringbone patterns, have some fun with discrete math and clock arithmetic:

gmic                             \
   -input 20,20,1,3              \
   -fill. 205,127,63,91,31,16,8  \
   -resize. 150,150,[-1],[-1],1  \
   -output. herr_00.png

gmic                                           \
   depth='{2^16-1}'                            \
   -input 17,17,1,3                            \
   -fill. 189,129,6,207,84,104,239,153,6,78,85 \
   -resize. 300,300,[-1],[-1],5                \
   -normalize 0,'{$depth}'                     \
   -output. herr_01.png

Tweak the numbers following the -fill command. Add a few. Take a few away. The game is to choose pel streams with lengths that do not evenly divide the capacity of selected images, or which do not evenly divide various image dimensions. Pel streams with lengths equal to prime numbers furnish a never-ending source of delight, should it be your wont for that kind of herringbone pattern.

In these examples, we use -resize with different interpolation schemes to upscale and filter the herringbone pattern for additional effect. In the last example, we re-normalize the image immediately before output to ensure that the pels range within the limits imposed by the 16 bit/channel Portable Network Graphics format, since the cubic interpolation scheme we specified for -resize can carry values outside this range. See Images As Datasets elsewhere on this site for a discussion on why this step is sometimes necessary.

Images

The -fill command can populate one image with the contents of another. Instead of a stream of pels, we pass an image selection parameter in square brackets. The images in the selection becomes the fill source; Targets live in the selector decorator adorning the -fill command.

Beware: -fill is not a cheap version of -array. The command effectively flattens the source into a one dimensional stream of pels. The very first pel in the source stream corresponds to column zero, row zero, channel 0 of slice 0 of the target image. The -fill command then increments, first by columns, then by rows, then channels, then slices, following the same hierarchy it imposes on literal pel streams. Indeed, as far as -fill is concerned, an image is just a convenient way to specify a pel stream. Unless the source and target images happen to have identical dimensions, -fill will not preserve image coherency. For that, investigate -array.

This notwithstanding, we can still play herringbone games by building a herringbone pattern in stages:

small_herringbone.png big_herringbone.png	gmic \ -input 13,13,1,3 \ -fill. 65,71,36,154,27,84,39,233,146,82 \ -input 256,256,1,3 \ -fill. [-2] \ -output. big_herringbone.png \ -output.. small_herringbone.png

This example illustrates using -fill with a pel stream, making an initial small herringbone image using a hand-written pel stream, then making a large herringbone using the initial pattern as an image source.

Fill Pitfalls

Here is an example you don’t want to do – unless you are being devious.

face.png	face.png is a 72 X 72 pixel image of a happy face with a nice transparent background – it has red, green, blue and alpha channels.
face-times-four.png	When the selected target image has dimensions 144 X 144, the `-fill` command will not make a 2 X 2 array of four 72 X 72 happy faces. That is, instead of the image on the left, the commands: gmic \ -input face.png \ -input 144,144,1,4 \ -fill. [-2] \ -output funny.png
funny_faces.png	produces this strange image, which is probably not what we want. So how did we go off the rails?

We’ve heard it said that -fill is a low-level command. It just streams pels into an image. That's really all it does. Essentially -fill turns source images into pel streams. Thus, the four channel, 72 X 72 happy face becomes a linear stream of 82,944 pels. Having flattened its source image, -fill then populates its selected 144 X 144 X 1 X 4 target image using the precise same hierarchy we stated before: it completely fills the first slice of the target before commencing the second, and within a slice, it completely fills channel zero before commencing on channel one, and so on.

In particular, the -fill command took the first 72 pel row of channel 0 of the source image and filled just the left half of the first row of the target image, which is, at 144 pels, twice as long. The second 72 pel row of the source image now fills the right half of the first target row. And so on. That explains the pair of flattened happy faces: the left half of the target image consists of the even rows; the right half contains the odd rows of the source image.

One quarter of the way toward filling the target’s channel zero, the source channel zero taps out. No matter. The -fill command just begins streaming in the source’s second channel, copying those second channel pels into the target’s still not entirely filled first channel. Since the dimensions of the target image are exactly twice that of the source, with the area of each target channel exactly four times that of the source, the entire happy face source image exactly fills just channel zero of the target. The source has been entirely used up. Still no matter. The -fill command goes back to the beginning of its source and commences to fill channel one of the target image through a second pass of the source.

At the end of the day, the source image is copied four times, each time filling just one channel of the target. That makes the four channels of the target image exact duplicates – hence, the target image is a study in gray of flattened happy faces. Could be they are not that happy.

The -fill command is not a blitter – a block copier. It does not think in rectangles, just streams. To do what you really want to do, try -array instead:

gmic -input face.png -array. 2,2,2 -output. face-times-four.png

Math Expressions

The third way to stream pels is really quite different. The source – well, there is hardly a source in the literal sense. The source is what an algorithm encoded in a math expression produces when given a pel. The math expression is a pel processor written in G’MIC’s math scripting language.

-fill, along with its cousin -input, both furnish hooks for pel processors, enabling fine-grained computations at the pel level. When the G’MIC command processor discovers that -fill has been given a math expression as an argument, it suspends its own processing and invokes the math expression parser. That math expression now sources a virtual stream of pels. It is invoked once for every pel in the target image and the expression can do whatever it wants to produce that pel - accessing in any manner the current or neighboring pels in the target image, pels in other images on the list, or none of these sources at all, as in the first example following, which simply pulls random numbers out of the ether without regard to the content of the target or any other image on the list. The final computation of the math expression, whatever that result may be, becomes the new value of the pel, an implicit assignment.

Here are a few examples:

colornoise.png (exerpt)	gmic -input 3840,2160,298,3 -fill. u -normalize. 0,'{2^8-1}' -split. z -output noisy.mov The pel processor is one character, ‘`u`’, which is a math expression hooked into a random number generator. It’s invoked 7,415,193,600 times, once for each pel in the image. On each invocation, it produces a random value between the limits zero and one inclusive. Since this is the last computation this processor performs (as well as its first), G’MIC implicitly assigns the resulting random value to the pel. A one character script, and you’ve written a generator for an Ultra High Definition video stream that, for ten seconds, emulates a 1967 color television set tuned between stations. I bet if you thought about it for any length of time, you’d wind up writing a much more involved script for this very same FX, assuming you wanted such a thing.
bendexter.png	gmic -input 300,300,1,1 -fill. x==y Draws a one pixel bend dexter. Pels on the diagonal have equal column and row coordinates, so, for those elements, the math expression returns 1 (True), otherwise it returns 0 (False). That is, the pel is set with the result of the comparison. Without further commands, the dynamic range of the image is from zero to one, which may not be what you want. `x` and `y` are two of many predefined variables available to math expressions, these providing the current pel's column (`x`) and row (`y`) coordinates. Similarly, there are current slice (`z`) and channel (`c`) coordinates that are predefined, while `i` furnishes the value of the pel addressed by these four coordinates. There are many other such variables – too many to list here. Mathematical Expressions of the GMIC Handbook furnishes a complete listing.
bendsinister.png	gmic -input 300,300,1,1 -fill. x==h-y Draws a one pixel bend sinister by means almost identical to the previous example, with just a variation in logic. `h` and its companion `w` are more predefined variables which tell the pel processor about the height and width of the current image.
bendboth.png	gmic -input 300,300,1,1 -fill. x==y;x==h-y The expectation may be to draw bend dexter and bend sinister by combining the two previous means. To this end, the math expression has two statements, separated by a semicolon. Combining the logic statements, one may suppose, combines the image. Alas, the results do not bear this reasoning out. G'MIC implicitly assigns the results of the last executed statement in the script. The last executed statement is the comparison of the current column coordinate to the height of the image, less the current row coordinate. That means sinister wins over dexter every time. Better luck next season.
x-eleven-crimson.png	gmic \ -input 300,300,1,3 \ -fill. '[220,20,60]' \ -output x-11-crimson.png Yet Another Way To Fill An Image With X-11 Crimson. With newer releases of G'MIC (2.0.x) math expressions acquired Vectors and Matrices. Math expressions that return vectors, such as the minimalist one here, set all of the pels of the pixel in one fell swoop. The number of iterations drop as well. Math expressions which return vectors equal to the channel count iterate over the entire image once – not once per channel. They are true 'pixel processors.' If the length of the return vector does not equal the channel count of the image, G'MIC will fill missing channels with '0' or ignore excess channels, as the case may be. Write vectors as comma separated items enclosed in square brackets. With vectors come a small host of vector and matix functions. G'MIC treats complex numbers as two element vectors, thus 3.5 + 1.27j looks like `[3.5, 1.27]` and `cabs([3.5, 1.27])` finds the hypotenuse (i.e., 'complex absolute value'): 3.7233, or there abouts, so if you have been writing `sqrt(a^2+b^2)` a lot, you may find `cabs([a,b])` a terser equivalent. Skip forward to the closing example of filling images with Newton-Raphson basins of attraction for more working examples of vectors, and vectors as complex numbers.
	gmic \ -input 256,125,1,1 \ -fill. x \ -normalize. 0,255 \ -output left2right_ramp.png Obtain ramps by the simple expedient of citing a current coordinate, which G'MIC implicitly assigns to the current pel. The low column coordinates begin at zero and end with w – 1 – dark to light. The dynamic range of the image depends on image geometry, which rarely reconciles with the dynamic range of image formats. Pel processors return values to images exactly as computed, no more, no less. It is on your watch to ensure that the image works within the native dynamic range of your chosen output format. In many cases, the results of simple expedients can be set to rights by the judicious application of `-normalize` See also Images as Datasets for the usual cautionaries.
	gmic \ -input 256,125,1,1 \ -fill. y \ -normalize. 0,255 \ -output top2bottom_ramp.png Obtain a top to bottom ramp by the same expedience as the prior example, just substituting the current row coordinate, `y`, for the current column coordinate, `x`.
	gmic \ ang=32.6 \ -input 256,256,1,1 \ -fill. *rad=2pi$ang/360;xcos(rad)+ysin(rad)*\ -normalize. 0,255 \ -output. anyangle_ramp.png Ramps at arbitrary angles may be obtained by utilizing other math expression features. A pel processor may have multiple statements, these separated by semi-colons. In the last example, the first statement converts its argument, the desired orientation of the ramp, from degrees to radians. The second statement references the converted argument, along with current column and row coordinates, to produce a pel value. The math parser then assigns the evaluation of this second statement to the pel addressed by the current values of `x`, `y`, `c`, and `z`. So far, we've seen just predefined variables, such as `pi`, the ratio of the circumference of a circle to its diameter. In this expression, we have defined our own, harnessing two species of variables, in fact. The native variant is `rad,` created in the math expression and then referenced in the math expression. Being native to the math expression, It does not take a `$` prefix; adding one may or may not raise a syntax error. The difficulty (see below) is with how the command processor pre-processes math expressions, substituting the contents of command level variables appearing in math expressions for the variables themselves. Since an inadvertent `$` prefix on a math expression variable most likely references a never-before seen command level variable, G'MIC completes the substitution by supplying the default values for newly-minted command level variables: empty strings, which in numerical contexts are zero. So-substituted, the resulting math expression is syntactically correct, more often than not, but likely not behaving as you want. The scope of these math expression variables is limited just to the math expression itself – they furnish a local, non-persistent scratchpad that disappears when the math expression ceases execution, and their values do not persist from one invocation of the pel processor to the next. The second variant is an import from the command environment. Command level variables like `ang`, above, can be referenced in the math expression when adorned with their `$` prefix. In fact, behind the scenes, the command processor performs a substitution pass on math expressions, looking for variables with `$` prefixes and substituting for these their current values. Consequently, command level variables turn into numeric constants before the math expression is invoked. Command level variables, then, do persist between invocations of pel processors, but pel processors have no real knowledge of their existence, because command level variables manifest themselves as numeric constants within math expressions, not in their original form.
turb.png turb_xor.png	gmic \ -input 256,256,1,2 \ -turbulence 100,9,4 \ -input [0],[0],1,1 \ -fill. 'xor(i(#-2,x,y,0,0),i(#-2,x,y,0,1))' \ -normalize. 0,255 \ -output. turb_xor.png This math expression harnesses two of the many builtin functions available to such. The builtin `xor(a,b)` performs a bit-wise exclusive OR operation on its arguments, rounding them down to the nearest integer, if necessary. `i(x,y,z,c,interpolation_type, boundary_conditions)` samples the image at a specified set of coordinates, which need not be integral and need not be on the image. The coordinates may have fractional parts, implying sub-pixel sampling. In this case, the `interpolation_type` flag chooses between nearest neighbor and linear interpolation. The `boundary_condition` flag furnishes policy when the sampling coordinates are off-image. Images Have Edges, Now What? reviews the various types, but briefly: Dirichlet (0) means all off-image samples are black, by definition. Neumann (1) extends the values of the images out to infinity; off-image samples reflect the values of those pixels on the nearest edge. Cyclic (2) treat images as torii; off-image samples wrap around to the far edge. Mirror (3) furnishes off-image samples by reflecting the off-image position back into the interior of the image. For terseness, arguments may be omitted, counting from the right, so `i(x,y)` references a pel at column `x` and row `y`. Omitted, slice (`z`) and channel (`c`) default to zero, nearest neighbor interpolation and Dirichlet conventions at the boundary. Math expressions are not limited to their operand image. Any predefined variable may be adorned with a `#<image selector>` suffix to reference any other image on the list. The numbering conventions follow those of the command processor's selector. Positive indices count from the beginning of the list; negative indices count from the end of the list. On a list with four images, `iM#1` refers to the maximum intensity of the second image on the list, counting from the beginning (the image list follows the zero-indexing convention). `iM#-3` references the same image, counting from the end of the list. This convention extends to the image sampling function `i(...)` and `j(...)` where the initial argument can be a selector. That is the key facet of this example, which samples each pel of a two channel image other than the operand, containing synthesized turbulence, and XOR's them together, producing the corresponding pel of the operand image.

This has been kind of a Whitman Sampler of G'MIC math expressions. Since this article is about the -fill command, it is not a propitious place to lay out a full-blown math expression tutorial – for such would overwhelm the -fill command itself. The definitive and complete, if terse, reference to math expressions remains the Mathematical Expressions section in the GMIC Handbook. Everything's there, sans examples and artfully illustrative errors. The key takeaway here is that the -fill command lets you hook up a pel or pixel level processor to your image, so you can operate on the raster itself.

New and Improved Processors

In the run-up to G‘MIC 2.0, released last November, 2016, the features and capabilities of the math processor grew considerably. The examples presented thus far, one might say, are ‘classic’ G‘MIC math expressions – capable of a series of operations but lacking looping constructs, dynamic variables, being somewhat sensitive to whitespace and limited to returning just scalars. In the first-generation math processor, an RGBA pixel could be visited by a pel processor as often as four times, since the processor could alter only one pel at a time.

Most of these limitations have been lifted. Math expression can now utilize control structures, query other images in the pipeline, support vector data types as well as scalars, so that an entire pixel can be a return value, not just one pel. While it does not have a direct bearing on semantic content, math expressions can now also be written with a liberal amount of white space, so that the structure of the math expression is clear.

To illustrate some of the newer features, we've written a G‘MIC command that plots the basins of attraction that arise from Newton-Raphson root polishing. This computation is inherently iterative, as the only practical way to find the roots of higher degree polynomials is to start with a trial value and, exploiting certain properties of polynomials and their first derivatives, iteratively refine intermediary trials. Iteration was not something that could be built into a math expression until recently.

The Newton-Raphson iteration:

z_{i + 1} = z_{i} - \frac{f (z_{i})}{f^{'} (z_{i})}

says we may find a better trial value for a root λ_k of polynomial f(z), z_i+1, by subtracting a “correction factor” from a current trial value, z_i. We obtain the correction factor by evaluating the polynomial and its first derivative at z_i and taking the ratio of the evaluation and derivative. Graphically, we extend the tangent of the function at the evaluation point to the parameter axis, viz:

In this illustration, f(z) is almost linear in the neighborhood of the root and it is not hard to imagine how extending the tangent from a pretty good guess gives us a new guess even closer to the root. Clearly, evaluating again at z_i+1, and extending a new tangent at that point will walk us ever closer to the root λ_k and it will only be a matter of taste to decide, at some point, that the trial value is close enough.

But what is a pretty good guess? One way to answer that question is to draw a map: -newtonraphson, listed below, when given a set of roots, paints such a map on the complex plane, using a grid of values (say, the position of pixels) as initial trial values. It performs a series of Newton-Raphson iterations on each grid point until an intermediary is really close to one of the given roots. When an iteration finishes, -newtonraphson reports the result by plotting a two-channel pixel in a new image. Channel zero is an index uniquely identifying one of the given roots. Channel one holds the iterations it took to arrive at the root from the initial grid point. -newtonraphson plots the pixel in the same position as the grid point, essentially reporting how that point behaves under Newton-Raphson. Stand back from this image, squint, and we should see the aggregate of grid points which converge to particular roots – the basins of attraction.

So a basin of attraction is a set of complex numbers which converge to a particular root under a Newton-Raphson iteration. It happens that the division of the complex plane into the various basins of attraction for a polynomial form striking fractal boundaries. In fact, the fractal boundaries separating the basins of attraction underlie why Newton-Raphson root polishing can sometimes be ill-behaved. Pick an initial trial in the midsts of the fractal territories, and convergence, at least initially, might not be that convergent in the early part. As we will see later, the initial trial may not converge at all.

If you are not familiar with how these basins of attraction can be imaged, read the Wikipedia article on Newton fractals, and, in particular, see Simon Tatham’ page, Fractals derived from Newton-Raphson Iteration. We ~~stole~~ borrowed quite a few pointers from Mr. Tatham's page.

To begin, include this G‘MIC command in your Mac OSX, Linux: $HOME/.gmic file, or, for Microsoft Windows your %APP_DATA%/<user>.gmic file, where <user> stands in for your login name:

#@gmic
#@cli newtonraphson : [roots],_center_x,center_y,_mag,_maxiter
#@cli : Draw on the selected images the basins of attraction of the polynomial
#@cli : roots given in the argument image using Newton-Raphson root
#@cli : polishing. Argument image is N,1,1,2 - an N vector for N, complex,
#@cli : number of roots.  Channel zero contains real, and channel one
#@cli : imaginary components. Default values: 'center_x = 0', 'center_y = 0',
#@cli : 'magnification = 1', 'maximum_iterations = 100'
newtonraphson : -check ${-is_image_arg\ $1} -skip ${2=0},${3=0},${4=1},${5=100}
   -verbose -
   sig=s
   -if {$!<2}
      sig=
   -endif
   -verbose +
   -echo[^-1] "Use Newton-Raphson root polishing on image"$sig$?" to paint basins of attraction associated with the roots given by $1."
   -verbose -
   ppu=255.5
   ox=$2
   oy=$3
   mag=$4
   lim=$5
   cx={$ox-w/(2*$ppu*$mag)}
   cy={$oy+h/(2*$ppu*$mag)}
   err=1e-10
   -pass$1 1
   -repeat {$!-1}
      -local[$>,-1]
         -fill[0] "Q  = vectorw#-1();
                   ic = 0;
                   zz = ["$cx"+x/("$ppu"*"$mag"),"$cy"-y/("$ppu"*"$mag)"];
                   dowhile(
                            acc = [0,0];
                            for(
                                  k=0, k < w#-1, k++,
                                  Q[k] = cabs(zz-I(#-1,k,0));
                                  acc += [1,0]//(zz-I(#-1,k,0))
                               );
                            dif  = [1,0]//acc;ic+=1;
                            zz   = zz - dif,
                            cabs(dif) > "$err" && ic < "$lim"
                          );
                  if( 
                      ic < "$lim",
                      [argmin(Q),ic],
                      [w# - 1,"$lim"]
                    )"
      -endlocal
   -done
   -rm.
   -verbose +

Implementation Notes

If you look at the math expression which we pass to -fill, you might wonder how it is we're doing Newton-Raphson. We're in the business of plotting basins of attraction, not finding roots – in fact, we start from known roots. Since we have the roots, we can rephrase the Newton-Raphson process to suit our convenience. And it is convenient for us to write the polynomial in factored form, each factor being a difference between the trial value, z_i, and a member of the set of given roots, {λ_n}:

$f (z_{i}) = \prod_{k = 0}^{n - 1} z_{i} - λ_{k}$
$f^{'} (z_{i}) = \sum_{j = 0}^{n - 1} \prod_{k = 0}^{n - 1, j \neq k} (z_{i} - λ_{k})$
$\frac{f^{'} (z_{i})}{f (z_{i})} = \sum_{k = 0}^{n - 1} \frac{1}{(z_{i} - λ_{k})}$
$z_{i + 1} = z_{i} - \frac{f (z_{i})}{f^{'} (z_{i})}$ → $z_{i + 1} = z_{i} - \frac{1}{\sum_{k = 0}^{n - 1} \frac{1}{(z_{i} - λ_{k})}}$

In this form, the denominator of the correction factor has an interesting structure. It is a sum of reciprocals. Should, at some iteration i, the trial value z_i approaches any one of the roots λ_k, the difference between the trial and that root commences to vanish. That drives the corresponding reciprocal to a very large value. The sum of reciprocal differences itself grows very large, being dominated by that one exploding term, so that the correction that is applied to the trial value for the next iteration becomes vanishingly small. On that condition, we exit the dowhile() loop – as we have found a pretty near enough approximation to a root.

Upon that exit, we prudently check if the overall iteration count has exceeded the limit set in $lim. If it does, we return a failure vector. The failure vector uses the root count, n, as an index, which is not assigned to any root. These are indexed from 0 to n – 1. For an iteration count, the failure vector contains the value of $lim. Downstream rendering can mark pixels with these two special channel values however they see fit.

Otherwise, we return a normal success vector. To correctly identify which root the trial value approximates, we keep off to one side a check vector, Q, with components |(z_i – λ_k )|. Applying argmin() to this vector returns the index of the smallest component in the vector, identifying the root to which z_i converges. This index becomes the first component of the return vector. The number of iterations is recorded in the expression variable ic and becomes the second component of the return vector.

Why would this expression fail to converge? Our game relies on extending tangents to the parameter axis, then walking straight up until we hit the curve, our new evaluation point. Is it possible that this new evaluation point could be the same as an earlier trial value? Oh, yes, indeed, in a significant class of polynomials, there are honey-pot traps for those taking the Newton-Raphson route: trial values z_i transform to z_{i – k} - a prior trial value tested k steps before the present trial. In that case, we are in a circular iteration of k steps – and we will never get out.

This maladay is manifested by so-called connected Fatou sets that can also appear in our maps. Entwined among the basins of attraction of every degree n polynomial, f_n(z) is a Fatou set of one lower degree. It is typically disconnected, the set forming Cantor dust., but can be connected. That is, the set is bounded by a connected Julia set which cleaves the complex plane into disjoint sets. One set, the trial values that are parts of honey-pot loops, form the interiors of the Fatou sets. Trial values within these sets never converge on a root, they just boomerang around. With disconnected Fatou sets, convergence is almost always possible; it is vanishingly unlikely to choose a dust mote of Cantor dust and bounce around forever.

Failure of convergence because we're gooed to a honey-pot is not catastrophic for our purposes. We are painters of basins of attraction, and, by extension, basins of non-attraction. If a trial value runs to the weeds, we just paint it another pretty color, an additional bit for the rikka.

Rendering Basins

The image that -newtonraphson produces is not strikingly beautiful. It is intended to be rendered. By design, we separate computation and presentation into distinct tasks. A basic rendering script might employ -map to associate colors to the root-assignment channel, then alter this coloration as a function of the iteration values. Here's a use of -newtonraphson that goes along those lines:

$ gmic                   \
   wp='{2^16-1}'         \
   -input 300,300,1,2    \
   -input '(1,-0.5,-0.5^0.0,0.866025403784,-0.866025403784)' \
   -newtonraphson[-2] [-1],-0.75,0,3.5                       \
   -rm.                  \
   -split. c             \
   -input '(150,250,200^30,260,150^180,7,220)' \
   -map[-3] [-1]         \
   -rm.                  \
   -mul. -1              \
   -normalize. 0,1       \
   -pow. 0.75            \
   -mul[-2,-1]           \
   -normalize. 0,'{$wp}' \
   -output. nr_00.png

Here, we are depicting the basins of attraction for the cubic polynomial z³ – 1 which has three roots in the complex plane at $1 + 0 j$ , $\frac{- 1}{2} + \sqrt{\frac{3}{4}} j$ and $\frac{- 1}{2} - \sqrt{\frac{3}{4}} j$ .

And this is the result:

$Newton-Raphson fractal$

Newton-Raphson fractal

New Math Expression Features

Creating a map of how a grid of values behaves under a Newton-Raphson iteration was not possible with earlier versions of G'MIC; the basic "do until these trial values converges to a root", which lies at the heart of -newtonraphson, requires a proximity test at the end of each iteration. Now the language offers dowhile(...) and for(...) for such open-ended, "until such-and-so prevails" testing. dowhile(...) executes an iteration, then tests; for(...) tests, then executes an iteration. whiledo(...) is a convenience variant of for(...) with an argument list that mirrors dowhile(...).

With these, math expressions may perform conditional iteration, giving rise to more elaborate algorithms. To accommodate longer scripts, the math expression parser is more flexible with white space, so that expressions may be formatted in a way to illuminate their structure and workflow.

Expanded vector and matrix datatypes turn math expressions into true pixel processors, as math expressions which return vectors can set the entire pixel in one atomic assignment. This reduces iteration counts as it is not necessary to iterate over individual channels.

Control Loops

The math expression control loop constructs are implemented as functions which take expressions as arguments and return the result of the last executing statement, or, in the case of for(...), the Not_A_Number token for the possibility that an initial terminating condition is positive at the outset, so that none of the other expression arguments execute. for(...) tests and executes; dowhile(...) executes and tests. whiledo(...) is a convenience function that, like for(...), tests and executes, but has an argument list that mirrors dowhile(...), and, in certain circumstances, expresses logic more clearly and is easier to set up, having two arguments instead of four. Like for(...), whiledo(...) could return NAN if its condition expression initially terminates without ever executing the procedure body.

Whitespace

If we take a deep breath and squint sideways, we see that the form of the -fill command is no difference in quality than what has gone on before:

-fill[<targets>] “<math_expression>”

which can be rendered at a slightly higher resolution:

-fill[<targets>] “<math_partial_expression>”$<command_variable>“<math_partial_expression…>”

Here, a double-quotation marks alert the math expression parser that white space may be used to "pretty-format" the math expression. Historically, G'MIC math expressions had all the charm of Perl Regular Expressions, which, from a casual glance, had the appearance of those final sputter of characters before communication lines drop. Double quoted math expressions may be broken into several lines, indented, and, with care, have white space inserted around operator tokens.

Double-quoted math expressions require special precautions when mixing with command level variables, as the sequence $<variable_name> will be protected from substitution when such appears within double quotes. Mixing in command level variables entails furnishing a closing quote to the partial math expression immediately preceding the command level variable, then an opening quote immediately following the command level variable. One can easily, and mistakenly, invert this context and think of the command line variable being quoted instead. When performing its substitutions, the command parser concatenates the substituted variables and the math sub-expressions, making one complete math expression, which is then given to the math expression parser.

Vectors

Lists have been with G'MIC math expressions for some time, which supported basic, vector-like operations. G'MIC supports vectors more formally with the 2.0 release.

Specify vectors by bracketing lists with square brackets '[' and ']'.

gmic foo='{[3,4,9,2,1]}'
    [gmic]-0./ Start G'MIC interpreter.
    [gmic]-0./ Set local variable foo='3,4,9,2,1'.
    [gmic]-0./ End G'MIC interpreter.

Or use the vector(...) constructor:

gmic foo='{vector(3,4,9,2,1)}'
[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ Set local variable foo='3,4,9,2,1'.
[gmic]-0./ End G'MIC interpreter.

Appending a number to the constructor name creates vectors of specific length:

gmic foo='{vector7(3,4,9)}'
[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ Set local variable foo='3,4,9,3,4,9,3'.
[gmic]-0./ End G'MIC interpreter.

You do not have to provide an argument list equal to the specified length of the vector. The math expression parser iterates over the argument list as many times as necessary to automatically complete the vector. Here, we specified a vector length of seven elements, but provided only three. The math expression parser iterated over the argument list three times (creating nine elements), then truncating the 2 excess elements in final iteration, making a vector of the desired length of seven.

The number following the constructor name can be one of the predefined positive environmental variables such as those describing image topology. Observe that in -newtonraphson: we set the check vector, Q to the length of the argument image:

```
Q  = vectorw#-1()
```

Note the expression w#-1, appended to the constructor name where a literal number could have been written. With this notation, the length of the check vector can be however long it needs to be to accommodate the number of roots contained in the argument image.

Note too that the argument list is empty in this example. In this case, the math expression parser creates a vector with zero-valued elements.

The math expression parser performs element-wise operations with vectors:

gmic five_sixes='{vector(1,2,3,4,5)+vector(5,4,3,2,1)}'
[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ Set local variable five_sixes='6,6,6,6,6'.
[gmic]-0./ End G'MIC interpreter.

gmic zero_cross='{vector(1,2,3,4,5)-vector(5,4,3,2,1)}'
[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ Set local variable zero_cross='-4,-2,0,2,4'.
[gmic]-0./ End G'MIC interpreter.

gmic up_slope='{vector(1,2,3,4,5)/vector(5,4,3,2,1)}'
[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ Set local variable up_slope='0.20000000000000001,0.5,1,2,5'.
[gmic]-0./ End G'MIC interpreter.

gmic peak='{vector(1,2,3,4,5)*vector(5,4,3,2,1)}'
[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ Set local variable peak='5,8,9,8,5'.
[gmic]-0./ End G'MIC interpreter.

The math expression parser throws an exception if vectors in these kind of pairwise operations do not have equal lengths.

Scalars may offset or scale vectors:

gmic below_zero='{vector(1,2,3,4,5)-6}'
[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ Set local variable below_zero='-5,-4,-3,-2,-1'.
[gmic]-0./ End G'MIC interpreter.

gmic triple_thensome='{vector(1,2,3,4,5)*3.5}' 
[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ Set local variable triple_thensome='3.5,7,10.5,14,17.5'.
[gmic]-0./ End G'MIC interpreter.

You can address specific vector scalar components with square bracket operators.

gmic twocube='{cubes=vector(1,2,3,4,5)^3;cubes[1]}'
[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ Set local variable twocube='8'.
[gmic]-0./ End G'MIC interpreter.

You can address a contiguous sub-vector with a comma operator inside square brackets

gmic three_middle_cubes='{cubes=vector(1,2,3,4,5)^3;cubes[1,3]}'[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ Set local variable three_middle_cubes='8,27,64'.
[gmic]-0./ End G'MIC interpreter.

Two-element vectors can be interpreted as complex numbers. a + bj = [a, b]. Use complex operators ** (complex multiplication), ^^ (complex power) and // (complex division) to distinguish complex multiplication, power or division from their element-wise counterparts. element-wise addition and substraction are equivalent for both vector2 and complex objects. -newtonraphson employed cabs([a, b]) to compute the complex absolute value; cconj([a, b]) returns the complex conjugate a – bj. cexp([a, b]) returns ε^a+bj.

As noted above, pixels are vector-like. I(x, y, z, interpolation_type, boundary_conditions) samples the image at column x, row y and slice z, returning a vector with elements corresponding to the channels of the image, the first element corresponding to channel zero. Like its pel counterpart, i(x, y, z, c), I(x, y, z) can sample between pixels and takes fractional arguments in support of this. To this end, optional fourth and fifth flag parameters set sampling behavior. The default, interpolation_type is nearest_neighbor and off-image pixels are black (Dirichlet), both these being requested by zero (0) flags. An interpolation_type of 1 induces linear interpolation of the sample point between neighboring pixels. A boundary_condition of 1 extends the values of edge pixels out to infinity. Under this policy, an off-image sample point is that of the nearest edge pixel (Neumann). A boundary flag of (2, Cyclic) presumes that the image has the manifold of a doughnut, or, more formally, a torus. An off-image sample point wraps around past the far edge and draws data from the corresponding modulo point on the other side of the image. Finally, a boundary flag of (3, Mirror) reflects an off-image sample point off the edge back into an image.

Postscript

The G'MIC commands -fill and -input both provide bridges from the command level, where commands operate on images, and the much more fine-grained level of operations on pels and pixels. Even in the earliest versions of G'MIC the math expression parser furnished a limited means of pixel- and pel-level computations, but the scope was akin to scratch-pad computations without a great deal of algorithmic elaboration possible. Perhaps you may have had an exotic convolution kernel in mind. Functions i(x, y, z, c) and j(dx, dy, dz, dc) helped with sampling around a central pixel, and a decent repertoire of mathematics functions were on hand to implement your ideas. That said, G'MIC as an exploratory tool of algorithms operating on regularly sampled pixel spaces needed to grow in at least two realms: a physical realm, to avail itself of as many core a computer might provide, and an algorithmic realm, to avail itself of the increasingly complex ideas that a graphics filter designer may commit to code. To expand in these two realms, the G'MIC math expression parser underwent about a two year evolution. Somewhere midway through that evolution, David Tschumperlé posted to his blog Image processing made easier with a powerful math expression evaluator, which relates many of the motivations he had for undertaking the evolution. Nearly two years on, it is still a good read for those new to this “deep” side of G'MIC.

The limiting factor now is the dearth of good working examples of pel and pixel processors, presented in a nice, overarching framework. These will be forthcoming here in the sweet bye-and-bye.

Garry Osgood