snibgo's ImageMagick pages

GrowCut segmentation

Images can be segmented (broken into pieces, clustered) according to the similarity of adjacent pixels. The GrowCut method grows clusters from pre-defined seeds according to the strength (or weighting) of cluster-membership.

• Basic algorithm
• Arguments
• Working practice
• The algorithm
• 3-D objects
• Animation
• Visualising the strength
• Automatic swiping
• Large images
• Future possibilities
• Scripts

The module entirely ignores the alpha channel of the input image. It ignores the general "-channel" setting (but has its own private "channel" setting).

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

For an alternative segmentation method, see Canny edge detection: Segmenting.

Basic algorithm

The process takes two inputs and creates one output. The first input is an image to be segmented; the second is a "swipe" image, a mask of seed points for clusters. The algorithm grows clusters until no more changes are made.

Segmenting this image ...	toes.png
... from this swipe ...	seg_toes_swipe.png
... makes this result: %IMDEV%convert ^ toes.png ^ seg_toes_swipe.png ^ -virtual-pixel Mirror ^ -process growcut ^ -background White ^ -layers flatten ^ gc_gc1.png Resulting pixels are exactly red or green, but with an alpha representing the strength. Flattening against white ensures uncertain areas show up.
Blurring the input image improves the result. %IMDEV%convert ^ toes.png ^ -blur 0x3 ^ seg_toes_swipe.png ^ -virtual-pixel Mirror ^ -process growcut ^ -background White ^ -layers flatten ^ gc_gc1b.png

The actual colours used for the swipe are not important. I use colours that bear some relation to the scene, but that is just my convention.

The output image has each pixel taking its colour from one of the swipes. The alpha channel represents the "strength" of the result, or how certain the algorithm is that the pixel is in that cluster. Every alpha will be greater than zero, but this might be rounded to zero in the final output image file.

We refine the result by painting more swipes in Gimp. (For this web page, I put them in a second layer.)

This is how it looks in Gimp:
%IMDEV%convert ^ toes.png ^ ( seg_toes_swipe.png ^ seg_toes_swipe2.png ^ -composite ^ ) ^ -virtual-pixel Mirror ^ -process 'growcut' ^ gc_gc2.png
Repeat, but with a blur. %IMDEV%convert ^ toes.png ^ -blur 0x3 ^ ( seg_toes_swipe.png ^ seg_toes_swipe2.png ^ -composite ^ ) ^ -virtual-pixel Mirror ^ -process 'growcut' ^ gc_gc2b.png
We can use any of the swipe colours for an opacity mask. %IM%convert ^ toes.png ^ ( gc_gc2b.png ^ -alpha off ^ -fill Black +opaque Red ^ -fill White -opaque Red ^ ) ^ -compose CopyOpacity -composite ^ gc_gc3b.png

In the examples on this page, swipes are fully opaque. They could be partially transparent, so they would act merely as hints.

Arguments

The process can take some arguments:

Samples using options:

%IMDEV%convert ^ toes.png ^ ( seg_toes_swipe.png ^ seg_toes_swipe2.png ^ -composite ^ ) ^ -virtual-pixel Mirror ^ -process 'growcut verbose' ^ gc_gc2v1.png growcut options: verbose numImages 2 threshold 0 OpenViews... SetPixPack... BasicAlgorithm... iter_num=45
%IMDEV%convert ^ toes.png ^ ( seg_toes_swipe.png ^ seg_toes_swipe2.png ^ -composite ^ ) ^ -virtual-pixel Mirror ^ -process 'growcut n_connected 4 verbose' ^ gc_gc2v2.png growcut options: n_connected 4 verbose numImages 2 threshold 0 OpenViews... SetPixPack... BasicAlgorithm... iter_num=59
%IMDEV%convert ^ toes.png ^ ( seg_toes_swipe.png ^ seg_toes_swipe2.png ^ -composite ^ ) ^ -virtual-pixel Mirror ^ -process 'growcut max_iterations 10 verbose2' ^ gc_gc2v3.png growcut options: max_iterations 10 verbose2 numImages 2 threshold 0 OpenViews... SetPixPack... BasicAlgorithm... num_ch = 4554 num_ch = 5287 num_ch = 6130 num_ch = 6895 num_ch = 7643 num_ch = 8335 num_ch = 8962 num_ch = 9675 num_ch = 10064 num_ch = 10055 iter_num=10

Working practice

In Gimp I create a file seg_toes.xcf from toes.png as a background layer, and create a second transparent layer. I name this layer "swipe", and draw a couple of coloured swipes. I extract the layers with extrXcfLayers (see Gimp and IM: From Gimp XCF to IM):

call %PICTBAT%extrXcfLayers seg_toes.xcf

I typically have Gimp open in one window, and a command window to run the command, and the output displayed in another window. If Explorer is open for the output directory, it notifies Microsoft Photo Viewer when its image changes. After editing the swipe, I press ctrl-S in Gimp to save the XCF file, then run the command that extracts the swipe mask from the XCF and runs convert, and the result appears in the viewer.

For this web page, I use two layers for swipes so I can show the effect of adding more swipes. Normally, I have all swipes in just one Gimp layer.

The algorithm

References:

• "GrowCut" - Interactive Multi-Label N-D Image Segmentation By Cellular Automata, Vladimir Vezhnevets and Vadim Konouchine, 2005?
• GrowCut Segmentation In Matlab, Shawn Lankton Online, 6 March 2008.

Regions are grown from the seed points, which are defined by swipe colours. The algorithm assigns pixels to one of the swipe colours. If there are (n) different swipe colours, there will be exactly (n) clusters.

The algorithm works like this:

1. Each pixel has a label (from the colour of the corresponding pixel in the swipe mask) and a strength (scale 0 to 1, from the alpha channel of the swipe). The strength is a form of weighting.
2. For each pixel, examine the 8 (or 4) neighbours. For each neighbour, calculate the colour-difference (scale 0 to 1) from the pixel. If the difference multiplied by the neighbour's strength is greater than the pixel's strength, then the pixel is labelled from the neighbour, and its strength becomes the neighbour's strength multiplied by the difference.

   if (strength(p) < diff(p,n) * strength(n) {
     new_label(p) = label(n);
     new_strength(p) = diff(p,n) * strength(n);
   }

3. Step two is repeated until it makes no more changes.

In the algorithm and implementation in the references above, a pixel takes its new label and strength from the neighbour that was examined last.

Some consequences from the algorithm:

• If the algorithm starts with no labels, there will be no clusters.
• If the algorithm starts with (n) different labels, all with strength = 1.0, it will end with (n) clusters.
• When the algorithm ends, each pixel will be in a cluster.
• Each pixel may be in a cluster before the algorithm ends. As the algorithm continues, strengths and labels of pixels may change.
• Clusters may not be connected. If the clusters are connected in the swipe mask, they will be connected in the result.
• A pixel with a strength=1 is guaranteed to not change cluster.
• The algorithm isn't fast. The time is proportional to the image size multiplied by the number of iterations. More swipes will reduce the time.

Perhaps the speed could be improved by first working on a smaller version of the image, then taking the result as the new swipe, after making pixels near a cluster's boundary transparent.

Vladimir Vezhnevets's paper proposes two additional rules for smoother boundaries. Each rule has a threshold. Defines a pixel's "enemies" as the neighbours (up to 8) that are of a different label to the pixel. E(p) is the number of enemies of pixel p.

1. If E(p) ≥ T₁, it doesn't attack any neighbour.
2. If E(p) ≥ T₂, it is occupied by its weakest neighbour.

Rule 2 is a simple modification of the main algorithm. Rule 1 isn't so simple. If implemented as a modification of the main algorithm, it requires that the neighbours of a pixel's neighbours are examined; very slow. Better if implemented as a separate pass prior to the main algorithm.

I haven't implemented either of those rules. Instead, the label and strength of a pixel is taken from its strongest neighbour (instead of from the last neighbour examined). This rule creates smoother boundaries, and also makes the algorithm converge more quickly.

(If you want to use the last neighbour examined, remove the code prev_diff < diff_str from the if test, around line 580.)

3-D objects

Stripping out the lightness is useful for segmentation, especially when the objects are rounded, the light is directional, and there are shadows. This is what the image looks like with constant lightness:

%IMDEV%convert ^ toes.png ^ -colorspace Lab ^ -channel R -evaluate set 50%% +channel ^ -colorspace sRGB ^ gc_no_light.png

Visually, it is hard to read this image. The blades of grass have visually merged together.

We can tell the module to segment the L*a*b* image using all three channels, or to ignore the first channel (L*) in comparisons, using just channels a* and b*.

%IMDEV%convert ^ toes.png ^ -colorspace Lab ^ seg_toes_swipe.png ^ -virtual-pixel Mirror ^ -process growcut ^ -background White ^ -layers flatten ^ -set colorspace sRGB ^ gc_lab.png
%IMDEV%convert ^ toes.png ^ -colorspace Lab ^ seg_toes_swipe.png ^ -virtual-pixel Mirror ^ -process 'growcut channel ab' ^ -background White ^ -layers flatten ^ -set colorspace sRGB ^ gc_ab1.png

Animation

The write option writes each iteration to an image file. The filename should include a printf-style escape, such as %03 (for Windows, double the %.) This always creates one file per image.

Each image is written as soon as it is created.

%IMDEV%convert ^ toes.png ^ seg_toes_swipe.png ^ -virtual-pixel Mirror ^ -process 'growcut write gc_f_%%06d.png' ^ NULL: %IM%convert ^ gc_f_*.png ^ NULL: ^ toes.png ^ -compose DstOver -layers composite ^ -layers optimize ^ gc_anim.gif del gc_f_*.png

The process module could be changed to keep all the frames in memory.

Visualising the strength

We can visualise the strengths generated by the GrowCut process, where lighter = stronger.

Extract the generated alpha. %IMDEV%convert ^ toes.png ^ seg_toes_swipe.png ^ -virtual-pixel Mirror ^ -process growcut ^ -alpha extract ^ -auto-level ^ gc_ax1.png

Greater strengths (greater certainty) is lighter; lower strength (lower certainty) is darker.

The darkest pixel is at the bottom edge, coordinate (158,228), near the red/green boundary. The corresponding pixel in toes.png has the colour #B4A09608A315. We can add a pixel of this colour in this position to the swipe mask.

Add a third point to the swipe. %IMDEV%convert ^ toes.png ^ ( seg_toes_swipe.png ^ -fill #B4A09608A315 -draw "point 158,228" ^ ) ^ -virtual-pixel Mirror ^ -process growcut ^ -write gc_ap1.png ^ -alpha extract ^ -auto-level ^ gc_apx1.png

Declaring a swipe at the most uncertain point is logically dubious because uncertain areas occur at boundaries, where we wouldn't expect to find a seed point for a cluster. But we might take each found cluster, find the centre of each, and use that for a swipe. For "centre", we have the choice of centroid (which isn't guaranteed to be within the cluster, or a morphology gradient (which is guaranteed to be within the cluster, and at the greatest distance from a cluster edge).

%IMDEV%convert ^
  gc_ap1.png ^
  -alpha off -unique-colors txt: 

# ImageMagick pixel enumeration: 3,1,4294967295,srgb
0,0: (4.29497e+09,0,0)  #FF0000  red
1,0: (0,4.29497e+09,0)  #00FF00  lime
2,0: (3.03174e+09,2.50961e+09,2.72857e+09)  #B495A2  srgb(180,149,162)

For example, to find the centre of the red cluster:

%IM%convert ^ gc_ap1.png ^ -alpha off ^ -fill black +opaque Red -fill white -opaque Red ^ -bordercolor black -border 1 -morphology Distance Euclidean:7 -shave 1x1 ^ -auto-level ^ gc_mph1.png

We consider the first white point to be the centre. We could loop through all the clusters, finding the centre of each.

Automatic swiping

The process module needs a swipe mask. This can be created manually, or automatically generated from an image. One of the techniques in Details, details can be used, for example finding a gradient magnitude, and selecting a few points with low gradients.

call %PICTBAT%slopeMag toes.png gc_toes_sm.png %IM%convert gc_toes_sm.png -negate gc_toes_sm.png set nlDEBUG=1 set nlPOINTSIZE= set nlSTROKEWIDTH= call %PICTBAT%nLightest gc_toes_sm.png 5 50 >gc_nl.lis set nlDEBUG= 236,88 40,0 22,174 211,209 118,100

From these coordinates, we create a list of draw commands of the colours from toes.png, and create the swipe mask.

call %PICTBAT%getDrawPoints toes.png gc_nl.lis gc_nl_draw.txt

fill #A76E8704840F point 236,88 
fill #5BA363E53FDB point 40,0 
fill #61C96B685E65 point 22,174 
fill #23D11FBD1F6D point 211,209 
fill #AA0D90408E30 point 118,100 

%IM%convert ^ toes.png ^ -alpha transparent ^ -alpha background ^ -draw @gc_nl_draw.txt ^ gc_auto_swipe.png -alpha background reduces the size of the file. The isolated colour pixels are difficult to see against the transparent background; we can make them more visible. %IM%convert ^ gc_auto_swipe.png ^ -blur 5x1000 -alpha off ^ gc_as_b.png

We use this automatic swipe for GrowCut:

Apply the automatic swipe. %IMDEV%convert ^ toes.png ^ -colorspace Lab ^ gc_auto_swipe.png ^ -virtual-pixel Mirror ^ -process growcut ^ -set colorspace sRGB ^ gc_as1.png Compose over the toes image: %IM%convert ^ toes.png ^ gc_as1.png ^ -composite ^ gc_as1_t.png

The grass is noisier than the toes, so the boundary between the two is within the grass. As above, blurring the image before GrowCut removes some grass detail and improves the result.

Blur before applying the automatic swipe. %IMDEV%convert ^ toes.png ^ -colorspace Lab ^ -blur 0x3 ^ gc_auto_swipe.png ^ -virtual-pixel Mirror ^ -process growcut ^ -set colorspace sRGB ^ gc_as2.png Compose over the toes image: %IM%convert ^ toes.png ^ gc_as2.png ^ -composite ^ gc_as2_t.png
A heavier blur before applying the automatic swipe. %IMDEV%convert ^ toes.png ^ -colorspace Lab ^ -blur 0x6 ^ gc_auto_swipe.png ^ -virtual-pixel Mirror ^ -process growcut ^ -set colorspace sRGB ^ gc_as3.png Compose over the toes image: %IM%convert ^ toes.png ^ gc_as3.png ^ -composite ^ gc_as3_t.png

If more that is known about the image (for example, it might be a photograph of an object against a roughly constant background that extends to the edges), then automatic swiping can be more specific.

Large images

Running time is proportional to the number of pixels in the image multiplied by the number of iterations. The number of iterations is proportional to the linear dimensions of the image.

For a 35 megapixel image, my computer takes about 4 seconds per iteration.

Future possibilities

Divider mask

A separate mask could be used to enforce a separation of clusters. This would be used when pixels from different but adjacent objects happen to be similar. The mask would be opaque black where a pixel is not to be compared with neigbours. Put it another way, that pixel would be considered at the maximum colour-difference from all its neighbours.

Aliasings

The algorithm assigns every pixel to exactly one cluster. In real photographs, pixels at the edge of objects are coloured by both the object and the background (they are anti-aliased). We could blur the result, but this is too simplistic. For each segment, we could mark outlying edge pixels, and outliers adjacent to those.

Fuzz

The process module ignores the current fuzz setting. It might use the fuzz setting, by subtracting this from the found colour-differences. This would have a similar effect to blurring the input image (but would be much faster).

Cluster mean

The algorithm compares each pixel with the colour of each neighbour. Instead, it might compare with the current average colour of the cluster to which the neighbour belongs. As this average would be updated as the algorithm proceeds, the result will depend on the order of the pixels. For example, starting at bottom-right would give a different result to starting from top-left. I have experimented with this, and don't like the results.

More usefully, an option might be provided to re-colour each segment to the average of the pixels in that cluster.

Scripts

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

growcut.c

/*
    Reference: im.snibgo.com/growcut.htm

    Updated:
      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"

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


#define FIRST_LINE_SHORTCUT 1

#if MAGICKCORE_QUANTUM_DEPTH==8
#define FIRST_LINE_SHORTCUT 0
#endif


typedef enum {
  a_scan_merge = 0,
  a_grow_cut = 1
} algorithmT;

typedef enum {
  direct = 0,
  clust_mean = 1
} comp_methodT;

typedef struct {
  MagickBooleanType
    do_verbose,
    do_verbose2,
    do_write_frames;

  int
    max_iter;

  float
    threshold, // like fuzz, eg 20%. Could come from "-fuzz".
    sq_threshold; // from 0.0 to 1.0

  unsigned int
    channels; // flags

  int
    num_chan;

  int
    precision,
    n_connected;

  int
    neighbour_offset[8];

  //ssize_t
  //  current_y;

  char
    write_filename[MaxTextExtent];

  // Input images
  Image
    *input_image,
    *swipe_image,
    *divide_image;

  // Output images
  Image
    *cluster_of, // R channel is cluster number of this pixel; B and G unused.
    *swipe_of;

  CacheView
    *input_view,
    *swipe_view,
    *divide_view;

  //const VIEW_PIX_PTR
  //  *input_pp3,
  //  *divide_pp;
} growcutT;


/*
We could output the unique colours of the cluster averages,
possible filtered by minimum size or the N largest clusters.
*/

static void usage (void)
{
  printf ("Usage: -process 'growcut [OPTION]...'\n");
  printf ("Segments an image.\n");
  printf ("\n");
  printf ("  i,  max_iterations N    maximum number of iterations (0=no limit)\n");
  printf ("  n,  n_connected N       connected (4 or 8)\n");
  printf ("  c,  channel string      which channels to use in comparisons\n");
  printf ("  v,  verbose             write text information to stdout\n");
  printf ("  v2, verbose2            write more 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 MagickBooleanType menu(
  const int argc,
  const char **argv,
  growcutT * ssm
)
// Returns MagickTrue if okay.
{
  int
    i;

  MagickBooleanType
    status;

  status = MagickTrue;

  ssm->max_iter = 0;
  ssm->n_connected = 8;
  ssm->channels = (CH_R | CH_G | CH_B);
  ssm->do_write_frames = MagickFalse;
  ssm->write_filename[0] = '\0';
  ssm->do_verbose = ssm->do_verbose2 = MagickFalse;

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

    if (IsArg (pa, "i", "max_iterations")==MagickTrue) {
      i++;
      if (!isdigit ((int)*argv[i]))
      {
        fprintf (stderr, "Non-numeric argument to 'max_iterations'.\n");
        status = MagickFalse;
      }
      else ssm->max_iter = (comp_methodT)atoi(argv[i]);
    } else if (IsArg (pa, "n", "n_connected")==MagickTrue) {
      i++;
      if (!isdigit ((int)*argv[i]))
      {
        fprintf (stderr, "Non-numeric argument to 'n_connected'.\n");
        status = MagickFalse;
      }
      else ssm->n_connected = (comp_methodT)atoi(argv[i]);
    } else if (IsArg (pa, "c", "channel")==MagickTrue) {
      i++;
      ssm->channels = 0;
      const char * p = argv[i];
      while (*p) {
        switch (toupper ((int)*p)) {
          case 'R':
            ssm->channels |= CH_R;
            break;
          case 'G':
            ssm->channels |= CH_G;
            break;
          case 'B':
            ssm->channels |= CH_B;
            break;
          case 'L':
            ssm->channels |= CH_R;
            break;
          case 'A':
            ssm->channels |= CH_G;
            break;
          default:
            fprintf (stderr, "Invalid 'channels' [%s]\n", argv[i]);
            status = MagickFalse;
        }
        p++;
      }
    } else if (IsArg (pa, "w", "write")==MagickTrue) {
      ssm->do_write_frames = MagickTrue;
      i++;
      CopyMagickString (ssm->write_filename, argv[i], MaxTextExtent);
      if (!*ssm->write_filename) {
        fprintf (stderr, "Invalid 'write' [%s]\n", argv[i]);
        status = MagickFalse;
      }
    } else if (IsArg (pa, "v", "verbose")==MagickTrue) {
      ssm->do_verbose = MagickTrue;
    } else if (IsArg (pa, "v2", "verbose2")==MagickTrue) {
      ssm->do_verbose = MagickTrue;
      ssm->do_verbose2 = MagickTrue;
    } else {
      fprintf (stderr, "growcut: ERROR: unknown option [%s]\n", pa);
      status = MagickFalse;
    }
  }

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

  if (ssm->do_verbose) {
    fprintf (stderr, "growcut options: ");

    if (ssm->max_iter) fprintf (stderr, " max_iterations %i", ssm->max_iter);

    if (ssm->n_connected != 8) fprintf (stderr, " n_connected %i", ssm->n_connected);

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

    if (ssm->do_verbose2) fprintf (stderr, " verbose2");
    else {
      if (ssm->do_verbose) fprintf (stderr, " verbose");
    }
    fprintf (stderr, "\n");
  }

  if (status == MagickFalse)
    usage ();

  return (status);
}


static MagickBooleanType InitSsm (
  Image **images,
  growcutT *pssm,
  ExceptionInfo *exception)
{
  int
   numImages,
   i;


  pssm->precision = GetMagickPrecision();

  pssm->threshold = 0;
  pssm->num_chan = 3;
  pssm->input_image = (Image *)NULL;
  pssm->swipe_image = (Image *)NULL;
  pssm->divide_image = (Image *)NULL;

  numImages = (int)GetImageListLength(*images);

  if (pssm->do_verbose) {
    fprintf (stderr, "  numImages %i\n", numImages);
  }

  if (numImages < 2) {
    fprintf (stderr, "GrowCut needs 2 or 3 images\n");
    return -1;
  }

  pssm->input_image = GetFirstImageInList(*images);

  if (pssm->input_image == (Image *)NULL) {
    fprintf (stderr, "Bad first image\n");
    return -1;
  }

  if (numImages >= 2) {
    pssm->swipe_image = GetNextImageInList(pssm->input_image);
    if (pssm->swipe_image == (Image *)NULL) {
      fprintf (stderr, "Bad swipe image\n");
      return -1;
    }
    if (pssm->swipe_image->columns != pssm->input_image->columns ||
        pssm->swipe_image->rows    != pssm->input_image->rows) {
      fprintf (stderr, "Swipe image dimensions don't match input image\n");
      return -1;
    }
  }

  if (numImages >= 3) {
    pssm->divide_image = GetNextImageInList(pssm->swipe_image);
    if (pssm->divide_image == (Image *)NULL) {
      fprintf (stderr, "Bad divide image\n");
      return -1;
    }
    if (pssm->divide_image->columns != pssm->input_image->columns ||
        pssm->divide_image->rows    != pssm->input_image->rows) {
      fprintf (stderr, "Divide image dimensions don't match input image\n");
      return -1;
    }
  }

  pssm->threshold = pssm->input_image->fuzz / QuantumRange;
  pssm->sq_threshold = pssm->threshold * pssm->threshold;

  if (pssm->do_verbose) {
    fprintf (stderr, "  threshold %.*g\n", pssm->precision, pssm->threshold);
    //fprintf (stderr, "sq_threshold %.*g\n", pssm->precision, pssm->sq_threshold);
  }

  // First four are 4-connected:
  // Neighbours are numbered:
  //   4 0 5
  //   2 . 3
  //   6 1 7

  {
    int n_x[] = {0,0,-1,+1, -1,+1,-1,+1};
    int n_y[] = {-1,+1,0,0, -1,-1,+1,+1};

    for (i = 0; i < 8; i++) {
      pssm->neighbour_offset[i] = 
        (n_y[i] + 1) * (pssm->input_image->columns + 2) + (n_x[i] + 1);

      //printf ("offset[%i]=%i\n", i, pssm->neighbour_offset[i]);
    }
  }

  return MagickTrue;
}



/*
   In the mask,
   a _black_ pixel prevents any "<= threshold" test at that pixel from succeeding.
*/

static void OpenViews (growcutT *pssm, ExceptionInfo *exception)
{
  pssm->input_view=AcquireVirtualCacheView (pssm->input_image, exception);
  if (pssm->swipe_image)  pssm->swipe_view=AcquireVirtualCacheView (pssm->swipe_image, exception);
  if (pssm->divide_image) pssm->divide_view=AcquireVirtualCacheView (pssm->divide_image, exception);
}

static void CloseViews (growcutT *pssm)
{
  if (pssm->divide_image) pssm->divide_view = DestroyCacheView (pssm->divide_view);
  if (pssm->swipe_image)  pssm->swipe_view = DestroyCacheView (pssm->swipe_view);
  pssm->input_view = DestroyCacheView (pssm->input_view);
}

static MagickBooleanType SyncViews (growcutT *pssm, ExceptionInfo *exception)
{
  return MagickTrue;
}

static MagickBooleanType SetPixPack (growcutT *pssm, ExceptionInfo *exception)
{
  //if (pssm->divide_image) {
  //  pssm->divide_pp=GetCacheViewVirtualPixels(pssm->divide_view,0,0,pssm->divide_image->columns,1,exception);
  //  if (pssm->divide_pp == (const VIEW_PIX_PTR *) NULL) return MagickFalse;
  //}
  return MagickTrue;
}

static MagickBooleanType SetThisRow (growcutT *pssm, ssize_t y, ExceptionInfo *exception)
{
  //pssm->input_pp3=GetCacheViewVirtualPixels(pssm->input_view,-1,y-1,pssm->input_image->columns+2,3,exception);
  //if (pssm->input_pp3 == (const VIEW_PIX_PTR *) NULL) return MagickFalse;

  //pssm->current_y = y;

  return MagickTrue;
}

static inline double PixPacketDiff (
  growcutT *pssm,
  const VIEW_PIX_PTR *p0,
  const VIEW_PIX_PTR *p1)
// Returns the difference, range 0.0 to 1.0.
{
  double
    d,
    diff;

  diff = 0;
  if (pssm->channels & CH_R) {
    d = (GET_PIXEL_RED(pssm->input_image, p1)   - GET_PIXEL_RED(pssm->input_image, p0)) / QuantumRange;
    diff = d * d;
  }
  if (pssm->channels & CH_G) {
    d = (GET_PIXEL_GREEN(pssm->input_image, p1) - GET_PIXEL_GREEN(pssm->input_image, p0)) / QuantumRange;
    diff += d * d;
  }
  if (pssm->channels & CH_B) {
    d = (GET_PIXEL_BLUE(pssm->input_image, p1)  - GET_PIXEL_BLUE(pssm->input_image, p0)) / QuantumRange;
    diff += d * d;
  }

  return sqrt (diff / pssm->num_chan);
}



/*---
static void chkpix (const char * msg, const VIEW_PIX_PTR *p)
{
  if (GetPixelRed (p) < 0)   printf ("%s Red < 0 ", msg);
  if (GetPixelGreen (p) < 0) printf ("%s Green < 0 ", msg);
  if (GetPixelBlue (p) < 0)  printf ("%s Blue < 0 ", msg);
}
---*/



static MagickStatusType WriteGcFrame (
  growcutT *pssm,
  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, pssm->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 inline double StrengthOf (const Image * image, const VIEW_PIX_PTR * p)
{
  double s = GET_PIXEL_ALPHA (image, p) / (double)QuantumRange;

  return (s <= 0) ? 0 : s;
}

/*---
static inline MagickBooleanType IsRgbaDiff(const VIEW_PIX_PTR * p0, const VIEW_PIX_PTR * p1)
{
  return (
    GetPixelRed (p0) != GET_PIXEL_RED(image,p1) ||
    GetPixelGreen (p0) != GET_PIXEL_GREEN(image,p1) ||
    GetPixelBlue (p0) != GET_PIXEL_BLUE(image,p1)
  );
}
---*/

static MagickBooleanType GrowCutIter (
  growcutT *pssm,
  Image **copy_swipe,
  ExceptionInfo *exception)
// Does one iteration of GrowCut.
// Replaces copy_swipe with next version of image.
// Returns whether a change was made.
{
  MagickBooleanType
    changed;

  int
    num_ch,
    cols_plus_3;

  Image
    *clone_swipe_out_image,
    *new_list;

  CacheView
    *in_view,
    *clone_swipe_out_view;

  ssize_t
    y;

  MagickBooleanType
    status = MagickTrue;


  in_view = AcquireVirtualCacheView(*copy_swipe,exception);

  // Clone the swipe image, same size, copied pixels:
  //
  clone_swipe_out_image=CloneImage(*copy_swipe, 0, 0, MagickTrue, exception);
  if (clone_swipe_out_image == (Image *) NULL)
    return(MagickFalse); // FIXME: raise error

  clone_swipe_out_view = AcquireAuthenticCacheView(clone_swipe_out_image,exception);


  changed = MagickFalse;

  num_ch = 0;

  cols_plus_3 = pssm->input_image->columns + 3;

#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static,4) shared(status) \
    MAGICK_THREADS(clone_swipe_out_image,clone_swipe_out_image,clone_swipe_out_image->rows,1)
#endif
  for (y=0; y < (ssize_t) pssm->input_image->rows; y++)
  {
    ssize_t
      x;

    int
      i;

    double
      strength_this,
      strength_neighbour,
      diff_str,
      prev_diff;

    const VIEW_PIX_PTR
      *pp_this_swipe,
      *pp_neighbour_swipe,
      *pp_this_pix,
      *pp_neighbour;

    const VIEW_PIX_PTR
      *pp_in3;

    VIEW_PIX_PTR
      *q_clone_swipe,
      *qq_out;

    const VIEW_PIX_PTR
      *input_pp3;

    if (!status) continue;
    //status = MagickTrue;

    if (!SetThisRow (pssm, y, exception)) {
      status = MagickFalse;
      continue;
      //return MagickFalse;
    }

    pp_in3 = GetCacheViewVirtualPixels(in_view,-1,y-1,(*copy_swipe)->columns+2,3,exception);
    if (!pp_in3) {
      status=MagickFalse;
      continue;
    }

    q_clone_swipe = GetCacheViewAuthenticPixels(clone_swipe_out_view,0,y,clone_swipe_out_image->columns,1,exception);
    if (!q_clone_swipe) {
      status=MagickFalse;
      continue;
    }

    input_pp3=GetCacheViewVirtualPixels(pssm->input_view,-1,y-1,pssm->input_image->columns+2,3,exception);
    if (!input_pp3) {
      status=MagickFalse;
      continue;
    }

    for (x=0; x < (ssize_t) pssm->input_image->columns; x++)
    {
      pp_this_swipe = pp_in3 + cols_plus_3 + x;
      strength_this = StrengthOf (*copy_swipe, pp_this_swipe);

      if (strength_this < 1.0) { // If pixel can't be changed, don't bother checking.
        pp_this_pix = input_pp3 + cols_plus_3 + x;
        qq_out = q_clone_swipe + x;
        prev_diff = 0;
        for (i = 0; i < pssm->n_connected; i++) {
          pp_neighbour_swipe = pp_in3 + pssm->neighbour_offset[i] + x;
          strength_neighbour = StrengthOf (*copy_swipe, pp_neighbour_swipe);

          // This test isn't needed, but improves performance x2.
          if (strength_neighbour > 0) {
            pp_neighbour = input_pp3 + pssm->neighbour_offset[i] + x;
            //diff = 1 - PixPacketDiff (pssm, pp_this_pix, pp_neighbour);
            diff_str = (1 - PixPacketDiff (pssm, pp_this_pix, pp_neighbour))
                     * strength_neighbour;
            if (    prev_diff < diff_str
                 && strength_this + 1.0e-7 < diff_str
               )
            {
              prev_diff = diff_str;
              SET_PIXEL_RED   (clone_swipe_out_image, GET_PIXEL_RED   (*copy_swipe, pp_neighbour_swipe), qq_out);
              SET_PIXEL_GREEN (clone_swipe_out_image, GET_PIXEL_GREEN (*copy_swipe, pp_neighbour_swipe), qq_out);
              SET_PIXEL_BLUE  (clone_swipe_out_image, GET_PIXEL_BLUE  (*copy_swipe, pp_neighbour_swipe), qq_out);
              SET_PIXEL_ALPHA (clone_swipe_out_image, diff_str * QuantumRange, qq_out); // FIXME: PLUS_HALF
              num_ch++;
              changed = MagickTrue;
            }
          }
        }
      }
    }
    if (!SyncCacheViewAuthenticPixels(clone_swipe_out_view,exception))
      status=MagickFalse;

    if (!SyncViews (pssm, exception)) {
      status=MagickFalse;
      continue;
    }

    //if (!status) break;
  }

  if (!status) return MagickFalse;

  clone_swipe_out_view=DestroyCacheView(clone_swipe_out_view);
  in_view=DestroyCacheView(in_view);

  DestroyImageList(*copy_swipe);

  new_list = NewImageList();
  AppendImageToList(&new_list, clone_swipe_out_image);

  *copy_swipe = GetFirstImageInList(new_list);

  if (pssm->do_verbose2) {
    fprintf (stderr, "  num_ch = %i\n", num_ch);
  }

  return changed;
}

static Image * GrowCut (
  growcutT *pssm,
  ExceptionInfo *exception)
{
  int
    iter_num;

  iter_num = 0;

  Image
    *copy_swipe;

  // Clone the swipe image, same size, copied pixels:
  //
  copy_swipe=CloneImage(pssm->swipe_image, 0, 0, MagickTrue, exception);
  if (copy_swipe == (Image *) NULL)
    return(MagickFalse); // FIXME: raise error

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

  while (GrowCutIter (pssm, &copy_swipe, exception)) {
    if (pssm->do_write_frames)
      if (!WriteGcFrame (pssm, copy_swipe, iter_num, exception)) break;

    iter_num++;
    if (pssm->max_iter && iter_num >= pssm->max_iter) break;
  }

  if (pssm->do_verbose) {
    fprintf (stderr, "iter_num=%i\n", iter_num);
  }

  return copy_swipe;
}


ModuleExport size_t growcutImage(Image **images,const int argc,
  const char **argv,ExceptionInfo *exception)
{
  Image
    *copy_image,
    *new_list;

  growcutT
    ssm;

  MagickBooleanType
    status;

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

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

  status = InitSsm (images, &ssm, exception);

  if (status != MagickTrue) {
    fprintf (stderr, "InitSsm failed\n");
    return (-1);
  }

  if (ssm.do_verbose) fprintf (stderr, "OpenViews...\n");
  OpenViews      (&ssm, exception);
  if (ssm.do_verbose) fprintf (stderr, "SetPixPack...\n");
  if (!SetPixPack     (&ssm, exception)) {
    fprintf (stderr, "SetPixPack failed");
    return (-1);
  }

  if (ssm.do_verbose) fprintf (stderr, "BasicAlgorithm...\n");

  copy_image = GrowCut (&ssm, exception);

  if (copy_image == (Image *)NULL) {
    fprintf (stderr, "BasicAlgorithm failed\n");
    return (-1);
  }

  CloseViews (&ssm);

  DestroyImageList(*images);

  new_list = NewImageList();

  AppendImageToList(&new_list, copy_image);

  *images = GetFirstImageInList(new_list);

  if (!status) return (-1);

  return(MagickImageFilterSignature);
}

%IM%identify -version

Version: ImageMagick 6.9.5-3 Q16 x86 2016-07-22 http://www.imagemagick.org
Copyright: Copyright (C) 1999-2015 ImageMagick Studio LLC
License: http://www.imagemagick.org/script/license.php
Visual C++: 180040629
Features: Cipher DPC Modules OpenMP 
Delegates (built-in): bzlib cairo flif freetype jng jp2 jpeg lcms lqr openexr pangocairo png ps rsvg tiff webp xml zlib

%IMDEV%identify -version

Version: ImageMagick 6.9.3-7 Q32 x86_64 2018-04-03 http://www.imagemagick.org
Copyright: Copyright (C) 1999-2016 ImageMagick Studio LLC
License: http://www.imagemagick.org/script/license.php
Features: Cipher DPC HDRI Modules OpenMP 
Delegates (built-in): bzlib cairo fftw fontconfig freetype fpx jbig jng jpeg lcms ltdl lzma pangocairo png rsvg tiff webp wmf x xml zlib

Page version v1.0 14-September-2014.

Page created 07-Apr-2018 15:28:06.

Copyright © 2018 Alan Gibson.