Filling holes in priority order

Filling pixels in priority order improves quality.

Description
Arguments
Examples
Photographs
Use which module?
Cleanup
Possible improvements
Acknowledgments
Source code

The general technique is also known as infilling or inpainting.

This is a process module. The source code is provided below: fillholespri.c. For details of integrating this with IM, see my page Process modules.

Description

This process module, fillholespri, is very similar to fillholes. Both modules fill holes from the edges inwards by selecting a square of pixels around each edge pixel and searching for this as a subimage in the remainder of the image. When the best exemplar is found, pixels are copied to that area of the hole. It might copy just the central pixel, or the entire window, or a small sub-window.

The simpler fillholes module fills all the edge pixels in a hole before moving inwards and filling the lext layer -- it proceeds like the rings of an onion. See the description at Fill holes.

This fillholespri module prioritises the edge pixels. It processes the pixel with the highest priority, searching and copying exemplar pixels as above. After each copy, some pixels are no longer on an edge (because they have been filled) and some pixels that were inside the hole have become on the edge, so it re-calculates the priority of affected pixels. Then it moves on to fill the pixel with the next highest priority.

The pixels at the edges of holes are assigned priority, not the holes themselves.

The priority order is the only functional difference between the modules. They both work from hole-edge inwards, and they both copy exemplars from the input image, with the same search algorithm. The priority module will complete the important parts of the onion ring, and will continue filling from there, before returning to the unimportant parts of the outer ring.

What defines the importance of pixels at hole edges? The product of three factors:

The confidence we have that an opaque pixel is correct. This is set to 1.0 for pixels that are initially opaque, and is multipled by a factor less than 1.0 for each generation.
The energy of an opaque pixel.
The number of opaque neighbours that an edge pixel has, divided by 8.

For edge pixels, confidence is inherited as the lowest confidence of any neighbouring opaque pixel multiplied by the factor, and energy is the highest energy of any neighbour.

And what is the energy of an opaque pixel? It is the highest absolute slope magnitude with any of its opaque neighbouring pixels. In turn, the slope is the difference between colour channel values divided by distance to the neighbour, which is 1.0 where the neighbours share an edge, or sqrt(2) where they don't.

The module works but isn't yet fully stable. For example, when initially reading I define an edge pixel to be a transparent pixel that is 8-connected to an opaque pixel. Then, when processing pixels, it seems better to employ a 4-connected test. I'm not yet sure about such design questions.

fillholespri needs more memory than fillholes, approximately 50 bytes per pixel of the source image, so a 35 megapixel image needs about 1.75GB extra memory. Optimisations could be added to reduce this, eg making the priority queue dynamic.

For the same images and the same options, fillholespri will be slower than fillholes. It will do the same number of searches (which account for the majority of the processing time), but has greatly increased housekeeping per pixel. For example, after copying a block of pixels (which could be 1x1), all those pixels must be removed from the priority queue, the pixels on the edge of the block must have energy re-calculated, and pixels immediately outside the block may have become edge pixels so need to have their priority re-calculated and to be inserted into the queue or to change their position in the queue.

To increase the speed, we need to reduce the number of searches or the size of each search. fillholespri should give higher quality results than fillholes, so there is more opportunity for choosing options that may lower quality but improve performance.

Arguments

For the options, and their descriptions, see the Fill holes page.

Both modules take the same options, and the options have the same defaults.

The auto_limit_search option is on by default. This records the bounding box of the first N matches, where N is the initial number of edge pixels. Subsequent searches are made only within that box.

This can give unfilled pixels, eg if holes with high-priority edges have solutions in one part of the image, but holes with low-priority edges (eg blank sky) have solutions in another part. The module attempts to correct for this. (It realises that applying the auto_limit_search, after the ordinary lsr, would reduce the search space to zero, so doesn't apply the auto_limit_search.)

Examples

This page shows the same examples, with the same options, as those on the Fill holes page.

Make a test image, with a hole. %IMG7%magick ^ -size 100x100 xc:lime ^ -stroke Red -draw "line 0,0 99,99" ^ -stroke Blue -strokewidth 3 -draw "line 0,20 69,99" ^ ( +clone -fill White -colorize 100 ^ -fill Black -draw "Rectangle 30,20 79,69" ^ -threshold 50%% ^ ) ^ -alpha off ^ -compose CopyOpacity -composite ^ fhp_src1.png
Fill any holes. %IM7DEV%magick ^ fhp_src1.png ^ -virtual-pixel None ^ -process 'fillholespri' ^ fhp_src1_fh.png Unlike `fillholes`, the result with default options is perfect.
As previous, but with homogeneity check off. %IM7DEV%magick ^ fhp_src1.png ^ -virtual-pixel None ^ -process 'fillholes hc off' ^ fhp_src1_fhnhc.png This makes the result worse.
A larger radius cures the problem with the blue line. %IM7DEV%magick ^ fhp_src1.png ^ -virtual-pixel None ^ -process 'fillholespri hc off window_radius 2' ^ fhp_src1_fh2.png
We can make an animation of the infill process. %IM7DEV%magick ^ fhp_src1.png ^ -virtual-pixel None ^ -process 'fillholespri wr 2 write2 fhp_x_fr_%%06d.png' ^ null: %IMG7%magick ^ fhp_x_fr_*.png ^ -layers optimize ^ fhp_x.gif

The animation illustrates the priority order.

A more complex graphic, with five holes.

%IMG7%magick ^
  -size 180x180 xc:lime ^
  -fill Black ^
  -draw "translate  60,60  circle 0,0 0,25" ^
  -draw "translate  60,120 circle 0,0 0,25" ^
  -draw "translate 120,60  circle 0,0 0,25" ^
  -draw "translate 120,120 circle 0,0 0,25" ^
  ( +clone ^
    -fill White -colorize 100 ^
    -fill Black ^
    +antialias ^
    -draw "translate  30,30  circle 0,0 0,20" ^
    -draw "translate 150,30  circle 0,0 0,22" ^
    -draw "translate  30,150 circle 0,0 0,24" ^
    -draw "translate 150,150 circle 0,0 0,26" ^
    -draw "rectangle 60,60 119,119" ^
  ) ^
  -alpha off -compose CopyOpacity -composite ^
  fhp_src2.png

%IM7DEV%magick ^
  fhp_src2.png ^
  -virtual-pixel None ^
  -process ^
    'fillholespri' ^
  fhp_src2a_fh.png

The window radius is too small.

The smallest radius to get the desired result is 4.
(fillholes needs at least 6.)

%IM7DEV%magick ^
  fhp_src2.png ^
  -virtual-pixel None ^
  -process ^
    'fillholespri window_radius 6' ^
  fhp_src2_fh.png

The result looks like exactly what we would expect.

The module's verbose option writes some text information to stderr:

%IM7DEV%magick ^
  fhp_src2.png ^
  -virtual-pixel None ^
  -process 'fillholespri window_radius 6 verbose' ^
  null: >fhp_verb.lis 2^>^&1

fillhole options:   window_radius 6  limit_search_radius 0  auto_lsr on  search entire  hom_chk 0.1  refine on  similarity_threshold 0  copy onepixel  copy_radius 0  verbose
pfh->WindowRadPix = 6
pfh->Match.limSrchRad = 0x0
pfh->CopyRadPix = 0
pfh->CompWind.WorstCompare = 1.61611e+27
pfh->sqDim = 13
pfh->sqCopyDim = 1
pfh->Match.RandSearchesI = 0
pfh->Match.SkipNumI = 0
pfh->Match.simThresholdSq = 0
pfh->Match.dissimThresholdSq = 1.61611e+27
Alloc 180x180
Populate
End populate

Photographs

A reduced-size photo for tests. set fhp_PHSRC=fhp_phsrc.png if not exist %fhp_PHSRC% %IMG7%magick ^ %PICTLIB%20130117\AGA_1135.jpg ^ -resize 400x400 ^ +depth ^ %fhp_PHSRC%
Add holes, and fill them, with default settings. %IM7DEV%magick ^ fh_ph0.png ^ -process 'fillholespri' ^ fhp_ph0_fh.png call StopWatch 0 00:07:27 StopWatch shows the elapsed time in `days hours:minutes:seconds`. The result looks horrible because the window_radius is too small.
Increase the `window_radius`. %IM7DEV%magick ^ fh_ph0.png ^ -process 'fillholespri wr 5' ^ fhp_ph1_fh.png call StopWatch 0 00:35:41 The visual result is good but performance is bad.
Copy the entire `window_radius`. %IM7DEV%magick ^ fh_ph0.png ^ -process 'fillholespri wr 5 copy window' ^ fhp_ph1a_fh.png call StopWatch 0 00:01:24 This has decreased the quality.
Copy just the central part of the window. %IM7DEV%magick ^ fh_ph0.png ^ -process 'fillholespri wr 5 copy window cr 3' ^ fhp_ph1b_fh.png call StopWatch 0 00:02:44 Quality has improved.
Limit the search radius. %IM7DEV%magick ^ fh_ph0.png ^ -process 'fillholespri wr 5 lsr 10%%' ^ fhp_ph2_fh.png call StopWatch 0 00:00:44 Limiting the search radius improved performance by a factor of about 30. It has also improved the infill on the left of the tower.
As previous but with random-search. %IM7DEV%magick ^ -seed 1234 ^ fh_ph0.png ^ -process 'fillholespri wr 5 lsr 10%% search random' ^ fhp_ph2r_fh.png 0 00:00:04 Very quick.
As previous but with skip-search. %IM7DEV%magick ^ fh_ph0.png ^ -process 'fillholespri wr 5 lsr 10%% search skip' ^ fhp_ph2k_fh.png 0 00:00:05 Also very quick.
Another example, with noise. When holes are close together, we need a small `window_radius`. %IM7DEV%magick ^ fh_ph3.png ^ -process 'fillholespri lsr 10%%' ^ fhp_ph3_fh.png call StopWatch 0 00:01:04 The output contains noise, most evident in the flagpole. But the result is impressive.
This has left some holes unfilled. We run `fillholespri` again, with default options. %IM7DEV%magick ^ fhp_ph3_fh.png ^ -process 'fillholespri' ^ fhp_ph3x_fh.png call StopWatch 0 00:00:12

In the current algorithm, transparent pixels within a search window create a "not found" condition, so we need small window_radius when holes are close together. But large holes need a large window_radius.

When an image has widely-spaced large holes and close-spaced small holes, there is no ideal window_radius. The module can be run multiple times, with increasing limit_srch_radius or window_radius.

Increment limit_srch_radius.

 %IM7DEV%magick ^
  fh_ph4.png ^
  -process 'fillholespri wr 1 lsr 1 als off' ^
  -process 'fillholespri wr 1 lsr 2 als off' ^
  -process 'fillholespri wr 1 lsr 3 als off' ^
  -process 'fillholespri wr 1 lsr 4 als off' ^
  -process 'fillholespri wr 1 lsr 5 als off' ^
  -process 'fillholespri wr 1 lsr 6 als off' ^
  -process 'fillholespri wr 1 lsr 7 als off' ^
  fhp_ph4_fh.png

call StopWatch

0 00:00:27

Very fast, but the edges are noisy.

Increment window_radius.

 %IM7DEV%magick ^
  fh_ph4.png ^
  -process 'fillholespri wr 1 lsr 10%%' ^
  -process 'fillholespri wr 2 lsr 10%%' ^
  -process 'fillholespri wr 3 lsr 10%%' ^
  -process 'fillholespri wr 4 lsr 10%%' ^
  -process 'fillholespri wr 5 lsr 10%%' ^
  -process 'fillholespri wr 6 lsr 10%%' ^
  -process 'fillholespri wr 7 lsr 10%%' ^
  fhp_ph5_fh.png

call StopWatch

0 00:01:05

Slower, but the result is less noisy.

A "hole" doesn't need to be surrounded by image pixels. We can extend an image, for example to the right and bottom:

%IM7DEV%magick ^
  %fhp_PHSRC% ^
  -background None ^
  -size 10x1 xc:None +append +repage ^
  -size 1x10 xc:None -append +repage ^
  -process 'fillholespri wr 5 lsr 10%%' ^
  fhp_ph6_fh.png

call StopWatch

0 00:00:42

An example with a more "natural" image:

toes.png.

With Gimp, erase (make transparent) all non-grass pixels.

toes_holed.png.

The erasing process does not need to be precise, and can overflow into the grass.

The hole covers most of the image, so limit_search_radius would leave many pixels unfilled. We might search the entire image, or use limit_search_radius but repeat the process until no unfilled pixels remain.

Search the entire image. %IM7DEV%magick ^ toes_holed.png ^ -process 'fillholespri wr 5' ^ fhp_th1_fh.png call StopWatch 0 01:05:49 Visually quite good. Horribly slow.
Limit the search radius, and repeat. %IM7DEV%magick ^ toes_holed.png ^ -process 'fillholespri wr 5 lsr 10%%' ^ -process 'fillholespri wr 5 lsr 10%%' ^ -process 'fillholespri wr 5 lsr 10%%' ^ -process 'fillholespri wr 5 lsr 10%%' ^ -process 'fillholespri wr 5 lsr 10%%' ^ -process 'fillholespri wr 5 lsr 10%%' ^ -process 'fillholespri wr 5 lsr 10%%' ^ fhp_th2_fh.png call StopWatch 0 00:01:30 Quicker but looks horrible, with many repeated patterns.
Search entire image, window mode. %IM7DEV%magick ^ toes_holed.png ^ -process 'fillholespri wr 5 copy window' ^ fhp_th3_fh.png call StopWatch 0 00:01:48 Even quicker and looks better.
Search entire image, window mode, with reduced copy radius. %IM7DEV%magick ^ toes_holed.png ^ -process 'fillholespri wr 5 cr 3 copy window' ^ fhp_th3a_fh.png call StopWatch 0 00:03:54 Even quicker and looks better.
Limit the search radius, window mode, and repeat. %IM7DEV%magick ^ toes_holed.png ^ -process 'fillholespri wr 5 lsr 10%% copy window' ^ -process 'fillholespri wr 5 lsr 10%% copy window' ^ -process 'fillholespri wr 5 lsr 10%% copy window' ^ -process 'fillholespri wr 5 lsr 10%% copy window' ^ -process 'fillholespri wr 5 lsr 10%% copy window' ^ -process 'fillholespri wr 5 lsr 10%% copy window' ^ -process 'fillholespri wr 5 lsr 10%% copy window' ^ fhp_th4_fh.png call StopWatch 0 00:00:04 Very quick but looks horrible.
Search entire image with a larger window, window mode. %IM7DEV%magick ^ toes_holed.png ^ -process 'fillholespri wr 10 copy window' ^ fhp_th5_fh.png call StopWatch 0 00:01:24 The long grass blade has been replicated five times.
As previous, with homogeneity check off. %IM7DEV%magick ^ toes_holed.png ^ -process 'fillholes wr 10 copy window hc off' ^ fhp_th6_fh.png 0 00:02:28
As previous, but search random. The seed makes it repeatable. %IM7DEV%magick ^ -seed 1234 ^ toes_holed.png ^ -process 'fillholespri wr 10 search random copy window' ^ fhp_thr1_fh.png 0 00:00:02 Very fast, but blocky result.
As previous, but reduce the copy radius. %IM7DEV%magick ^ -seed 1234 ^ toes_holed.png ^ -process 'fillholespri wr 10 search random copy window cr 1' ^ fhp_thr2_fh.png 0 00:01:00 Yuck.
Instead, a smaller search window. %IM7DEV%magick ^ -seed 1234 ^ toes_holed.png ^ -process 'fillholespri wr 5 search random copy window' ^ fhp_thr3_fh.png 0 00:00:01 Better but not great.
As previous but with ten times as many random searches. %IM7DEV%magick ^ -seed 1234 ^ toes_holed.png ^ -process 'fillholespri wr 5 search random rand_searches 1000%% copy window' ^ fhp_thr4_fh.png 0 00:00:05 Acceptable.
As previous, but search skip. %IM7DEV%magick ^ toes_holed.png ^ -process 'fillholespri wr 10 search skip copy window' ^ fhp_ths1_fh.png 0 00:00:03 The result is credible, aside from being over-sharp.
As previous, but repeated, with dissimilarity threshold. %IM7DEV%magick ^ toes_holed.png ^ -process 'fillholes wr 10 search skip dt 0.05 copy window' ^ -process 'fillholes wr 5 search skip dt 0.05 copy window' ^ -process 'fillholes wr 2 search skip copy window' ^ fhp_ths2_fh.png 0 00:00:31 Slower and not so good.

The process has introduced an artificial sharpness within the filled area. We can blur a result, placing it behind the original image so only the filled area is blurred.

%IMG7%magick ^
  fhp_th3_fh.png -blur 0x1 ^
  toes_holed.png ^
  -compose SrcOver -composite ^
  fhp_th3_fhp_bl.png

cmd /c exit /B 0

It is now very hard to identify which pixels have been filled, even at increased magnification. Of course, blurring does not remove the repeated pattern where fillholes has become trapped.

Use which module?

Of the two process modules, fillholes and fillholespri, neither is better than the other, in an absolute sense.

fillholespri is up to 25% slower. However, it currently contains debugging code, so will improve.
fillholespri generally gives better results.
fillholespri takes more memory, 48 bytes per image pixel.
fillholespri is far more complex, increasing the chance of bugs and the cost of maintenance.

Cleanup

We don't need to keep the frame files, so delete them.

del fhp_x_fr_*.png
del fhp_x_fr2_*.png

Possible improvements

See my main Filling holes page.

Acknowledgments

My priority algorithm is strongly based on:

Object Removal by Exemplar-Based Inpainting Criminisi, Pérez and Toyama, 2004

That paper doesn't give full implementation details, and it is likely that my implementation differs from theirs.

For other acknowledgements and references, see my main Filling holes page.

Source code

For convenience, this .c code and .bat are also available in a zip file. See Zipped BAT files.

fillholespri.c

/*
    Reference: http://im.snibgo.com/fillholes.htm

    Created 21-Nov-2015

    Updated:
      17-July-2016 AutoRepeat.
      3-April-2018 for v7.0.7-28
*/


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <ctype.h>
#include "vsn_defines.h"

#include "fillholescommon.inc"

#define DEBUG 0
#define VERIFY 0

#define CONF_FACT 0.9

typedef double fhpValueT;

typedef enum {
  psOpaq,    // Input pixel was opaque, or has become opaque.
  psEdge,    // Currently pixel is transparent but has an opaque neighbour.
  psInner,   // Currently pixel is transparent and has only transparent neghbours.
  psCantFill // Pixel is transparent and can't be filled, so remains transparent.
} pixStateT;

// Valid state transitions:
//   Inner -> Edge -> Opaq
//   Inner -> Opaq
//   Edge -> CantFill

typedef struct {
  pixStateT
    pixState;

  fhpValueT
    confidence,
    energy,
    priority;
} fhPriPixT;

typedef struct {
  ssize_t
    x,
    y;
} pqNodeT;

typedef struct {
  ssize_t
    width,
    height;

  fhPriPixT 
    **fhPriPix;

  ssize_t
    pqNumNodes;

  ssize_t
    pqMaxNodes;

  pqNodeT
    *pqNodes;
} fhPriT;


//---------------------------------------------------------------------------
// Priority queue functions.
// This is a MAX queue: higher values will come out first.

#define LCHILD(x) 2 * x + 1
#define RCHILD(x) 2 * x + 2
#define PARENT(x) (x - 1) / 2

static fhpValueT inline pqDataOfNode (fhPriT * fhp, ssize_t NodeNum)
{
  ssize_t x = fhp->pqNodes[NodeNum].x;
  ssize_t y = fhp->pqNodes[NodeNum].y;
  return fhp->fhPriPix[y][x].priority;
}

static void pqVerify (fhPriT * fhp, MagickBooleanType verbose)
{
  ssize_t
    i;

  if (verbose) printf ("pqVerify pqNumNodes=%li\n", fhp->pqNumNodes);

  for (i=0; i < fhp->pqNumNodes; i++) {
    ssize_t lch = LCHILD(i);
    ssize_t rch = RCHILD(i);

    MagickBooleanType HasL = lch < fhp->pqNumNodes;
    MagickBooleanType HasR = rch < fhp->pqNumNodes;

    if (HasL && pqDataOfNode (fhp, lch) > pqDataOfNode (fhp, i))
      printf ("** Bad %li left\n", i);

    if (HasR && pqDataOfNode (fhp, rch) > pqDataOfNode (fhp, i))
      printf ("** Bad %li right\n", i);
  }
}

#if 0
static void pqDump (FillHoleT * fh, fhPriT * fhp)
{
  ssize_t
    i;

  printf ("pqNumNodes=%li\n", fhp->pqNumNodes);
  printf ("  # x,y, state, confidence, energy, priority\n");

  for (i=0; i < fhp->pqNumNodes; i++) {
    pqNodeT * pqn = &fhp->pqNodes[i];
    fhPriPixT * pn = &(fhp->fhPriPix[pqn->y][pqn->x]);
    char cState;
    if      (pn->pixState == psOpaq) cState = 'o';
    else if (pn->pixState == psEdge) cState = 'e';
    else cState = 'i';

    printf ("  %li  %li,%li %c %g %g %g\n",
      i, pqn->x, pqn->y, cState, pn->confidence, pn->energy, pn->priority);
  }

  pqVerify (fhp, fh->debug);
}
#endif

static void inline pqSwapNodes (fhPriT * fhp, ssize_t n1, ssize_t n2)
{
  ssize_t t = fhp->pqNodes[n1].x;
  fhp->pqNodes[n1].x = fhp->pqNodes[n2].x;
  fhp->pqNodes[n2].x = t;

  t = fhp->pqNodes[n1].y;
  fhp->pqNodes[n1].y = fhp->pqNodes[n2].y;
  fhp->pqNodes[n2].y = t;
}

static void pqBalanceNode (fhPriT * fhp, ssize_t NodeNum)
// "Bubble-down".
// Finds the smallest of this node and its children.
// If the node isn't the smallest of the three,
//   swaps data with that child and recursively balances that child.
{
  ssize_t largest = NodeNum;
  ssize_t lch = LCHILD(NodeNum);
  ssize_t rch = RCHILD(NodeNum);

  fhpValueT SmData = pqDataOfNode (fhp, NodeNum);

  if (lch < fhp->pqNumNodes) {
    fhpValueT ldist = pqDataOfNode (fhp, lch);
    if (ldist > SmData) {
      SmData = ldist;
      largest = lch;
    }
  }

  if (rch < fhp->pqNumNodes) {
    fhpValueT rdist = pqDataOfNode (fhp, rch);
    if (rdist > SmData) {
      SmData = rdist;
      largest = rch;
    }
  }

  if (largest != NodeNum) {
    pqSwapNodes (fhp, largest, NodeNum);
    pqBalanceNode (fhp, largest);
  }
}

static void pqInsertNew (fhPriT * fhp, ssize_t x, ssize_t y)
{
  //printf ("pqIns %li,%li ", x, y);

  if (fhp->pqNumNodes >= fhp->pqMaxNodes) {
    printf ("pq bust\n");
  }

  fhpValueT newPri = fhp->fhPriPix[y][x].priority;
  // Priority could be zero, eg for graphics files with flat colours.

  ssize_t n = fhp->pqNumNodes;

  // "Bubble up".
  // If this data is more than node's parent,
  //   drop parent data into this node
  //   and consider putting new data into that parent.
  while (n && newPri > pqDataOfNode (fhp, PARENT(n))) {
    fhp->pqNodes[n] = fhp->pqNodes[PARENT(n)];
    n = PARENT(n);
  }
  fhp->pqNodes[n].x = x;
  fhp->pqNodes[n].y = y;

  fhp->pqNumNodes++;

  //pqVerify (fhp, MagickTrue);
}

static ssize_t pqFindNode (
  fhPriT * fhp, ssize_t x, ssize_t y, fhpValueT val, ssize_t NodeStart)
// Returns -1 if not present.
// Note: recursive.
{
  if (fhp->pqNumNodes==0) return -1;

  pqNodeT * pqn = &fhp->pqNodes[NodeStart];

  if (pqn->x == x && pqn->y == y) return NodeStart;

  if (pqDataOfNode (fhp, NodeStart) < val) return -1;

  ssize_t lch = LCHILD(NodeStart);
  if (lch < fhp->pqNumNodes) {
    ssize_t nch = pqFindNode (fhp, x, y, val, lch);
    if (nch >= 0) return nch;
  }

  ssize_t rch = RCHILD(NodeStart);
  if (rch < fhp->pqNumNodes) {
    ssize_t nch = pqFindNode (fhp, x, y, val, rch);
    if (nch >= 0) return nch;
  }

  return -1;
}

static ssize_t pqDelNode (
  fhPriT * fhp, ssize_t x, ssize_t y)
// Returns -1 if not present.
{
  fhpValueT oldPri = fhp->fhPriPix[y][x].priority;

  ssize_t n = pqFindNode (fhp, x, y, oldPri, 0);

  if (n < 0) {
    printf ("** Can't find to del %li,%li %g\n", x, y, oldPri);
    return n;
  }

  // Move last to here, delete last, rebalance.

  pqNodeT * pqn = &fhp->pqNodes[n];
  pqNodeT * pqLast = &fhp->pqNodes[fhp->pqNumNodes-1];
  pqn->x = pqLast->x;
  pqn->y = pqLast->y;
  fhp->pqNumNodes--;

  fhpValueT newPri = fhp->fhPriPix[pqn->y][pqn->x].priority;

  // "Bubble up".
  // If this data is more than node's parent,
  //   swap with parent and iterate.
  ssize_t bn = n;
  while (bn && newPri > pqDataOfNode (fhp, PARENT(bn))) {
    pqSwapNodes (fhp, bn, PARENT(bn));
    bn = PARENT(bn);
  }

  // Bubble down from here.
  pqBalanceNode (fhp, n);

  // Bubble down from the top.
  // FIXME: do we need this?
  pqBalanceNode (fhp, 0);

#if VERIFY==1
  pqVerify (fhp, MagickFalse);
#endif

  return n;
}

/*==
static void pqFindMax (fhPriT * fhp, ssize_t *x, ssize_t *y)
{
  if (!fhp->pqNumNodes) printf ("** Bad: fm: empty");

  pqNodeT * pqn = &fhp->pqNodes[0];
  *x = pqn->x;
  *y = pqn->y;
}
==*/

static int pqRemoveMax (fhPriT * fhp, ssize_t *x, ssize_t *y)
// Returns 1 if okay, or 0 if no data.
{
  if (!fhp->pqNumNodes) {
    //printf ("** Bad: rm: empty");
    return 0;
  }

  pqNodeT * pqn = &fhp->pqNodes[0];
  *x = pqn->x;
  *y = pqn->y;

  // Put last into root, and rebalance.

  fhp->pqNodes[0] = fhp->pqNodes[--fhp->pqNumNodes];

  pqBalanceNode (fhp, 0);

  return 1;
}


static void UpdIfData (
  FillHoleT * fh,
  fhPriT * fhp, ssize_t x, ssize_t y, fhpValueT oldPri, fhpValueT newPri)
// If the pixel is in the queue, updates the priority.
{
  //assert(newPri >= oldPri);

  printf ("UpdIfData %li,%li: ", x, y);

  ssize_t n = pqFindNode (fhp, x, y, oldPri, 0);

  if (n < 0) {
    // Not in queue.
    return;
  } else {
    printf ("%li\n", n);
  }

  fhp->fhPriPix[y][x].priority = newPri;

  // "Bubble up".
  // If this data is more than node's parent,
  //   swap with parent and iterate.
  while (n && newPri > pqDataOfNode (fhp, PARENT(n))) {
    pqSwapNodes (fhp, n, PARENT(n));
    n = PARENT(n);
  }

  // Bubble down from the top.
  // FIXME: do we need this?
  pqBalanceNode (fhp, 0);

#if VERIFY==1
  pqVerify (fhp, fh->do_verbose);

  if (fh->debug==MagickTrue) {
    ssize_t n = pqFindNode (fhp, x, y, newPri, 0);
    printf ("  UpdIfData: n=%li\n", n);
    assert (n >= 0);
  }
#endif
}

static void UpdInsData (
  FillHoleT * fh,
  fhPriT * fhp, ssize_t x, ssize_t y, fhpValueT oldPri, fhpValueT newPri)
// Updates the priority of a pixel, or inserts if not already present.
{
  ssize_t n = pqFindNode (fhp, x, y, oldPri, 0);

  if (fh->debug==MagickTrue) printf ("UpdInsData %li,%li: %li\n", x, y, n);

  fhp->fhPriPix[y][x].priority = newPri;

  if (n < 0) {
    pqInsertNew (fhp, x, y);
  } else {
    if (fh->debug==MagickTrue) printf ("UpdInsData %g -> %g\n", oldPri, newPri);

    //assert(newPri >= oldPri);

    if (newPri >= oldPri) {
      // "Bubble up".
      // If this data is more than node's parent,
      //   swap with parent and iterate.
      while (n && newPri > pqDataOfNode (fhp, PARENT(n))) {
        pqSwapNodes (fhp, n, PARENT(n));
        n = PARENT(n);
      }
    } else {
      pqBalanceNode (fhp, n);
    }

    // Bubble down from the top.
    // FIXME: do we need this?
    pqBalanceNode (fhp, 0);

#if VERIFY==1
    pqVerify (fhp, fh->debug);
#endif
  }

  if (fh->debug==MagickTrue) {
    ssize_t n = pqFindNode (fhp, x, y, newPri, 0);
    printf ("  UpdInsData: n=%li\n", n);
    assert (n >= 0);
  }
}


//---------------------------------------------------------------------------

/* Energy of a pixel is defined as the maximum gradient between that pixel
   and an 8-connected neighbour.

   The gradient is delta(channel)/distance,
   where distance is 1 for adjacent NSEW or sqrt(2) for NW, NE, SW, SE.

   Knowing dRed, dGreen and dBlue, do we use the max of these, or RMS or what? Max, I think.

   Neighbours that are fully transparent are ignored.

   Neighbours that are partially transparent???


   NO!!! We need the energy of a transparent (unknown) pixel, as defined by its neighbours.
   Bugger.

   Maybe we could take "this pixel" to be the average of is neighbours, then calc as above.

   Or, we say above definition is for opaque pixels.
   For edge pixels, energy is the maximum energy of its opaque neighbours.
*/

/*
  We initially put only the edge pixels in the queue.
  The others have zero priority, so no need to put them in until we know the priority.
*/

#if 0
static void DumpPriPix (
  fhPriT * fhp
)
{
  ssize_t x, y;

  printf ("x,y, state, confidence, energy, priority\n");
  for (y = 0; y < fhp->height; y++) {
    for (x = 0; x < fhp->width; x++) {
      fhPriPixT * pn = &(fhp->fhPriPix[y][x]);
      printf ("%li,%li %i %g %g %g\n",
       x, y, pn->pixState, pn->confidence, pn->energy, pn->priority);
    }
  }
}
#endif


#define CALC_CE \
  if (pno->pixState == psOpaq) { \
    if (pno->confidence > 0 && confid > pno->confidence * CONF_FACT) confid = pno->confidence * CONF_FACT; \
    if (energy < pno->energy) energy = pno->energy; \
    num_opaq++; \
  }

static fhpValueT inline CalcPixConfEn (
  FillHoleT * fh,
  fhPriT * fhp,
  ssize_t x,
  ssize_t y
// Calculate pixel confidence and energy.
// Returns calculated priority (but doesn't set it).
)
{
  return 0;
}


static void inline CalcOnePixEdgeData (
  FillHoleT * fh,
  fhPriT * fhp,
  ssize_t x,
  ssize_t y
)
// Confidence of an edge pixel is the lowest confidence of 8-connected opaque neighbours,
//   multiplied by a factor.
// Energy of an edge pixel is the highest energy of 8-connected opaque neighbours.
// Priority is confidence * energy * num_opaq/8.
{
  if (fh->debug==MagickTrue) printf ("      coped %li,%li", x, y);

  fhPriPixT * pno;

  fhPriPixT * pn = &(fhp->fhPriPix[y][x]);

  fhpValueT oldPri = pn->priority;  // In case we need to update queue.

  pn->priority = 0;


  fhpValueT confid = 99;
  fhpValueT energy = 0;
  int num_opaq = 0;

  if (y > 0) {
    if (x > 0) {
      pno = &(fhp->fhPriPix[y-1][x-1]);
      CALC_CE;
    }
    fhPriPixT * pno = &(fhp->fhPriPix[y-1][x]);
    CALC_CE;
    if (x < fhp->width-1) {
      pno = &(fhp->fhPriPix[y-1][x+1]);
      CALC_CE;
    }
  }

  if (x > 0) {
    pno = &(fhp->fhPriPix[y][x-1]);
    CALC_CE;
  }
  if (x < fhp->width-1) {
    pno = &(fhp->fhPriPix[y][x+1]);
    CALC_CE;
  }

  if (y < fhp->height-1) {
    if (x > 0) {
      pno = &(fhp->fhPriPix[y+1][x-1]);
      CALC_CE;
    }
    pno = &(fhp->fhPriPix[y+1][x]);
    CALC_CE;
    if (x < fhp->width-1) {
      pno = &(fhp->fhPriPix[y+1][x+1]);
      CALC_CE;
    }
  }

  assert (num_opaq > 0);
  assert (confid > 0);

  pn->confidence = confid;
  pn->energy = energy;

  if (pn->energy==0) pn->energy=0.00001; // FIXME?

  // If this is edge pixel, we need to update queue.

  fhpValueT newPri = pn->confidence * pn->energy * num_opaq/8.0;

  if (fh->debug==MagickTrue) printf (" co=%g en=%g no=%i  ", pn->confidence, pn->energy, num_opaq);

  if (fh->debug==MagickTrue) printf (" coped %g->%g \n", oldPri, newPri);

  if (pn->pixState == psEdge) {
    if (fh->debug==MagickTrue) printf (" edge\n");
    UpdInsData (fh, fhp, x, y, oldPri, newPri);
  } else {
    if (fh->debug==MagickTrue) printf (" notEdge\n");
    //pn->priority = newPri;
    // Pixel may be in queue, even though no longer an edge.
    assert (1==0);
    if (oldPri != newPri) UpdIfData (fh, fhp, x, y, oldPri, newPri);
  }
}


static void CalcEdgeData (
  FillHoleT * fh,
  fhPriT * fhp
)
// Calculates confidence, energy and priority of all edge pixels,
// and insert them into queue.
{
  ssize_t
    x, y;

  for (y = 0; y < fhp->height; y++) {
    for (x = 0; x < fhp->width; x++) {
      fhPriPixT * pn = &(fhp->fhPriPix[y][x]);
      if (pn->pixState == psEdge) {
        if (fh->debug==MagickTrue) printf ("ced %li,%li ", x, y);
        CalcOnePixEdgeData (fh, fhp, x, y);
      }
    }
  }

#if VERIFY==1
  pqVerify (fhp, fh->debug);
#endif
}


static fhpValueT CalcEnergy (
  const Image * image,
  const VIEW_PIX_PTR * p1,
  const VIEW_PIX_PTR * p2,
  fhpValueT dist)
{
  fhpValueT
    d,
    dMax;

  dMax = fabs(GET_PIXEL_RED(image,p1) - GET_PIXEL_RED(image,p2));
  d = fabs(GET_PIXEL_GREEN(image,p1) - GET_PIXEL_GREEN(image,p2));
  if (dMax < d) dMax = d;
  d = fabs(GET_PIXEL_BLUE(image,p1) - GET_PIXEL_BLUE(image,p2));
  if (dMax < d) dMax = d;

  dMax = dMax / (QuantumRange * dist) * (GET_PIXEL_ALPHA(image,p2) / QuantumRange);

  return dMax;
}

/*===
static void CalcEdgePixelData (
  FillHoleT * fh,
  fhPriT * fhp)
// Edge pixels are exactly those in the priority queue.
{
  int i;

  for (i=0; i < fhp->pqNumNodes; i++) {
    pqNodeT * pqn = &fhp->pqNodes[i];
    CalcOnePixEdgeData (fh, fhp, pqn->x, pqn->y);
  }
}
===*/


static MagickBooleanType Initialise (const Image *image,
  FillHoleT * fh,
  fhPriT * fhp,
  ExceptionInfo *exception
)
{
  ssize_t
    x, y, xMult;

  MagickBooleanType
    status = MagickTrue;

  CacheView
    *image_view;

  fhp->width = image->columns;
  fhp->height = image->rows;

  // Allocate memory

  if (fh->do_verbose) printf ("Alloc %lix%li\n", fhp->width, fhp->height);

  fhp->pqMaxNodes = fhp->height * fhp->width;
  fhp->pqNodes = (pqNodeT *) AcquireQuantumMemory(fhp->pqMaxNodes, sizeof(pqNodeT));
  if (fhp->pqNodes == (pqNodeT *) NULL) {
    return MagickFalse;
  }

  fhp->pqNumNodes = 0;

  fhp->fhPriPix = (fhPriPixT **) AcquireQuantumMemory(fhp->height, sizeof(*fhp->fhPriPix));
  if (fhp->fhPriPix == (fhPriPixT **) NULL) {
    RelinquishMagickMemory(fhp->pqNodes);
    return MagickFalse;
  }
  for (y = 0; y < fhp->height; y++) {
    fhp->fhPriPix[y] = (fhPriPixT *) AcquireQuantumMemory(fhp->width, sizeof(**fhp->fhPriPix));
    if (fhp->fhPriPix[y] == (fhPriPixT *) NULL) break;
  }
  if (y < fhp->height) {
    for (y--; y >= 0; y--) {
      if (fhp->fhPriPix[y] != (fhPriPixT *) NULL)
        fhp->fhPriPix[y] = (fhPriPixT *) RelinquishMagickMemory(fhp->fhPriPix[y]);
    }
    fhp->fhPriPix = (fhPriPixT **) RelinquishMagickMemory(fhp->fhPriPix);
    RelinquishMagickMemory(fhp->pqNodes);
    return MagickFalse;
  }

  xMult = Inc_ViewPixPtr (image);

  // Populate values

  InitAutoLimit (fh);

  if (fh->do_verbose) printf ("Populate\n");

  image_view = AcquireVirtualCacheView(image,exception);

#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp parallel for schedule(static,4) shared(status) \
      MAGICK_THREADS(image,image,image->rows,1)
#endif
  for (y = 0; y < fhp->height; y++) {
    VIEW_PIX_PTR const
      *p, *pxy, *pup, *pdn;

    if (status == MagickFalse)
      continue;

    // We use virtual pixels for energy calculation.

    p=GetCacheViewVirtualPixels(image_view,-1,y-1,image->columns+2,3,exception);
    if (p == (const VIEW_PIX_PTR *) NULL)
      status=MagickFalse;

    // 31 Jan 2015: FIXME: Eek. Following were wrong for v7?
    pup = p + xMult;
    pxy = p + xMult * (image->columns + 3);
    pdn = p + xMult * (2*image->columns + 5);

    for (x = 0; x < fhp->width; x++) {
      fhpValueT
        e;

      //if (fh->debug==MagickTrue) printf ("Pop %i,%i ", (int)x, (int)y);

      fhPriPixT * pn = &(fhp->fhPriPix[y][x]);
      pn->confidence = GET_PIXEL_ALPHA(image,pxy)/QuantumRange;

      pn->energy = CalcEnergy (image, pxy, pup-xMult, M_SQRT2);
      e = CalcEnergy (image, pxy, pup, 1);
      if (pn->energy < e) pn->energy = e;
      e = CalcEnergy (image, pxy, pup+xMult, M_SQRT2);
      if (pn->energy < e) pn->energy = e;

      e = CalcEnergy (image, pxy, pxy-xMult, 1);
      if (pn->energy < e) pn->energy = e;
      e = CalcEnergy (image, pxy, pxy+xMult, 1);
      if (pn->energy < e) pn->energy = e;

      e = CalcEnergy (image, pxy, pdn-xMult, M_SQRT2);
      if (pn->energy < e) pn->energy = e;
      e = CalcEnergy (image, pxy, pdn, 1);
      if (pn->energy < e) pn->energy = e;
      e = CalcEnergy (image, pxy, pdn+xMult, M_SQRT2);
      if (pn->energy < e) pn->energy = e;

      //if (fh->debug==MagickTrue) printf (" e=%g\n", pn->energy);

      if (pn->confidence <= 0) {
        // FIXME? If any 8-connected neighbour is opaque, this is an edge.
        if (GET_PIXEL_ALPHA(image,pup-xMult)>0
         || GET_PIXEL_ALPHA(image,pup)>0
         || GET_PIXEL_ALPHA(image,pup+xMult)>0
         || GET_PIXEL_ALPHA(image,pxy-xMult)>0
         || GET_PIXEL_ALPHA(image,pxy+xMult)>0
         || GET_PIXEL_ALPHA(image,pdn-xMult)>0
         || GET_PIXEL_ALPHA(image,pdn)>0
         || GET_PIXEL_ALPHA(image,pdn+xMult)>0)
        {
          pn->pixState = psEdge;
          fh->Match.nLimitCnt++;
        } else {
          pn->pixState = psInner;
        }

      } else {
        pn->pixState = psOpaq;
      }
      if (pn->confidence == 0 ) pn->energy = 0;
      pn->priority = 0; // Until we know better.

      // We can't find data for edge pixels until we know data for all edge neighbours.

      // FIXME? following were wrong for v7.
      p   += xMult;
      pup += xMult;
      pxy += xMult;
      pdn += xMult;
    }
  }

  if (fh->debug==MagickTrue) printf ("fh->Match.nLimitCnt=%li\n", fh->Match.nLimitCnt);

  image_view=DestroyCacheView(image_view);

  pqVerify (fhp, fh->debug);

  CalcEdgeData (fh, fhp);

  //if (fh->debug==MagickTrue) pqDump (fh, fhp);
  pqVerify (fhp, fh->debug);
  //if (fh->debug==MagickTrue) DumpPriPix (fhp);

  if (fh->do_verbose) printf ("End populate\n");

  return (status);
}


static void deInitialise(
  FillHoleT * fh,
  fhPriT * fhp
)
{
  ssize_t
    y;

  for (y = 0; y < fhp->height; y++) {
    fhp->fhPriPix[y] = (fhPriPixT *) RelinquishMagickMemory(fhp->fhPriPix[y]);
  }
  fhp->fhPriPix = (fhPriPixT **) RelinquishMagickMemory(fhp->fhPriPix);

  RelinquishMagickMemory(fhp->pqNodes);
}

static MagickBooleanType inline IsStateOpaque (
  fhPriT * fhp,
  ssize_t x,
  ssize_t y)
{
  if (x >= 0 && x < fhp->width &&
      y >= 0 && y < fhp->height)
  {
    return fhp->fhPriPix[y][x].pixState == psOpaq;
  }
  return MagickFalse;
}

static void inline MakeEdgePix (
  FillHoleT * fh,
  fhPriT * fhp,
  ssize_t x,
  ssize_t y)
// If not outside image, and inner, set to edge and calc priority and add to queue.
// If it's already an edge, re-calc data and update queue,
{
  if (fh->debug==MagickTrue) printf ("  mep %li,%li\n", x, y);

  if (x >= 0 && x < fhp->width &&
      y >= 0 && y < fhp->height)
  {
    fhPriPixT * pn = &(fhp->fhPriPix[y][x]);
    if (pn->pixState == psInner) {
      // FIXME: only if 4-connected neighbours are opaque?

      if (
          IsStateOpaque (fhp, x,   y-1) ||
          IsStateOpaque (fhp, x-1, y)   ||
          IsStateOpaque (fhp, x+1, y)   ||
          IsStateOpaque (fhp, x,   y+1)
       /*=== ||
          IsStateOpaque (fhp, x-1, y-1) ||
          IsStateOpaque (fhp, x+1, y-1) ||
          IsStateOpaque (fhp, x-1, y+1) ||
          IsStateOpaque (fhp, x+1, y+1) ===*/
         )
      {
        if (fh->debug==MagickTrue) printf ("    mepI %li,%li\n", x, y);

        pn->pixState = psEdge;

        CalcOnePixEdgeData (fh, fhp, x, y);

        // FIXME: again, batch entry probably quicker.
      } else {
        if (fh->debug==MagickTrue) printf ("    mepU %li,%li no_opaque_neighbours", x, y);
      }
    } else if (pn->pixState == psEdge) {
      if (fh->debug==MagickTrue) printf ("    mepU %li,%li already_edge", x, y);

      CalcOnePixEdgeData (fh, fhp, x, y);
    }
  }
}


static MagickBooleanType inline IsPixelOpaque (
  Image * image,
  CacheView * view,
  ssize_t x,
  ssize_t y,
  ExceptionInfo *exception)
{
  const VIEW_PIX_PTR
    *src;

  src = GetCacheViewVirtualPixels(view,x,y,1,1,exception);
  if (src == (const VIEW_PIX_PTR *) NULL) {
    printf ("bad src");
    return MagickFalse;
  }

  return (GET_PIXEL_ALPHA (image, src) > 0);
}


static fhpValueT inline PixelEnergy (
  Image * image,
  CacheView * view,
  ssize_t x,
  ssize_t y,
  ExceptionInfo *exception)
{
  const VIEW_PIX_PTR
    *src,
    *this,
    *other;

  int
    i;

  src = GetCacheViewVirtualPixels(view,x-1,y-1,3,3,exception);
  if (src == (const VIEW_PIX_PTR *) NULL) { printf ("bad src"); return MagickFalse; }

  // FIXME: for v7

  ssize_t xMult = Inc_ViewPixPtr (image);

  other = src;
  this = src + xMult*4;
  fhpValueT e = 0;

  for (i = 0; i < 9; i++) {
    if (other != this && GET_PIXEL_ALPHA (image, other) > 0) {
      fhpValueT v = CalcEnergy (image, this, other,
        (i==0 || i==2 || i==6 || i==8) ? M_SQRT2 : 1);
      if (e < v) e = v;
    }
    other += Inc_ViewPixPtr (image);
  }

  return (e);
}


static MagickBooleanType ProcessPixels (
  Image *image,
  Image *new_image,
  FillHoleT * fh,
  fhPriT * fhp,
  MagickBooleanType *ChangedSome,
  MagickBooleanType *unfilled,
  ExceptionInfo *exception)
{
  Image
    *holed_image;

  CacheView
    *copy_inp_view,
    *new_view;

  ssize_t
    x, y;

  int
    GotOne,
    frameNum = 0;

  MagickBooleanType
    changedAny,
    status=MagickTrue;

  CopyWhereT
    cw;

  if (fh->copyWhat == copyWindow) {
    cw.wi = fh->sqCopyDim;
    cw.ht = fh->sqCopyDim;
  } else {
    cw.wi = 1;
    cw.ht = 1;
  }

  // Clone the image, same size, copied pixels:
  //
  holed_image=CloneImage(image, 0, 0, MagickTrue, exception);
  if (holed_image == (Image *) NULL)
    return MagickFalse;

  if (SetNoPalette (holed_image, exception) == MagickFalse)
    return MagickFalse;

  fh->CompWind.ref_image = image;
  fh->CompWind.ref_view = AcquireVirtualCacheView (image, exception);

  new_view        = AcquireAuthenticCacheView (new_image, exception);
  fh->CompWind.sub_image = new_image;
  fh->CompWind.sub_view = new_view;
  copy_inp_view   = CloneCacheView (fh->CompWind.ref_view);

  *ChangedSome = MagickFalse;
  changedAny = MagickFalse;
  *unfilled = MagickFalse;

//#if DEBUG==1
  if (fh->debug==MagickTrue) printf ("ProcessPixels: start do\n");
//#endif

  do {
    GotOne = pqRemoveMax (fhp, &x, &y);

    if (fh->debug==MagickTrue) printf ("GotOne? %i\n", GotOne);

    if (GotOne)
    {
      assert (fhp->fhPriPix[y][x].pixState == psEdge);

      if (fh->debug==MagickTrue) printf ("pp %li,%li\n", x, y);

      changedAny = MagickFalse;

      fh->Match.nLimitCnt--;

      status = MatchAndCopy (
        fh,
        new_image,
        new_view,
        copy_inp_view,
        x,
        y,
        &cw,
        &frameNum,
        unfilled,
        &changedAny,
        exception);

      if (!changedAny) {
        if (fh->debug==MagickTrue) printf ("  Can't change %li,%li r=%i\n", x, y, (int)fh->CompWind.Reason);

        fhp->fhPriPix[y][x].pixState = psCantFill;
      } else {
        ssize_t i, j;

        if (fh->debug==MagickTrue)
          printf (" doneMaC st=%i  %lix%li+%li+%li\n", status, cw.wi, cw.ht, cw.dstX, cw.dstY);

        *ChangedSome = MagickTrue;

        // We may not have changed them all.
        MagickBooleanType OpaqAny = MagickFalse;
        if (fh->debug==MagickTrue) printf ("toOpaq ");

        assert (cw.dstX >= 0 && cw.dstX+cw.wi <= fh->Width);
        assert (cw.dstY >= 0 && cw.dstY+cw.ht <= fh->Height);

        fhpValueT conf = 99;

        for (j=0; j < cw.ht; j++) {
          ssize_t adjY = cw.dstY+j;
          for (i=0; i < cw.wi; i++) {
            ssize_t adjX = cw.dstX+i;
            fhPriPixT * pn = &(fhp->fhPriPix[adjY][adjX]);
            if (pn->pixState != psOpaq) {
              if (IsPixelOpaque (new_image, new_view, adjX, adjY, exception)) {
                if (fh->debug==MagickTrue) printf (" %li,%li\n", adjX, adjY);
                if (pn->pixState == psEdge && (adjX != x || adjY != y)) {
                  if (pqDelNode (fhp, adjX, adjY) < 0) status=MagickFalse;
                  fh->Match.nLimitCnt--;
                }
                pn->pixState = psOpaq; // FIXME: or "filled"?
                if (pn->confidence > 0 && conf > pn->confidence) conf = pn->confidence;
                // We will also update energy, at least of pixels on edge of copied block.
                // We can only do this after we know which neighbours are opaque.
                OpaqAny = MagickTrue;
              }
            }
          }
        }

        if (fh->debug==MagickTrue) printf ("done toOpaq\n");

        if (conf == 99) conf = 0.001;
        else conf *= CONF_FACT;

        if (OpaqAny) {
          // Calc energies of pixels on edge of copied block.
          // FIXME: only for cw.wi>1 or cw.ht>1?
          // Possible optimisation: If a pixel is surrounded by opaque pixels,
          // we will never care what energy it has.
          if (fh->debug==MagickTrue) printf ("calcE\n");
          for (i=0; i < cw.wi; i++) {
            // Top
            fhPriPixT * pn = &(fhp->fhPriPix[cw.dstY][cw.dstX+i]);
            pn->energy = PixelEnergy (new_image, new_view, cw.dstX+i, cw.dstY, exception);
            if (pn->confidence==0) pn->confidence = conf;
            if (cw.ht > 1) {
              // Bottom
              pn = &(fhp->fhPriPix[cw.dstY+cw.ht-1][cw.dstX+i]);
              pn->energy = PixelEnergy (new_image, new_view, cw.dstX+i, cw.dstY+cw.ht-1, exception);
              if (pn->confidence==0) pn->confidence = conf;
            }
          }
          for (j=1; j < cw.ht-1; j++) {
            // Left
            fhPriPixT * pn = &(fhp->fhPriPix[cw.dstY+j][cw.dstX]);
            pn->energy = PixelEnergy (new_image, new_view, cw.dstX, cw.dstY+j, exception);
            if (pn->confidence==0) pn->confidence = conf;
            if (cw.wi > 1) {
              // Right
              pn = &(fhp->fhPriPix[cw.dstY+j][cw.dstX+cw.wi-1]);
              pn->energy = PixelEnergy (new_image, new_view, cw.dstX+cw.wi-1, cw.dstY+j, exception);
              if (pn->confidence==0) pn->confidence = conf;
            }
          }

          // Walk around border outside the copied pixels.
          // If not outside image, and inner, set to edge and calc priority and add to queue.
          // But if we haven't changed them all, some of these may not really be edges.
          if (fh->debug==MagickTrue) printf ("\nWalk border from %li,%li\n", cw.dstX-1, cw.dstY-1);
          for (i=0; i < cw.wi+2; i++) {
            MakeEdgePix (fh, fhp, cw.dstX-1+i, cw.dstY-1);
            MakeEdgePix (fh, fhp, cw.dstX-1+i, cw.dstY+cw.ht);
          }
          for (j=0; j < cw.ht; j++) {
            MakeEdgePix (fh, fhp, cw.dstX-1, cw.dstY+j);
            MakeEdgePix (fh, fhp, cw.dstX+cw.wi, cw.dstY+j);
          }
          if (fh->debug==MagickTrue) printf ("End walk border\n");
          // FIXME: should we also update data on pixels that are already edges, updating queue?
          // I think so.
        }
      }
      if (fh->debug==MagickTrue) printf ("\n");
    }

    //if (fh->Match.SetAutoLimit) printf ("fh->Match.nLimitCnt=%li\n", fh->Match.nLimitCnt);

    if (fh->Match.SetAutoLimit && fh->Match.nLimitCnt <= 0) UseAutoLimit (fh);

  } while (GotOne==1 && status==MagickTrue);

#if DEBUG==1
  printf ("ProcessPixels: done do\n");
#endif

  if (status==MagickFalse) printf ("** Bug: ProcessPixels: status==MagickFalse\n");

  if (*unfilled) printf ("** Some pixels are not filled\n");

  copy_inp_view   = DestroyCacheView (copy_inp_view);
  new_view        = DestroyCacheView (new_view);

  fh->CompWind.ref_view = DestroyCacheView (fh->CompWind.ref_view);

//  srch_holed_view = DestroyCacheView (srch_holed_view);
//  inp_view        = DestroyCacheView (inp_view);

  if (fh->debug==MagickTrue) printf ("unfilled=%i\n", *unfilled);

  return status;
}


// The next function corresponds in style to functions in transform.c
// It takes one image, and returns an image with filled holes.
//
static Image *fillholespri (
  Image *image,
  FillHoleT * fh,
  ExceptionInfo *exception)
{
  Image
    *new_image;

  fhPriT
    fhPri;

  MagickBooleanType
    ChangedSome,
    unfilled;

  ResolveImageParams (image, fh);

  if (fh->debug==MagickTrue) printf ("sizeof(fhPriPixT)=%li sizeof(pqNodeT)=%li\n",
    sizeof(fhPriPixT), sizeof(pqNodeT));

  if (fh->debug==MagickTrue) printf ("Initialise ...");
  if (Initialise (image, fh, &fhPri, exception) == MagickFalse)
    return (Image *) NULL;

  // Clone the image, same size, copied pixels:
  //
  new_image=CloneImage(image, 0, 0, MagickTrue, exception);
  if (new_image == (Image *) NULL)
    return(new_image);

  if (fh->debug==MagickTrue) printf ("Process pixels ...");
  ProcessPixels (image, new_image, fh, &fhPri, &ChangedSome, &unfilled, exception);

  if (fh->debug==MagickTrue) printf ("Deinitialise ...\n");
  deInitialise (fh, &fhPri);

  while (unfilled && ChangedSome && fh->AutoRepeat) {
    if (fh->do_verbose) printf ("fillholespri: doing it again.\n");

    Image *tmp_copy = CloneImage(new_image, 0, 0, MagickTrue, exception);
    if (tmp_copy == (Image *) NULL) return(tmp_copy);

    if (Initialise (tmp_copy, fh, &fhPri, exception) == MagickFalse)
      return (Image *) NULL;

    ProcessPixels (tmp_copy, new_image, fh, &fhPri, &ChangedSome, &unfilled, exception);

    tmp_copy = DestroyImage (tmp_copy);

    if (fh->debug==MagickTrue) printf ("Deinitialise ...\n");
    deInitialise (fh, &fhPri);

  }

  return (new_image);
}

ModuleExport size_t fillholespriImage (
  Image **images,
  const int argc,
  const char **argv,
  ExceptionInfo *exception)
{
  Image
    *image,
    *new_image;

  MagickBooleanType
    status;

  FillHoleT
    fh;

  assert(images != (Image **) NULL);
  assert(*images != (Image *) NULL);
  assert((*images)->signature == MAGICK_CORE_SIG);

  status = menu (argc, argv, &fh);
  if (status == MagickFalse)
    return (-1);

  InitRand (&fh);

  for (image=(*images); image != (Image *) NULL; image=GetNextImageInList(image))
  {
    new_image = fillholespri (image, &fh, exception);

    ReplaceImageInList(&image,new_image);
    *images=GetFirstImageInList(image);
  }

  DeInitRand (&fh);

  return(MagickImageFilterSignature);
}

fillholescommon.inc

This code is common to both fillholes.c and fillholespri.c.

/*
    Reference: http://im.snibgo.com/fillholes.htm

    Data structures and functions common to fillholes process modules.

    Created 21-Nov-2015
*/

#define CH_R 1
#define CH_G 2
#define CH_B 4

#include "compwind.h"
#include "match.h"


typedef enum {
  copyOnePixel,
  copyWindow
} CopyWhatT;

typedef struct {
  CompWindT
    CompWind;

  MatchT
    Match;

  /* Set by user. */
  MagickBooleanType
    do_verbose,
    AutoRepeat,
    CopyRad_IsPc,
    WindowRad_IsPc,
    debug;

  double
    WindowRad,
    CopyRad,
    ImgDiagSq;

  int
    WindowRadPix,
    CopyRadPix,
    write_frames;

  char
    write_filename[MaxTextExtent];

  CopyWhatT
    copyWhat;

  /* Calculated. */
  int
    sqDim,
    sqCopyDim;

  ssize_t
    Width,
    Height;

  RandomInfo
    *restrict random_info;

} FillHoleT;

typedef struct {
  ssize_t
    srcX,
    srcY,
    dstX,
    dstY,
    wi,
    ht;
} CopyWhereT;

#include "compwind.inc"
#include "match.inc"


static void usage (void)
{
  printf ("Usage: -process 'fillhole [OPTION]...'\n");
  printf ("Populates transparent pixels, InFill.\n");
  printf ("\n");
  printf ("  wr,  window_radius N        radius of search window, >= 1\n");
  printf ("  lsr, limit_search_radius N  limit radius from transparent pixel to search\n");
  printf ("                                for source, >= 0\n");
  printf ("                                default = 0 = no limit\n");
  printf ("  als,  auto_limit_search X   automatic limit search, on or off\n");
  printf ("                                default on\n");
  printf ("  hc,  hom_chk X              homogeneity check, X is off or a small number\n");
  printf ("                                default 0.1\n");
  printf ("  e,   search_to_edges X      search for matches to image edges, on or off\n");
  printf ("                                default on\n");
  printf ("  s,   search X               X=entire or random or skip\n");
  printf ("                                default entire\n");
  printf ("  rs,  rand_searches N        number of random searches (eg 100)\n");
  printf ("                                default 0\n");
  printf ("  sn,  skip_num N             number of searches to skip in each direction\n");
  printf ("                                (eg 10)\n");
  printf ("                                default 0\n");
  printf ("  ref, refine X               whether to refine random and skip searches,\n");
  printf ("                                on or off\n");
  printf ("                                default on\n");
  printf ("  st,  similarity_threshold N\n");
  printf ("                              stop searching when RMSE <= N (eg 0.01)\n");
  printf ("                                default 0\n");
  printf ("  dt,  dissimilarity_threshold N\n");
  printf ("                              don't copy if best RMSE >= N (eg 0.05)\n");
  printf ("                                default: no threshold\n");
  printf ("  cp,  copy X                 X=onepixel or window\n");
  printf ("  cr,  copy_radius N          radius of pixels to copy, >= 1, <= wr\n");
  printf ("                                default: 0 or wr\n");
  printf ("  a,   auto_repeat            if pixels changed but any unfilled, repeat\n");
  printf ("  w,   write filename         write frames to files\n");
  printf ("  w2,  write2 filename        write frames to files\n");
  printf ("  v,   verbose                 write text information to stdout\n");
  printf ("\n"); 
}


static MagickBooleanType IsArg (char * pa, char * ShtOpt, char * LongOpt)
{
  if ((LocaleCompare(pa, ShtOpt)==0) || (LocaleCompare(pa, LongOpt)==0))
    return MagickTrue;

  return MagickFalse;
}


static inline MagickBooleanType EndsPc (const char *s)
{
  char c = *(s+strlen(s)-1);
  if (c == '%' || c == 'c')
    return MagickTrue;
  return MagickFalse;
}


static MagickBooleanType menu(
  const int argc,
  const char **argv,
  FillHoleT * pfh
)
// Returns MagickTrue if okay.
{
  int
    i;

  MagickBooleanType
    status;

  status = MagickTrue;

  pfh->do_verbose = MagickFalse;
  pfh->WindowRad = 1.0;
  pfh->CopyRad = 0.0;
  pfh->WindowRad_IsPc = pfh->CopyRad_IsPc = MagickFalse;
  pfh->WindowRadPix = 1;
  pfh->CopyRadPix = 0;
  pfh->write_frames = 0;
  pfh->write_filename[0] = '\0';
  pfh->copyWhat = copyOnePixel;
  pfh->debug = MagickFalse;
  pfh->AutoRepeat = MagickFalse;

  pfh->Match.LimSrchRad = 0.0;
  SetDefaultMatch (&pfh->Match);

  pfh->CompWind.HomChkOn = MagickTrue;
  pfh->CompWind.HomChk = 0.1;
  pfh->CompWind.nCompares = 0;


  for (i=0; i < argc; i++) {
    char * pa = (char *)argv[i];

    if (IsArg (pa, "wr", "window_radius")==MagickTrue) {
      i++;
      if (!isdigit ((int)*argv[i]))
      {
        fprintf (stderr, "Non-numeric argument to 'window_radius'.\n");
        status = MagickFalse;
      }
      else {
        pfh->WindowRad = atof(argv[i]);
        pfh->WindowRad_IsPc = EndsPc (argv[i]);
        if (!pfh->WindowRad_IsPc)
          pfh->WindowRadPix = pfh->WindowRad + 0.5;
      }
      if (pfh->WindowRad <= 0) {
        fprintf (stderr, "Bad 'window_radius' value.\n");
        status = MagickFalse;
      }
    } else if (IsArg (pa, "cr", "copy_radius")==MagickTrue) {
      i++;
      if (!isdigit ((int)*argv[i]))
      {
        fprintf (stderr, "Non-numeric argument to 'copy_radius'.\n");
        status = MagickFalse;
      }
      else {
        pfh->CopyRad = atof(argv[i]);
        pfh->CopyRad_IsPc = EndsPc (argv[i]);
        if (!pfh->CopyRad_IsPc)
          pfh->CopyRadPix = pfh->CopyRad + 0.5;
      }
      if (pfh->CopyRad <= 0) {
        fprintf (stderr, "Bad 'copy_radius' value.\n");
        status = MagickFalse;
      }
    } else if (IsArg (pa, "lsr", "limit_search_radius")==MagickTrue) {
      // FIXME: also allow percentage of smaller image dimension.
      i++;
      if (!isdigit ((int)*argv[i])) {
        fprintf (stderr, "Non-numeric argument to 'limit_search_radius'.\n");
        status = MagickFalse;
      } else {
        pfh->Match.LimSrchRad = atof(argv[i]);
        pfh->Match.LimSrchRad_IsPc = EndsPc (argv[i]);
        if (!pfh->Match.LimSrchRad_IsPc)
          pfh->Match.limSrchRadX = pfh->Match.limSrchRadY = pfh->Match.LimSrchRad + 0.5;
      }
      if (pfh->Match.LimSrchRad < 0) {
        fprintf (stderr, "Bad 'limit_search_radius' value.\n");
        status = MagickFalse;
      }

/*---
    } else if (IsArg (pa, "c", "channel")==MagickTrue) {
      i++;
      pfh->channels = 0;
      const char * p = argv[i];
      while (*p) {
        switch (toupper ((int)*p)) {
          case 'R':
            pfh->channels |= CH_R;
            break;
          case 'G':
            pfh->channels |= CH_G;
            break;
          case 'B':
            pfh->channels |= CH_B;
            break;
          case 'L':
            pfh->channels |= CH_R;
            break;
          case 'A':
            pfh->channels |= CH_G;
            break;
          default:
            fprintf (stderr, "Invalid 'channels' [%s]\n", argv[i]);
            status = MagickFalse;
        }
        p++;
      }
---*/
    } else if (IsArg (pa, "als", "auto_limit_search")==MagickTrue) {
      i++;
      if (LocaleCompare(argv[i], "on")==0)
        pfh->Match.AutoLs = MagickTrue;
      else if (LocaleCompare(argv[i], "off")==0)
        pfh->Match.AutoLs = MagickFalse;
      else {
        fprintf (stderr, "Invalid 'auto_limit_search' [%s]\n", argv[i]);
        status = MagickFalse;
      }
    } else if (IsArg (pa, "ref", "refine")==MagickTrue) {
      i++;
      if (LocaleCompare(argv[i], "on")==0)
        pfh->Match.Refine = MagickTrue;
      else if (LocaleCompare(argv[i], "off")==0)
        pfh->Match.Refine = MagickFalse;
      else {
        fprintf (stderr, "Invalid 'refine' [%s]\n", argv[i]);
        status = MagickFalse;
      }
    } else if (IsArg (pa, "s", "search")==MagickTrue) {
      i++;
      if (LocaleCompare(argv[i], "entire")==0)
        pfh->Match.searchWhat = swEntire;
      else if (LocaleCompare(argv[i], "random")==0)
        pfh->Match.searchWhat = swRandom;
      else if (LocaleCompare(argv[i], "skip")==0)
        pfh->Match.searchWhat = swSkip;
      else {
        fprintf (stderr, "Invalid 'search' [%s]\n", argv[i]);
        status = MagickFalse;
      }
    } else if (IsArg (pa, "rs", "rand_searches")==MagickTrue) {
      i++;
      if (!isdigit ((int)*argv[i])) {
        fprintf (stderr, "Non-numeric argument to 'rand_searches'.\n");
        status = MagickFalse;
      } else {
        pfh->Match.RandSearchesF = atof (argv[i]);
        pfh->Match.RandSearches_IsPc = EndsPc (argv[i]);
        if (!pfh->Match.RandSearches_IsPc)
          pfh->Match.RandSearchesI = pfh->Match.RandSearchesF + 0.5;
      }
      if (pfh->Match.RandSearchesF <= 0) {
        fprintf (stderr, "Bad 'rand_searches' value.\n");
        status = MagickFalse;
      }
    } else if (IsArg (pa, "sn", "skip_num")==MagickTrue) {
      i++;
      if (!isdigit ((int)*argv[i])) {
        fprintf (stderr, "Non-numeric argument to 'skip_num'.\n");
        status = MagickFalse;
      } else {
        pfh->Match.SkipNumF = atof (argv[i]);
        pfh->Match.SkipNum_IsPc = EndsPc (argv[i]);
        if (!pfh->Match.SkipNum_IsPc)
          pfh->Match.SkipNumI = pfh->Match.SkipNumF + 0.5;
      }
      if (pfh->Match.SkipNumF <= 0) {
        fprintf (stderr, "Bad 'skip_num' value.\n");
        status = MagickFalse;
      }
    } else if (IsArg (pa, "hc", "hom_chk")==MagickTrue) {
      i++;
      if (LocaleCompare(argv[i], "off")==0) {
        pfh->CompWind.HomChkOn = MagickFalse;
      } else {
        pfh->CompWind.HomChkOn = MagickTrue;
        if (!isdigit ((int)*argv[i])) {
          fprintf (stderr, "'hom_chk' argument must be number or 'off'.\n");
          status = MagickFalse;
        } else {
          pfh->CompWind.HomChk = atof (argv[i]);
        }
      }
    } else if (IsArg (pa, "e", "search_to_edges")==MagickTrue) {
      i++;
      if (LocaleCompare(argv[i], "on")==0)
        pfh->Match.SearchToEdges = MagickTrue;
      else if (LocaleCompare(argv[i], "off")==0)
        pfh->Match.SearchToEdges = MagickFalse;
      else {
        fprintf (stderr, "Invalid 'search_to_edges' [%s]\n", argv[i]);
        status = MagickFalse;
      }
    } else if (IsArg (pa, "a", "auto_repeat")==MagickTrue) {
      pfh->AutoRepeat = MagickTrue;
    } else if (IsArg (pa, "cp", "copy")==MagickTrue) {
      i++;
      if (LocaleCompare(argv[i], "onepixel")==0)
        pfh->copyWhat = copyOnePixel;
      else if (LocaleCompare(argv[i], "window")==0)
        pfh->copyWhat = copyWindow;
      else {
        fprintf (stderr, "Invalid 'copy' [%s]\n", argv[i]);
        status = MagickFalse;
      }
    } else if (IsArg (pa, "st", "similarity_threshold")==MagickTrue) {
      i++;
      pfh->Match.simThreshold = atof (argv[i]);
    } else if (IsArg (pa, "dt", "dissimilarity_threshold")==MagickTrue) {
      i++;
      pfh->Match.DissimOn = MagickTrue;
      pfh->Match.dissimThreshold = atof (argv[i]);
    } else if (IsArg (pa, "w", "write")==MagickTrue) {
      pfh->write_frames = 1;
      i++;
      CopyMagickString (pfh->write_filename, argv[i], MaxTextExtent);
      if (!*pfh->write_filename) {
        fprintf (stderr, "Invalid 'write' [%s]\n", argv[i]);
        status = MagickFalse;
      }
    } else if (IsArg (pa, "w2", "write2")==MagickTrue) {
      pfh->write_frames = 2;
      i++;
      CopyMagickString (pfh->write_filename, argv[i], MaxTextExtent);
      if (!*pfh->write_filename) {
        fprintf (stderr, "Invalid 'write' [%s]\n", argv[i]);
        status = MagickFalse;
      }
    } else if (IsArg (pa, "v", "verbose")==MagickTrue) {
      pfh->do_verbose = MagickTrue;
    } else if (IsArg (pa, "d", "debug")==MagickTrue) {
      pfh->debug = MagickTrue;
    } else {
      fprintf (stderr, "fillhole: ERROR: unknown option [%s]\n", pa);
      status = MagickFalse;
    }
  }

/*
  pfh->num_chan = 0;
  if (pfh->channels & CH_R) pfh->num_chan++;
  if (pfh->channels & CH_G) pfh->num_chan++;
  if (pfh->channels & CH_B) pfh->num_chan++;
  if (pfh->channels == 0) {
    fprintf (stderr, "No channels\n");
    status = MagickFalse;
  }
*/


  if (pfh->CopyRadPix==0 || pfh->CopyRadPix > pfh->WindowRadPix) {
    pfh->CopyRadPix = pfh->WindowRadPix;
  }
  if (pfh->copyWhat==copyOnePixel) pfh->CopyRadPix=0;

  pfh->Match.simThresholdSq = pfh->Match.simThreshold * pfh->Match.simThreshold;
  pfh->Match.dissimThresholdSq = pfh->Match.dissimThreshold * pfh->Match.dissimThreshold;

  pfh->CompWind.win_wi = pfh->sqDim;
  pfh->CompWind.win_ht = pfh->sqDim;

  if (pfh->do_verbose) {
    fprintf (stderr, "fillhole options: ");

    fprintf (stderr, "  window_radius %g", pfh->WindowRad);
    if (pfh->WindowRad_IsPc) fprintf (stderr, "%c", '%');

    fprintf (stderr, "  limit_search_radius %g", pfh->Match.LimSrchRad);
    if (pfh->Match.LimSrchRad_IsPc) fprintf (stderr, "%c", '%');

    fprintf (stderr, "  auto_lsr %s", pfh->Match.AutoLs ? "on" : "off");

    fprintf (stderr, "  search");
    switch (pfh->Match.searchWhat) {
      case swEntire: fprintf (stderr, " entire"); break;
      case swSkip:
        fprintf (stderr, " skip  skip_num %g", pfh->Match.SkipNumF);
        if (pfh->Match.SkipNum_IsPc) fprintf (stderr, "%c", '%');
        break;
      case swRandom:
        fprintf (stderr, " random  rand_searches %g", pfh->Match.RandSearchesF);
        if (pfh->Match.RandSearches_IsPc) fprintf (stderr, "%c", '%');
        break;
      default: fprintf (stderr, " ??");
    }

    fprintf (stderr, "  hom_chk ");
    if (pfh->CompWind.HomChkOn) {
      fprintf (stderr, "%g", pfh->CompWind.HomChk);
    } else {
      fprintf (stderr, "off");
    }

/*-
    if (pfh->max_iter) fprintf (stderr, " max_iterations %i", pfh->max_iter);

    if (pfh->channels != (CH_R | CH_G | CH_B)) {
      fprintf (stderr, " channels ");
      if (pfh->channels & CH_R) fprintf (stderr, "R");
      if (pfh->channels & CH_G) fprintf (stderr, "G");
      if (pfh->channels & CH_B) fprintf (stderr, "B");
    }
-*/

    fprintf (stderr, "  refine %s", pfh->Match.Refine ? "on" : "off");

    if (!pfh->Match.SearchToEdges) fprintf (stderr, "  search_to_edges off");

    fprintf (stderr, "  similarity_threshold %g", pfh->Match.simThreshold);

    if (pfh->Match.DissimOn == MagickTrue)
      fprintf (stderr, "  dissimilarity_threshold %g", pfh->Match.dissimThreshold);

    if (pfh->AutoRepeat) fprintf (stderr, "  auto_repeat");

    fprintf (stderr, "  copy ");
    if (pfh->copyWhat==copyOnePixel) fprintf (stderr, "onepixel");
    else if (pfh->copyWhat==copyWindow) fprintf (stderr, "window");
    else fprintf (stderr, "??");

    if (pfh->CopyRadPix != pfh->WindowRadPix) {
      fprintf (stderr, "  copy_radius %g", pfh->CopyRad);
      if (pfh->CopyRad_IsPc) fprintf (stderr, "%c", '%');
    }

    if (pfh->debug) fprintf (stderr, "  debug");

    if (pfh->do_verbose) fprintf (stderr, "  verbose");

    fprintf (stderr, "\n");
  }

  if (status == MagickFalse)
    usage ();

  return (status);
}



static MagickStatusType WriteFrame (
  FillHoleT *pfh,
  Image * img,
  int frame_num,
  ExceptionInfo *exception)
{
  ImageInfo
    *ii;

  MagickStatusType
    status;

  Image
    *copy_img;

  ii = AcquireImageInfo ();

  copy_img = CloneImage(img, 0, 0, MagickTrue, exception);
  if (copy_img == (Image *) NULL)
    return(MagickFalse); // FIXME: raise error

  copy_img->scene = frame_num;

  CopyMagickString (copy_img->filename, pfh->write_filename, MaxTextExtent);

#if IMV6OR7==6
  status = WriteImage (ii, copy_img);
#else
  status = WriteImage (ii, copy_img, exception);
#endif

  DestroyImageList(copy_img);

  ii = DestroyImageInfo (ii);

  return status;
}


/*===
static double CompareWindow (
  FillHoleT * pfh,
  CacheView * inp_view,
  CacheView * subimage_view,
  ExceptionInfo *exception,
  ssize_t hx, ssize_t hy,
  ssize_t sx, ssize_t sy)
// Compares subimage window with hole starting at top-left=hx,hy
// with window in inp_view (which may contain hole) starting at top-left=sx, sy.
// Returns positive number MSE; closer to zero is closer match.
// If sx,sy has same number or more transparent pixels as hx,hy,
//   or either window is entirely transparent,
//   returns very large number.
{
  ssize_t
    xy;

  const VIEW_PIX_PTR
    *sxy, *sxySave,
    *hxy, *hxySave;

  VIEW_PIX_PTR
    minSub,
    maxSub;

  double
    score = 0;

  int
    cntNonTrans = 0;

  //if (pfh->debug==MagickTrue) printf ("CompareWindow %i %i  %i %i  ", (int)hx, (int)hy, (int)sx, (int)sy);

  if (sx==hx && sy==hy)
    return pfh->CompWind.WorstCompare;

  sxy = GetCacheViewVirtualPixels(inp_view,sx,sy,pfh->sqDim,pfh->sqDim,exception);
  if (sxy == (const VIEW_PIX_PTR *) NULL) {
    printf ("bad sxy");
    return pfh->CompWind.WorstCompare;
  }
  // If central pixel in sxy is transparent, forget about it.
  if (GetPixelAlpha (sxy + pfh->WindowRadPix*(pfh->sqDim+1)) <= 0)
    return pfh->CompWind.WorstCompare;

  hxy = GetCacheViewVirtualPixels(subimage_view,hx,hy,pfh->sqDim,pfh->sqDim,exception);
  if (hxy == (const VIEW_PIX_PTR *) NULL) {
    printf ("bad hxy");
    return pfh->CompWind.WorstCompare;
  }

  sxySave = sxy;
  hxySave = hxy;
  instantPp (&minSub);
  instantPp (&maxSub);

  MagickBooleanType InitMinMax = MagickFalse;
  for (xy=0; xy < pfh->sqDim * pfh->sqDim; xy++) {
    double sa = GetPixelAlpha (sxy);
    double ha = GetPixelAlpha (hxy);

    if (sa <= 0 && ha > 0) return pfh->CompWind.WorstCompare;

    if (sa > 0 && ha > 0) {
      double d = (GetPixelRed(sxy)-GetPixelRed(hxy)) / QuantumRange;
      score += d*d;
      d = (GetPixelGreen(sxy)-GetPixelGreen(hxy)) / QuantumRange;
      score += d*d;
      d = (GetPixelBlue(sxy)-GetPixelBlue(hxy)) / QuantumRange;
      score += d*d;
      cntNonTrans++;
    }

    if (pfh->CompWind.HomChkOn && ha > 0) {
      if (!InitMinMax) {
        minSub = *hxy;
        maxSub = *hxy;
        InitMinMax = MagickTrue;
      } else {
        // FIXME: or could record limits of sxy instead of hxy?
        if (GetPixelRed(&minSub) > GetPixelRed(hxy)) SetPixelRed(&minSub, GetPixelRed(hxy));
        if (GetPixelGreen(&minSub) > GetPixelGreen(hxy)) SetPixelGreen(&minSub, GetPixelGreen(hxy));
        if (GetPixelBlue(&minSub) > GetPixelBlue(hxy)) SetPixelBlue(&minSub, GetPixelBlue(hxy));
        if (GetPixelRed(&maxSub) < GetPixelRed(hxy)) SetPixelRed(&maxSub, GetPixelRed(hxy));
        if (GetPixelGreen(&maxSub) < GetPixelGreen(hxy)) SetPixelGreen(&maxSub, GetPixelGreen(hxy));
        if (GetPixelBlue(&maxSub) < GetPixelBlue(hxy)) SetPixelBlue(&maxSub, GetPixelBlue(hxy));
      }
    }

    sxy++;
    hxy++;
  }

  if (cntNonTrans == 0) return pfh->CompWind.WorstCompare;

  if (pfh->CompWind.HomChkOn) {
    sxy = sxySave;
    hxy = hxySave;
    int outside = 0;

    VIEW_PIX_PTR tppMin, tppMax;

    double hcp1 = pfh->CompWind.HomChk + 1;

    SetPixelRed   (&tppMin, hcp1*GetPixelRed(&minSub)   - pfh->CompWind.HomChk*GetPixelRed(&maxSub));
    SetPixelGreen (&tppMin, hcp1*GetPixelGreen(&minSub) - pfh->CompWind.HomChk*GetPixelGreen(&maxSub));
    SetPixelBlue  (&tppMin, hcp1*GetPixelBlue(&minSub)  - pfh->CompWind.HomChk*GetPixelBlue(&maxSub));

    SetPixelRed   (&tppMax, hcp1*GetPixelRed(&maxSub)   - pfh->CompWind.HomChk*GetPixelRed(&minSub));
    SetPixelGreen (&tppMax, hcp1*GetPixelGreen(&maxSub) - pfh->CompWind.HomChk*GetPixelGreen(&minSub));
    SetPixelBlue  (&tppMax, hcp1*GetPixelBlue(&maxSub)  - pfh->CompWind.HomChk*GetPixelBlue(&minSub));

    for (xy=0; xy < pfh->sqDim * pfh->sqDim; xy++) {
      double sa = GetPixelAlpha (sxy);
      double ha = GetPixelAlpha (hxy);

      if (sa > 0 && ha <= 0) {
        if      (GetPixelRed (sxy) < GetPixelRed (&tppMin)) outside++;
        else if (GetPixelRed (sxy) > GetPixelRed (&tppMax)) outside++;

        if      (GetPixelGreen (sxy) < GetPixelGreen (&tppMin)) outside++;
        else if (GetPixelGreen (sxy) > GetPixelGreen (&tppMax)) outside++;

        if      (GetPixelBlue (sxy) < GetPixelBlue (&tppMin)) outside++;
        else if (GetPixelBlue (sxy) > GetPixelBlue (&tppMax)) outside++;
      }
    }

    if (outside > 0) return pfh->CompWind.WorstCompare;
    // FIXME? Seems a bit harsh. Maybe make score just a bit worse.
  }

  score /= (cntNonTrans);

  return score;
}
===*/

static void inline CopyOnePix (
  FillHoleT * pfh,
  Image * image,
  CacheView * new_view,
  CacheView * inp_view,
  CopyWhereT *cw,
  ExceptionInfo *exception)
// Copies one pixel from inp_view to new_view.
// Adjusts cw coords.
{
  const VIEW_PIX_PTR
    *src;

  VIEW_PIX_PTR
    *dst;

  cw->srcX += pfh->WindowRadPix;
  cw->srcY += pfh->WindowRadPix;
  cw->dstX += pfh->WindowRadPix;
  cw->dstY += pfh->WindowRadPix;

  if (pfh->debug==MagickTrue) printf ("CopyOnePix %li,%li => %li,%li\n",
    cw->srcX, cw->srcY,
    cw->dstX, cw->dstY);

  src = GetCacheViewVirtualPixels(
    inp_view,cw->srcX,cw->srcY,1,1,exception);
  if (src == (const VIEW_PIX_PTR *) NULL) printf ("bad src");

  dst = GetCacheViewAuthenticPixels(
    new_view,cw->dstX,cw->dstY,1,1,exception);
  if (dst == (const VIEW_PIX_PTR *) NULL) printf ("bad dst");

  SET_PIXEL_RED   (image, GET_PIXEL_RED   (image, src), dst);
  SET_PIXEL_GREEN (image, GET_PIXEL_GREEN (image, src), dst);
  SET_PIXEL_BLUE  (image, GET_PIXEL_BLUE  (image, src), dst);
  SET_PIXEL_ALPHA (image, QuantumRange, dst);
}


static void inline CopyWindow (
  FillHoleT * pfh,
  Image * image,
  CacheView * new_view,
  CacheView * inp_view,
  CopyWhereT *cw,
  ExceptionInfo *exception)
// Copies non-transparent pixels from inp_view to replace transparent pixels in new_view.
// If CopyRad < WindowRad, or if out of bounds, the cw coords will be adjusted.
{
  const VIEW_PIX_PTR
    *src;

  VIEW_PIX_PTR
    *dst;

  int
    ij;

  // dst pixels could be outside authentic
  // We need to adjust GCV parameters so we are not outside authentic.

  // FIXME: Check somewhere for case of sqDim <= image dimensions.

  if (pfh->CopyRadPix < pfh->WindowRadPix) {
    int dRad = pfh->WindowRadPix - pfh->CopyRadPix;
    cw->srcX += dRad;
    cw->srcY += dRad;
    cw->dstX += dRad;
    cw->dstY += dRad;
  }

  int offsX = 0;
  int offsY = 0;

  if (cw->dstX < 0) offsX = -cw->dstX;
  if (cw->dstY < 0) offsY = -cw->dstY;

  cw->wi = pfh->sqCopyDim - offsX;
  cw->ht = pfh->sqCopyDim - offsY;

  cw->srcX += offsX;
  cw->srcY += offsY;

  cw->dstX += offsX;
  cw->dstY += offsY;

  // FIXME: chk next 22-Nov-2015:
  if (cw->wi + cw->dstX > pfh->Width)  cw->wi = pfh->Width - cw->dstX;
  if (cw->ht + cw->dstY > pfh->Height) cw->ht = pfh->Height - cw->dstY;

  if (pfh->debug==MagickTrue) printf ("CW offsXY=%i,%i dstXY=%li,%li wi=%li ht=%li\n", offsX, offsY, cw->dstX, cw->dstY, cw->wi, cw->ht);

  assert (cw->wi > 0 && cw->wi <= pfh->sqCopyDim);
  assert (cw->ht > 0 && cw->ht <= pfh->sqCopyDim);

  src = GetCacheViewVirtualPixels(inp_view,cw->srcX,cw->srcY,cw->wi,cw->ht,exception);
  if (src == (const VIEW_PIX_PTR *) NULL) { printf ("bad src"); return; }

  dst = GetCacheViewAuthenticPixels(new_view,cw->dstX,cw->dstY,cw->wi,cw->ht,exception);
  if (dst == (const VIEW_PIX_PTR *) NULL) { printf ("bad dst"); return; }

  // FIXME: "image" in following may be wrong.
  for (ij=0; ij < cw->wi*cw->ht; ij++) {
    if (GET_PIXEL_ALPHA (image, src) > 0 && GET_PIXEL_ALPHA (image, dst) == 0) {
      SET_PIXEL_RED   (image, GET_PIXEL_RED (image, src), dst);
      SET_PIXEL_GREEN (image, GET_PIXEL_GREEN (image, src), dst);
      SET_PIXEL_BLUE  (image, GET_PIXEL_BLUE (image, src), dst);
      SET_PIXEL_ALPHA (image, QuantumRange, dst);
    }
    src += Inc_ViewPixPtr (image);
    dst += Inc_ViewPixPtr (image);
  }

  if (pfh->debug==MagickTrue) printf ("doneCW\n");
}


static void ResolveImageParams (
  const Image *image,
  FillHoleT * pfh
)
{
  pfh->Width = image->columns;
  pfh->Height = image->rows;

  pfh->Match.ref_columns = image->columns;
  pfh->Match.ref_rows = image->rows;

  pfh->ImgDiagSq = image->rows * image->rows + image->columns * image->columns;

  if (pfh->WindowRad_IsPc)
    pfh->WindowRadPix = (0.5 + pfh->WindowRad/100.0
                         * (pfh->Width<pfh->Height ? pfh->Width:pfh->Height));

  if (pfh->WindowRadPix < 1) pfh->WindowRadPix = 1;

  pfh->Match.matchRadX = pfh->Match.matchRadY = pfh->WindowRadPix;

  if (pfh->CopyRad_IsPc)
    pfh->CopyRadPix = (0.5 + pfh->CopyRad/100.0
                         * (pfh->Width<pfh->Height ? pfh->Width:pfh->Height));

  if (pfh->CopyRadPix < 0) pfh->CopyRadPix = 0;

  if (pfh->Match.LimSrchRad_IsPc)
    pfh->Match.limSrchRadX = pfh->Match.limSrchRadY = (0.5 + pfh->Match.LimSrchRad/100.0
                         * (pfh->Width<pfh->Height ? pfh->Width:pfh->Height));

  pfh->sqDim = 2 * pfh->WindowRadPix + 1;
  pfh->sqCopyDim = 2 * pfh->CopyRadPix + 1;

  pfh->CompWind.win_wi = pfh->sqDim;
  pfh->CompWind.win_ht = pfh->sqDim;

//  pfh->WorstCompare = (4*QuantumRange*QuantumRange * pfh->sqDim * pfh->sqDim);
//  pfh->WorstCompare *= (pfh->ImgDiagSq * 2);

  pfh->CompWind.WorstCompare = (4*QuantumRange*QuantumRange * pfh->sqDim * pfh->sqDim);
  pfh->CompWind.WorstCompare *= (pfh->ImgDiagSq * 2);

  if (pfh->Match.DissimOn == MagickFalse) {
    pfh->Match.dissimThreshold = pfh->CompWind.WorstCompare;
    pfh->Match.dissimThresholdSq = pfh->CompWind.WorstCompare;
  }

  pfh->CompWind.AllowEquCoord = MagickFalse;

  if (pfh->do_verbose) {
    fprintf (stderr, "pfh->WindowRadPix = %i\n", pfh->WindowRadPix);
    fprintf (stderr, "pfh->Match.limSrchRad = %lix%li\n",
      pfh->Match.limSrchRadX, pfh->Match.limSrchRadY);
    fprintf (stderr, "pfh->CopyRadPix = %i\n", pfh->CopyRadPix);
    fprintf (stderr, "pfh->CompWind.WorstCompare = %g\n", pfh->CompWind.WorstCompare);
    fprintf (stderr, "pfh->sqDim = %i\n", pfh->sqDim);
    fprintf (stderr, "pfh->sqCopyDim = %i\n", pfh->sqCopyDim);
    fprintf (stderr, "pfh->Match.RandSearchesI = %i\n", pfh->Match.RandSearchesI);
    fprintf (stderr, "pfh->Match.SkipNumI = %i\n", pfh->Match.SkipNumI);
    fprintf (stderr, "pfh->Match.simThresholdSq = %g\n", pfh->Match.simThresholdSq);
    fprintf (stderr, "pfh->Match.dissimThresholdSq = %g\n", pfh->Match.dissimThresholdSq);
  }
}


static MagickBooleanType MatchAndCopy (
  FillHoleT * pfh,
  Image * new_image,
  CacheView *new_view,        // destination for changes
  CacheView *copy_inp_view,   // source for changes
  ssize_t x,
  ssize_t y,
  CopyWhereT *cpwh,
  int *frameNum,
  MagickBooleanType *unfilled,
  MagickBooleanType *changedAny,
  ExceptionInfo *exception)
// x,y is central pixel of window to be matched.
// Returns true if okay, or false if serious problem.
{
  MagickBooleanType
    status = MagickTrue;

  if (pfh->debug==MagickTrue) printf ("MaC %li,%li \n", x, y);

  MatchT * pm = &pfh->Match;
  CompWindT * cmpwin = &pfh->CompWind;

  Match (x, y, x, y, pm, cmpwin, pfh->random_info, exception);

  if (pfh->debug==MagickTrue) printf ("Mac Best ref xy=%li,%li %g  from sub xy %li,%li\n",
      pm->bestX, pm->bestY, pm->bestScore, cmpwin->subX, cmpwin->subY);

  if (pm->bestScore < pm->dissimThresholdSq) {

    // Copy functions expect these are the top-left of the search windows.
    cpwh->srcX = pm->bestX - pm->matchRadX;
    cpwh->srcY = pm->bestY - pm->matchRadY;
    cpwh->dstX = cmpwin->subX;
    cpwh->dstY = cmpwin->subY;

    // Note: cpwh->wi and cpwh->ht is the window for copy, not search.

    if (pfh->debug==MagickTrue) printf (" MaC cpwh: %li,%li %lix%li+%li+%li\n",
      cpwh->srcX, cpwh->srcY, cpwh->wi, cpwh->ht, cpwh->dstX, cpwh->dstY);

    if (pfh->copyWhat == copyOnePixel)
      CopyOnePix (pfh, new_image, new_view, copy_inp_view, cpwh, exception);
    else
      CopyWindow (pfh, new_image, new_view, copy_inp_view, cpwh, exception);

    if (SyncCacheViewAuthenticPixels (new_view,exception) == MagickFalse)
      status=MagickFalse;

    *changedAny = MagickTrue;
    if (pfh->write_frames == 2) {
      WriteFrame (pfh, new_image, *frameNum, exception);
      (*frameNum)++;
    }
  } else {
    *unfilled = MagickTrue;
  }
  if (pfh->debug==MagickTrue) printf ("doneMac\n");
  return (status);
}

static void InitRand (FillHoleT * pfh)
{
  pfh->random_info=AcquireRandomInfo();

  // There seems to be a problem: the first few values show coherency,
  // so skip over them.

  int i;
  for (i=0; i < 20; i++) {
    GetPseudoRandomValue(pfh->random_info);
  }
}

static void DeInitRand (FillHoleT * pfh)
{
  pfh->random_info=DestroyRandomInfo(pfh->random_info);
}

static void InitAutoLimit (FillHoleT * pfh)
{
  pfh->Match.SetAutoLimit = MagickTrue;
  pfh->Match.UseAutoLimit = MagickFalse;

  pfh->Match.limLeft = pfh->Width;
  pfh->Match.limTop = pfh->Height;
  pfh->Match.limRight = 0;
  pfh->Match.limBot = 0;

  pfh->Match.nLimitCnt = 0;
}

static void UseAutoLimit (FillHoleT * pfh)
{
  pfh->Match.SetAutoLimit = MagickFalse;
  pfh->Match.UseAutoLimit = pfh->Match.AutoLs;

  if (pfh->debug==MagickTrue)
    printf ("UseAutoLimit: LTRB %li,%li %li,%li\n",
      pfh->Match.limLeft, pfh->Match.limTop, pfh->Match.limRight, pfh->Match.limBot);
}

All images on this page were created by the commands shown.

To improve internet download speeds, some images may have been automatically converted (by ImageMagick, of course) from PNG or TIFF or MIFF to JPG.

My usual version of IM is:

%IMG7%magick -version

Version: ImageMagick 7.1.1-15 Q16-HDRI x64 a0a5f3d:20230730 https://imagemagick.org
Copyright: (C) 1999 ImageMagick Studio LLC
License: https://imagemagick.org/script/license.php
Features: Cipher DPC HDRI OpenCL OpenMP(2.0) 
Delegates (built-in): bzlib cairo freetype gslib heic jng jp2 jpeg jxl lcms lqr lzma openexr pangocairo png ps raqm raw rsvg tiff webp xml zip zlib
Compiler: Visual Studio 2022 (193532217)

%IM7DEV%magick -version

Version: ImageMagick 7.1.1-13 (Beta) Q32-HDRI x86_64 a8de149e1:20230703 https://imagemagick.org
Copyright: (C) 1999 ImageMagick Studio LLC
License: https://imagemagick.org/script/license.php
Features: Cipher DPC HDRI Modules OpenCL OpenMP(4.5) 
Delegates (built-in): bzlib cairo fftw fontconfig freetype heic jbig jng jpeg lcms ltdl lzma pangocairo png raqm raw rsvg tiff webp wmf x xml zip zlib
Compiler: gcc (11.3)

This customised development version is:

%IM7DEV%magick identify -version

Version: ImageMagick 7.1.1-13 (Beta) Q32-HDRI x86_64 a8de149e1:20230703 https://imagemagick.org
Copyright: (C) 1999 ImageMagick Studio LLC
License: https://imagemagick.org/script/license.php
Features: Cipher DPC HDRI Modules OpenCL OpenMP(4.5) 
Delegates (built-in): bzlib cairo fftw fontconfig freetype heic jbig jng jpeg lcms ltdl lzma pangocairo png raqm raw rsvg tiff webp wmf x xml zip zlib
Compiler: gcc (11.3)

Source file for this web page is fillholespri.h1. To re-create this web page, execute "procH1 fillholespri".

This page, including the images, is my copyright. Anyone is permitted to use or adapt any of the code, scripts or images for any purpose, including commercial use.

Anyone is permitted to re-publish this page, but only for non-commercial use.

Anyone is permitted to link to this page, including for commercial use.

Page version v1.0 29-Nov-2015.

Page created 29-Sep-2023 06:39:43.